以前、bedtoolsの一つであるintersectBedの使い方についてご紹介しましたが、今回はcoverageBedについてご紹介します。
マッピング結果がどのくらいゲノム全体をカバーできているか知りたい時、coverageBedを使うとカバレッジが簡単に計算できます。
■用意するもの
・マッピング結果のBEDファイル
__--マッピング領域の染色体名・スタート・エンド
BEDは1行に最低3列、オプションで12列まで書くことができるフォーマットですが、coverageBedでは3列のBEDしか受け付けないようです。
列数が3より多い場合は、以下のコマンドで列数を3列だけにしておきます。
$ awk '{print $1"¥t"$2"¥t"$3}' original_map.bed > map.bed
今回、マッピング結果のBEDファイルとして以下のようなファイルを作成しました。
chr1__0__100
chr2__0__100
・ゲノムのBEDファイル
__--各染色体の染色体名・スタート・エンド
ゲノムのBEDファイルは、ゲノムのfastaファイル（ここではgenome.faという名前のものを使いました）から以下のコマンドで作ることができます。
染色体の長さが分かれば手作業で作ってもよいです。
$ samtools faidx genome.fa
__→genome.fa.faiというファイルができます
$ awk '{print $1"¥t0¥t"$2}' genome.fa.fai > genome.bed
__→以下のような内容のgenome.bedができます
chr1__0__249250621
chr2__0__243199373
__:
chrY__0__59373566
chrM__0__16571
以上で準備は完了です。次回でカバレッジの計算を行います。
