NGSデータ解析マシンのスペックによるデータ解析時間の違い

あけましておめでとうございます。今年もよろしくお願いします。

　NGSのデータ解析で時折聞かれるのが、解析マシンのスペックはどれくらいあれば十分かというような質問である。一般的にはメモリリッチで、I/Oが高速で、CPUのマルチスレッドに対応したマシンが適してると思われるが、どのくらいが必要十分な性能で、どれくらいがオーバースペックになるのだろうか？NGSのツールは並列処理に対応しているものも多いが、スレッド数を増やしても頭打ちになる数値は意外に小さい。メモリも16GB程度でかなりのリード数をアセンブルできるケースもあれば、歯が立たないケースもある。また、いつもは数分で終わる処理が、複数のジョブを走らせてしまい、HDDのI/Oが足枷となり数倍の時間がかかってしまうこともある。先入観があり、意外な思い込みがあるかもしれない。シミュレーションのデータ解析を複数ハードウエアで実施して調べてみることにする。

　テストには、ゲノムサイズが100-Mb程度で、nの核相世代がある紅藻Chondrus crispusを用いる。100メガ 、1ギガ以上の動植物ゲノムは大型のクラスタシステムでde novoアセンブルするのが一般的なので、100メガというサイズは単体のLinuxマシン、MAC 、PCで扱う最大規模のデータサイズと言える。このような単体マシンには負荷が大きいシミュレーションデータの解析時間を複数ハードで比較することで、NGSのデータ解析にはどのようなスペックが必要十分なのか推察する。

流れは再現できるようにしています。興味があればテストしてみてください（注　ループは各マシン４回ずつ行っていますが、１回ループでも１日近く掛かります）。

使用したハードウエア

１、Apple Macbook pro retina15.4インチ  2015上位 (リンク)

• 　CPU: Core i7 2.5 GHz/4コア8スレッド
• 　memory: DDR3L 1600 MHz 合計16 GB
• 　GPU: Intel Iris Pro Graphics/AMD Radeon R9 M370X
• 　system: drive: M.2 SSD：512GB
• 　OS: OSX sierra

２、Apple mac pro 2012 early （リンク）

• 　CPU: xeon (Westmere EP) 2.66 GHz/6コア x 2 (合計24スレッド)
• 　memory: DDR3 1,333 MHz 8GB x 8 合計64 GB
• 　GPU: GTX 660Ti
• 　system: drive: SATA3 SSD：240 GB x 2 (RAID0のストライピング)
• 　OS: OSX sierra

３、Apple mac pro 2012 early （リンク）

• 　CPU: xeon (Westmere EP) 3.46 GHz/6コア x 2 (合計24スレッド)
• 　memory: DDR3 1,333 MHz 8GB x 12 合計96 GB
• 　GPU: GTX 680
• 　system: drive: SATA3 SSD：512 GB x 3 (RAID0のストライピング)
• 　OS: OSX sierra

４、Dell  T7610（リンク）

• 　CPU: xeon E5 2697 v2  2.7 GHz/12コア x 2 (合計48スレッド)
• 　memory: DDR3 1,333 MHz 8GB x 4 合計32 GB
• 　GPU: quadro600
• 　system: drive: SATA3 SSD：240 GB
• 　OS: ubuntu14.04

５、自作ブレードサーバ

• 　CPU: xeon E5 2650 v3  2.5 GHz/10コア x 2 (合計40スレッド)
• 　memory: DDR4 1,600 MHz 16GB x 12 合計192 GB
• 　GPU: quadro600
• 　system: drive: SATA3 7200rpm HDD：3 TB
• 　OS: cent OS6

簡単な解説

