目的

同一または異なる抗体によるChIP-Seqデータを比較する手段の一つとして、片方のデータで検出されたピーク領域をqueryとして、query領域周辺にマップされたリード数を比較する方法がある。
ここでは、bedtoolsを用いてピーク領域におけるリード数をカウントし、Rでカウントする方法を説明する。

方法

0. 用意するもの

ピーク領域を定義するファイル（BED）
リードのアラインメントファイル（BAM or BED）
bedtools のインストール
Rのインストール

1. ピーク領域におけるリード数をカウントする

ピーク領域（query）: peak_0.bed
アラインメントファイル: mapped_0.bam, mapped_1.bam
queryで定義された範囲から±100 bp の範囲のリード数をカウントする（awk コマンドの -100と+100）

awk 'BEGIN{OFS="\t"}{print $1, $2 - 100, $3 + 100}' peak_0.bed | bedtools intersect -c -a stdin -b mapped_0.bam -bed > result_query_0.txt
awk 'BEGIN{OFS="\t"}{print $1, $2 - 100, $3 + 100}' peak_0.bed | bedtools intersect -c -a stdin -b mapped_1.bam -bed > result_query_1.txt

これにより、それぞれのBAMにおいて各ピーク領域に重なるリード数がカウントされる

2. ピーク領域におけるリード数を散布図を描いて比較する

ピーク座標をIDとする
IDを元に２つのデータフレームをマージする
散布図を作成する

d0 = read.table("result_query_0.txt")
paste(d0[,1], d0[,2], sep = ":")
colnames(d0)[ncol(d0)] = "tag0"

d1 = read.table("result_query_1.txt")
d1 = cbind(id = paste(d1[,1], d1[,2], sep = ":"), d1)
colnames(d1)[ncol(d1)] = "tag1"

d = merge(d0, d1, all = T)

plot(d$tag0, d$tag1)

資料

BED format: http://bedtools.readthedocs.org/en/latest/content/general-usage.html
bedtools: http://bedtools.readthedocs.org/en/latest/index.html
Rのmerge関数: http://cse.naro.affrc.go.jp/takezawa/r-tips/r/43.html