macでインフォマティクス

macでインフォマティクス

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

SPAdesアセンブラ

2018 タイトル修正、tips追加

2019 tweet追加、誤字修正ンストールバージョン3.13.1に更新、spades3.14に修正

2020 インストール追記、リンク追加、新しい論文引用、ツイート追記

2021 1/14 3.15にhelpを更新、ツイート追記、6/10 更新

2022/08/08 インストールバージョン更新

2023/01/31 追記、04/23 誤字修正

2024/02/03 追記

 

 ref.1

 人体から海洋までほとんどの環境のバクテリアは研究所でクローン化できないため、既存のNGS(Next Generation Sequencing)技術を使用してシーケンスを決定することはできない。これは、Human Microbiome Project(HMP)(論文より Gill et al、2006)から抗生物質発見(Li and Vederas、2009)までのさまざまなプロジェクトの主要なボトルネックを表している。例えば、HMPの重要な問題は、バクテリアが互いにどのように相互作用するかである。これらの相互作用は、しばしば、他のバクテリアとのコミュニケーションまたはそれらを殺すために産生される様々なペプチドによって行われる。しかしながら、質量分析(このような研究のための重要な技術)は、かなり完全なプロテオームの知識を必要とするため、ヒトマイクロバイオームのpeptidomics研究は今や限られている。一方、新しいペプチド抗生物質の研究は、Non-Ribosomal Peptide Syntetases(NRPS)をコードする遺伝子のDNA配列決定(Sieber and Marahiel、2005)により大きく恩恵を受けるが、既存のメタゲノミクスアプローチは、通常このような例外的に長い遺伝子(60,000ヌクレオチド)をアセンブリすることができない。

 HMPと新しい抗生物質の発見は、Single-Cell Sequencing(SCS)によって革命を起こす多くのプロジェクトの2つの例にすぎない。 SCSの実験的(Ishoey et al、2008; Navin et al、2011; Islam et al、2011)および計算(Chitsaz et al、2011)の両局面における最近の改良により、1細胞からバクテリアゲノムを決定する可能性が開かれた。特に、Chitsaz et al(2011)は、SCSで生物の代謝を推定するのに十分な多数の遺伝子を捕捉できることを示した。多くのアプリケーション(プロテオミクスや抗生物質の発見を含む)では、大多数の遺伝子を捕捉することは、完全なゲノムと同じくらい有用である。

 現在、MDA(Multiple Displacement Amplification)はシングルセル増幅の主要技術である(Dean et al、2001)。しかし、MDAは極端な増幅バイアス(異なる領域間でカバレッジの差が桁違いに大きい)をもたらし、キメラリードとキメラペアリードが生じ、これが次のアセンブリを複雑にする。既存のアセンブラはこれらの問題に対処できないことから、Rodrigue et al. (2009) は、SCSが直面している課題は実験的ではなく計算的なものになっていると指摘した。最近の論文(Marcy et al、2007; Woyke et al、2010; Youssef et al、2011; Blainey et al、2011; Grindberg et al、2011)は、SCSにおけるフラグメントアセンブリの課題を示している。

 Chitsaz et al (2011)はE+V-SC assemblerを導入し、modified EULER-SR とmodified velvet を組み合わせ、SCSデータのフラグメントアセンブリの大幅な改善を達成した。しかし、Chitsaz et alの共著者である我々(この論文の著者ら)(G.T.およびP.A.P.)は、 (2011)では、SCSの可能性を最大限に活用するために、(既存のツールを変更するだけでなく)アルゴリズム設計を変更する必要があることに気付いた。

 著者らは、SPAdesアセンブラを発表し、SCSと標準(マルチセル)バクテリアの両方のデータセットのために、数多くの新しいアルゴリズムソリューションを導入し、最先端のアセンブラを改良した。フラグメントアセンブリは、そのk-mersのセットから文字列を再構成する問題として、しばしば抽象化される。この抽象化は、多くのフラグメントアセンブリアルゴリズムの基礎となるアセンブリへのde Bruijnアプローチに自然につながる。しかしながら、NGSデータのより高度な抽象化は、文字列中の距離≒dのk-mersのペア(以下、k-bimersと呼ぶ)から文字列を再構築する問題を考慮する。残念ながら、前者の抽象化には単純なアルゴリズムがある一方で、後者の抽象化の分析は、主にde Bruijn graphの後処理ヒューリスティックになっている。ペアリード(以下、bireads2と呼ぶ)の解析のための多くのヒューリスティックが提案されているが(Pevzner et al、2001; Pevzner and Tang、2001; Zerbino and Birney、2008; Butler et al、2008; Simpson et al、2009 ; Chaisson et al、2009; Li et al、2010)、bireads2の利用はアセンブリのステージの中でおそらく最も不十分な段階であると考えられる。 Medvedev et al(2011a)は、後者の抽象化に適した新しい手法であり、標準的なDe Bruijnグラフより重要な利点を持つPaired de Bruijnグラフ(PDBG)を導入した。しかし、PDBGは理論的には実用的ではなくむしろ理論的に導入されており、主にリード間を固定距離とした非現実的なケースを対象としていた。

 IduryとWaterman(1995)がフラグメントアセンブリのためにde Bruijn graphを導入したとき、多くの人々はサンガーのリードのエラー率が高いため、このアプローチを非実用的と見なしていた(ref.3)。Pevzner et al (2001)は、このボトルネックを解消するためにエラー訂正手順を導入した。同様に、PDBGsは、NGSに特徴的なbiread(and thus k-bimer)距離の変動のために実用的でないように見えるかもしれない。 SPAdesは、k-bimersの正確な距離を明らかにするk-bimer調整を導入し、PDBGに触発されたペアのアセンブリグラフを導入することによって、このボトルネックに対処する。特に、SPAdesはペアリードを利用することができる。 E + V-SCは、リードを使用したが、キメラリードペアのレベルの上昇によって生じるミスアセンブリを回避するために、ペアエンドのpairingを無視した。

 読者はPevzner et al(2004 pubmed)で紹介されたA-Bruijn graphの概念に精通していると仮定する 。この論文のDe Bruijn graph、PDBG、およびその他のgraphは、A-Bruijn graphの特別なケースである。 SPAdesは、最初のde Bruijnグラフを作成するためだけにk-mersを使用し、後でそれらについて「忘れる」という意味で普遍的なA-Bruijnアセンブラである。それ以降のステージでは、k-mersによってラベル付けする必要がないgraphに対してのみgraph理論操作を実行する。操作はgraphトポロジーカバレッジ、およびシーケンス長に基づいているが、シーケンスそのものには基づいていない。最後の段階で、コンセンサスDNA配列が復元される。著者らは、同じフレームワーク内でA-Bruijn graphのいくつかのバリエーション(例えば、Bruijn graphの対と多重化)を実装し、これらのgraphが有用であることが証明されている他のアプリケーションに適用する汎用アセンブラを設計した(Bandeira et al、 2007、2008; Pham and Pevzner、2010)。

 セクション2では、SPAdesのステージの概要を説明する。セクション3では、de Bruijn graph、Multisized de Bruijn graph、およびPDBGを定義する。セクション4-6は、SPAdesのさまざまな段階をカバーしている。セクション7では、シングルセルおよび培養された大腸菌(E.coli)データセットについて、SPAdesおよび他のアセンブラベンチマークする。セクション8では、アセンブリgraphの構築とデータ構造に関する追加の詳細を示す。第9章では、ペアアセンブリgraphを構築する詳細な例を示す。セクション10では、ユニバーサルアセンブラの概念についてさらに説明する。