• １、macbook proは、CPU、メモリ、GPU、ストレージのバランスが非常によく考えられている。M.2 SSDはSATA SSDよりランダムアクセスも速く、CPUのアーキテクチャも比較対象中では最も新しい。ただしメモリは最大16GBしかなく、冷却系もノートPCのため貧弱で、冷却が追いつかなければCPU、SSDの性能が低下する可能性がある（価格comにレビュー書いてます（リンク））。
•  ２、MAC PRO 2012はCPUコア数は多いが、CPUのアーキテクチャは比較対象では最も古い。macbook proよりメモリを4倍多く積んでいる。
• ３、２と同型機だが、CPUをMAC PRO 2012標準の2.66GHzから3.4GHzまでクロックアップしたX5690に換装し、メモリ容量は96GBまで増設、システムディスクをPCIスロットに繋ぎ3台でストライピングするなどの改良を施している。旧mac proではほぼ最速構成。
• ４、xeon E5のv2のワークステーション。似たモデルにHPのz820がある。E5 2697v2は現mac proの最上位モデルと同じCPUで、それが2-wayで搭載されている。CPUの世代は２−３のmac proより２世代分新しく、現行mac proと同じモデルとなる（ただし現行mac proはなぜか1-way）。CPUは高性能だが、メモリ容量は32GBしかない。
• 5、xeon E5のv3を2-wayで積んだ一般的なブレードサーバ。自分で組んだ。この中では唯一DDR4メモリになる。メモリ容量は192GBで、NGSでまことしやかに重要と言われるメモリは比較対象中で圧倒的に多い。ただしシステムドライブもデータドライブも7200rpmのHDDである。

ベンチマーク時のハードウエア設定

部屋の温度15-20度。他のジョブは実行しない。再起動して10分後に測定開始。mac proのファンは最大回転速度に設定。macbook proは電源につなぎ、USBファンで風を当てながら実行。

プログラムのバージョン

• samtools 1.3.1
• bwa 0.7.17
• art 2.5.7
• GATK 3.8
• SOAPdenovo 2.0.4
• spades 3.11
• afterQC 0.9.6

--準備--　

Chondrus_crispusのゲノムをensembl plantからダウンロードする。

NGSのツールは一般にNを想定していないので、sedで除いてから実験を開始する。100kbp以下のcontigも除く。

#Ensembl plantからマスクなしのゲノムをダウンロード。
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-37/fasta/chondrus_crispus/dna/Chondrus_crispus.ASM35022v2.dna.toplevel.fa.gz

#解凍
gzip -dv Chondrus_crispus.ASM35022v2.dna.toplevel.fa.gz

#Nを除き、空行をけずり、80文字改行。
sed -e 's/N//g' Chondrus_crispus.ASM35022v2.dna.toplevel.fa |grep -v "^s*#" |grep -v "^s*$" |fold -w 80 > C_crispus_change1.fa

#長さで降順ソート。100-kb以下の配列があれば捨てる。(seqkit)
seqkit sort -r -L 100000 --quiet -l C_crispus_change1.fa > C_crispus_change2.fa

#念のためembossのseqretで修復(embossはbrewで導入可能)
seqret

#出力名はC_crispus_custom.faとした。

以下のようなスクリプトsimulation.shを用意し、４回ループで実行して、各コマンド処理時間をtimeコマンドで計測する。

simulation.sh

#!/bin/sh 
#benchmark 
#注1 timeの結果は標準エラー出力なので2>となる。 
#注2 wgsimのランダム発生器はマイループ異なる数値になるようにした。 


f='C_crispus_custom.fa' 
thread='8' 　 #CPUスレッド数 #macbook proは8、他は24とした。
cov='100' 　　#カバレッジ 
insert='500' #インサートサイズ 
size='150' 　 #リード長 
stdev='50' 　 #SD

#index作成 
bwa index -a is $f #2GB以下のsmall genomeはモードis、 prefixはfastaファイル名
samtools faidx $f 
picard CreateSequenceDictionary R=C_crispus_custom.fasta O=C_crispus_custom.dict 

#java1.8
export JAVA_HOME=$(/usr/libexec/java_home -v 1.8) #java1.8にする。 

