kuroの覚え書き

96の個人的覚え書き

塩基単位でのread depthを求める

以前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

このあたり参照