以前BEDToolsのcoverageBedでdepthを各塩基ごとに求める方法を書いた。
$ coverageBed -a genome.bed -b sample.bed -d
というのが基本書式だ。
genome.bedはreference.fasta.faiから
$ awk '{print $1"\t0\t"$2}' reference.fasta.fai > genome.bed
という感じで作成しておく。
さらに、sample.bedはbamにそのまま置き換えてもちゃんと機能するらしいということで
$ coverageBed -a genome.bed -b sample.bam -d > depth.bed
としてやれば良さそうな感じ。
実際それでやってみると
1 0 30427671 1 0 1 0 30427671 2 0 1 0 30427671 3 0 1 0 30427671 4 0 1 0 30427671 5 0 1 0 30427671 6 0 1 0 30427671 7 0 1 0 30427671 8 0 1 0 30427671 9 0 1 0 30427671 10 0 ・・・・ 1 0 30427671 28106 1 1 0 30427671 28107 1 1 0 30427671 28108 1 1 0 30427671 28109 1 1 0 30427671 28110 1 1 0 30427671 28111 1 1 0 30427671 28112 1 1 0 30427671 28113 1 1 0 30427671 28114 1 1 0 30427671 28115 1 ・・・・
こんな感じのファイルが出力されてくる
いわゆるbed3にpositionとdepthが付加された感じかな
ちなみにこのbamファイルで最初に出てくるgenome positionが28106
これからdepth>0のpositionだけを抽出するなら
$ cat depth.bed | grep .'\t'0$ -v > depth1over.bed
などとすると
1 0 30427671 28106 1 1 0 30427671 28107 1 1 0 30427671 28108 1 1 0 30427671 28109 1 1 0 30427671 28110 1 1 0 30427671 28111 1 1 0 30427671 28112 1 1 0 30427671 28113 1 1 0 30427671 28114 1 1 0 30427671 28115 1 ・・・・
というふうに0カウントな行が取り除かれる。
さて、coverageBedについての情報をもうちょっと集めてみると、bamファイルから直接処理するなら
$ coverageBed -abam sample.bam -b genome.bed -d > depth.bed
とするようにと書かれていたりする。
これを実際にやってみると
1 28105 28191 NB501315:33:HGYFWBGXY:1:23206:2970:18726 60 +28105 28191 0,0,0 1 86, 0, 1 1 1 28105 28191 NB501315:33:HGYFWBGXY:1:23206:2970:18726 60 +28105 28191 0,0,0 1 86, 0, 2 1 1 28105 28191 NB501315:33:HGYFWBGXY:1:23206:2970:18726 60 +28105 28191 0,0,0 1 86, 0, 3 1 1 28105 28191 NB501315:33:HGYFWBGXY:1:23206:2970:18726 60 +28105 28191 0,0,0 1 86, 0, 4 1 1 28105 28191 NB501315:33:HGYFWBGXY:1:23206:2970:18726 60 +28105 28191 0,0,0 1 86, 0, 5 1 1 28105 28191 NB501315:33:HGYFWBGXY:1:23206:2970:18726 60 +28105 28191 0,0,0 1 86, 0, 6 1 1 28105 28191 NB501315:33:HGYFWBGXY:1:23206:2970:18726 60 +28105 28191 0,0,0 1 86, 0, 7 1 1 28105 28191 NB501315:33:HGYFWBGXY:1:23206:2970:18726 60 +28105 28191 0,0,0 1 86, 0, 8 1 1 28105 28191 NB501315:33:HGYFWBGXY:1:23206:2970:18726 60 +28105 28191 0,0,0 1 86, 0, 9 1 1 28105 28191 NB501315:33:HGYFWBGXY:1:23206:2970:18726 60 +28105 28191 0,0,0 1 86, 0, 10 1
こんな感じの結果が帰ってくる。
どうも、-aの方にsampleが来ることで、readごとにまとめてカウントすることになってしまうようで、これじゃちょっと違う気がする。bed12にpositionとdepthがくっつく形だな。
なお-abamを普通に-aのままbamを与えても、全く同じ結果になるので -abamというオプションをあえて使う意味はなくなっているのかもしれない。
bedtoolsの使い方をいろいろ調べ、試してみているうちに、なかなか良い使い方に出くわした。
$ bedtools intersect -a sample.bam -b background.bam -wa -v > true_candidate.bam
こういう使い方をすると、sample.bamからbackground.bamと重複する領域を持つリードを除いたbamを作ることができるようだ。これは使える。(とおもう)
https://bi.biopapyrus.jp/rnaseq/mapping/format-convert.html
http://kazumaxneo.hatenablog.com/entry/2017/07/23/213505
http://staffblog.amelieff.jp/entry/2012/12/13/112308
https://cell-innovation.nig.ac.jp/surfers/BEDtools.html
このあたり参照