for i in 1 2 3 4 
 do
 pyenv global 2.7.9 #python2.7環境で実行
 echo ${i}_time_start >> running_time


 #1 illumina Hiseq2500 simulation. insert size 500bp, SD50 paired-end.coverage x100
 echo =================art================= >> running_time
 time(art_illumina -ss HS25 -i $f -p -l $size -f $cov -m $insert -s $stdev -o HS25 2> /dev/null) 2>> running_time
 rm *aln
 fastqF=HS251.fq
 fastqF=HS251.fq
 echo >> running_time


 #2 30> quality trimming, and 5-prime 5bp cutted.
 echo =================after.py================= >> running_time
 time(after.py -1 $fastqF -2 $fastqR -q 30 -s 20 -f 5 --no_overlap 2> /dev/null) 2>> running_time
 mv good/* .
 rm -rf bad/ good/
 fastqFQC=${fastqF%fq}good.fq
 fastqRQC=${fastqR%fq}good.fq
 echo >> running_time


 #3 Pacbio long read simulation
 echo =================pbsim================= >> running_time
 time(pbsim --prefix pbsim --data-type CLR --depth 5 --model_qc data/model_qc_clr $f 2> /dev/null) 2>> running_time
 cat pbsim*fastq > long_reads.fastq
 rm pbsim*
 long=pbsim_long_reads.fastq
 echo >> running_time


 #4 read statistics
 echo =================seqkit================= >> running_time
 (time seqkit stats $fastqFQC $fastqRQC $long 1>> reads_summary) 2>> running_time
 echo >> running_time


 #5 bwa mem,sort,index, and merge. (sortでsamからダイレクトに.bamを出力できるが、gatkでエラーになるので、viewにパイプで渡す)。
 echo =================bwa_mem_and_samtools_view_sort_index================= >> running_time
 #ショートリード
 echo short read >> running_time
 time(bwa mem -R "@RG	ID:X	LB:Y	SM:Z	PL:ILLUMINA" -t $thread $f $fastqFQC $fastqRQC 2>/dev/null |samtools view -S -b -o pair.bam - && samtools sort -@ $thread pair.bam > pair_sorted.bam && samtools index pair_sorted.bam) 2>> running_time
 #ロングリード
 echo long read >> running_time
 time(bwa mem -R "@RG	ID:X	LB:Y	SM:Z	PL:PACBIO" -x pacbio -t $thread $f $long 2>/dev/null |samtools view -S -b -o long.bam - && samtools sort -@ $thread long.bam > long_sorted.bam && samtools index long_sorted.bam) 2>> running_time
 #マージ time(samtools merge -@ $thread -f merged.bam long_sorted.bam pair_sorted.bam 2> /dev/null && samtools index merged.bam) 2>> running_time
 echo >> running_time


 #6 GATK haplotypecaller(short read)
 echo =================gatk haplotypecaller================= >> running_time
 time(gatk -T HaplotypeCaller -R $f -I pair_sorted.bam -o raw_variants.vcf 2>/dev/null) 2>> running_time
 echo >> running_time


 #7 de novo assembly1 SOAPdenovo2(63mer)
 echo =================SOAPdenovo2================= >> running_time
 time(SOAPdenovo-63mer all -p $thread -s SOAPdenovo.config -o assembly 2>/dev/null) 2>> running_time
 mkdir SOAPdenovo$i 
 cp assembly.scafSeq SOAPdenovo$i/SOAP_scaffolds.fasta 
 rm assembly*
 echo >> running_time


 #8 de novo assembly2 spades short + long(k=auto)
 time(spades.py -t $thread --careful -1 $fastqFQC -2 $fastqRQC --pacbio $long -o spades_assembly) 2>> running_time
 echo >> running_time
 

 #9 quast assembly evaluation
 echo =================quast================= >> running_time
 time(quast SOAPdenovo$i/SOAP_scaffolds.fasta spades_hybrid_assembly$i/scaffolds.fasta -R $f -o assembly_performance$i -t $thread) 2>> running_time
 echo >> running_time


 rm *bam*
 echo $i simulation end

done

 pyenvでpythonバージョンを切り替えます。上のコードの通り実行するならpyenvを導入して"pyenv install 2.7.9"と"pyenv install 3.4.3"を実行してから、おこなって下さい。導入していないツールがあれば、brewで導入しておいてください。

brew install samtools bwa gatk spades unicycler soapdenovo 

pbsimはソースをダウンロードしてビルドします（解説リンク）。 seqkitはバイナリをダウンロードしてパスを通します（解説リンク ）。

実行方法

chmod u+x simulation.sh
bash simulation.sh

簡単に流れを書いておきます。

１、artのHiseq2500のエラー設定で、Chondrus crispusのゲノムサイズのx100リードを発生。

　↓

２、afterQCでQ30>トリミング

　↓

３、Pacbioの設定でロングリードをChondrus crispusのゲノムサイズのx5リードを発生。今回はエラートリミングはしない。

　↓

４、リード数をカウント。特に他の作業と関連性はない。ルーチンワークの１つとして取り込んだ。

　↓

５A、bwa mem でショートリードをアライメント。ポジションソートしたbamにする。

　↓

５B、bwa mem でロングリードをアライメント。ポジションソートしたbamにする。

　↓

５C、AとBのbamをマージ。

　↓

６、GATKでSNVとsmall indelをコール。ショートリードのbamを使う（indel周辺のリアライメントはここでは行わない）。

　↓

７、SOAPdenovo2でショートリードをde novo assembly。(GAP closureはなし)

　↓

８、spadesでショートリードとロングリードをhybridモードでde novo assembly。

　↓

９、quastでアセンブリを評価。

テスト回数４回 (n=4)

グラフのバーはいずれも短い方が高速です。テストによって横軸の単位は秒(s)、分(m)、時間(h)に変えています。

１、artを使い、ショートリードをx100カバレッジ発生させるのに要した時間。

このコマンドで9.3GBのfastqが２つできます。artは計算の並列化には対応していません。CPUの世代が新しいほど早いという傾向がありそうですが、それにしてはmacbook proだけかなりの時間がかかっています。大量のリードを発生させるならxeonのサーバーで行った方が良さそうです。

２、afterQCでクオリティの分析とトリミングに要した時間。

afterQCも計算の並列化には対応していません。CPUのアークテクチャが新しく、クロックも高い方が高速という結果となりました。mac pro2.66GHzでは2時間以上かかってしまいますが、xeonマシンやmacbook proに切り替えれば１時間でクオリティチェックとトリミングが完了します。

３、pbsimを使い、ロングリードをx5カバレッジ発生させるのに要した時間。

pbsimも計算の並列化には対応していません。artとは異なり、処理時間に大きな差はないようです。macbook proが一番早く終わります。

４、seqkitでリードの分析に要した時間。

ショートリード2800万x2とロングリード15万x1を分析しています。seqkit は並列化に対応してますが、それよりCPUのアーキテクチャの世代の方が聞いているように見えます。がしかし、数秒の差なので、この程度なら処理時間に大きな差はないと言えそうです。

５A、ショートリードx2をbwa memでアライメントして、bamに変換してソートし、indexをつける処理に要したトータルの時間。

以外にもmac proが一番高速という結果になりました。xeon E5 2697v2の48threadでもmac pro 2.66GHzに追いつけていません。ただ、今回bwaはパイプでbwa memとviewとsortのトータルタイムを計ってしまっています。複数の要素が関係しているので考察は困難です。面倒でも別々にデータを取った方が良かったかもしれません。

それでも結果をまとめるなら、bwaはmac OS系の処理時間が早い、と言えそうです。

2019 1/11追記

私のケアレスミスでsamtols viewとindexなどがマルチスレッドになっていませんでした。そのため、特定のフローだけシングルスレッドで動作したことになります。このことに起因してCPUのクロックが高いmac proが早く終わった可能性があります。

ぷっぷくさん、ご指摘ありがとうございます。

５B、ロングリードをbwa memでアライメントして、bamに変換してソートし、indexをつける処理に要したトータルの時間。

ショートリードと同様にmac proが一番高速という結果になりました。

７、SOAPdenovo2でショートリードx2 (10GB x 2) をアセンブルするのに要した時間 (k=63)。

de novoアセンブルはノートPCには厳しいのが分かります。mac proよりサーバーマシン２台の方が早く終わり、10GBx2のデータを１時間でアセンブルできています（ゲノムがシンプルなのが大きい）。threadを48まで増やしても差は出ていません。macbook proで時間がかかった理由は、スレッド数が少ないこと、メモリ不足、などが考えられます（swapのログは取ってません）。無事アセンブルできたのは意外でした。 

８、spadesでショートリードとロングリードをハイブリッドアセンブルするのに要した時間。

spadesはまずエラー訂正を行い、それからk-merを変えて、アセンブルを繰り返します。ピークメモリ使用量は30GBくらい使っていました。

処理はサーバーマシンが圧倒的に高速でした。mac pro3.4GHzで14時間かかる処理をxeon E5 2697v2サーバー（32GBメモリしかないにもかかわらず）はわずか2時間超で終えます。一方、macbook proはエラー訂正でメモリ不足でspadesがクラッシュしました。 エラー訂正なしでランすれば可能かもしれませんが、mac proでも14時間程度かかる処理をmacbook proでするのはおすすめしません(頻繁なswapがあればSSDを痛める)。

しかしmac proとxeonサーバーでここまで差がついた原因ははっきりしません。実験前は、3.4GHz 96GBメモリmac pro とxeon E5 2697v2 32GBメモリサーバは僅差と予想していましたが、スレッド数を揃えてもサーバーの方が3倍高速にアセンブルを完了しました。スレッド数を倍にすると6倍〜7倍の差が出ています。spadesは内部で大きく４つの処理を行いますが、そのうちのいずれかのプロセスでxeon E5v2サーバーは劇的に高速化しているのかもしれません。しかし、CPUの世代が効いていると仮定するとxeon E5 v3サーバは遅すぎます。データが変わると、大きく結果が変わってくる可能性もあります。

quastとGATKですが、サーバーの方がエラーを吐いたので結果は載せていません。

  まとめ

 もう少しCPUの世代、クロックにリニアに比例した結果になるかと予想していましたが、ジョブにより結果は大きく変動しました。spadesのようにxeonのe5世代が圧倒的に高速な処理もあれば、bwaのようにmac proがxeon e5世代より早く処理できるジョブもありました。また、メモリが16GBしかないmacbook proも思った以上に戦えているのも印象的でした。次世代macbook proはメモリ32GBやCPU6コアのモデルが期待されますが、そうなればNGSの解析マシンとしてもっとバリバリ使えそうですね (既定路線であって時間の問題)。

話を戻します。今回の結果から、いずれのジョブでも常に高速な万能マシンを選ぶのは簡単でないことが分かりました。もちろんもっと周波数が高くアーキテクチャも改良されコア数も増えたCPUや (e.g., 7980xe)、高クロック大容量メモリに切り替えれば安定して高速な結果を得られるのでしょうが、かなりの費用が発生します。それよりも使用スレッドの検討や (下の追試参照)、スワップが頻繁に発生しているならメモリ増設を行う方が費用効果は高いでしょう。また、NGSに限らずインフォマティクスのツールは、年々、より精度が高くよりメモリ効率が高いツールが発表され続けています。テストして良好なパフォーマンスを示すなら、切り替えも検討するべきでしよう。

追試

bwa memベンチマーク

bwaは複数コマンドのトータルタイムを取っており、解釈が難しかったので、10mのゲノムを使い、スレッドを変えながらアライメント時間を測ってみた(bwa memだけの時間)。

使用マシン#4

#４、Dell  T7610（リンク）

• 　CPU: xeon E5 2697 v2  2.7 GHz/12コア x 2 (合計48スレッド)
• 　memory: DDR3 1,333 MHz 8GB x 4 合計32 GB
• 　GPU: quadro600
• 　system: drive: SATA3 SSD：240 GB
• 　OS: ubuntu14.04

結果

スレッド数が20〜30で完全に頭打ちになっている。48コアだと逆に速度が低下した。並列化が進みすぎて個々のCPUを効率的に使えていない可能性がある。

追記　もしくはハイパースレッディングの論理コア使用で頭打ちになっているのかもしれない（検証には、よりたくさんの物理コアが1wayで搭載された１ノードのサーバーマシンが必要）。

追記

目視で確認する限り、リニアに処理時間が短くなっているのは、せいぜい2コア、または4コアまでか。

注１

今回は32GBしかメモリを載せていないlinuxサーバーでも10GBx2リードのDe brujinのgraphを使ったアセンブルをできてますが、ピークメモリ使用量は30GBを超えていました。サーバーのswap領域を"swapon -s"コマンドで確認したところ、33GBほど確保されていました。スワップが起きてもメモリ要求量がおよそ65GBになるまでランできる計算になりますが、スワップがおきると処理速度が大幅に低下します。また、一般にリアルデータはシミュレーターで再現できないエラーを含んでおり、k-merの多様性はシミュレーションよりリアルデータのほうが大きくなります。10GBx2のリアルデータを32GB程度のメモリのマシンで解析するのは難しいでしょう。特にpolypoidのゲノムや、純化されてないゲノム、メタゲノムだったりすると困難です。いちおうswap領域のサイズは変更できるので（変更例）、十分割り当てればランはできるのでしょうが、それよりは素直にメモリ増設を検討したほうが幸せになれると思います（またはメモリ使用量が少ない方法論を探す）。

参考

Trinity CPU Usage