macでインフォマティクス

macでインフォマティクス

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

高効率なプロテインレベルのアセンブリツール PLASS

 2019 6/25 twitter追記

2024/04/03 引用追加

2024/11/02インストール手順修正、help更新

 

 メタゲノム研究の主な制限は、ショートリードの大部分(土壌で80% - 90%[1])を、遺伝子およびタンパク質配列の予測を可能にするのに十分な長さの連続した配列(contigs)にアセンブリすることができないことである。
 低存在量のゲノムはアセンブリが困難であり、アセンブリされていないリードは遺伝的多様性の大きな不均衡と恐らくはさらに大きな割合の生物学的新規性を含み、それはその後の分析でほとんど失われる。
 この損失を減らし、リファレンスゲノムへの依存を少なくするために、遺伝子中心のアプローチが開発されてきた。 1つの環境からの何百ものサンプルがプールされ、コンティグ内の遺伝子が予測され、95%の同一性で遺伝子カタログにまとめられられる[ref.2-5]。各サンプル中の遺伝子存在量は、リードをリファレンス遺伝子クラスターにマッピングすることによって見出される。このようにして、メタゲノムサンプルの機能的および分類学的構成、ならびにそれらの環境パラメータへの依存性を調べることができる。また、ゲノムベースの分析はビニングによって可能になる。
 メタゲノムのショートリードの最先端のアセンブラ[ref.7 - 9]は、de-Bruijnグラフのパスとしてコンティグを見つける。このグラフは、リード内の各k-merワードと、リード内で連続して発生するk-mer間のエッジに対するノードを持っている。メタゲノムのデータでは、de-Bruijnアセンブラは、感度と選択性のトレードオフが制限されている。偽のエッジでグラフが爆発するのを避けるために、k-merは長くspecificでなければならない。しかし、集団内の多様性が高く、重複したリードがSNPによるミスマッチを含むことが多い場合、長いk-merは感度を欠く。どのk値が選択されようとも、k-merは種間で保存されているゲノム領域では十分に特異的であるには短すぎ、集団内の多様性が高い領域では十分に敏感であるには長すぎる。このジレンマはアセンブリを断片化し短くする。

 微生物集団のほとんどのSNPはコードされたタンパク質配列中のアミノ酸変化をもたらさない(補足図1)。したがって、ORFome [ref.10]およびSFA-SPA [ref.11]は、ヌクレオチドの代わりにタンパク質をアセンブリすることを提案した。しかし、それらは大規模なメタゲノムに対して実行するには遅すぎる。

 タンパク質ははるかに少なく短いリピートしか有さない。さらに、非常に類似したタンパク質配列間のキメラアセンブリ(例えば、97%以上の配列同一性)は、ゲノム中でどの遺伝子が一緒に存在するかについての誤った結論をもたらさないため、DNAベースのアセンブリでの問題を回避できる。そのため、タンパク質レベルのアセンブリでは、DNA配列ではキメラアセンブリのリスクを抑え、アセンブリできない領域のアセンブリも可能になることでカバー率が向上する。
 Plassはこれらの利点を活用してアセンブリプロセスを大幅に簡素化する。Plassはグラフフリーのgreedy iterative assembly strategyを使い、all-versus-all のオーバーラップを線形時間で計算する。これにより、単一サーバー上で巨大なリードセットをオーバーラップベースでアセンブリすることが可能になる。最も重要なことは、Plassはショートk-merのみに頼るのではなくフルアラインメントのオーバーラップを計算することによって、de-Bruijnアセンブラの特異性と感度のトレードオフの限界を克服し、複雑なメタゲノムサンプルにおいて数倍以上のタンパク質配列を回収する。
 Plassはランダムディスクアクセスを回避するためにタンパク質配列をメインメモリに格納する必要がある。そのため、入力のリードから変換されたアミノ酸ごとに1バイトのメインメモリが必要で、2〜3億の2×150 bpリードのアセンブリには〜500 GBのRAMが必要になる。これと比較すると、オーバーラップグラフのアセンブラのメモリ要件とランタイムは、処理データと超直線的に比例する。したがって、Plassは、オーバーラップグラフアセンブラの高いspecificity- sensitivityと、de-Bruijnグラフアセンブラの線形のランタイムスケーリング性能および線形のメモリスケーリング性能を組み合わせたものである。(以下略)

 Plass(Protein-Level ASSembler)は、ショートシークエンスリードを6フレームで翻訳し、タンパク質レベルでアセンブリするソフトウェアである。 Plassの主な目的は、複雑なメタゲノムデータセットアセンブリである。 Plassは土壌メタゲノムにてMegahitより10倍多いタンパク質残基をアセンブリする。 Plassは、C ++で実装され、LinuxおよびmacOSで利用可能なGPLライセンスオープンソースソフトウェアである。 ソフトウェアは複数のコア上で効率的に動作するように設計されている。 著者らは2つの冗長性フィルタリングされたリファレンスタンパク質カタログ、 640の土壌試料(SRC)からの20億の配列、および775の海洋真核生物メタトランスクリプトーム(MERC)からの292ミリオンの配列をアセンブリし、最大のタンパク質配列の無料コレクションを構築した。

 



