bedtoolsによる異なるChIP-Seqデータのタグ数の比較

  • 2
    Like
  • 2
    Comment
More than 1 year has passed since last update.

目的

同一または異なる抗体による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を元に2つのデータフレームをマージする
  • 散布図を作成する
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