目的
同一または異なる抗体による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