インストール

ubuntu18.04でテストした。

Github

#bioconda 
mamba install -c conda-forge -c bioconda plass -y

#linuxバイナリ
 # static build with AVX2 (fastest)
 wget https://mmseqs.com/plass/plass-linux-avx2.tar.gz; tar xvfz plass-linux-avx2.tar.gz; export PATH=$(pwd)/plass/bin/:$PATH
 # static build with SSE4.1
 wget https://mmseqs.com/plass/plass-linux-sse41.tar.gz; tar xvfz plass-linux-sse41.tar.gz; export PATH=$(pwd)/plass/bin/:$PATH
 # universal build with macOS (Intel or Apple Silicon)
 wget https://mmseqs.com/plass/plass-osx-universal.tar.gz; tar xvfz plass-osx-universal.tar.gz; export PATH=$(pwd)/plass/bin/:$PATH

> plass

$ plass 

Protein Level Assembler.

 

Plass Version: 5.cf8933

© Martin Steinegger (martin.steinegger@mpibpc.mpg.de)

 

usage: plass <command> [<args>]

 

Main workflows for database input/output

  assemble          Assemble protein sequences by iterative greedy overlap assembly

 

> plass assemble

usage: plass assemble <i:fast(a|q)File[.gz]> | <i:fastqFile1_1[.gz] <i:fastqFile1_2[.gz] ... <i:fastqFileN_1[.gz] <i:fastqFileN_2[.gz]> <o:fastaFile> <tmpDir> [options]

 By Martin Steinegger <martin.steinegger@mpibpc.mpg.de> 

options: prefilter:                      

 --alph-size TWIN                 Alphabet size (range 2-21) [13]

 --spaced-kmer-mode INT           0: use consecutive positions in k-mers; 1: use spaced k-mers [0]

 --spaced-kmer-pattern STR        User-specified spaced k-mer pattern

 --mask INT                       Mask sequences in k-mer stage: 0: w/o low complexity masking, 1: with low complexity masking [0]

 --mask-lower-case INT            Lowercase letters will be excluded from k-mer search 0: include region, 1: exclude region [0]

 -k INT                           k-mer length (0: automatically set to optimum) [14]

 --split-memory-limit BYTE        Set max memory per split. E.g. 800B, 5K, 10M, 1G. Default (0) to all available system memory [0]

 --add-self-matches BOOL          Artificially add entries of queries with themselves (for clustering) [0]

align:                          

 --min-seq-id FLOAT               Overlap sequence identity threshold [0.0, 1.0] [0.900]

 --cov-mode INT                   0: coverage of query and target

                                  1: coverage of target

                                  2: coverage of query

                                  3: target seq. length has to be at least x% of query length

                                  4: query seq. length has to be at least x% of target length

                                  5: short seq. needs to be at least x% of the other seq. length [0]

 -c FLOAT                         List matches above this fraction of aligned (covered) residues (see --cov-mode) [0.000]

 --wrapped-scoring BOOL           Double the (nucleotide) query sequence during the scoring process to allow wrapped diagonal scoring around end and start [0]

 -e DOUBLE                        Extend sequences if the E-value is below [0.0, inf] [1.000E-05]

 -a BOOL                          Add backtrace string (convert to alignments with mmseqs convertalis module) [0]

 --min-aln-len INT                Minimum alignment length (range 0-INT_MAX) [0]

 --seq-id-mode INT                0: alignment length 1: shorter, 2: longer sequence [0]