セクション2

以下では、SCSで特に面倒な問題を扱うSPAdesの4つの段階について説明する。

(1)ステージ1(アセンブリgraph構築)はすべてのNGSアセンブラによって扱われ、しばしばde Bruijn graph簡略化(EULER / Velvetにおけるバルジ/バブル除去)と呼ばれる。この論文では、複数のk-merを使ったde Bruijn graphを使用し、新しいバルジ/チップ除去アルゴリズムを実装し、キメラを検出して除去し、biread情報を距離ヒストグラムに集約し、実行されたgraph操作をバックトラックすることを可能にするアセンブリグラフ構築の新しいアプローチを提案する。
(2)ステージ2(k-bimer調整)は、アセンブリgraphの距離ヒストグラムとパスのジョイント解析を使用して、ゲノム内のk-mer間の正確な距離推定値(アセンブリgraph内のエッジ)を導き出す。
(3)ステージ3は、PDBGアプローチからインスパイアされたペアのアセンブリgraphを構築する。
(4)ステージ4(コンティグ構築)は、サンガーシーケンス時代に十分に研究された(Ewing et al、1998)。 NGSプロジェクトは通常、高いカバレッジを特徴とするため、NGSアセンブラはかなり正確なコンティグを生成する(ただしSCSの精度は低下する)。 SPAdesは、コンティグのDNA配列を構築し、グラフの簡略化(論文 8.6節参照)を目的としてコンティグへリードをマッピングする。
これまでの研究では、さまざまなアセンブラをエラー訂正ツールと組み合わせることで性能が向上することが実証されている(Pevzner et al、2010; Kelley et al、2010; Ilie et al、2010; Gnerre et al、2011)。しかし、 (Quake [Kelley et al、2010]などの)ツールは、暗黙的にほぼ均一なカバレッジを仮定しているため、シングルセルのデータでは機能しない。 Chitsaz et al。(2011)はVelvet-SCとEULERのエラー訂正(Chaisson and Pevzner、2008)を結びつけ、E + V-SCツールとなった。この論文では、SPAdesは、アセンブリ前にエラー修正と品質トリミングのためにHammer(Medvedev et al、2011b)(SCSを対象としています)によるエラー修正を行う。

 

