Palmsonntagmorgen

NGSデータを使った解析と、その周辺。

HISAT-StringTie-Ballgown を試してみよう その2

前回の記事の続きです。

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に格納されたものを確認
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)

# 発現変動遺伝子の推定(transcript単位)
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)

ここではおそらく性別間で差がある遺伝子構造を可視化しているのですが、差がない遺伝子を指定していまっているのでよくわかりません。

おわりに

最後ミスったためしまりがないことになってしまいましたが、ともかく一通りの作業ができるようになりました。
新規転写物を探索したい時には有効なワークフローなので、機会があれば利用されると良いと思います。