kmermatcher:                    

 --kmer-per-seq INT               k-mers per sequence [60]

 --kmer-per-seq-scale TWIN        Scale k-mer per sequence based on sequence length as kmer-per-seq val + scale x seqlen [nucl:0.200,aa:0.000]

 --adjust-kmer-len BOOL           Adjust k-mer length based on specificity (only for nucleotides) [0]

 --hash-shift INT                 Shift k-mer hash initialization [67]

 --include-only-extendable BOOL   Include only extendable [1]

 --ignore-multi-kmer BOOL         Skip k-mers occurring multiple times (>=2) [1]

profile:                        

 --num-iterations INT             Number of assembly iterations [1, inf] [12]

misc:                           

 --dbtype INT                     Database type 0: auto, 1: amino acid 2: nucleotides [0]

 --shuffle BOOL                   Shuffle input database [1]

 --createdb-mode INT              Createdb mode 0: copy data, 1: soft link data and write new index (works only with single line fasta/q) [0]

 --id-offset INT                  Numeric ids in index file are offset by this value [0]

 --rescore-mode INT               Rescore diagonals with:

                                  0: Hamming distance

                                  1: local alignment (score only)

                                  2: local alignment

                                  3: global alignment

                                  4: longest alignment fulfilling window quality criterion [3]

 --min-length INT                 Minimum codon number in open reading frames [45]

 --max-length INT                 Maximum codon number in open reading frames [32734]

 --max-gaps INT                   Maximum number of codons with gaps or unknown residues before an open reading frame is rejected [2147483647]

 --contig-start-mode INT          Contig start can be 0: incomplete, 1: complete, 2: both [2]

 --contig-end-mode INT            Contig end can be 0: incomplete, 1: complete, 2: both [2]

 --orf-start-mode INT             Orf fragment can be 0: from start to stop, 1: from any to stop, 2: from last encountered start to stop (no start in the middle) [1]

 --forward-frames STR             Comma-separated list of frames on the forward strand to be extracted [1,2,3]

 --reverse-frames STR             Comma-separated list of frames on the reverse strand to be extracted [1,2,3]

 --translation-table INT          1) CANONICAL, 2) VERT_MITOCHONDRIAL, 3) YEAST_MITOCHONDRIAL, 4) MOLD_MITOCHONDRIAL, 5) INVERT_MITOCHONDRIAL, 6) CILIATE

                                  9) FLATWORM_MITOCHONDRIAL, 10) EUPLOTID, 11) PROKARYOTE, 12) ALT_YEAST, 13) ASCIDIAN_MITOCHONDRIAL, 14) ALT_FLATWORM_MITOCHONDRIAL

                                  15) BLEPHARISMA, 16) CHLOROPHYCEAN_MITOCHONDRIAL, 21) TREMATODE_MITOCHONDRIAL, 22) SCENEDESMUS_MITOCHONDRIAL

                                  23) THRAUSTOCHYTRIUM_MITOCHONDRIAL, 24) PTEROBRANCHIA_MITOCHONDRIAL, 25) GRACILIBACTERIA, 26) PACHYSOLEN, 27) KARYORELICT, 28) CONDYLOSTOMA

                                   29) MESODINIUM, 30) PERTRICH, 31) BLASTOCRITHIDIA [1]

 --translate INT                  Translate ORF to amino acid [0]

 --use-all-table-starts BOOL      Use all alternatives for a start codon in the genetic table, if false - only ATG (AUG) [0]

 --keep-target BOOL               Keep target sequences [1]

 --protein-filter-threshold FLOAT filter proteins lower than threshold [0.0,1.0] [0.200]

 --filter-proteins INT            filter proteins by a neural network [0,1] [1]