オリジナルのSPAdesは以下のようなフローでアセンブリを実行する。

f:id:kazumaxneo:20180614103714j:plain

プレゼンテーションPDFより転載。 de Brujin graph構築にはショートリード情報が用いられ、ロングリード情報はエラーの多寡に関わらず使われない。

 

 

公式ページ

http://cab.spbu.ru/software/spades/

3.12マニュアル

http://cab.spbu.ru/files/release3.12.0/manual.html 

 

New Frontiers of Genome Assembly with SPAdes 3.1


ppt1

f:id:kazumaxneo:20180613212251j:plain

ppt2

f:id:kazumaxneo:20180613212329j:plain

 

インストール

mac os 10.14でコンパイル済みのプログラムをダウンロードしてテストした。

Github

 バイナリをダウンロードする。Linuxは上のマニュアル参照(バイナリが用意されている)。

#3.15.5 linux binary
wget http://cab.spbu.ru/files/release3.15.5/SPAdes-3.15.5-Linux.tar.gz
tar -xzf SPAdes-3.15.5-Linux.tar.gz
cd SPAdes-3.15.5-Linux/bin/
export PATH=$PATH:$PWD

#darwin(mac) binary
curl http://cab.spbu.ru/files/release3.15.5/SPAdes-3.15.5-Darwin.tar.gz -o SPAdes-3.15.5-Darwin.tar.gz
tar -zxf SPAdes-3.15.5-Darwin.tar.gz
cd SPAdes-3.15.5-Darwin/bin/
export PATH=$PATH:$PWD

#bioconda(link)
conda install -c bioconda -y sppades

#brew
brew install spades

#help (v3.15.0)

./spades.py 

# ./spades.py 

SPAdes genome assembler v3.15.0

 

Usage: spades.py [options] -o <output_dir>

 

Basic options:

  -o <output_dir>             directory to store all the resulting files (required)

  --isolate                   this flag is highly recommended for high-coverage isolate and multi-cell data

  --sc                        this flag is required for MDA (single-cell) data

  --meta                      this flag is required for metagenomic data

  --bio                       this flag is required for biosyntheticSPAdes mode

  --rna                       this flag is required for RNA-Seq data

  --plasmid                   runs plasmidSPAdes pipeline for plasmid detection

  --metaviral                 runs metaviralSPAdes pipeline for virus detection

  --metaplasmid               runs metaplasmidSPAdes pipeline for plasmid detection in metagenomic datasets (equivalent for --meta --plasmid)

  --rnaviral                  this flag enables virus assembly module from RNA-Seq data

  --iontorrent                this flag is required for IonTorrent data

  --test                      runs SPAdes on toy dataset

  -h, --help                  prints this usage message

  -v, --version               prints version

 

