アメリエフの技術ブログ

Amelieff Staff Blog

coverageBedの使い方(2)

coverageBedの使い方(1)のつづきです。

■カバレッジ計算
マッピング結果とゲノム、2つのBEDが用意できたら、以下のコマンドでカバレッジを計算します。

$ coverageBed -a map.bed -b genome.bed

以下のような結果が出力されます。

chr1___0__249250621__1__100__249250621__0.0000004
chr10__0__135534747__0____0__135534747__0.0000000
__:
chr2___0__243199373__1__100__243199373__0.0000004
chr20__0___63025520__0____0___63025520__0.0000000
__:

出力項目は、左から
1列目「genome.bedの1列目」
2列目「genome.bedの2列目」
3列目「genome.bedの3列目」
4列目「genome.bedにオーバーラップしたmap.bedの領域数」
5列目「カバレッジが0より大きいgenome.bedの塩基数」
6列目「genome.bedの領域長(=3列目と同じ値)」
7列目「map.bedによりオーバーラップされたgenome.bedの割合」
です。

7列目を見れば、染色体ごとのマッピング率がわかります。

ゲノム全体のマッピング率は、以下のコマンドで計算できます。
(今回のテストデータでは、6.46059e-08になりました)

$ coverageBed -a map.bed -b genome.bed | awk '{a+=$3;b+=$5} END {print b/a}' -

-bで指定するBEDを特定の領域(Exomeのキャプチャ領域など)にすれば、その領域に対するカバレッジを計算できますし、
以下のようなゲノム領域を10kb刻みで区切ったBEDにすれば、10kbごとのカバレッジを計算することができます。

chr1____0__100
chr1__100__200
chr1__200__300
__:

また、-aの代わりに-abamを使うと、BAMファイルを指定することもできます。
$ coverageBed -abam map.bam -b genome.bed