common:                         

 --compressed INT                 Write compressed output [0]

 -v INT                           Verbosity level: 0: quiet, 1: +errors, 2: +warnings, 3: +info [3]

 --sub-mat TWIN                   Substitution matrix file [nucl:nucleotide.out,aa:blosum62.out]

 --max-seq-len INT                Maximum sequence length [65535]

 --threads INT                    Number of CPU-cores used (all by default) [128]

 --db-load-mode INT               Database preload mode 0: auto, 1: fread, 2: mmap, 3: mmap+touch [0]

 --delete-tmp-inc INT             Delete temporary files incremental [0,1] [1]

 --remove-tmp-files BOOL          Delete temporary files [0]

 --mpi-runner STR                 Use MPI on compute cluster with this MPI command (e.g. "mpirun -np 42")

expert:                         

 --write-lookup INT               write .lookup file containing mapping from internal id, fasta id and file number [1]

 --filter-hits BOOL               Filter hits by seq.id. and coverage [0]

 --sort-results INT               Sort results: 0: no sorting, 1: sort by E-value (Alignment) or seq.id. (Hamming) [0]

 --create-lookup INT              Create database lookup file (can be very large) [0]

 

references:

 - Steinegger M, Mirdita M, Soding J: Protein-level assembly increases protein sequence recovery from metagenomic samples manyfold. Nature Methods, 16(7), 603-606 (2019)

> plass nuclassemble #最近のバージョンには存在しない

$ plass nuclassemble

plass nuclassemble:

Extends sequence to the left and right using ungapped alignments.

 

Please cite:

Steinegger, M. Mirdita, M., & Soding, J. Protein-level assembly increases protein sequence recovery from metagenomic samples manyfold. biorxiv, https://doi.org/10.1101/386110 (2018)

© Martin Steinegger <martin.steinegger@mpibpc.mpg.de> 

 

Usage: <i:fast(a|q)File[.gz]> | <i:fastqFile1_1[.gz] ... <i:fastqFileN_1[.gz] <i:fastqFile1_2[.gz] ... <i:fastqFileN_2[.gz]> <o:fastaFile> <tmpDir> [options]

 

align options              default   description [value range]

  -e                       0.000     Extend sequences if the E-value is below [0.0, inf]         

  --min-seq-id             0.900     Overlap sequence identity threshold [0.0, 1.0]              

 

profile options            default   description [value range]

  --num-iterations         12        Number of assembly iterations [1, inf]                      

 

common options             default   description [value range]

  --threads                56        number of cores used for the computation (uses all cores by default)

  -v                       3         verbosity level: 0=nothing, 1: +errors, 2: +warnings, 3: +info

 

An extended list of options can be obtained by calling 'plass nuclassemble -h'.

2 Database paths are required

dockerイメージも用意されている。

 

テストラン

git clone https://github.com/soedinglab/plass.git
cd plass/examples/

 ペアエンドfastq。ラン時は作業ディレクトリを指定する。

plass assemble reads_1.fastq.gz reads_2.fastq.gz assembly.fa tmp

 シングルエンドfastq

 plass assemble reads_1.fastq.gz assembly.fa tmp

 

感想

HIseqやNovaseqクラスの複数のメタゲノムデータセットアセンブリには、準スパコンクラスのマシンとスケーリング性能が極めて高いde novoアセンブラの両方が必要なのかと考えてましたが(Extreme Scale De Novo Metagenome Assemblyの"C. Performance Results"を参照)、Preprintを読むと、PLASSは、メモリ効率の良いmegahitがsegmentation errorを起こすようなデータセットでも(プロテインの)アセンブリが可能なようです。アセンブリされたプロテインのトータルサイズも、プロテインの結果はMEGAHITの結果よりずっと長く、特に複雑性の高いデータセットで差が顕著になっています。プロテインデータベースを構築してスクリーニングする用途にも使えそうです。

引用

Protein-level assembly increases protein sequence recovery from
metagenomic samples manyfold
Martin Steinegger, Milot Mirdita, Johannes Söding

bioRxiv preprint first posted online Aug. 7, 2018

 

追記

Protein-level assembly increases protein sequence recovery from metagenomic samples manyfold
Martin Steinegger, Milot Mirdita & Johannes Söding 
Nature Methods volume 16, pages603–606 (2019)

 

Strain-resolved de-novo metagenomic assembly of viral genomes and microbial 16S rRNAs
Annika Jochheim,  Florian A. Jochheim, Alexandra Kolodyazhnaya, Étienne Morice, Martin Steinegger,  Johannes Söding

bioRxiv, Posted March 29, 2024.

 

関連