２つのサンプルから得られたピークセットがどのくらい重なるのか調べたい！という時の方法です。
今回はBEDtoolsを使うやり方と、拙作のcompare_bsを使うやり方を紹介します。
ピークデータのダウンロード
ピークはBED形式であれば何でもよいのですが、ここでは例としてUCSCからH1細胞のCjunとMaxのピークファイルをダウンロードします。
(.narrowPeak という拡張子のファイルですが、これらはMACS2で得られたピークファイルで、BED形式として扱えます。)
$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak.gz $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak.gz
解凍後、ピーク数を確認してみましょう。
$ gunzip *.gz $ wc -l wgEncode* 2148 wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak 11129 wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak 13277 total
Cjunのピークはやや少ないようです。
bedtools
bedtoolsのintersectを使うと、共通するピーク領域を出力してくれます。その後数をカウントすればピーク数が判明します。
（Bedtoolsの詳細はこちら: BEDtoolsワンライナー覚書 - Palmsonntagmorgen）
$ bedtools intersect \
-a wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak \
-b wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak \
> and.bed
$ wc -l and.bed
377 and.bed
Cjunの2148ピークのうち、Maxと重なるものは377ピークになりました。
bedtools intersectは２つのピークが1塩基でも重なっている場合にoverlapと判定します。
また、サンプルAの１つのピークに大してサンプルBの３つのピークが重なる場合、出力されるピーク数は３つになります。コマンドに -wa オプションをつけると、1つのピークとして出力されます。
intersectで出力される領域については、以下に図解があります。
intersect — bedtools 2.29.2 documentation
compare_bs
compare_bsは以下の記事で説明したChIPseqTools の中にあります。
GitHubからプログラムをダウンロード・インストール - Palmsonntagmorgen
compare_bsを使うと以下のようになります。
$ compare_bs -and \ -1 wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak \ -2 wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak \ > and.bed $ head and.bed #file1: wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak #file2: wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak #num1: 2148 num2: 11129 num1_overlap: 373 (17.4%) num1_notoverlap: 1775 (82.6%) num2_overlap: 378 (3.4%)num2_notoverlap: 10751 (96.6%) #peakwidth total1: 523816 bp peakwidth total2: 3534874 bp overlappeaks total: 70413 bp (13.44% / 1.99%) chromosome start end summit length enrich intensity chr10 80917404 80917527 80917465 123 0.00 149.81 chr22 39570815 39571059 39570937 244 0.00 142.53 chr18 9643685 9643754 9643719 69 0.00 133.02 chr6 20811494 20811738 20811616 244 0.00 132.72 chr1 205263766 205264010 205263888 244 0.00 116.91
compare_bsの出力では、冒頭にピーク数のスタッツが表示されます。
この場合、Cjunのピーク数2148、Maxのピーク数11129で、CjunのうちMaxと重なるものが373ピーク、MaxのうちCjunと重なるものが378ピークとわかります。
(CjunとMaxでそれぞれ重なったピーク数がずれているのは、ピークの重なりが1対1対応でないことに起因します。つまり１つのピークに対して複数の相手方ピークが重なる可能性があるということです。)
その下の行の #peakwidth total は、ピークの数でなく延べピーク幅で重なりをカウントしたものです。
compare_bs はデフォルトでは「サンプル１のピークのうち、サンプル２と重ならないもの」を出力し、-and オプションをつけると重なるものを出力します。
ですので、出力されるピークは bedtools intersect の -waオプションと同じく、サンプル1のピーク幅になります。
重なりの定義もbedtools と同じく、1塩基でも重なればoverlapと判定します。
-l 10000 とオプションをつけると、10kbpまで離れたピークをoverlapとして判定します。
$ ~/git/ChIPseqTools/bin/compare_bs -and -l 10000 \ -1 wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak \ -2 wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak \ -nobs #file1: wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak #file2: wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak #num1: 2148 num2: 11129 num1_overlap: 689 (32.1%) num1_notoverlap: 1459 (67.9%) num2_overlap: 787 (7.1%)num2_notoverlap: 10342 (92.9%) #peakwidth total1: 523816 bp peakwidth total2: 3534874 bp overlappeaks total: 6727767 bp (1284.38% / 190.33%)
ここでは-nobs オプションをつけて、ピークリストを出力せずに冒頭のスタッツのみを表示しています。
10kbpまで重なりの閾値をゆるめたことで、少し重なりが増えました。
（注：overlapped peakwidthは10kbp extensionも考慮してoverlapを出すので、パーセンテージは少しおかしくなってしまいます）
まとめ
bedtools, compare_bsで得られる結果を図にまとめると以下のようになります。
得られるピークサイトはほぼ同じですが、compare_bsの優位点としては、statsを出してくれることとサンプル１と２のoverlap peak数をそれぞれ出してくれること、peak widthでもstatsを出してくれることがあります。
ピーク比較の時にはピークの幅は一定であると暗黙に仮定してしまいがちですが、幅の分布はかなり広がっているのが普通で、
１つの大きなピークが細切れになってカウントされていたりすることもあるので、ピーク数をカウントしてものを言いたい時には実際に得られたピークを確認する作業を忘れないようにしましょう。