Input data:

  --12 <filename>             file with interlaced forward and reverse paired-end reads

  -1 <filename>               file with forward paired-end reads

  -2 <filename>               file with reverse paired-end reads

  -s <filename>               file with unpaired reads

  --merged <filename>         file with merged forward and reverse paired-end reads

  --pe-12 <#> <filename>      file with interlaced reads for paired-end library number <#>.

                              Older deprecated syntax is -pe<#>-12 <filename>

  --pe-1 <#> <filename>       file with forward reads for paired-end library number <#>.

                              Older deprecated syntax is -pe<#>-1 <filename>

  --pe-2 <#> <filename>       file with reverse reads for paired-end library number <#>.

                              Older deprecated syntax is -pe<#>-2 <filename>

  --pe-s <#> <filename>       file with unpaired reads for paired-end library number <#>.

                              Older deprecated syntax is -pe<#>-s <filename>

  --pe-m <#> <filename>       file with merged reads for paired-end library number <#>.

                              Older deprecated syntax is -pe<#>-m <filename>

  --pe-or <#> <or>            orientation of reads for paired-end library number <#> 

                              (<or> = fr, rf, ff).

                              Older deprecated syntax is -pe<#>-<or>

  --s <#> <filename>          file with unpaired reads for single reads library number <#>.

                              Older deprecated syntax is --s<#> <filename>

  --mp-12 <#> <filename>      file with interlaced reads for mate-pair library number <#>.

                              Older deprecated syntax is -mp<#>-12 <filename>

  --mp-1 <#> <filename>       file with forward reads for mate-pair library number <#>.

                              Older deprecated syntax is -mp<#>-1 <filename>

  --mp-2 <#> <filename>       file with reverse reads for mate-pair library number <#>.

                              Older deprecated syntax is -mp<#>-2 <filename>

  --mp-s <#> <filename>       file with unpaired reads for mate-pair library number <#>.

                              Older deprecated syntax is -mp<#>-s <filename>

  --mp-or <#> <or>            orientation of reads for mate-pair library number <#> 

                              (<or> = fr, rf, ff).

                              Older deprecated syntax is -mp<#>-<or>

  --hqmp-12 <#> <filename>    file with interlaced reads for high-quality mate-pair library number <#>.

                              Older deprecated syntax is -hqmp<#>-12 <filename>

  --hqmp-1 <#> <filename>     file with forward reads for high-quality mate-pair library number <#>.

                              Older deprecated syntax is -hqmp<#>-1 <filename>

  --hqmp-2 <#> <filename>     file with reverse reads for high-quality mate-pair library number <#>.

                              Older deprecated syntax is -hqmp<#>-2 <filename>

  --hqmp-s <#> <filename>     file with unpaired reads for high-quality mate-pair library number <#>.

                              Older deprecated syntax is -hqmp<#>-s <filename>

  --hqmp-or <#> <or>          orientation of reads for high-quality mate-pair library number <#> 

                              (<or> = fr, rf, ff).

                              Older deprecated syntax is -hqmp<#>-<or>

  --sanger <filename>         file with Sanger reads

  --pacbio <filename>         file with PacBio reads

  --nanopore <filename>       file with Nanopore reads

  --trusted-contigs <filename>

                              file with trusted contigs

  --untrusted-contigs <filename>

                              file with untrusted contigs

 

Pipeline options:

  --only-error-correction     runs only read error correction (without assembling)

  --only-assembler            runs only assembling (without read error correction)

  --careful                   tries to reduce number of mismatches and short indels

  --checkpoints <last or all>

                              save intermediate check-points ('last', 'all')

  --continue                  continue run from the last available check-point (only -o should be specified)

  --restart-from <cp>         restart run with updated options and from the specified check-point

                              ('ec', 'as', 'k<int>', 'mc', 'last')

  --disable-gzip-output       forces error correction not to compress the corrected reads

  --disable-rr                disables repeat resolution stage of assembling

 

Advanced options:

  --dataset <filename>        file with dataset description in YAML format

  -t <int>, --threads <int>   number of threads. [default: 16]

  -m <int>, --memory <int>    RAM limit for SPAdes in Gb (terminates if exceeded). [default: 250]

  --tmp-dir <dirname>         directory for temporary files. [default: <output_dir>/tmp]

  -k <int> [<int> ...]        list of k-mer sizes (must be odd and less than 128)

                              [default: 'auto']

  --cov-cutoff <float>        coverage cutoff value (a positive float number, or 'auto', or 'off')

                              [default: 'off']

  --phred-offset <33 or 64>   PHRED quality offset in the input reads (33 or 64),

                              [default: auto-detect]

  --custom-hmms <dirname>     directory with custom hmms that replace default ones,

                              [default: None]

 

