2021 1/21, 1/30 文章修正
HTSのシークエンシングリードのデータ解析では、巨大なテキストファイルを読み込み、何らかの計算を行なって結果をファイル保存します。これを繰り返して最終的に生物学的な洞察に繋げるわけですが、この繰り返しで、大きなファイルの読み出し、書き出しが繰り返し発生することになります。これにはどれくらいの時間がかかっているのでしょうか？マッピングを例に、ファイルの書き出しが計算時間の中でどれくらいの比率を占めるのか調べてみました。
方法
１、データについて
A、スモールデータ
ショートリードのシミュレーターARTを使って酵母ゲノム（GCF_000002945.2）の30x カバレッジのペアエンドfastq（150-bp x 2）をシミュレートした。得られたfastq.gz（1.1GB x 2）をtimeの"real"で計測しながらGCF_000002945.2アセンブリにマッピングした。マッピングにはminimap2を-x srのパラメータプリセットで使用した。minimap2はcondaの仮想環境を作って最新版のバイナリをインストールした。
conda create -n benchmark python=3.9 -y
conda activate benchmark
conda install -c bioconda -y minimap2 samtools
B、ラージデータ
手元にあった海洋生物ゲノムのショートリードのペアエンドfastqを使用した。リード数は68,379,301x2（非圧縮だと19.7GBx２のファイルサイズ）。このfastqを該当するゲノムアセンブリにマッピングした。他の条件は１と同じ。 このデータのSAM出力ファイルは、samtools sortの読み込みにも使用した。
３、時間の計測方法
timeコマンドを使った。 timeコマンドの計測はそのままでは扱いにくいので、awkにパイプしてrealの時間だけ保存した。
for i in `seq 1 5`;
do
printf "$i\tTR3990x\t" >> bench;
time -p (minimap2 -ax sr -t 16 ref.fa R1.fq.gz R2.fq.gz 1> /dev/null 2> log) 2> time;
cat time |awk '/real/'|awk -F" " '{print $2}' >> bench;
done
上のコードを実行すると、下のようなランタイム計測結果が得られる。
> cat bench
1 TR3990x 4.88
2 TR3990x 4.90
3 TR3990x 4.89
4 TR3990x 4.93
5 TR3990x 4.93
（単位は秒）
後半の計測では、timeの代わりに自動でコマンドのベンチマークを取ってくれるhyperfineを使った。
４、ハードウエア
TR 3990xのCPUを積んだ計算機を使った。物理メモリは128GB。
結果
１、まずはminimap2を使ってfastqをゲノムにマッピングし、結果を/dev/nullに捨てた時とディスクに書き出した時の２条件で、どれくらいランタイムが異なるか調べてみます。方法に記載したスモールデータを使います。微妙な時間変化があるのかどうか調べるため、50回ループ実行します。
同一ディスクに書き出す場合
#タイムをbench1として保存
for i in `seq 1 50`;
do
printf "$i\tTR3990x\t" >> bench1;
time -p (minimap2 -ax sr -t 16 GCF_000002945.fna R1.fq.gz R2.fq.gz 1> out.sam 2> log) 2> time;
cat time |awk '/real/'|awk -F" " '{print $2}' >> bench1;
rm out.sam time log;
done
結果
> cat bench1
1 TR3990x 4.92
2 TR3990x 4.94
3 TR3990x 4.94
4 TR3990x 4.97
5 TR3990x 4.94
6 TR3990x 4.93
7 TR3990x 4.94
8 TR3990x 4.99
9 TR3990x 4.98
10 TR3990x 4.94
11 TR3990x 4.97
12 TR3990x 4.96
13 TR3990x 4.98
14 TR3990x 4.99
15 TR3990x 4.96
16 TR3990x 5.02
17 TR3990x 4.97
18 TR3990x 4.96
19 TR3990x 4.96
20 TR3990x 4.94
21 TR3990x 4.94
22 TR3990x 4.97
23 TR3990x 4.96
24 TR3990x 5.00
25 TR3990x 4.94
26 TR3990x 4.97
27 TR3990x 4.95
28 TR3990x 4.97
29 TR3990x 4.96
30 TR3990x 4.97
31 TR3990x 4.97
32 TR3990x 4.99
33 TR3990x 4.99
34 TR3990x 4.97
35 TR3990x 4.97
36 TR3990x 4.96
37 TR3990x 4.91
38 TR3990x 4.95
39 TR3990x 4.99
40 TR3990x 4.97
41 TR3990x 4.99
42 TR3990x 4.97
43 TR3990x 4.99
44 TR3990x 4.98
45 TR3990x 4.98
46 TR3990x 4.96
47 TR3990x 4.96
48 TR3990x 4.99
49 TR3990x 5.00
50 TR3990x 4.94
（単位は秒）
次に/dev/nullに捨てる場合。”1> /dev/null”としています。
#タイムはbench2として保存
for i in `seq 1 50`;
do
printf "$i\tTR3990x\t" >> bench2;
time -p (minimap2 -ax sr -t 16 GCF_000002945.fna R1.fq.gz R2.fq.gz 1> /dev/null 2> log) 2> time;
cat time |awk '/real/'|awk -F" " '{print $2}' >> bench2;
rm time log;
done
結果
> cat bench2
1 TR3990x 4.79
2 TR3990x 4.88
3 TR3990x 4.84
4 TR3990x 4.82
5 TR3990x 4.90
6 TR3990x 4.85
7 TR3990x 4.88
8 TR3990x 4.89
9 TR3990x 4.90
10 TR3990x 4.89
11 TR3990x 4.85
12 TR3990x 4.86
13 TR3990x 4.87
14 TR3990x 4.88
15 TR3990x 4.88
16 TR3990x 4.88
17 TR3990x 4.87
18 TR3990x 4.87
19 TR3990x 4.87
20 TR3990x 4.87
21 TR3990x 4.90
22 TR3990x 4.91
23 TR3990x 4.92
24 TR3990x 4.88
25 TR3990x 4.87
26 TR3990x 4.88
27 TR3990x 4.90
28 TR3990x 4.86
29 TR3990x 4.87
30 TR3990x 4.87
31 TR3990x 4.88
32 TR3990x 4.90
33 TR3990x 4.89
34 TR3990x 4.91
35 TR3990x 4.93
36 TR3990x 4.91
37 TR3990x 4.89
38 TR3990x 4.89
39 TR3990x 4.87
40 TR3990x 4.87
41 TR3990x 4.90
42 TR3990x 4.92
43 TR3990x 4.87
44 TR3990x 4.88
45 TR3990x 4.89
46 TR3990x 4.90
47 TR3990x 4.91
48 TR3990x 4.91
49 TR3990x 4.89
50 TR3990x 4.92
50回の結果をボックスプロットにまとめます。
図は縦軸がタイムです。０スタートでないことに注意して下さい。ランタイムの平均は、/dev/nullに捨てるbench2が4.87秒（S.D 0.02）、samファイル保存するbench1が4.97秒（S.D 0.03）でした。有意差検定は行なっていませんが、わずかにランタイムの違いが出ていそうです。
このベンチマークのbench1では、SAMの書き出しを読み出しと同一のディスクで行なっています。そのため、読み書きのオーバーヘッドが生じて、これも差の一因になった可能性があります。もし本当にオーバーヘッドが生じているなら、読み出しと書き出しを物理的に異なるディスクで行うことで、書き出し時間自体は短くできそうです。
２、この仮説が妥当かどうか検証するため、別の物理ディスクにファイル書き出しするテストを行います。微妙な差を判断しやすくするため、非圧縮で19.7GBx２もある大きめのfastqを使います（方法のラージデータ）。
計測結果をまとめるのが手間なので、以後の計測にはhyperfineを利用します。コマンドを指定回数だけ実行してランタイムの算術平均・中央値と分散、標準偏差（S.D）をまとめてくれる便利なコマンドです（紹介）。このhyperfineを使って、/dev/nullに捨てる時、読み出しと同一の物理ディスクにsamファイルを書き出す時、異なる物理ディスクにsamファイルを書き出す時のマッピングのランタイムを比較します。
hyperfineデフォルトの10回反復を実行
#/dev/null
hyperfine 'minimap2 -ax sr -t 24 ref.fasta R1.fq.gz R2.fq.gz > /dev/null'
#同一ディスクに書き出し
hyperfine 'minimap2 -ax sr -t 24 ref.fasta R1.fq.gz R2.fq.gz > temp.sam'
#異なるディスクに書き出し
hyperfine 'minimap2 -ax sr -t 24 ref.fasta R1.fq.gz R2.fq.gz > /media/kazu/6TB/temp.sam'
結果です。
- /dev/null => 419.9 s ± 0.4 s
- 同一ディスク（M.2 SSD）へ書き出し => 617.2 s ± 74.4 s
- 異なるディスク（SATA3接続6TB HDD）へ書き出し => 449.4 s ± 62.9 s
/dev/null書き出しと比べると、異なるHDDへの書き出しは1.07倍、同一ディスクへの書き出しは1.47倍のランタイムになりました。予想と傾向は一致しており、異なるディスクへの書き出しが同一ディスク書き出しより顕著にランタイムは短くなるという結果です（ただし分散は大きい）。
私が知る限り、minimap2はfastq/fastaを一定のリード数のチャンクの単位でメモリに読み込み、マッピングと書き出しのサイクルを繰り返す事でマッピングを進めます（"-K"で指定するメモリサイズ次第でチャンクサイズは変わります）。標準エラー出力を見ているとこれを観察できます。
m$ minimap2 -K 10M -t 4 -ax sr ref.fasta R1.fastq.gz R2.fastq.gz > 1.sam
[M::mm_idx_gen::0.123*1.25] collected minimizers
[M::mm_idx_gen::0.137*1.53] sorted minimizers
[M::main::0.137*1.53] loaded/built the index for 8 target sequence(s)
[M::mm_mapopt_update::0.137*1.53] mid_occ = 1000
[M::mm_idx_stat] kmer size: 21; skip: 11; is_hpc: 0; #seq: 8
[M::mm_idx_stat::0.144*1.51] distinct minimizers: 641305 (98.57% are singletons); average occurrences: 1.028; average spacing: 5.996
[M::worker_pipeline::0.871*3.64] mapped 39906 sequences
[M::worker_pipeline::1.399*3.93] mapped 38970 sequences
[M::worker_pipeline::1.933*4.04] mapped 38994 sequences
[M::worker_pipeline::2.464*4.11] mapped 39140 sequences
[M::worker_pipeline::3.001*4.16] mapped 39412 sequences
[M::worker_pipeline::3.547*4.20] mapped 39476 sequences
[M::worker_pipeline::4.113*4.22] mapped 40010 sequences
[M::worker_pipeline::4.654*4.24] mapped 39460 sequences
[M::worker_pipeline::5.287*4.25] mapped 39500 sequences
[M::worker_pipeline::5.868*4.26] mapped 39476 sequences
[M::worker_pipeline::6.440*4.27] mapped 39412 sequences
[M::worker_pipeline::7.019*4.27] mapped 39742 sequences
[M::worker_pipeline::7.475*4.21] mapped 35244 sequences
[M::main] Version: 2.14-r883
minimap2標準エラー出力の例（-K 10M）
minimap2は非常に高速にアラインしていくため、極めて短いタイミングで読み出し、書き出しのサイクルが繰り返されていきます。同一の物理ディスクに書き出すことで、読み書きの繰り返しによるオーバーヘッドが生じたり、シリコンディスクへ長時間強い負荷がかかって予想が難しい確率的な性能低下が発生しているのでしょうか（だからS.Dが大きい？）。
３、原因ははっきりしませんが、２つのテストから、大きなfastqのマッピングでは、読み出しとは異なるディスクにファイルを書き出した方が、マッピングタイムは大幅に短縮されるという結果が得られました（環境やプログラム次第で結果は大きく変わる可能性あり）。では、samからcoordinate_sorted.bamを作る時も同じ傾向が出るのでしょうか？次はこれを調べてみます。
２のsamファイル出力（およそ50GB）を利用します。samtoolsのsortコマンドでsamからbamに変換してソートします。
#/dev/null
hyperfine 'samtools sort -O BAM -@ 8 temp.sam > /dev/null'
#同一ディスクに書き出し
hyperfine 'samtools sort -O BAM -@ 8 temp.sam > sort.bam'
#異なるディスクに書き出し
hyperfine 'samtools sort -O BAM -@ 8 temp.sam > /media/kazu/6TB/temp.bam'
”-O BAM”は指定なしでもsort.bamができます（samtool v1.11使用)。ディスクに負荷がかかりやすいように、並列数は少し大きめの”-@ 8”を指定しました。
結果です。
- /dev/null => 200.4 s ± 48.4 s
- 同一ディスク（M.2 SSD）へ書き出し => 280.7 s ± 21.3 s
- 異なるディスク（SATA3接続6TB HDD）へ書き出し =>477.9 s ± 20.2 s
今度は２と反対の傾向、即ち、異なるディスク（HDD）へ書き出しした方が、同一の物理ディスク（SSD）に保存するより顕著に長い時間がかかるという結果でした。
なぜでしょうか？これはsamtools sortコマンドの処理の仕方が原因と考えられます。少し長くなりますが、自分が知る限り、samtools sortはデフォルトでプログラムが利用可能なメモリサイズのチャンクごとにsamをメモリに読み込み、ソートしてtemporalなsorted.bamを出力し、最後にそれらを結合します（最後の小さめのチャンクはインメモリからダイレクト）。プログラムのデフォルトでは-m768Mなので、768MBのSAMファイルまで読み込みバイナリ変換&ソートします。この挙動のため、バクテリアなどの小さなfastqを扱う場合は、temporalなsorted.bamは１つも出力されず最終のbamが出力されることも多くありますが、真核生物の巨大なsamファイルをバイナリ変換&ソートすると、ラン進行に伴って、作業パスに多くのtemporalなsorted.bamが書き出されていくのを観察する事ができます（大きい時には100以上のファイルが一時的に生じる）。このsamtools sortコマンドは、"-@ <INT>"というオプションも持っています。-@ <INT>は並列処理を行うためのフラグとして機能します。例えば”-m 500M -@ 4”のオプションを立てた時は、４つ並列でそれぞれ500MBまでSAMのデータをメモリに収納し、バイナリ変換&ソートして sort.bamに出力します。これをsamファイルを全て読み出し切るまで繰り返します（４つずつtemporalなsorted.bamが書き出され、最後に全ファイルが結合される）。この一旦書き出されるtemp.bamは、一定レベルで圧縮したバイナリファイルであり（*2）、temp.bamのファイルサイズは500MBよりずっと小さくなりますが、それでも大きなファイルを並列にディスクに書き込んでいるわけです。
samtools sortを小さなメモリサイズで実行すると、これを観察することができます。
samtools sortのtemp.bam出力の例（-m 10M）。最後に結合される。これらのtemp.bamはマッピングのサブセットそのものである。IGVに読み込めば、ゲノム全体にリードがまばらにマッピングされているのを確認できる。
説明が長くなりましたが、上のランで異なるディスクのランタイムが増えたのは、HDDをファイル書き出し先に指定した上に、並列数を8まで増やしてしまい、HDDへの書き出しが遅くなってしまったのが原因と推測されます（下の参考参照。今回のHDDはAHCIモードではない）。
このように書くともっともらしく聞こえるかもしれませんが、本当でしょうか？別の高速なディスク、ここではSSDの物理ディスクに書き出してこれを検証してみます。
#異なるSSDのディスクに書き出し
hyperfine 'samtools sort -O BAM -@ 8 temp.sam > /media/kazu/500GBSSD/temp.bam'
- 異なるディスク（NVMe 500GB SSD）へ書き出し => 201.6 s ± 0.5 s
/dev/nullと同じくらいのランタイムになりました。これなら上記の仮説は概ね正しいと言えそうです。しかしこの結果を見ると、外部SSDをファイル書き出し先に指定した場合はファイル書き出しはプログラムのボトルネックになってないですね。高速なファイル書き出し先が確保できる時は、並列数をより増やしてリニアにランタイムを短縮できるかもしれませんね。
まとめ
大小のマッピングファイルの書き出しを例に、ファイルの保存がランタイム全体に与える影響の大きさを示しました。samtools sortのボトルネックはファイル書き出しになりやすく（複数ファイルのシーケンシャルライト）、一方でminimap2のボトルネックは、プログラムは非常に高速にも関わらず、マッピングプロセスつまりプログラムそのものにあり、ファイルの書き出し（単独のファイルのシーケンシャルライト）はボトルネックにはなりにくいということだと解釈しています。とは言え、同一ディスクで読み書きすると、ボトルネックが読み書きに移ってくるということです。
今回のテストでは、プログラムに応じてランタイムを短縮するための処方箋は異なる、という割と当たり前の結論が得られました。しかし意外と意識していない事かもしれません。コマンドを実行する時の参考にして下さい。パイプ処理のタイム検証もしたかったのですが、解釈が難しくなるため今回は控えました。
2021 2/19 追記
久しぶりに読み返してみると、テストの動機と結論が繋がってなくて”おかしな記事”になっていることに気づきました。言いたかったのは、ファイルの読み書きを考慮してコマンドを実行しないとランタイムに大きな差が出てしまうケースがあるということ。そればかりか、ストレージに大きな負担をかけて故障リスクを上げてしまう可能性もあるのではないかということでした。
引用
Minimap2: pairwise alignment for nucleotide sequences.
Li H
Bioinformatics. 2018 May 10.
https://github.com/sharkdp/hyperfine
*1
sambamba sortでは９段階で圧縮レベルを選べる。レベル0で非圧縮バイナリー出力。レベル9で最大圧縮率だが計算負荷は最大。I/O負荷が減るメリットとトレードオフ。
*2
並列処理なので、メモリ要求量は最大指定数（-m）x並列数（-@）の倍数で増えることに注意が必要です。
また、並列処理するとディスクにはかなりの負担が掛かる点にも注意が必要です。-@ 16なら、temporalなsort.bamの書き出しでは、当然ながら16ファイル同時にシーケンシャルライトを実行していることになります。100MB〜数百MBくらいのファイルを16ファイル同時にコピーするのを想像してみて下さい。このジョブと並行して何らかの作業を行うなら、I/Oが律速になるだろうことは容易に想像できます。
参考
