RNA-seqその６、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)
結果はこんな感じ
ちなみに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
