読者です 読者をやめる 読者になる 読者になる

アメリエフのブログ

Amelieff Staff Blog

coverageBedの使い方(1)

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

以上で準備は完了です。次回でカバレッジの計算を行います。