macでインフォマティクス

macでインフォマティクス

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

SPAdesとUnicyclerでlarge k-merを使う part2 (SPAdesのテスト)

benchmark assembly tips

2019 12/8 誤字修正

 

127以上のk-merを使うために、SPAdesとUnicyclerをビルドし直した（リンク）。今回は、実際に127以上のk-mer値でアセンブリを行い、アセンブリ性能がどのように変化するか簡単にテストした結果を書く。 Real dataの傾向が知りたいので、GAGE-BのMiseqでシーケンシングされたVibrio cholerae CO 0132データセット(250bp x 2、x100)を使う（GAGE-B paperの表１参照）。

HP:  https://ccb.jhu.edu/gage_b/datasets/index.html

f:id:kazumaxneo:20181215114845j:plain

 テストOS

  • ubuntu 18.04、Anaconda3.5.0

 

結果

Vibrio cholerae CO 0132の250bp x 2のペアエンドfastqをHPからダウンロードした。

 

１、merged.fqを作成

bbmerge.sh in1=reads1.fq in2=reads2.fq out=merged.fq outu1=unmerged1.fq outu2=unmerged2.fq ihist=ihist.txt

#リード長を分析
seqkit stats *fq.gz

file                 format  type  num_seqs      sum_len  min_len  avg_len  max_len

GAGEB_reads_1.fq.gz  FASTQ   DNA    796,812  199,999,812      251      251      251

GAGEB_reads_2.fq.gz  FASTQ   DNA    796,812  199,999,812      251      251      251

merged.fq.gz         FASTQ   DNA    755,374  149,383,904       39    197.8      493

mergedの平均サイズが短い。あまり長いライブラリは作られていない。

 

２、fastpによるクオリティチェック（出力指定なしだとfastpはQT前後の図だけ出力する）

fastp -i GAGEB_reads_1.fq.gz -I GAGEB_reads_2.fq.gz

GAGEB_reads_1 (before)

f:id:kazumaxneo:20181215120025j:plain

GAGEB_reads_2 (before)

f:id:kazumaxneo:20181215120021j:plain

いつものように、R1は3'末端1bpだけクオリティが低い。R2はより悪く、200bpあたりでQ30を下回っている。エラー率がアセンブリへ与える影響を調べるためには、トリミングの有り・無しでアセンブリ結果を比較すべきだが、今回は時間の関係でQTあり/無しの比較は行わない。トリミングなしで進め、k-merの選択による違いだけ評価する。

 

３、spadesパッケージのBayesHammer(link)を使ってerror correction

spades.py -1 GAGEB_reads_1.fq.gz -2 GAGEB_reads_2.fq.gz --merged.fq.gz --only-error-correction -o error_correction


