Palmsonntagmorgen

NGSデータを使った解析と、その周辺。

BEDtoolsワンライナー覚書

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 d1.bed d2.bed d3.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
...

1~3列目はa.bedと同一です。4列目はそのサイトと重なるb.bed内のサイトの数、5列目はサイト長に対するb.bedがoverlapしている領域の長さ、6列目はサイト長、7列目は割合 (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