mecobalamin’s diary

人間万事塞翁が馬

RNA-seqその6、featureCountsでリードカウント

前回作ったbamファイルを使って発現量を調べる
mecobalamin.hatenablog.com

使うコマンドはfeatureCount

こちらのサイトを参考にした
featureCount で bam と gtf からリードカウントを得る | Tips for NGS Data Analysis
https://ncrna.jp/blog/item/387-featurecount-2018

featureCountで検索かけるとRを使った解析ソフトウェア群の
Bioconductorが見つかる
Bioconductor - Home
BiocManagerというパッケージマネージャーもある。
(biocLiteはバージョンの古いマネージャー)

featureCountはRsubreadに含まれるが
こちらのwslでは動かなかった

なのでコマンドライン版のfeatureCountを使う
The Subread package
ここからダウンロードできる。
最新の1.6.3と、1.5.0-p2はこちらの環境では動かなかった
1.6.3は実行後エラーが出てすぐに止まった
1.5.0-p2はエラーメッセージが表示され続けた

1.5.3は動いたのでこのバージョンを使う
ファイルを展開後makeして$HOME/binにリンクを張った
WEHI Bioinformatics - Subread package

31 May 2019追記
最新版1.6.4が公開済みだった
アーカイブをダウンロードして展開
バイナリが含まれていたので$HOME/binにリンクを張った
ちゃんと動いた
追記ここまで

gtfファイルとbamファイルを読み込ませて実行する
結果はcounts.txtに出力される

$ featureCounts -p -O -T 4 -t exon -g gene_id -a ref_GRCh38.p12_top_level.gtf -o counts.txt hogehoge.bam

6 March 2019追記

  • pと-Oオプションを追加
  • pはペアエンド、-Oは重なりがあるリードがあるときに指定するらしい。

Assign reads to all their overlapping meta-features (or features if -f is specified).
featureCounts | 各遺伝子にマッピングされたリード数を計数
WEHI Bioinformatics - featureCounts
追記ここまで

このcounts.txtをRに読み込ませて解析するが、
とりあえずヒストグラムを作る

ここからはRで作業する
counts.txtのあるディレクトリに作業ディレクトリを移動し、
read.table()でファイルを読み込ませる
一行目のコメントとヘッダーを飛ばして、
タブ区切りにする

どれがカウントデータかよくわからないけど
colnames()で確認すると

[1] "Geneid"              "Chr"                 "Start"               "End"                 "Strand"             
[6] "Length"              "hogehoge.bam"

なのでおそらく最後のカラム

試しにコピー数が3から30のだけを取り出して
ヒストグラムにする

setwd("/path/to/counts.txt")
fc <- read.table("counts.txt", header = TRUE, skip = 1, sep="\t")
d <- subset(fc, fc$hogehoge.bam < 30)
d <- subset(d, d$hogehoge.bam > 3)
h <- hist(d$hogehoge.bam)

結果はこんな感じ
f:id:mecobalamin:20190310000502p:plain

ちなみにpng形式で保存するには
Rに表示されているグラフを消してから
png()のコマンドでファイル名とサイズを宣言
hist()を実行、そしてdev.off()でファイルを閉じる

png("histogram.png", width = 600, height = 450)
h <- hist(d$hogehoge.bam, main = "", xlab= "hogehoge.bam")
dev.off()

10 March 2019修正
dev.off()が抜けてたので追加
必要ならpng()の前にもdev.off()を実行して

null device 
          1

が表示されてから上記のコマンドを実行する

hにはhistogramのパラメータが入っている
$breaksはhist()の引数breaksを指定することで
変更可能でbin幅が変えられる

> h
$breaks
 [1]  4  6  8 10 12 14 16 18 20 22 24 26 28 30

$counts
 [1] 1114  358  250  196  154  134   87  104   93   64   54   51   34

$density
 [1] 0.206832529 0.066468622 0.046416636 0.036390642 0.028592648 0.024879317 0.016152989
 [8] 0.019309320 0.017266988 0.011882659 0.010025993 0.009468994 0.006312662

$mids
 [1]  5  7  9 11 13 15 17 19 21 23 25 27 29

$xname
[1] "d$hogehoge.bam"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

次は目的の遺伝子がいくつあるか
グラフ化したい

11 March 2019追記
HeLaのbamファイルとhg38のgtfファイルを使ってfeatureCountsを試した
hg38のfastaとgtfはリンク先を参照
mecobalamin.hatenablog.com
HeLaのsraファイルはSRR6799791で、
hg38のシークエンスデータにhisat2を使ってマッピングした
できたファイルがHeLa.bam

featureCountsを実行する

$ featureCounts -p -B -T 4 -O -t exon -g gene_id -a rename_UCSC.hg38.gtf -o counts.txt HeLa.bam

できたcounts.txtから試しにEGFRの数を抜き出してみる
https://www.ncbi.nlm.nih.gov/gene/1956
上皮成長因子受容体 - Wikipedia

EGFRにはスプライシングバリアントがある
それぞれのaccession numberでカウントをcounts.txtから抜き出す
Rを使う

acc.egfr <- c("NM_005228", "NM_201282", "NM_201283", "NM_201284", "NM_001346897", "NM_001346898", "NM_001346899", "NM_001346900", "NM_001346941")
iso.egfr <- c("a", "b", "c", "d", "e", "f", "g", "h", "i")
var.egfr <- t(rbind(acc.egfr, iso.egfr))

fc <- read.table("counts.txt", header = TRUE, skip = 1, sep="\t")

egfr <- c()
for(i in var.egfr[,1]){egfr <- rbind(egfr, subset(fc, grepl(i, fc$Geneid)))}
cbind(iso.egfr, egfr[c(1,7)])

結果はこうなったけど
これをそのままmRNAのコピー数としていいのか。。。?

       iso.egfr         Geneid            HeLa.bam
134652        a    NM_005228.4                 981
134654        b    NM_201282.1                 254
134653        c    NM_201283.1                 190
134655        d    NM_201284.1                 255
134648        e NM_001346897.1                 494
134649        f NM_001346898.1                 494
134650        g NM_001346899.1                 981
134657        h NM_001346900.1                 963
134651        i NM_001346941.1                 945