chromapを試す（その２）

高速にマッピングできるchromap*1 の続きです。今回はHi-Cデータに対してマッピングし、BWAと比較します。

前回の記事はこちら

rnakato.hatenablog.jp

Hi-Cファイルの取得

今回はHaarhuisらのCell論文*2 のデータを使いました。ヒトHap1細胞でのHi-Cデータです。

fastq-dump --split-files --gzip SRR5266584

BWAでマッピング

BWAを使ってfastqデータをマッピングします。マッピングのパラメータは 4D Nucleome project *3 を参考にしています。用いるCPU数は32です。

$ time bwa mem -t 32 -SP5M bwa-indexes/UCSC-hg38 \
   SRR5266584_1.fastq.gz SRR5266584_2.fastq.gz \
   > mapped.bwa.sort.bam
... (中略)
[M::mem_pestat] analyzing insert size distribution for orientation RR...
[M::mem_pestat] (25, 50, 75) percentile: (1778, 3622, 6305)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 15359)
[M::mem_pestat] mean and std.dev: (4152.72, 2755.39)
[M::mem_pestat] low and high boundaries for proper pairs: (1, 19886)
[M::mem_process_seqs] Processed 604710 reads in 302.831 CPU sec, 9.743 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t 32 -SP5M bwa-indexes/UCSC-hg38 SRR5266584_1.fastq.gz SRR5266584_2.fastq.gz
[main] Real time: 1874.166 sec; CPU: 56598.868 sec

________________________________________________________
Executed in   31.24 mins   fish           external
   usr time  928.66 mins  595.00 micros  928.66 mins
   sys time   14.65 mins  231.00 micros   14.65 mins

32コアで31分程度で終わりました。実際、マッピングのみであればHi-Cであってもそこまで遅くはありません。

実際のHi-C解析では、以下のようにpairtools を用いてその後のマップリードのソートや冗長リードの除去などに更に時間が必要になります。コマンドの詳細は割愛しますが、ここでは72分かかっています。31分 + 72分 = 103分を使っています。

$ time pairtools parse --min-mapq 30 \
   --walks-policy 5unique --max-inter-align-gap 30 \
   --nproc-in 8 --nproc-out 8 --chroms-path genometable.hg38.txt mapped.bwa.bam \
 | pairtools sort --nproc 8 \
 | pairtools dedup --nproc-in 8 --nproc-out 8 --mark-dups \
 > mapped.bwa.pairsam

real    72m54.795s
user    79m7.940s
sys     17m2.803s

chromapでマッピング

chromapでマッピングしてみます。Hi-Cデータのため、--preset hicオプションをつけています。これにより、出力形式が4DNで用いられているpairs形式になります。また、--remove-pcr-duplicatesオプションで冗長リードの除去ステップを追加しています。

なお、chromapではマッピングの時にもゲノムデータ (hg38.fa)を指定する必要があります。

$ time chromap --preset hic --remove-pcr-duplicates \
   -t 32 -x chromap-indexes/UCSC-hg38 -r hg38.fa \
   -1 SRR5266584_1.fastq.gz -2 SRR5266584_2.fastq.gz \
   -o mapped.chromap.pairs
Preset parameters for Hi-C are used.
Start to map reads.
Parameters: error threshold: 4, min-num-seeds: 2, max-seed-frequency: 500,1000, max-num-best-mappings: 1, max-insert-size: 1000, MAPQ-threshold: 1, min-read-length: 30, bc-error-threshold: 1, bc-probability-threshold: 0.90
Number of threads: 32
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.
Allow split alignment.
Output mappings in pairs format.
Reference file: hg38.fa
Index file: chromap-indexes/UCSC-hg38
1th read 1 file: SRR5266584_1.fastq.gz
1th read 2 file: SRR5266584_2.fastq.gz
Output file: mapped.chromap.pairs
Loaded all sequences successfully in 28.22s, 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 112.23s.
Mapped 500000 read pairs in 5.55s.
（中略）
Mapped 500000 read pairs in 1.56s.
Mapped 100349 read pairs in 0.56s.
Mapped all reads in 304.95s.
Number of reads: 120200698.
Number of mapped reads: 108983560.
Number of uniquely mapped reads: 100363412.
Number of reads have multi-mappings: 8620148.
Number of candidates: 4829748453.
Number of mappings: 108983560.
Number of uni-mappings: 100363412.
Number of multi-mappings: 8620148.
Sorted, deduped and outputed mappings in 57.58s.
# uni-mappings: 53267577, # multi-mappings: 868789, total: 54136366.
Number of output mappings (passed filters): 47034561
Total time: 508.85s.

________________________________________________________
Executed in  511.87 secs   fish           external
   usr time  100.27 mins  475.00 micros  100.27 mins
   sys time    1.93 mins  234.00 micros    1.93 mins

なんと、32コアで8分あまりでマッピング終了しました。BWAでは冗長リードの除去が別途pairtools で必要になることを考えると、これは相当速いです。

ただしchromapは制限酵素サイトを考慮するオプションを持っていないので、フラグメントのフィルタリングに別途pairtoolsなどを利用する必要があります*4。手元で試したところでは、chromapの出力をpairtoolsに入れるとエラーになるため、更なる調査が必要そうです。Micro-Cの場合は制限酵素サイトを考慮する必要がないので、この出力ファイルをそのまま利用することが可能です。