bedtoolsはBEDフォーマットのファイルを扱うのに便利なツール群です。
今回はその中の1つ、intersectBedについてご紹介します。
intersectBedを使うと、複数BED間で重複している領域を簡単に抽出することができます。
テストデータとして、2つのBEDファイルを用意しました。
==> A.bed <==
chr1 10 30 A_1
chr1 40 50 A_2
chr1 70 80 A_3
==> B.bed <==
chr1 10 15 B_1
chr1 20 25 B_2
chr1 30 40 B_3
chr1 60 90 B_4
A.bedとB.bedの位置関係は次のようになります。
◆-aで指定したBED中の要素のうち、-bで指定したBED中の要素にオーバーラップするものだけを抽出するには、-waオプションをつけます。
$ intersectBed -a A.bed -b B.bed -wa
chr1 10 30 A_1
chr1 10 30 A_1
chr1 70 80 A_3
→A_1とA_3だけが出力されました。
A_1が2回出力されるのは、2回オーバーラップしているためです。
◆-uオプションをつけると、重複した結果がマージされます。
$ intersectBed -a A.bed -b B.bed -wa -u
chr1 10 30 A_1
chr1 70 80 A_3
→A_1の出力が1回だけになりました。
◆-waの代わりに-wbを使うと、bの情報も出力されます。
$ intersectBed -a A.bed -b B.bed -wb
chr1 10 15 A_1 chr1 10 15 B_1
chr1 20 25 A_1 chr1 20 25 B_2
chr1 70 80 A_3 chr1 60 90 B_4
→左4列は「a中のオーバーラップ領域」、
右4列は「オーバーラップしたbの情報」です。
◆-fオプションを使うと、オーバーラップの割合を指定できます。
例えば、全長の3割以上がbにオーバーラップされるものだけを抽出するには、-f 0.3と指定します。
$ intersectBed -a A.bed -b B.bed -wa -f 0.3
chr1 70 80 A_3
→A_3だけが出力されました。
A_1（20bp）は、B_1によって5bp、B_2によって5bpで計10bpオーバーラップしているのですが、出力されません。おそらく積算ではなく、1要素によるオーバーラップ領域長で見ているのだと思われます。
◆-vオプションをつけると、オーバーラップ「しなかった」ものが出力されます。
$ intersectBed -a A.bed -b B.bed -wa -v
chr1 40 50 A_2
◆-cオプションをつけると、オーバーラップしたbの数も出力されます。
$ intersectBed -a A.bed -b B.bed -wa -c
chr1 10 30 A_1 2
chr1 40 50 A_2 0
chr1 70 80 A_3 1
◆-waoオプションをつけると、オーバーラップしたbの情報と、オーバーラップした塩基数が出力されます。
$ intersectBed -a A.bed -b B.bed -wao
chr1 10 30 A_1 chr1 10 15 B_1 5
chr1 10 30 A_1 chr1 20 25 B_2 5
chr1 40 50 A_2 . -1 -1 . 0
chr1 70 80 A_3 chr1 60 90 B_4 10
以上、A.bedを主体にした結果を出しましたが、Bを主体にした結果を出したい場合は、-a B.bed -b A.bedのように、-aと-bの指定をひっくり返せばOKです。
また、今回のテストBEDにはストランド情報がないので実行していませんが、-sオプションで「同じ向きに限定」、-Sで「逆の向きに限定」することができます。
すごいのはこのintersectBed、BAM・VCF・GFFにも使えることです（※BAMの場合は-aの代わりに-abamでファイルを指定してください）。
マッピング結果への簡単なアノテーション付けくらいならこれだけで行えてしまうくらいなので、まだ使ったことがない方はぜひお試しください。
