Ballgown は主に StringTie で定量した遺伝子（トランスクリプト）発現量を利用して、発現変動解析を行う R/Bioconductor のパッケージである。Ballgown を利用するには、StringTie で定量した結果をすべて同じフォルダにまとめておと便利である。例えば、StringTie の結果をすべて ballgown フォルダにまとめてあったとして、Ballgown では次のようにデータを読み込むことができる。なお、デフォルトでは、各々の結果のフォルダ名がサンプル名として読み込まれるため、必要であれば pData 関数で書き換えることも可能である。
b <- ballgown(dataDir = 'ballgown', samplePattern = 'SRR')
pData(b) <- data.frame(id = c('ctrl_1', 'ctrl_2', 'ctrl_3', 'cold_1', 'cold_2', 'cold_3'),
group = c('ctrl', 'ctrl', 'ctrl', 'cold', 'cold', 'cold'))転写産物の発現量 FPKM を取得したい場合は texpr 関数を利用する。
head(texpr(b))
## FPKM.SRR1792924 FPKM.SRR1792925 FPKM.SRR1792926 FPKM.SRR1792927
## 1 4.125814 3.844283 4.593955 4.123420
## 2 2.796588 2.798470 3.924363 3.454246
## 3 0.534306 0.282695 0.668188 0.766483
## 4 0.725205 1.565021 1.152122 1.377866
## 5 0.619181 0.664552 0.453376 0.846209
## 6 1.171823 1.892226 0.772719 0.946749
## FPKM.SRR1792928 FPKM.SRR1792929
## 1 3.694850 3.848086
## 2 3.479400 3.589532
## 3 0.444744 0.635831
## 4 1.048007 0.984064
## 5 0.533180 0.154207
## 6 1.959617 1.625064遺伝子ごとの発現量を取得する場合は gexpr 関数を使用する。
head(gexpr(b))
## FPKM.SRR1792924 FPKM.SRR1792925 FPKM.SRR1792926 FPKM.SRR1792927
## AT1G01030 1.25360604 1.22689400 1.29295800 1.31943000
## AT1G01110 1.66181793 1.92547586 1.80818711 1.90685011
## AT1G01130 2.27273100 4.40685400 3.13476900 3.65375300
## AT1G01150 0.06257279 0.04701473 0.08466496 0.08402109
## AT1G01270 0.00000000 0.00000000 0.00000000 0.00000000
## AT1G01280 0.00000000 0.00000000 0.00000000 0.00000000
## FPKM.SRR1792928 FPKM.SRR1792929
## AT1G01030 1.48861482 1.3497740
## AT1G01110 1.79535574 1.6217058
## AT1G01130 3.83300400 3.7215440
## AT1G01150 0.08348642 0.0971909
## AT1G01270 0.00000000 0.0000000
## AT1G01280 0.00000000 0.0000000解析に先立ち、各サンプルの FPKM の分布状況を確認する。
fig_countdist <- as.data.frame(texpr(b)) %>%
pivot_longer(everything(), names_to = 'library') %>%
mutate(value = value + 1) %>%
ggplot(aes(x = value)) +
geom_histogram(bins = 50) +
scale_x_log10() +
facet_wrap(~ library, nrow = 2)
print(fig_countdist)
次に、発現量量の分散がある程度大きい遺伝子を抽出して、主成分分析を行ってみる。この際に、FPKM のスケールが大きいので対数化しておく。
b_highexp <- subset(b, 'rowVars(texpr(b)) > 1', genomesubset = TRUE)
tx_highexp <- texpr(b_highexp)
x_pca <- prcomp(t(log10(tx_highexp + 1)), scale = FALSE)
fig_pca <- as.data.frame(x_pca$x) %>%
mutate(libname = rownames(x_pca$x)) %>%
ggplot(aes(x = PC1, y = PC2, label = libname)) +
geom_label()
print(fig_pca)
サンプル同士の関係を確認してから、発現変動遺伝子（トランスクリプト）を検出する。検定は Ballgown の statetest 関数を使用する。比較する際のグループ情報は covariate オプションで指定する。このオプションに pData 関数に代入した実験デザインの列名を指定する。
det_table <- stattest(b_highexp, feature = 'transcript', covariate = 'group', getFC = TRUE, meas = 'FPKM')
det_table <- data.frame(gene = geneNames(b_highexp),
gene_id = geneIDs(b_highexp),
log_mean = log10(rowMeans(texpr(b_highexp))),
det_table)
head(det_table)
## gene gene_id log_mean feature id fc pva qval
## <chr> <chr> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 30 . MSTRG.9 1.26729358 transcript 30 0.8769386 0.455128749 0.7547076
## 34 . MSTRG.11 2.07281850 transcript 34 1.2155198 0.651011430 0.8625935
## 35 . MSTRG.11 0.68589443 transcript 35 0.6007765 0.780760395 0.9188453
## 36 . MSTRG.11 1.34838370 transcript 36 0.1875270 0.689512569 0.8807716
## 37 . MSTRG.11 -0.09451057 transcript 37 2.4095289 0.636693122 0.8543992
## 40 . MSTRG.12 1.61330432 transcript 40 1.2792621 0.009098968 0.1227936MA プロットを描く場合は、平均と foldchange を検定結果から抽出して描けばよい。発現量に有意差のあるトランスクリプトを赤色（#d62728）の点で示した。
fig_ma <- ggplot(det_table, aes(x = log_mean, y = log2(fc), color = qval < 0.1)) +
geom_point() +
scale_color_manual(values = c('#666666', '#d62728'))
print(fig_ma)