アメリエフの技術ブログ

Amelieff Staff Blog

intersectBedの使い方

bedtoolsBEDフォーマットのファイルを扱うのに便利なツール群です。

今回はその中の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の位置関係は次のようになります。

f:id:amelieff:20190731121314p:plain




-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でファイルを指定してください)。
マッピング結果への簡単なアノテーション付けくらいならこれだけで行えてしまうくらいなので、まだ使ったことがない方はぜひお試しください。