featureCount で bam と gtf からリードカウントを得る

featureCountは、bamからの高速なリードカウントツール。 Rやコマンドラインから実行できる。コマンドライン版は Subreadに含まれている。R版は、Rsubreadに実装されている。

featureについて

featureCount は feature ごとにリードをカウントするツールである。feture とは、リファレンス配列上のある区間のこと。その feature をまとめたセットをmeta-featureと呼んでいる。 RNA-seq の場合は、feature は exon であるし、meta-feature は gene や transcript になる。

コマンドライン版 Subread

インストール

wget http://downloads.sourceforge.net/project/subread/subread-1.5.0-p2/subread-1.5.0-p2-source.tar.gz
cd subread-1.5.0-p2-source/src
make -f Makefile.Linux   # 正常に終了すると ../bin に実行ファイルができる
cd ../../
cp -a subread-1.5.0-p2-source $HOME/opt/local
export PATH=$HOME/opt/local/subread-1.5.0-p2-source/bin:$PATH

featureCount でリードカウントを得る

まず遺伝子アノテーションファイルを得る。

wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_mouse/release_M9/gencode.vM9.annotation.gtf.gz
gzip -d gencode.vM9.annotation.gtf.gz

ローカルディレクトリに gencode.vM9.annotation.gtf ができる。

次に、featureCountを実行する。

featureCounts -T 8 -t exon -g gene_id -a gencode.vM9.annotation.gtf -o counts.txt hoge01.bam hoge02.bam

-TはCPUコア数。
-t は feature type を指定し、これにマッチした GTF file の行だけがカウントされる。デフォルトは exon になっている。
-a はアノテーションファイルを指定する。
-gは meta-feature を指定する。デフォルトは、gene_idになっている。つまり gene-levelのカウントが得られる。例の GENCODE の GTFの場合は、Ensembl Gene ID (= gene_id)になる。exon_idを指定すれば、exon-levelのカウントが得られる。
-oは出力されるカウントデータのファイル名。最後にカウントの対象となる bam file を指定する。1つでも複数でも指定できる。複数指定した場合は、1つのファイルが生成され、列ごとにそれぞれの bam のカウントが表示される。

R版 Rsubread

コマンドライン版だけ使うなら、このステップは必要ない。

インストール

Rを起動する。

sudo R

以下を入力。

library("BiocInstaller")
biocLite("Rsubread")

カウントしてみる。以下の例は、gencode.vM9.annotation.gtf に応じて、exonごとにリードカウントしている。

library("Rsubread")
library("edgeR")

gtf = "gencode.vM9.annotation.gtf"

# results/hisat2/*/*.sort.bam にある bam についてカウントする
bams = list.files(path = "results/hisat2/", pattern = "*.sort.bam$", recursive = TRUE)
bams = paste0("results/hisat2/", bams)

x = featureCounts(
  bams,
  annot.ext           = gtf,
  isGTFAnnotationFile = TRUE,
  GTF.featureType     = "exon",
  GTF.attrType        = "exon_id",
  useMetaFeatures     = FALSE,
  allowMultiOverlap   = TRUE,
  nthreads            = 8
)

# このあと edgeR で正規化、解析する場合は、DGEList オブジェクトにする
y = DGEList(counts = x$counts, genes = x$annotation)

2016/04/19