macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

pblat: マルチスレッドに対応したblat

2022/12/29 コマンド追記

2023/01/02 追記

 

 Blat [論文より ref.1 link]は、DNA、RNAおよびタンパク質配列をリファレンスゲノムにマッピングするように設計された配列アラインメントツールである。これは一般に、リファレンスゲノム内の配列の検索、closely relatedな種のゲノムからの相同配列の発見、mRNA配列からのエキソン - イントロンジャンクションの同定、および遺伝子構造の決定、ならびにゲノムおよびトランスクリプトーム配列の構築およびアノテーションを助けるために用いられる[ref.2]。 BWA [ref.3]やBowtie [ref.4]のような多くの高速シーケンスアライナーは、ハイスループットシーケンシングによって生成されたショートシーケンスリードマッピングするために開発されたが、ロングリードやギャップが豊富なシーケンスや不連続からのスプライスシーケンスのマッピングはできない。対照的に、blatは高感度と高精度を両立しており、そのようなアプリケーションに理想的なツールである[ref.67]。

 しかし、ハイスループットシークエンシングプロジェクトによって生成されるシークエンスの量が増えるにつれて、blatは大規模解析や定期的に更新されるアノテーションに必要な速度要件を満たすことができなくっている。例えば、脊椎動物の全トランスクリプトーム配列をリファレンスゲノムにマッピングするために使用される場合、blatを終了するまでに数日かかるだろう。これは、blatアルゴリズムがシングルスレッドであるため、最新のマルチコアプロセッサを最大限に活用できないためである。 GNU parallel [ref.8]ツールを使用して、1台以上のコンピューターを使用してblatの複数のインスタンスを並列に実行することはできる。しかしながら、各blatのプロセスは、全リファレンスゲノムのコピーをロードし、そしてゲノムのインデックスを構築してメモリに格納することになり、これは、複数のブラットプロセスが同時に実行される場合、従来のコンピュータの利用可能な物理メモリを超える可能性がある。

 これらの制限を克服するために、pblat(parallel blat)を提示する。これは、マルチスレッドおよびクラスタコンピューティングのサポートを実装することによって、blatのアライメントを高速化するように機能する。pblatでは、すべてのスレッドがリファレンスゲノム全体とインデックスの同じメモリコピーを共有する。そのため、pblatはblatとほぼ同じ量のメモリを使用する。実行時間は使用されるスレッドの数とともに減少し、pblatの出力結果はblatのそれと同じになる。 pblatのクラスタバージョン(pblat-cluster)はMPI(Message Passing Interface)をサポートするコンピュータクラスタ上で動作する。これにより、blatの実行時間を数日から数分に短縮することができる。

 

f:id:kazumaxneo:20190119101205p:plain

Fig.1Performance evaluation of pblat.  論文より転載。

 

インストール 

mac os10.14でビルドし、動作確認した。

pblat can run on Linux and Mac OS.

本体 Github

git clone https://github.com/icebert/pblat.git
cd  pblat/
make

#conda
mamba install -c bioconda pblat -y

> ./pblat

$ ./pblat 

pblat - BLAT with parallel supports v. 35 fast sequence search command line tool

 

usage:

   pblat database query [-ooc=11.ooc] output.psl

where:

   database and query are each a .fa file

   -ooc=11.ooc tells the program to load over-occurring 11-mers from

               and external file.  This will increase the speed

               by a factor of 40 in many cases, but is not required

   output.psl is where to put the output.

 