#出力をカレントにコピー
mv error_correction/corrected/*gz .

#rename
mv GAGEB_reads_1.00.0_0.cor.fastq.gz GAGEB_1_cor.fq.gz
mv GAGEB_reads_2.00.0_0.cor.fastq.gz GAGEB_2_cor.fq.gz
mv merged.00.0_1.cor.fastq.gz merged_cor.fq.gz

 

４、de novo assembly

spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21 -o spades_k21 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31 -o spades_k21-31 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41 -o spades_k21-41 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51 -o spades_k21-51 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61 -o spades_k21-61 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71 -o spades_k21-71 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81 -o spades_k21-81 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91 -o spades_k21-91 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,101 -o spades_k21-101 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111 -o spades_k21-111 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121 -o spades_k21-121 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131 -o spades_k21-131 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141 -o spades_k21-141 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151 -o spades_k21-151 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161 -o spades_k21-161 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171 -o spades_k21-171 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181 -o spades_k21-181 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191 -o spades_k21-191 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191,201 -o spades_k21-201 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191,201,211 -o spades_k21-211 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191,201,211,221 -o spades_k21-221 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191,201,211,221,231 -o spades_k21-231 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191,201,211,221,231,241, -o spades_k21-241 &&
spades.py -1 GAGEB_1_cor.fq.gz -2 GAGEB_2_cor.fq.gz --merged merged_cor.fq.gz --careful -t 40 --only-assembler -k 21,31,41,51,61,71,81,91,111,121,131,141,151,161,171,181,191,201,211,221,231,241,249 -o spades_k21-249

 （xeon E5 v3 x2サーバーでトータル10時間ほどかかった）

 

５、評価

quast k21.scaffolds.fasta k21-41_scaffolds.fasta ...（省略） \
-t 24 -o quast_result -1 GAGEB_reads_1.fq.gz -2 GAGEB_reads_2.fq.gz \
-R GCF_000006745.1_ASM674v1_genomic.fna

k21は極端に悪い。その後、k-merサイズが増えるにつれ、アセンブリの連続性は徐々によくなる。k100以上になるとほぼ変化がない。

f:id:kazumaxneo:20181215114735j:plain

 Quast report1

リファレンスがあるので、それをgold standardとして、アセンブリ精度に関する分析も行う。

f:id:kazumaxneo:20181215114748j:plain

 Quast report2

f:id:kazumaxneo:20181215114655j:plain

icarus.html

 

 Quast report2の重要な項目のみ、表にまとめ直した(*2)。

f:id:kazumaxneo:20181217095430p:plain

  1.  scaffolds数はk-mer値が増えるにつれ減少
  2.  10kb以上のscaffolds数もk-mer値が増えるにつれ減少
  3.   NG50はk-mer値180付近で最大になり、それ以上k-merを増やすと減少
  4.   genome fraction（カバー率）はk-mer値180-200付近で最大になり、それ以上k-merを増やすと0.2%ほど減少 (10000-bp前後)
  5. ミスアセンブリはk-mer値が小さい方が少ない。k-mer値が190を越えるとミスアセンブリ数は10を越える。その大半がrelocations。
  6. ミスマッチ、indel数は大差なし

k-mer値が増えるにつれscaffolds数が減少しているが、NG50もk-merが増えると減少している。アセンブリの連続性が改善されてscaffolds数が減少しているなら、直感的にはNG50も増えるように思える。最初は理由がわからなかったが、アセンブリのGFAファイルをbandageで可視化してみると、何が起きているのか理解できてきた。6に続く。
 
 

６、Bandageを使いGFAファイルを可視化する。端末から呼び出し、SVG出力。

K21

Bandage image spades_k21/assembly_graph.fastg k21.svg

f:id:kazumaxneo:20181215130005j:plain

 K21-71

Bandage image spades_k21-71/assembly_graph.fastg k21-71.svg

f:id:kazumaxneo:20181215125714j:plain

 K21-101

Bandage image spades_k21-101/assembly_graph.fastg k21-101.svg

f:id:kazumaxneo:20181217124848j:plain

K21-131

Bandage image spades_k21-131/assembly_graph.fastg k21-131.svg

f:id:kazumaxneo:20181215125620j:plain

K21-171

Bandage image spades_k21-171/assembly_graph.fastg k21-171.svg

f:id:kazumaxneo:20181215125543j:plain

K21-191

Bandage image spades_k21-191/assembly_graph.fastg k21-191.svg

f:id:kazumaxneo:20181215125914j:plain

K21-231

Bandage image spades_k21-231/assembly_graph.fastg k21-231.svg

f:id:kazumaxneo:20181215125435j:plain
　Vibrio choleraeのゲノムはchromosome２つだけからなる。k-merが大きくなるについてgraphはシンプルになるが、k-mer190くらいからはorphanなcontigが増えてくる。k-merが小さいうちはこのような断片化は起きなかったことも加えて考えると、断片化は、k-merピークが得られるか得られないのボーダーギリギリまでk-merの値を増やすことで、k-merカバレッジが足りない領域ができて発生したと推測される。

　この仮説は、先ほどの表で、scaffolds数の減少とNG50減少が同時に起こったことも説明できる。すなわち、k-mer値を増やすと、より長いunitigを作り出せるが、同時にk-merカバレッジが不足する領域が確率的に発生し、そこでアセンブリが切れてしまう。k-merを減らしすぎると、断片化は防げ、計算も早くなるが、今度は短いリピートの解消もできなくなる。例えば、上のデータでは、k-merが100から200くらいがバランスがいいように感じられるが、これは、リード長はいくつか、どれだけのカバレッジが得られたか、どのくらいのエラー率か、どのくらい複雑なゲノムか（*3）、またどれだけ計算リソースを使えるか、などによって変化しうる（GAGE-Bの結論とだいたい同じ）。de Bruijn Graphアセンブリでより長いcontigを作る上でのベストなk-merサイズ選択の難しさを示している（*4）。

 

結論

一度考えましたが考察が甘かったので破棄しました。もう少し考えて追記します。

 

 

 

 以上、１例報告であまり参考になるようなものでなかったかもしれませんが、簡単なテストでした。

 

今回使用したツール


引用

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.

 

参考

*1

バーコードを付けた相乗りシーケンスなら、１サンプルあたりのコストは大きく抑えられるが、そうでないなら、ロングリードシーケンスのコストは、バクテリア向けには高い。

 

*2

上では20merずつ変化させたが、ほんとは10merずつ増やしている。見えないかもしれませんが、全データを載せておきます。

f:id:kazumaxneo:20181217022850p:plain 

 こちらはGAGE-B 当時のアセンブリstatistics。

f:id:kazumaxneo:20181216205517j:plain

 

*3

極端な話、全くリピートがないスモールゲノムなら、短いk-merでも完全なアセンブリが可能になる。

 

*4

https://github.com/rrwick/Unicycler

f:id:kazumaxneo:20190109161333j:plain

UnicyclerのGithub HPより転載