chromapは高速・高精度にリードをゲノムにマップするツールです*1。 既存ツールのbwaやbowtie2よりも遥かに高速だとのことなので、ここではChIP-seqデータを用いて使用法及び速度・精度面を確かめてみたいと思います。
chromapについて
chromap は minimap2*2の後継です。NGS解析では通常リードがどこにどのくらいマップされたかのカウントだけを見るため、１塩基単位での詳細なアラインメントは必要なく、その部分を省くlight-weight alignmentという方法を用いることで、bwaやbowtie2よりも数倍高速にマッピングを達成しています（ざっくりした説明）。 minimap2は50bp以上のseed領域の一致を必要とするため、scATAC-seqのようにバーコード配列を含むリードのマッピングには不向きであるという弱点がありました。chromapはその弱点を克服し、アラインメントだけでなく冗長リードの除去やアダプター配列の同定までを高速かつ包括的に行います。
ChIP-seqデータはシングルエンドかつアダプター配列も含まれていないのでchromapの真価は十分に発揮されないのですが、簡単のため（あと私が一番慣れているため）フィージビリティスタディとしてここでは利用します。
インストール
condaでインストールするのが簡単です。
conda install -c bioconda chromap
indexの作成
ゲノム配列を入力にしてindexを作成します（ここではhg38.fa）。-i がindex作成のためのオプションですね。
ゲノムファイルの作成方法はこちらを参照してください。
chromap -i -r hg38.fa -o chromap-UCSC-hg38
ChIP-seqファイルの取得
なんでもよいのですが、ここではPol2 ChIP-seq (K562)のデータをSRAよりダウンロードしました。SRR227359.fastq.gzという名前で保存されます。
fastq-dump --gzip SRR227359
マッピング
chromapでhg38にマッピングしてみましょう。わかりやすさのため、打ち込むコマンドの前には$マークをつけています。
chromapではマッピングの時にもゲノムデータを指定する必要があります。用いるCPU数は24です。
$ time chromap --preset chip -t 24 \ $ -x chromap-UCSC-hg38 -r hg38.fa \ $ -1 SRR227359.fastq.gz -o SRR227359.chromap.bed Preset parameters for ChIP-seq are used. Start to map reads. Parameters: error threshold: 8, min-num-seeds: 2, max-seed-frequency: 500,1000, max-num-best-mappings: 1, max-insert-size: 2000, MAPQ-threshold: 30, min-read-length: 30, bc-error-threshold: 1, bc-probability-threshold: 0.90 Number of threads: 24 Analyze bulk data. Won't try to remove adapters on 3'. Will remove PCR duplicates after mapping. Will remove PCR duplicates at bulk level. Won't allocate multi-mappings after mapping. Only output unique mappings after mapping. Only output mappings of which barcodes are in whitelist. Output mappings in BED/BEDPE format. Reference file: hg38.fa Index file: chromap-UCSC-hg38 1th read 1 file: SRR227359.fastq.gz Output file: SRR227359.chromap.bed Loaded all sequences successfully in 50.49s, number of sequences: 25, number of bases: 3088286401. Kmer size: 17, window size: 7. Lookup table size: 392946789, occurrence table size: 445474231. Loaded index successfully in 174.56s. Loaded 500000 reads in 1.01s. （中略） Loaded 314483 reads in 0.53s. Mapped in 0.56s. No more reads. Mapped in 0.18s. Mapped all reads in 30.17s. Number of reads: 21814483. Number of mapped reads: 12508635. Number of uniquely mapped reads: 11368866. Number of reads have multi-mappings: 1139769. Number of candidates: 227477151. Number of mappings: 12508635. Number of uni-mappings: 11368866. Number of multi-mappings: 1139769. Sorted, deduped and outputed mappings in 8.99s. # uni-mappings: 10833419, # multi-mappings: 1020796, total: 11854215. Number of output mappings (passed filters): 9211499 Total time: 86.82s. real 1m32.913s user 4m12.370s sys 0m27.555s
--preset chipオプションでChIP-seqに合わせたパラメータセットになります。 このオプションでは出力はBED形式になります。- 途中のログを見ると、冗長なリードは除去され、multi readを除いた unique readsのみが出力されているようです。
- マッピング終了後にスタッツが表示されていますね（"Number of ..."の行）。入力リード数は21814483、マップされたリードは12508635、unique readsは11368866となっています。
Number of output mappings (passed filters): 9211499となっていますが、調べた感じではpassed filtersはクオリティ（q-value）が30を超えたリードのことを指すようです。最終的に出力されるリード数がこれになります。
chromapでは1m32sという実行時間になりました。１回しか計測していないのでざっくりですが、だいぶ速いと言えそうです。
比較のため、bowtie2とBWAでもマッピングしてみましょう。index作成は省略します。
$ time bowtie2 -x bowtie2-UCSC-hg38 -p 24 SRR227359.fastq.gz > SRR227359.bowtie2.sam
21814483 reads; of these:
21814483 (100.00%) were unpaired; of these:
8003629 (36.69%) aligned 0 times
9880708 (45.29%) aligned exactly 1 time
3930146 (18.02%) aligned >1 times
63.31% overall alignment rate
real 1m50.454s
user 36m56.341s
sys 5m40.437s
$ time bwa mem -t24 bwa-UCSC-hg38 SRR227359.fastq.gz > SRR227359.bwa.sam
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 4705884 sequences (240000084 bp)...
[M::process] read 4705884 sequences (240000084 bp)...
[M::mem_process_seqs] Processed 4705884 reads in 1092.230 CPU sec, 47.042 real sec
[M::process] read 4705884 sequences (240000084 bp)...
[M::mem_process_seqs] Processed 4705884 reads in 1236.997 CPU sec, 51.901 real sec
[M::process] read 4705884 sequences (240000084 bp)...
[M::mem_process_seqs] Processed 4705884 reads in 1199.590 CPU sec, 50.164 real sec
[M::process] read 2990947 sequences (152538297 bp)...
[M::mem_process_seqs] Processed 4705884 reads in 1202.052 CPU sec, 50.408 real sec
[M::mem_process_seqs] Processed 2990947 reads in 708.905 CPU sec, 29.680 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t24 /work/Database/bwa-indexes/UCSC-hg38 SRR227359.fastq.gz
[main] Real time: 243.498 sec; CPU: 5451.335 sec
real 4m3.885s
user 88m52.498s
sys 1m59.045s
Bowtie2は1m50s, BWAは4m3sになりました。Bowtie2も速いですね。BWAだけちょっと遅いですが、もともとこんなだったでしょうか？
ChIP-seqデータを入力にした場合、既存法と比べてそこまで劇的に高速化しているという感じではありません。 しかし内部的に冗長リードの除去までやってくれていることを考えると（通常はSAMtoolsやPicardを使って除去する必要があるため）、chromapは高速かつ簡便であると言えそうです。
出力をSAMにしてみる
chromap ではChIP-seqデータに対してBED形式の出力が推奨されています。--SAMオプションをつけることでSAM形式で出力可能ですが、低速になるため推奨しないとのことです（マニュアル参照）。
どのくらい遅くなるのかやってみましょう。
$ time chromap --preset chip -t 24 \ $ -x chromap-UCSC-hg38 -r hg38.fa \ $ -1 SRR227359.fastq.gz --SAM -o SRR227359.chromap.sam Preset parameters for ChIP-seq are used. Start to map reads. （中略） Mapped all reads in 62.61s. Number of reads: 21814483. Number of mapped reads: 12508635. Number of uniquely mapped reads: 11368866. Number of reads have multi-mappings: 1139769. Number of candidates: 227477151. Number of mappings: 12508635. Number of uni-mappings: 11368866. Number of multi-mappings: 1139769. Sorted, deduped and outputed mappings in 31.94s. # uni-mappings: 10859912, # multi-mappings: 1023343, total: 11883255. Number of output mappings (passed filters): 9236028 Total time: 109.61s. real 2m4.987s user 6m24.888s sys 0m48.419s
2m4sなので、ちょっと遅くなりましたが、十分許容範囲内かと思います。BEDで出力した時と微妙にスタッツが変化しているのが謎です。
chromapは結果を標準出力に出力することができないため、BAMやCRAMファイルにするためには一旦SAM形式で出力する必要があります。 また、SAMで遅くなる理由は主にI/O速度であると思われるので、リード数が多いほど書き込みのオーバーヘッド時間が増加すると想像されます。
その他のオプション: 冗長リードの除去をオフにしたい場合
--preset chip オプションに--remove-pcr-duplicatesオプションが含まれていますので、冗長リードのフィルタをオフにしたい場合は --preset chip オプションを使わないようにします。
$ time chromap -l 2000 --low-mem --BED -t 24 \ $ -x chromap-UCSC-hg38 -r hg38.fa \ $ -1 SRR227359.fastq.gz --SAM -o SRR227359.chromap.include_dup.sam Start to map reads. Parameters: error threshold: 8, min-num-seeds: 2, max-seed-frequency: 500,1000, max-num-best-mappings: 1, max-insert-size: 2000, MAPQ-threshold: 30, min-read-length: 30, bc-error-threshold: 1, bc-probability-threshold: 0.90 Number of threads: 24 Analyze bulk data. Won't try to remove adapters on 3'. Won't remove PCR duplicates after mapping. Will remove PCR duplicates at bulk level. Won't allocate multi-mappings after mapping. Only output unique mappings after mapping. Only output mappings of which barcodes are in whitelist. （中略） Mapped all reads in 74.75s. Number of reads: 21814483. Number of mapped reads: 12508635. Number of uniquely mapped reads: 11368866. Number of reads have multi-mappings: 1139769. Number of candidates: 227477151. Number of mappings: 12508635. Number of uni-mappings: 11368866. Number of multi-mappings: 1139769. Sorted and outputed mappings in 42.73s. # uni-mappings: 11368866, # multi-mappings: 1139769, total: 12508635. Number of output mappings (passed filters): 9685899 Total time: 144.08s. real 2m29.591s user 8m41.938s sys 1m0.581s
その他のオプション: クオリティによるフィルタリングをオフにしたい場合
chromapはリードのクオリティで自動的にフィルタが入るようになっていますが、（そういう需要はあまりないと思いますが）フィルタをオフにしたい場合は-q 0をつければOKです。
$ time chromap -l 2000 --low-mem -q 0 --BED -t 24 \ $ -x chromap-UCSC-hg38 -r hg38.fa \ $ -1 SRR227359.fastq.gz --SAM -o SRR227359.chromap.include_dup.q0.sam Mapped all reads in 78.28s. Number of reads: 21814483. Number of mapped reads: 12508635. Number of uniquely mapped reads: 11368866. Number of reads have multi-mappings: 1139769. Number of candidates: 227477151. Number of mappings: 12508635. Number of uni-mappings: 11368866. Number of multi-mappings: 1139769. Sorted and outputed mappings in 47.17s. # uni-mappings: 11368866, # multi-mappings: 1139769, total: 12508635. Number of output mappings (passed filters): 12508635 Total time: 188.37s. real 3m11.691s user 8m59.719s sys 1m8.940s
マップデータの可視化（100kbp bin）
chromapの論文では、ピーク領域についてはbowtie2とほぼ同一の結果が得られることが示されています。 ここでは得られたマップデータのリード分布が同一かどうかをおおまかに確かめます。 染色体レベルでの大まかなリード分布を調べるため、100kbp binでサンプル間比較をしてみましょう。
データをbigWigに変換するためにDROMPAplus*3に含まれるparse2wig+を使います。 bowtie2とBWAの出力ファイルには冗長リードが含まれていますが、parse2wig+は内部で冗長リードを除去するため、得られるBigWigはどのサンプルでも冗長リードは除去されています。 parse2wig+の使い方詳細はこちらを参照してください。
gt=genome_table.hg38.txt bin=100000 # 100kbp parse2wig+ -i SRR227359.chromap.bed -o SRR227359.chromap.bed --gt $gt --binsize $bin -n GR -f TAGALIGN parse2wig+ -i SRR227359.chromap.sam -o SRR227359.chromap.sam --gt $gt --binsize $bin -n GR parse2wig+ -i SRR227359.bowtie2.sam -o SRR227359.bowtie2.sam --gt $gt --binsize $bin -n GR parse2wig+ -i SRR227359.chromap.include_dup.sam -o SRR227359.chromap.include_dup.sam --gt $gt --binsize $bin -n GR parse2wig+ -i SRR227359.bwa.sam -o SRR227359.bwa.sam --gt $gt --binsize $bin -n GR
1行目ではBED形式のマップデータを読み込むために-f TAGALIGNオプションをつけています。
マップデータのクオリティスタッツにはほとんど違いはありませんでしたが、chromapが最もマップ率が低い、すなわち厳しい閾値を用いている感じでした（結果は割愛）。
なおBowtie2, BWAの出力ファイルにはunmappedも含めた全リードが含まれているのに対し、chromapの出力ファイルにはmapped readしか含まれていないようです。
DROMPA+のGVコマンドを使って可視化します。結果はGV.pdfに保存されます。
gt=genome_table.hg38.txt ideogram=DROMPAplus/data/ideogram/hg38.tsv drompa+ GV --gt $gt \ -i parse2wigdir+/SRR227359.chromap.bed.100000.bw,parse2wigdir+/SRR227359.bowtie2.sam.100000.bw,chromap.bed \ -i parse2wigdir+/SRR227359.chromap.sam.100000.bw,parse2wigdir+/SRR227359.bowtie2.sam.100000.bw,chromap.sam \ -i parse2wigdir+/SRR227359.chromap.include_dup.sam.100000.bw,parse2wigdir+/SRR227359.bowtie2.sam.100000.bw,chromap.include_dup.sam \ -i parse2wigdir+/SRR227359.bwa.sam.100000.bw,parse2wigdir+/SRR227359.bowtie2.sam.100000.bw,bwa.sam \ -o GV --showctag 1 --showitag 1 --ideogram $ideogram
結果は以下のようになりました。
緑のラインがChIP, 青のラインがInput、赤で示されているラインがChIP/Input ratioです。 ここではbowtie2のマップデータをInputにし、左に示されている４つのサンプルをChIPとしてratioを計算していますので、ratioはbowtie2のマップデータに対する増減を示しています。
結果を見ると、染色体腕部ではCNVによるリード数のばらつきが認められますが、ratioではきちんとキャンセルされており、ツール間でほぼ同数のリード数が得られていることがわかります。 一方セントロメア周辺部では、chromapのマップリード数が極端に少ないという結果になりました。 chromapとparse2wig+で冗長リードの定義が異なるせいかと思い、chromapで冗長リードを除去しないオプションで得たマップファイルでも試してみましたが（３段目）、やはりセントロメア周辺のリード数が少なくなっています。 セントロメア周辺はhighly repetitive regionsなので、chromapのキャッシュを用いたアラインメントはリピート領域に弱いか、あるいはchromapが内部的にリピート領域に対して厳しいスコアリングをしていると思われます。通常のChIP-seq解析ではリピート領域は対象にしませんのでどのツールを使っても変わらないと思いますが、H3K9me3のようにセントロメア周辺に特に濃縮するようなサンプルの場合は検討が必要そうです。なお、BWAは逆にリピート領域により多量に張り付くという結果になりました。
（2022/3/23追記）
ひょっとしてクオリティフィルタでリピート領域が除外されているのかと思い、-q 0オプションを追加して得られたマップデータについても同様に可視化してみましたが、若干リード量が回復したものの同じにはならず、やはりchromapではリピート領域のリード数は少ないままでした。やはりlight-weight alignmentはリピートに弱いのかもしれません。
まとめ
chromapについて、ChIP-seqデータを用いて簡単に試してみました。出力形式をSAMにしてもそこまで遅くはならないので、SAM出力からBAM/CRAMに変換するようにしておくと、他のツールとの互換性が高まると思います。
次回はHi-Cデータに対してchromapを試してみたいと思います。
*1:H. Zhang et al., Fast alignment and preprocessing of chromatin profiles with Chromap, Nature Communications, 2021, DOI 10.1038/s41467-021-26865-w
*2:H. Li, Minimap2: pairwise alignment for nucleotide sequences, Bioinformatics, 2018, DOI 10.1093/bioinformatics/bty191
*3:R. Nakato, T. Sakata, Methods for ChIP-seq analysis: A practical workflow and advanced applications, Methods, 2020, DOI 10.1016/j.ymeth.2020.03.005