options:

   -t=type     Database type.  Type is one of:

                 dna - DNA sequence

                 prot - protein sequence

                 dnax - DNA sequence translated in six frames to protein

               The default is dna

   -q=type     Query type.  Type is one of:

                 dna - DNA sequence

                 rna - RNA sequence

                 prot - protein sequence

                 dnax - DNA sequence translated in six frames to protein

                 rnax - DNA sequence translated in three frames to protein

               The default is dna

   -prot       Synonymous with -t=prot -q=prot

   -ooc=N.ooc  Use overused tile file N.ooc.  N should correspond to 

               the tileSize

   -threads=N  number of threads to run

   -tileSize=N sets the size of match that triggers an alignment.  

               Usually between 8 and 12

               Default is 11 for DNA and 5 for protein.

   -stepSize=N spacing between tiles. Default is tileSize.

   -oneOff=N   If set to 1 this allows one mismatch in tile and still

               triggers an alignments.  Default is 0.

   -minMatch=N sets the number of tile matches.  Usually set from 2 to 4

               Default is 2 for nucleotide, 1 for protein.

   -minScore=N sets minimum score.  This is the matches minus the 

               mismatches minus some sort of gap penalty.  Default is 30

   -minIdentity=N Sets minimum sequence identity (in percent).  Default is

               90 for nucleotide searches, 25 for protein or translated

               protein searches.

   -maxGap=N   sets the size of maximum gap between tiles in a clump.  Usually

               set from 0 to 3.  Default is 2. Only relevent for minMatch > 1.

   -noHead     suppress .psl header (so it's just a tab-separated file)

   -makeOoc=N.ooc Make overused tile file. Target needs to be complete genome.

   -repMatch=N sets the number of repetitions of a tile allowed before

               it is marked as overused.  Typically this is 256 for tileSize

               12, 1024 for tile size 11, 4096 for tile size 10.

               Default is 1024.  Typically only comes into play with makeOoc.

               Also affected by stepSize. When stepSize is halved repMatch is

               doubled to compensate.

   -mask=type  Mask out repeats.  Alignments won't be started in masked region

               but may extend through it in nucleotide searches.  Masked areas

               are ignored entirely in protein or translated searches. Types are

                 lower - mask out lower cased sequence

                 upper - mask out upper cased sequence

                 out   - mask according to database.out RepeatMasker .out file

                 file.out - mask database according to RepeatMasker file.out

   -qMask=type Mask out repeats in query sequence.  Similar to -mask above but

               for query rather than target sequence.

   -repeats=type Type is same as mask types above.  Repeat bases will not be

               masked in any way, but matches in repeat areas will be reported

               separately from matches in other areas in the psl output.

   -minRepDivergence=NN - minimum percent divergence of repeats to allow 

               them to be unmasked.  Default is 15.  Only relevant for 

               masking using RepeatMasker .out files.

   -dots=N     Output dot every N sequences to show program's progress

   -trimT      Trim leading poly-T

   -noTrimA    Don't trim trailing poly-A

   -trimHardA  Remove poly-A tail from qSize as well as alignments in 

               psl output

   -fastMap    Run for fast DNA/DNA remapping - not allowing introns, 

               requiring high %ID. Query sizes must not exceed 5000.

   -out=type   Controls output file format.  Type is one of:

                   psl - Default.  Tab separated format, no sequence

                   pslx - Tab separated format with sequence

                   axt - blastz-associated axt format

                   maf - multiz-associated maf format

                   sim4 - similar to sim4 format

                   wublast - similar to wublast format

                   blast - similar to NCBI blast format

                   blast8- NCBI blast tabular format

                   blast9 - NCBI blast tabular format with comments

   -fine       For high quality mRNAs look harder for small initial and

               terminal exons.  Not recommended for ESTs

   -maxIntron=N  Sets maximum intron size. Default is 750000

   -extendThroughN - Allows extension of alignment through large blocks of N's

 

実行方法

リファレンスFASTAと入力のFASTAファイル database.fasta、query.fassta、outputの順番に指定する(fastqには対応していない)。

pblat -threads=12 database.fa query.fasta output.psl
  • -threads=<N>   number of threads to run

> head output

f:id:kazumaxneo:20190119104004p:plain

 

 

タンパク質配列をクエリとし、blatでターゲットとなるゲノムを検索する(参考)。

(データベースがゲノムで、そこにクエリのproteinをアラインする)

blat -t=dnax -q=prot genome.fasta protein.faa output.psl
  • -maxIntron=N    Sets maximum intron size. Default is 750000
  • -minScore=N    sets minimum score.  This is the matches minus the  mismatches minus some sort of gap penalty.  Default is 30

  • -minIdentity=N    Sets minimum sequence identity (in percent).  Default is 90 for nucleotide searches, 25 for protein or translated protein searches.

  •  -noTrimA    Don't trim trailing poly-A
  • -extendThroughN    Allows extension of alignment through large blocks of N's
  • -out=type    Controls output file format.  Type is one of:

                       psl - Default.  Tab separated format, no sequence

                       pslx - Tab separated format with sequence

                       axt - blastz-associated axt format

                       maf - multiz-associated maf format

                       sim4 - similar to sim4 format

                       wublast - similar to wublast format

                       blast - similar to NCBI blast format

                       blast8- NCBI blast tabular format

                       blast9 - NCBI blast tabular format with comments

 

追記

blast8形式で出力。minimum score90、minimum identity 60%、スレッド数20、NNN部分でアラインメントを止めない、最大イントロン長10000。

pblat -maxIntron=10000 -out=blast8 -noTrimA -extendThroughN -minScore=90 -minIdentity=60 -threads=20 -t=dnax -q=prot genome.fasta closely_related_species_proteome.faa blatout.blast8

 

BLATで近似アラインメントを行なってから、イントロン/エクソン境界をより正確に検出するGenewiseなどを使ってアラインメントをリファインすることがある。

 

samに変換するには以下を参考にしてください。

Question: Conversion Of Blat Output To Sam/Bam

https://www.biostars.org/p/4957/

 

アノテーションに使う場合、こちらが参考になると思います。

Using BLAT to Find Sequence Similarity in Closely Related Genomes

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4101998/(PMC)

引用

pblat: a multithread blat algorithm speeding up aligning sequences to genomes

Meng Wan, Lei Kong

BMC Bioinformatics 2019 20:28