2サンプル間ピーク比較

２つのサンプルから得られたピークセットがどのくらい重なるのか調べたい！という時の方法です。

今回は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で得られる結果を図にまとめると以下のようになります。

f:id:rnakato:20180806141533p:plain

得られるピークサイトはほぼ同じですが、compare_bsの優位点としては、statsを出してくれることとサンプル１と２のoverlap peak数をそれぞれ出してくれること、peak widthでもstatsを出してくれることがあります。
ピーク比較の時にはピークの幅は一定であると暗黙に仮定してしまいがちですが、幅の分布はかなり広がっているのが普通で、１つの大きなピークが細切れになってカウントされていたりすることもあるので、ピーク数をカウントしてものを言いたい時には実際に得られたピークを確認する作業を忘れないようにしましょう。