BEDtoolsの作者は開発熱心なので、できることがどんどん増えているような気がします。
手元のバージョンはv2.27.1です。
前準備
多くのコマンドはsorted BEDを要求しますので、事前に以下のコマンドで全てのBEDをソートしておくとストレスがないかと思います。
$ sort -k1,1 -k2,2n in.bed > in.sorted.bed
更に、いくつかのコマンドにはGenome tableファイルが必要になります。 Genome tableの作成は以下の記事を参照してください。
genome tableを作成する - Palmsonntagmorgen
a.bedに含まれる領域のうち、b.bedと重なる領域"のみ"を出力
$ bedtools intersect -a a.bed -b b.bed > a_and_b.bed
-bには複数のBEDファイルを指定することが可能です。
$ bedtools intersect -a a.bed -b b1.bed b2.bed b3.bed > a_and_b.bed
a.bedに含まれる領域のうち、b.bedと重なる領域を出力
$ bedtools intersect -wa -a a.bed -b b.bed > a_and_b.bed
intersectで出力される領域については、以下に図解があります。
intersect — bedtools 2.27.0 documentation
a.bedに含まれる領域のうち、b.bedと重ならない領域を出力
$ bedtools intersect -v -a a.bed -b b.bed > a_not_b.bed
a.bedの各領域がb.bedにどの程度カバーされるかを調べる
$ bedtools coverage -a a.bed -b b.bed > a.coverage.bed
出力は以下のようになります。
chr1 0 100 3 30 100 0.3000000 chr1 100 200 1 100 100 1.0000000 chr2 0 100 0 0 100 0.0000000 ...
１~3列目はa.bedと同一です。４列目はそのサイトと重なるb.bed内のサイトの数、５列目はサイト長に対するb.bedがoverlapしている領域の長さ、６列目はサイト長、７列目は割合 (5列目/6列目)です。
BEDに含まれる領域を両側に100bp伸ばす
$ bedtools slop -i in.bed -g genome_table -b 100 > in.extend.bed
片側のみに伸長したい場合は-l, -rオプションを使います。
BEDに含まれる領域のうちoverlapするものをマージ
$ bedtools merge -i in.bed > in.merged.bed
-dオプションでoverlapの基準をゆるめることができます。
BEDに含まれる領域を両側に1000bp伸ばし、重なったサイトはマージ
$ bedtools slop -i in.bed -g genome_table -b 1000 | bedtools merge -i - > in.extend.bed
3つのBEDファイルを結合し、重なるサイトはマージ
catするとsortされていないBEDが出力になるので、間にsortをかませています。
$ cat a.bed b.bed c.bed | sort -k1,1 -k2,2n | bedtools merge -i - > merged.bed
BEDに含まれて’いない’ゲノム領域を出力
$ bedtools complement -i in.bed -g genome_table > notin.bed
全ゲノム領域について、 BEDに含まれた領域にどの程度含まれているかをヒストグラムで出力
3つ以上のピークリストに対して、どの程度ピーク同士が重なっているかを調べる時などに使います。 BAMを入力にしてリードをカウントすることもできますが、遅いです。
$ bedtools genomecov -i in.bed -g genome_table > in.genomecov.bed
出力については以下の図解がわかりやすいです。
genomecov — bedtools 2.27.0 documentation
ゲノム内で100bpの領域を10000個ランダムに生成
ゲノムにはギャップ領域などがありますので、実際にはそれらを除外した領域を与えることが望ましいでしょう。
$ bedtools random -l 100 -n 10000 -g genome_table > random.bed
ランダムのため得られる結果は毎回変わりますが、-seedオプションを指定することで再現可能になります。
in.bedに含まれる領域をゲノム内でシャッフルして出力
randomではサイト長が固定ですが、こちらは実際に得られたピークデータを入力にするぶん、より実際に近いと思います。
Permutation testなどで使えそうですが、randomと同様、与えるゲノム領域を考慮する必要があります。
$ bedtools shuffle -i in.bed -g genome_table
シャッフルを同一染色体内に限定
$ bedtools shuffle -i in.bed -g genome_table -chrom
シャッフル後のサイトをinclude.bedに含まれる領域のみに限定
$ bedtools shuffle -i in.bed -g genome_table -incl include.bed
シャッフル後のサイトがexclude.bedに含まれる領域に含まれないようにする
$ bedtools shuffle -i in.bed -g genome_table -excl exclude.bed
シャッフル後のサイトが互いに重ならないようにする
$ bedtools shuffle -i in.bed -g genome_table -noOverlapping
BEDで指定されている領域の塩基配列を取得
モチーフ解析の時に必要になります。
$ bedtools getfasta -fi genome.fa -bed in.bed > in.fa
a.bedとb.bedの重なりをJaccard indexで計算
Jaccard indexとは、平たく言えばベン図の重なり部分の大きさを0~1の間のスコアで評価するものです。 便利そうで、あまり使う機会がありません。
$ bedtools jaccard -a a.bed -b b.bed
詳細は以下を参照してください。
jaccard — bedtools 2.27.0 documentation
複数のBAMファイルに対して、target_regions.bed に含まれる領域のマップリード数をカウント
$ bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed in.bed > in.count.bed
parallelコマンドがインストールされていれば、以下のコマンドでも可能です。
$ find *bam | parallel 'bedtools coverage -hist -abam {} -b target_regions.bed | grep ^all > {}.hist.all.txt'
これは以下のコマンドを実行したのと同じ意味になります。
$ bedtools coverage -hist -abam samp.01.bam -b target_regions.bed | grep ^all > samp.01.bam.hist.all.txt $ bedtools coverage -hist -abam samp.02.bam -b target_regions.bed | grep ^all > samp.02.bam.hist.all.txt $ bedtools coverage -hist -abam samp.03.bam -b target_regions.bed | grep ^all > samp.03.bam.hist.all.txt $ bedtools coverage -hist -abam samp.05.bam -b target_regions.bed | grep ^all > samp.05.bam.hist.all.txt $ bedtools coverage -hist -abam samp.06.bam -b target_regions.bed | grep ^all > samp.06.bam.hist.all.txt $ bedtools coverage -hist -abam samp.10.bam -b target_regions.bed | grep ^all > samp.10.bam.hist.all.txt
