シーケンスリードのカバレッジ(被覆率)を計算する

bedtools で計算する方法 (bam を bed形式に変換する)

bamからカバレッジを計算し、bedファイルを作成する。bedはテキストファイルなので、bamのサイズ(リード数)によっては、かなり巨大なファイルができる。遺伝子単位や限られたゲノム領域のみについて計算し、 自作のプログラムや表計算ソフトで、処理したいときになどに利用できる。全ゲノムレベルのカバレッジを計算し たい場合は、bigwigへの変換をお勧めする。

インストール

$ curl -O //github.com/arq5x/bedtools2/releases/download/v2.25.0/bedtools-2.25.0.tar.gz
$ tar zxvf bedtools-2.25.0.tar.gz
$ cd bedtools2/
$ make
$ cd ..
$ cp -a bedtools2 $HOME/opt/local

設定

$ emacs -nw ~/.zshenv
export PATH=$HOME/opt/local/bedtools/bin:$PATH
$ source ~/.zshenv

ここでは、bed ファイル (Pou5f1.bed) に記述された座標にマップされたリードのカバレッジを bam ファイルから(hoge.bam) 計算する。

$ coverageBed -d -abam hoge.bam -b Pou5f1.bed > Pou5f1.TruRNA-Seq01.coverage.bed

RSeQCで計算する方法 (bam を wig に変換する)

wig file はカバレッジのデータを保存するためによく使われるファイル形式である。ここでは、bam から wig を計算する方法を述べる。RNA-seqデータのQuality checkのツール集RSeQCに 含まれる、bam2wig.pyを利用する。

bam2wig.py -t 1000000000 -i hoge.bam -o hoge.wig -s hg38.chrom.sizes

wigToBigWig がインストールされており、パスが通っている場合は、自動的に bigwig を作ってくれる。

オプション-sは、染色体サイズの情報が含まれたテキストファイルを指定する。 ファイルは、UCSC Genome Browser からダウンロードできる。 http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes

また以下のコマンドで生成することもできる。

まず、コマンドのセットアップをする。wget か MySQLクライアントがインストールされている必要がある。

wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/fetchChromSizes # インストール
chmod +x fetchChromSizes

実行する。

fetchChromSizes mm10 > mm10.info

このコマンドは、http://hgdownload.cse.ucsc.edu/goldenPath/${DB}/bigZips/${DB}.chrom.sizes から wget コマンドでファイルをダウンロードするだけのもの。MySQLクライアントがインストールされている環境では、 genome-mysql.cse.ucsc.edu の MySQLサーバの ${DB} にある chromInfo テーブルの結果を表示するだけ。

wig を bigwig へ変換する

wig 形式はテキスト形式になっているが、このままではファイルサイズが巨大になる。そこで、バイナリ形式である、bigwig へ変換する。

wget //hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/wigToBigWig  # コマンドをインストール
wigToBigWig hoge.wig mm10.info hoge.bw

2016/05/16: bedtools-2.25.0へ修正