前回の記事の続きです。
HISAT-StringTie-Ballgown を試してみよう - Palmsonntagmorgen
発現量データ生成
前回の記事で発現量データを生成するのを忘れていたので、ここで生成します。stringtieを使います。
for prefix in ERR188044 ERR188104 ERR188234 ERR188245 ERR188257 ERR188273 ERR188337 ERR188383 ERR188401 ERR188428 ERR188454 ERR204916
do
stringtie -e -B -p 16 -G stringtie_merged.gtf \
-o ballgown/${prefix}/${prefix}_chrX.gtf \
${prefix}_chrX.bam
done
これでballgown/${prefix} というディレクトリが自動生成され、その中に発現量データが出力されます。
ディレクトリの中身を見てみましょう。
$ ls ballgown/ERR188044/
ERR188044_chrX.gtf e2t.ctab e_data.ctab i2t.ctab i_data.ctab t_data.ctab
gtfファイルと、いくつかのテキストファイルができています。
eはexon, iはid, tがtranscriptの情報を表しているようですね。
比較するサンプルの確認
サンプルの情報はchrX_data/geuvadis_phenodata.csv
に記述されています。
どんなサンプルなのか確認しましょう。
$ cat chrX_data/geuvadis_phenodata.csv
"ids","sex","population"
"ERR188044","male","YRI"
"ERR188104","male","YRI"
"ERR188234","female","YRI"
"ERR188245","female","GBR"
"ERR188257","male","GBR"
"ERR188273","female","YRI"
"ERR188337","female","GBR"
"ERR188383","male","GBR"
"ERR188401","male","GBR"
"ERR188428","female","GBR"
"ERR188454","male","YRI"
"ERR204916","female","YRI"
サンプルIDごとに、性別とpopulation情報が載っています。
解析内容は、population間のばらつきを除去しながら性別間の遺伝子発現変動を見る、というもののようです。
Ballgownで解析
BallgownはRで動きますので、同じディレクトリでRを起動しましょう。
以降のコマンドはR上で実行するものとします。Nature Protocolのコマンドから少し修正しています。
library(ballgown)
library(RSkittleBrewer)
library(genefilter)
library(dplyr)
library(devtools)
pheno_data = read.csv("chrX_data/geuvadis_phenodata.csv")
bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
bg_chrX
ballgown instance with 3520 transcripts and 12 samples
bg_chrX
という名前のballgownインスタンスが生成されました。
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)
results_transcripts = stattest(bg_chrX_filt,
feature="transcript",
covariate="sex",
adjustvars = c("population"),
getFC=TRUE,
meas="FPKM")
head(results_transcripts)
feature id fc pval qval
1 transcript 1 0.9347423 0.7043147 0.9373013
2 transcript 2 1.2006990 0.8716451 0.9723579
3 transcript 3 1.0173174 0.9896415 0.9990549
4 transcript 4 0.3814504 0.5200179 0.9192672
5 transcript 5 0.6014996 0.3216161 0.9192672
6 transcript 6 0.6540123 0.3264784 0.9192672
同様に、遺伝子単位の発現変動情報を抽出します。
results_genes = stattest(bg_chrX_filt,
feature="gene",
covariate="sex",
adjustvars = c("population"),
getFC=TRUE,
meas="FPKM")
head(results_genes)
feature id fc pval qval
1 gene MSTRG.1 1.0824443 0.6847873 0.9196130
2 gene MSTRG.10 0.9335310 0.8132130 0.9599962
3 gene MSTRG.100 0.6126762 0.5687947 0.8823732
4 gene MSTRG.1000 1.2044223 0.4500794 0.8360675
5 gene MSTRG.1001 1.3799856 0.6140598 0.8936018
6 gene MSTRG.1003 0.9087684 0.5909100 0.8867180
以下ではresults_transcripts
に遺伝子名の情報を付加します。
results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),
geneIDs=ballgown::geneIDs(bg_chrX_filt),
results_transcripts)
head(results_transcripts)
geneNames geneIDs geneNames.1 geneIDs.1 feature id fc
1 . MSTRG.2 . MSTRG.528 transcript 1668 0.030602170
2 PLCXD1 MSTRG.2 XIST MSTRG.528 transcript 1667 0.002990002
3 . MSTRG.2 . MSTRG.528 transcript 1666 0.015949119
4 . MSTRG.2 . MSTRG.528 transcript 1669 0.028359333
5 . MSTRG.3 TSIX MSTRG.527 transcript 1665 0.077438632
6 PLCXD1 MSTRG.2 . MSTRG.612 transcript 1862 7.335431029
pval qval
1 1.474474e-10 2.824911e-07
2 2.497711e-10 2.824911e-07
3 6.564345e-10 4.949516e-07
4 6.871521e-08 3.885845e-05
5 1.911435e-06 8.647331e-04
6 1.331176e-05 5.018534e-03
遺伝子名はid順にソートされていますが、発現変動遺伝子が見たいので、p値の小さい順にソートしなおします。
results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)
head(results_transcripts)
geneNames geneIDs feature id fc pval qval
1 . MSTRG.528 transcript 1668 0.030602170 1.474474e-10 2.824911e-07
2 XIST MSTRG.528 transcript 1667 0.002990002 2.497711e-10 2.824911e-07
3 . MSTRG.528 transcript 1666 0.015949119 6.564345e-10 4.949516e-07
4 . MSTRG.528 transcript 1669 0.028359333 6.871521e-08 3.885845e-05
5 TSIX MSTRG.527 transcript 1665 0.077438632 1.911435e-06 8.647331e-04
6 . MSTRG.612 transcript 1862 7.335431029 1.331176e-05 5.018534e-03
head(results_genes)
feature id fc pval qval
1 gene MSTRG.528 0.002409209 2.041411e-11 2.086323e-08
2 gene MSTRG.527 0.082689539 6.110906e-06 2.409310e-03
3 gene MSTRG.143 3.158100078 7.719104e-06 2.409310e-03
4 gene MSTRG.612 7.327271932 1.073484e-05 2.409310e-03
5 gene MSTRG.763 0.278924648 1.178723e-05 2.409310e-03
6 gene MSTRG.618 9.109365087 4.555636e-05 7.759766e-03
以下のコマンドを実行すると、得られたresults_transcripts
, results_genes
をテキストファイルに出力します。(ファイル名は任意)
write.csv(results_transcripts, "chrX_transcript_results.csv", row.names=FALSE)
write.csv(results_genes, "chrX_gene_results.csv", row.names=FALSE)
可視化
以下は可視化のコマンドです。得られた発現変動遺伝子についてさらに調査する時に用います。
tropical= c('darkorange', 'dodgerblue', 'hotpink', 'limegreen', 'yellow')
palette(tropical)
fpkm = texpr(bg_chrX,meas="FPKM")
fpkm = log2(fpkm+1)
boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')
ラベルが見切れていますが、サンプルごとの遺伝子発現量の分布が性別ごとに色分けされて箱ひげ図で表示されました。
縦軸のFPKMというのは遺伝子ごとにマップされたリード数を遺伝子長と全リード数で正規化したものです。
次に、特定の遺伝子について可視化してみましょう。
ここでは例として12番めの遺伝子について可視化します。
ballgown::transcriptNames(bg_chrX)[12]
12 "NR_027231"
ballgown::geneNames(bg_chrX)[12]
12 "LINC00685"
おや、、、何故だかNature Protocolの例と表示される遺伝子名が違います。順番が変わったのでしょうか。
ともかく表示してみましょう。
plot(fpkm[12,] ~ pheno_data$sex, border=c(1,2),
main=paste(ballgown::geneNames(bg_chrX)[12],' : ',
ballgown::transcriptNames(bg_chrX)[12]),
pch=19,
xlab="Sex",
ylab='log2(FPKM+1)')
points(fpkm[12,] ~ jitter(as.numeric(pheno_data$sex)),
col=as.numeric(pheno_data$sex))
サンプルごとの発現量のばらつきを性別で分けてプロットしているのですが、全てのサンプルで発現量が0のため、差が見えません。
この遺伝子ではない遺伝子が例になっていたはずなのですが。僕のスクリプトにどこか間違っているところを見つけたら教えてください。
plotTranscripts(ballgown::geneIDs(bg_chrX)[1729],
bg_chrX,
main=c('Gene XIST in sample ERR188234'),
sample=c('ERR188234'))
遺伝子の簡単なisoform情報も表示することができます。
plotMeans('MSTRG.56', bg_chrX_filt, groupvar="sex", legend=FALSE)
ここではおそらく性別間で差がある遺伝子構造を可視化しているのですが、差がない遺伝子を指定していまっているのでよくわかりません。
おわりに
最後ミスったためしまりがないことになってしまいましたが、ともかく一通りの作業ができるようになりました。
新規転写物を探索したい時には有効なワークフローなので、機会があれば利用されると良いと思います。