ほかにmetaspades.py、rnaspades.py、dipspades.py、plasmidspades.py、/truspades.pyがある。以下のような使い分けになる。

  • metaSPAdes – a pipeline for metagenomic data sets (see metaSPAdes options).
  • plasmidSPAdes – a pipeline for extracting and assembling plasmids from WGS data sets (see plasmidSPAdes options).
  • rnaSPAdes – a de novo transcriptome assembler from RNA-Seq data (see rnaSPAdes manual).
  • truSPAdes – a module for TruSeq barcode assembly (see truSPAdes manual).
  • dipSPAdes – a module for assembling highly polymorphic diploid genomes (deprecated, see dipSPAdes manual).

SPAdesは現在でも開発が続けられており、バージョンを変えると末端部位やカバレッジの低い領域の結果が揺れ動くことがあります。アップグレードする時は注意してください。

 

 ラン

ペアエンドリードfastqのアセンブリ。fastqはgzip圧縮にも非圧縮にも両方対応している。

spades.py -1 pair1.fq.gz -2 pair2.fq.gz -k auto --careful -t 12 -o output
  • -1     file with forward paired-end reads-1 <filename> file with forward paired-end reads
  • -2     file with reverse paired-end reads
  • -t     number of threads-t/--threads <int> number of threads [default: 16]
  • -k     comma-separated list of k-mer sizes (must be odd and-k <int,int,...> comma-separated list of k-mer sizes (must be odd and less than 128) [default: 'auto']
  • --careful     tries to reduce number of mismatches and short indels--careful tries to reduce number of mismatches and short indels
  • -o   directory to store all the resulting files (required)-o <output_dir> directory to store all the resulting files (required)

シングルエンドfastqは-s、インターレースfastqは--12で指定する。シングルセルシーケンスなどで複数fastqがある場合は、数値をつけて指定する。例えばインターレースfastqを複数指定するなら"--pe1-12 sample1.fq --pe2-12 sample2.fq --pe3-12 sample3.fq ..."。詳細はマニュアルの Specifying multiple librariesの節を参照。

 

ペアエンドリードfastqのアセンブリにマージしたfastqも加える。

spades.py -k 21,33,55,77,99,127 --careful -t 16 \
-1 pair1.fq -2 pair2.fq --merged merged.fq

 

ペアエンドリードfastqとONTのロングリードのハイブリッドアセンブリ

spades.py -1 pair1.fq.gz -2 pair2.fq.gz --nanopore long.fa -k auto --careful -t 12 -o output

 

ランが終わるとk-merごとにアセンブリデータフォルダができる。scaffolds.fastaが最終出力となる。correctedにはエラー訂正されたfq.gzも保存されている。

f:id:kazumaxneo:20180613230424j:plain

 

ペアエンドリードfastqとメイアペア(ロングインサートのペアエンド)のアセンブリ。ショートリードは2ライブラリ、メイトペアもインサートサイズの異なる2ライブラリ使うとする。

spades.py -k 21,33,55,77,99,127 --careful -t 16 \
--pe1-1 pair_sampleA_1.fq.gz --pe1-2 pair_sampleA_2.fq.gz \
--pe2-1 pair_sampleB_1.fq.gz --pe2-2 pair_sampleB_2.fq.gz \
--mp1-1 mate_sampleC1.fq.gz --mp1-2 sampleC2.fq.gz \
--mp2-1 mate_sampleD1.fq.gz --mp2-2 sampleD2.fq.gz \
-o spades_outdir


#mate-pairサンプル1はfrの向きのリードペア。"--mp-or <サンプル番号> リードの向き"、順に指定
spades.py -k 21,33,55,77,99,127 --careful -t 16 \
--pe1-1 pair_sampleA_1.fq.gz --pe1-2 pair_sampleA_2.fq.gz \
--pe2-1 pair_sampleB_1.fq.gz --pe2-2 pair_sampleB_2.fq.gz \
--mp1-1 mate_sampleC1.fq.gz --mp1-2 sampleC2.fq.gz \
--mp-or 1 fr \
-o spades_outdir
  • --mp-or <#> <orientation>.  orientation of reads for mate-pair library number <#>  (<orientation> = fr, rf, ff). 

     

 ケースバイケースですが、うまくパラメータを使えばパフォーマンスが上がる可能性があります。気づいたことを書いておきます。

  1. --carefulをつけると、付属のbwaを使ってアセンブリ後にミスマッチを修正するステップが挿入される。フラグを立てると、マッピングと修正の時間分だけランタイムは増えるが、スモールゲノムではフラグをつけることが推奨されている(defaultではオフ。エラーが増える可能性もゼロではないので要検討)。ラージゲノムでは非推奨。純粋な培養株のシーケンスを想定して修復するため、rnaSPADdesとmetaSPAdesにこのオプションはない。
  2. たとえハイクオリティなPacbioのロングリードでも、ロングリードだけでSPAdesを動かすことは非推奨になってます。
  3. --only-error-correctionをつけると、BayesHammerによるエラー訂正だけが実行されます。エラー訂正されたfastqが欲しいときに使えます。
  4. --only-assemblerをつけると、エラー訂正なしでアセンブリが実行されます。非推奨ですが、ほかのエラー訂正ツールで修復したfastqを使う時にも利用できます。
  5. 真菌くらい大きなゲノムになるとかなりメモリを食います。論文を読んで、ランできそうか確認してください。defaultでは-m --memoryは最大250GBで、越えるとランを停止します(メモリ最大値を増やしてコンパイルする)。
  6. 軽度のコンタミネーションは、--cov-cutoff <NUM>でカバレッジcut-offをつけることで抑制できます。アセンブリ後にカバレッジGC、taxonomy profileを分析して、異種生物ゲノムのcontigのコンタミカバレッジGCなどの配列特徴量によって区別できそうか調べることがよく行われます。
  7. unicyclerのようにk-merを増やすとパフォーマンスが上がることがある。unicyclerのdefault設定では、Miseqの300bpリードには27,47,63,77,89,99,107,115,121,127が使用される。
  8. 3.12で導入された--mergedのオプションを使い、オーバーラップのペアエンドをマージしたfastqを指定すると、少しだけパフォーマンスが上がることがある(リード長以上の大きなk値が使用可能になるメリットもある)。
  9. 3.14では"--isolate"の オプションが追加された。単離された菌のシーケンスデータで100前後のシークエンシングリードデプスがある場合はつけることが推奨されている(correct/がなくなり他のディレクトリ構造も大きく変わるので注意)*1。またbiosyntheticSPAdesの"--bio"オプション(link)、Hybrid transcriptome assemblyモードなどが追加された。
  10. "--isolate"オプションの実装は、最近のNGSのデータが極めてハイクオリティであり、十分に深くシークエンシングされている限り、単離された菌のアセンブリではエラーコレクションによる改善はほぼないことが背景にあるのだと推測される。エラーコレクションがないことで、SPAdesのランタイムは2倍かそれ以上に短くなる。
  11. SPADesは複数のk-merを使うIterative なDe Bruijn graph Assemblerで、小さなkmerサイズから始めて、得られたコンティグを仮想的なリードとして使用し、より長いkmerサイズで次のグラフを構築することで、反復的にグラフを構築する。小さいk-merはカバレッジが不足する領域やシークエンシングエラーが多い領域から配列を得るために重要になる。中程度のk-merは、小さいk-merで単純化されなかった領域を伸ばすために重要となる、大きいk-merは、中程度のk-merで単純化されなかった領域を伸ばすために重要となる。よって、短いk-merと大きいk-merだけ指定するより、小・中・大と適度な間隔でk値を指定した方がcontigの連続性は伸びる。
  12. 過剰なリードデプスについて。公開されているゲノム配列を調べる研究をしている際に、必須遺伝子の途中で切れていることが多い奇妙な細菌ゲノム配列を何個か見つけた。このようなゲノムのfastqをダウンロードすると、リード数が非常に多いという特徴があった(例えば14Gbなど)。リード数を100x程度に減らして自分の手でアセンブルし直すと必須遺伝子途中で切れる現象は見られなかった。これは何を示唆しているか?;SPAdesゲノムアセンブラに限らず、ウィルス以外のゲノムアセンブラの多くは適度なリードデプスを想定している。例えばバルクで読んでいるために、1%の割合でしか浸透していない変異もランダムにシークエンシングされているサンプルを仮定する。このサンプルに2000xのデプスがある時、この頻度の低い突然変異配列は平均して20xのk-mer配列として取り出される。これをグラフに使うと、アセンブリグラフをトラバースする時にアセンブルの分岐となり、結果として切れてしまう。おそらくこのような理由でコンティグの連続性が低くなる。ゲノムアセンブリに必要十分な、そして過剰でないリードデプスになっているか注意すること。
  13. 倍数体について。SPADesファミリーのアセンブラは、dipSPADes(開発終了)以外はハプロイドゲノムのアセンブラとなっています。小さな真核生物ゲノムアセンブリに使う時は、倍数体ゲノムに使っていないか注意してください。では倍数体ゲノムをdipSPADesでないSPAdesでアセンブルした場合、どのような結果になるのでしょうか?;相同染色体間のヘテロ接合度は染色体や領域によって割合に差があります。ヘテロ接合度の高さによって、その染色体領域が相同染色体間で区別してアセンブルされるかコンセンサスとして1つの配列にアセンブルされるか変わると言えます。どのコンティグが倍数体染色体配列間のコンセンサス配列になっていて、どのコンティグが区別されてアセンブルされているかは、シークエンシングリードのエラー率、ゲノムのヘテロ接合度やゲノムの複雑さに影響を受ける複雑な問題であり、単純な問題ではありません。よってこのような問題を想定していないSPAdesを使って倍数体ゲノムをアセンブリするとパフォーマンスが出ない可能性が高いといえます。既にSPAdesでアセンブルしているなら、この論文で使われているヘテロなバリアント頻度を推定する方法でコンティグを評価することを考えてください(例;マッピング => バリアントコール=> ヘテロなバリアントの単位長あたりの密度を計算)。このあたりの話題については、植物ならこちらの論文が参考になります(7. From haploid to diploid genome assemblyのセクション)。
  14. 上の13と関連して、現在所有しているドラフトゲノムアセンブリ配列について、相同染色体間の配列が区別されて組み立てられているのか、あるいはそうでないのか初期評価するなら、BUSCOの"duplicated busco"の値の高さを確認するのが有効です。2倍体の全ての相同染色体間の配列が区別されてアセンブルされているなら"100%"に(シングルコピーのBUSCO遺伝子が2コピーずつあるので)、逆に全ての領域が相同染色体間の配列のコンセンサスとして(=1個の重複がない配列として)表現されているなら"0%"に、混ざっているなら0~100%の間の値を取るはずです。ただし誤差はそれなりにあることに注意してください。
  15. 細菌ゲノムではSPAdesはできることを全部しているハイパフォーマンスなアセンブラですが、ある程度大きなゲノム(藻類や真菌など)のゲノムアセンブリに使用する場合、SPAdesではあまりパフォーマンスが出ない可能性があります。そのようなゲノムで且つショートリードしか利用できない場合、個人的にはMaSuRCaをお勧めします。昔は導入を含めて使いにくい部分がありましたが、継続的な開発により簡単に使えるようになり、性能もアップしています。ある程度複雑なゲノム配列のアセンブリででは、ショートリードのアセンブリでも、SPAdesより10倍ほど長い配列が得られることもあります(scaffolding後のサイズ、精度は考慮せず)。
  16. ロングリードonlyアセンブリが小さなプラスミドのアセンブリに失敗しやすいことが報告されています(論文)。この論文では、比較に使われたCanu、Flye、Miniasm、Ravenが10kb以下のプラスミドの再構成に失敗しやすい傾向があったのに対し、Unicyclerはプラスミド配列を100%回収することに成功したことから、細菌ゲノムアセンブリにおいてプラスミド回収の可能性を高めるために、Unicyclerを使用する事が推奨されています。unicyclerのアセンブリはspadesベースなので、spadesの--plasmidモードもunicyclerと同様に小さなプラスミドのアセンブリ性能は高い可能性があります(要検討)。上記の論文がIJSEMの機関ジャーナルから出版されたことは、この事が下流解析に深刻な影響を及ぼす可能性の大きさを憂慮してかもしれません。
  17. ハイブリッドアセンブリモードでも、SPAdesは基本的にイルミナデータのみを用いてDe Bruijnグラフを構築する。ショートリードからDe Bruijnグラフを構築してアセンブリグラフに変換後、バルジ、チップ、キメラエッジを除去する。その後、このグラフのギャップクローズとリピート解消のためにHybridSPAdesのEXSPANDERはロングリードを使用する。具体的には、ロングリードをアセンブリグラフにマッピングし、リードパスを生成する。それからロングリードのコンセンサスを使ってアセンブリグラフのギャップを埋める(論文)。イルミナベースのアセンブリが中心であることに変わりはないが、ロングリードはパスの解決に役立つ。特に複雑なグラフの場合に有用とされている。中心にあるのはショートリードからのde Bruijnグラフであり、よってショートリードが不足している場合、ロングリードが十分あってもHybridSPAdesのパフォーマンスが出ない可能性がある。
  18. ファージゲノムなどのアセンブリでは、高いレベルのバリエーションによりSPAdesでは連続性の高いアセンブリが得られないことが多い。その場合、metaSPAdesなどのメタゲノムアセンブラを検討する。あるいは、リードのランダムなサブサンプリング(seqkit sample <FLOAT>)でリード数を限定することを考える(12と同じ理由)。ファージゲノムなどのアセンブリでは、低レベルのバリエーションを除外するunicyclerアセンブラも有効(注;ウィルスゲノムのHaplotype resolvedアセンブリを焦点にしている時)。

 

 

Quest(アセンブリ評価)

http://quast.bioinf.spbau.ru

 

  • Unicycler

アセンブリ改善と時間短縮

大きなk-merを使うためビルドし直す

 

127-mer以上のk-merを指定する場合、この投稿などを参考にしてください。

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

 

サイズとカバレッジでフィルタリング

 

SPAdes ファミリー(サードのツールも含む)

 

引用

ref.1

SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing

Anton Bankevich, Sergey Nurk, Dmitry Antipov, Alexey A. Gurevich, Mikhail Dvorkin, Alexander S. Kulikov, Valery M. Lesin, Sergey I. Nikolenko, Son Pham, Andrey D. Prjibelski, Alexey V. Pyshkin, Alexander V. Sirotkin, Nikolay Vyahhi, Glenn Tesler, Max A. Alekseyev,Pavel A. Pevzner, 

J Comput Biol. 2012 May; 19(5): 455–477.

 

ref.2

ExSPAnder: a universal repeat resolver for DNA fragment assembly

Prjibelski AD, Vasilinetc I, Bankevich A, Gurevich A, Krivosheeva T, Nurk S, Pham S, Korobeynikov A, Lapidus A, Pevzner PA.

Bioinformatics. 2014 Jun 15;30(12):i293-301

 

ref.3

Assembling short reads from jumping libraries with large insert sizes

Vasilinetc I, Prjibelski AD, Gurevich A, Korobeynikov A, Pevzner PA.

Bioinformatics. 05 Oct 5;(0):6-8.

 

ref.4

hybridSPAdes: an algorithm for hybrid assembly of short and long reads

Antipov D, Korobeynikov A, McLean JS, Pevzner PA

Bioinformatics. 2016 Apr 1;32(7):1009-15.

 

ref.5

TruSPAdes: barcode assembly of TruSeq synthetic long reads

Bankevich A, Pevzner PA

Nat Methods. 2016 Mar;13(3):248-50.

 

ref.6

metaSPAdes: a new versatile metagenomic assembler

Nurk S, Meleshko D, Korobeynikov A, Pevzner PA

Genome Res. 2017 May;27(5):824-834.

 

ref.7

BayesHammer: Bayesian clustering for error correction in single-cell sequencing

Nikolenko SI, Korobeynikov AI, Alekseyev MA

BMC Genomics. 2013;14 Suppl 1:S7.

 

2020 10/6

Using SPAdes De Novo Assembler
Andrey Prjibelski Dmitry Antipov Dmitry Meleshko Alla Lapidus Anton Korobeynikov
First published: 19 June 2020, https://doi.org/10.1002/cpbi.102

 

参考

バクテリアのゲノムアセンブリ

Tychus: a whole genome sequencing pipeline for as- sembly, annotation and phylogenetics of bacterial genomes

https://www.biorxiv.org/content/biorxiv/early/2018/03/16/283101.full.pdf

 

アセンブリの進め方。

http://compgenomics2017.biology.gatech.edu/images/a/ad/Genome_Assembly_Final_Pipeline_and_Results.pdf

 

追記

metaspadesのメモリがout of memoryになってしまう場合、fastqを分割してerror correctionするのもありかもしれませんね。他にGCsplitなども使えそうです。

out of memoryになってしまう場合、digital normalizationも有効だと思います。 


 *1

--isolateをつけるとエラーコレクションステップが無くなりますが、このオプションが実装されたのは、単離ゲノムで100x程度あればエラーコレクションの有無は結果に影響しないと判断されたからなのかもしれません。シークエンシング精度が、エラーコレクションが必要だった初期のNGSより改善されていることも背景にあるのでしょう。エラーコレクションが無くなってランタイムが短くなるのは大きなアドバンテージです。

 

memo

AttributeError: module 'collections' has no attribute 'Hashable'エラーがでたが、issue#873 に従って、constructor.pyのコードを修正すると直った。