BAMから、depth10以上の領域を示すbedを作る方法【bedtoolsの使い方】

Fukuです。こんにちは。
NGS解析 で浮かんだ疑問を解決するコマンドを紹介します。

今日の目的

Whole Genome Sequenceや Whole Exome Sequenceにおいて
NGSシーケンサのアラインメントファイル（BAM, SAM）から、「Depthが10以上のゲノム領域」を示す領域ファイル（BED）を作成します。

BAM

↓

chr1    65838   65942
chr1    65948   65963
chr1    69470   69602
...

これを例えば、「Depthが10未満の領域は下流解析から除く」という目的に利用できます。（この目的なら、ほかの方法もありますが）

方法

bedtools を使います。 bedtools を使うと、さまざまな ゲノム配列に関する計算をすることができます。
BAM、BED、GFF / GTF、VCFなどの広く使用されているファイル形式の複数のファイルから、ゲノム領域を相互に交差、マージ、カウントなどをすることができます。
bedtools 公式ドキュメント

コマンド

bedtools genomecov -ibam example_sorted.bam -bg |  awk '$4 > 9 {print $0}'  | bedtools merge  > depth10.bed

• -ibam input.bam： 入力BAMファイルを指定
• -bg ：depthを1塩基単位ではなく、bedgraph形式で出力させる
• awkの中の >9 ： depthが9より大きい = 10以上 のデータだけを抽出すること
• awkの中の{print $0} ：上記条件に合う行の全体を出力する
• > リダイレクト。標準出力をファイルに書き出す。

3ステップの工程をパイプで繋いでいます。
以下、一つ一つ解説します。

① bedtools genomecov で各ポジションのdepthを計算する

参考：bedtools genomecov公式ドキュメント

$ bedtools genomecov -ibam example_sorted.bam -bg |  head                            
chr1    10020   10039   1
chr1    10039   10093   2
chr1    10093   10110   1
chr1    10829   10848   1
chr1    10904   10924   1
chr1    11194   11248   1
chr1    11248   11265   2
chr1    11265   11269   3
chr1    11269   11294   4
chr1    11294   11349   3

結果の１行目を見ると、chr1:10020-10039の領域はdepthが1 だと示されています。

② depthが◯より大きい領域を抽出する
(1) の結果から、depth すなわち4列目を基準にしてフィルタリングします。 linuxコマンド awk を使います。
ここでは 9より大きい、条件を設定していますが、数値を変えたい場合は 9を書き換えればOKです。

$ bedtools genomecov -ibam example_sorted.bam -bg |  awk '$4 > 9 {print $0}' | head  
chr1    65838   65841   10
chr1    65841   65863   11
chr1    65863   65865   12
chr1    65865   65870   13
chr1    65870   65875   14
chr1    65875   65876   13
chr1    65876   65877   12
chr1    65877   65884   13
chr1    65884   65886   14
chr1    65886   65890   13

たしかにdepthが10以上の行だけが抜きだされています。

③ bedtools merge で 接する/重なる領域をマージする

(2)の結果から、別々の行になっているけれど実は隣り合っている領域をマージします。

参考：bedtools merge公式ドキュメント

$ bedtools genomecov -ibam SRR10338722_sorted_rmdup.bam -bg |  awk '$4 > 9 {print $0}' | bedtools merge | head                
chr1    65838   65942
chr1    65948   65963
chr1    69470   69602
chr1    127231  127250
chr1    128343  128351
chr1    484311  484334
chr1    484949  484972
chr1    718209  718305
chr1    719890  719891
chr1    719905  720018

chr1:65838〜65942 は、(2)では10行以上にわたって分割されていましたが、bedtools mergeで1行にまとまりました。

bedtools のインストール方法

ダウンロード、解凍、コンパイルをします。

$ wget https://github.com/arq5x/bedtools2/releases/download/v2.29.1/bedtools-2.29.1.tar.gz
$ tar -zxvf bedtools-2.29.1.tar.gz
$ cd bedtools2
$ make

必要に応じて、パスを通します。
（$PATHが通っているところにシンボリックリンクを作成します。）

$ ln -s  /path/to/bedtools2/bedtools $HOME/bin/

おわりに

いかがでしたか？
ツールを使いこなして、やりたいデータ分析を実現していきましょう！