・・・・・¡

いろいろ

edgeRで3群間比較してみた

edgeRでの多群間発現変動解析について、少し実験してみました。

今回はこのような仮想的なデータを用意しました。G1、G2、G3という3つのグループを含むという設定です。今回用いるデータはかなりノイズが少ない現実のデータとは少し遠いシミュレーションデータを用いています。ただ、発現変動解析の挙動を把握するにはちょうど良いと思います。

発現パターン 遺伝子数
1 G1 = G2 = G3 16000
2 G1 < G2 = G3 500
3 G1 > G2 = G3 500
4 G2 < G3 = G1 500
5 G2 > G3 = G1 500
6 G3 < G1 = G2 500
7 G3 > G1 = G2 500
8 G1 < G2 < G3 500
9 G3 < G2 < G1 500

1行目の"G1 = G2 = G3"は非発現変動遺伝子nonDEGを表します。また、2,3行目はG1対他で発現変動、 4,5行目はG2対他で発現変動、 6,7行目はG3対他で発現変動している遺伝子を表します。最後に、8,9行目は全てのグループ間で発現変動している遺伝子を表します。これらの発現変動のパターンを図に示すと以下の様になります。また、図で表現するといろいろ面倒なので、下図中の矢印で示すようにそれぞれテキストでの表現を対応させます。この表現の仕方は、baySeqというRパッケージに倣っています。

f:id:menseki:20191126145039p:plain

table(truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
16000 1000 1000 1000 1000

下準備

今回はこの仮想的なデータを用いて、有名な発現変動解析用RパッケージedgeRを用いて3群間の比較を行っていきたいと思います。 まずは下準備として、DGEList関数によるedgeR用のオブジェクト生成と、calcNormFactors関数による正規化係数の計算を行います。

library(edgeR)
##> Loading required package: limma
group <- gl(3, 3, labels = c("G1", "G2", "G3"))
design <- model.matrix(~group)
y <- DGEList(counts = counts, group = gl(3, 3, labels = c("G1", "G2", "G3")))

y <- calcNormFactors(y)

続いてデータのばらつきの指標であるDispersionパラメータを推定します。

y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)

G1 対 G2 の exact test

まずは正確検定によるG1とG2の比較を行います。

et <- exactTest(y, pair = c("G1", "G2"))
res <- topTags(et, n = nrow(counts), sort.by = "none")
pred_g1g2 <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_g1g2, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 115 741 744 4 786
nonDEG 15885 259 256 996 214

edgeRではp値とともにFDRという値が得られるのですが、一般的にこの値を用いて帰無仮説を棄却するかを判定していくことが多いと思うので、ここではFDRが0.05以下の遺伝子を発現変動遺伝子とします。その結果、G1とG2の間で変動がある遺伝子がDEGと判定されていることがわかります。表の見方はRのtable関数のヘルプを参照してください。DEGと判定されているものを○で示すと以下の図のようになります。

f:id:menseki:20191126145042p:plain

G2 対 G3 の exact test

同様にしてG2対G3の正確検定もやってみます。

et <- exactTest(y, pair = c("G2", "G3"))
res <- topTags(et, n = nrow(counts), sort.by = "none")
pred_g2g3 <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_g2g3, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 125 1 725 781 784
nonDEG 15875 999 275 219 216

f:id:menseki:20191126145045p:plain

G3 対 G1 の exact test

さらに、G3対G1の正確検定もやってみます。

et <- exactTest(y, pair = c("G3", "G1"))
res <- topTags(et, n = nrow(counts), sort.by = "none")
pred_g3g1 <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_g3g1, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 132 766 8 755 1000
nonDEG 15868 234 992 245 0

ここで{G1},{G2},{G3}の遺伝子が100%DEGと判定されています。これは、今回の仮想データに含まれる発現パターンのうち、{G1},{G2},{G3}に対応するものはG1<G2<G3またはG3<G2<G1であり、G1とG3が発現量の差が最も大きいからだと考えられます。

f:id:menseki:20191126145049p:plain

exact testの全ての結果を統合

以上の結果をまとめて表示すると以下のようになります。

plyr::adply(table(pred_g1g2, pred_g2g3, pred_g3g1, truth), 1:3)
pred_g1g2 pred_g2g3 pred_g3g1 {G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG DEG DEG 0 0 2 1 608
nonDEG DEG DEG 12 0 4 677 176
DEG nonDEG DEG 16 654 2 3 178
nonDEG nonDEG DEG 104 112 0 74 38
DEG DEG nonDEG 14 1 628 0 0
nonDEG DEG nonDEG 99 0 91 103 0
DEG nonDEG nonDEG 85 86 112 0 0
nonDEG nonDEG nonDEG 15670 147 161 142 0

一番理解しやすいのは1行目の情報で、全ての比較でDEGと判定された遺伝子は{G1},{G2},{G3}を予測しているという結果を示しています。逆に最終行では全ての比較でnonDEGと表現された場合は{G1,G2,G3}を予測しているといえるでしょう。また、2行目のG1対G2ではnonDEGと予測され、その他の組み合わせではDEGと判定された遺伝子が{G3},{G1,G2}という発現パターンのDEGを予測しているという結果が示されています。3行目と5行目も同様です。残った4,6,7行目はDEGの判定の組み合わせが矛盾しているデータです。例えば4行目ではG1対G2とG2対G3両方で発現変動していないと判定されていてこの時点でG3とG1でも発現変動していないのが自然なのですが、G1とG3で発現変動していると判定しています。このことは、統計検定の非対称性に起因するのではないかというのが私の考えです。つまり、検定で帰無仮説が棄却されなかった、FDR>0.05のデータに関して、DEGではない(nonDEGである)といった積極的な結論は導けないということです。

GLMバージョンの結果も示しておきます。

GLMに基づくDispersionパラメータの推定とデータへの当てはめ

y <- estimateGLMCommonDisp(y,design)
y <- estimateGLMTrendedDisp(y,design)
y <- estimateGLMTagwiseDisp(y,design)
fit <- glmFit(y,design)

G1 対 G2 のLRT

lrt <- glmLRT(fit,coef=2)
res <- topTags(lrt, sort.by = "none", n = nrow(counts))
pred_coef2 <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_coef2, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 131 753 761 6 791
nonDEG 15869 247 239 994 209

f:id:menseki:20191126145052p:plain

G2 対 G3 でLRT

lrt <- glmLRT(fit,contrast = c(0, -1, 1))
res <- topTags(lrt, sort.by = "none", n = nrow(counts))
pred_cont <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_cont, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 144 1 742 791 795
nonDEG 15856 999 258 209 205

f:id:menseki:20191126145056p:plain

G3 対 G1 でLRT

lrt <- glmLRT(fit,coef=3)
res <- topTags(lrt, sort.by = "none", n = nrow(counts))
pred_coef3 <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_coef3, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 145 771 12 764 1000
nonDEG 15855 229 988 236 0

f:id:menseki:20191126145059p:plain

ANOVA的な方法

lrt <- glmLRT(fit,coef=2:3)
res <- topTags(lrt, sort.by = "none", n = nrow(counts))
pred_coef23 <- factor(ifelse(res@.Data[[1]]$FDR <= 0.05, "DEG", "nonDEG"))
table(pred_coef23, truth)
{G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG 229 835 822 831 1000
nonDEG 15771 165 178 169 0

f:id:menseki:20191126145103p:plain

LRTの全ての結果を統合

plyr::adply(table(pred_coef2, pred_cont, pred_coef3, truth), 1:3)
pred_coef2 pred_cont pred_coef3 {G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG DEG DEG 0 0 2 1 621
nonDEG DEG DEG 15 0 4 688 174
DEG nonDEG DEG 16 664 5 3 170
nonDEG nonDEG DEG 114 107 1 72 35
DEG DEG nonDEG 17 1 649 2 0
nonDEG DEG nonDEG 112 0 87 100 0
DEG nonDEG nonDEG 98 88 105 0 0
nonDEG nonDEG nonDEG 15628 140 147 134 0
plyr::adply(table(pred_coef2, pred_cont, pred_coef3, pred_coef23, truth), 1:4)
pred_coef2 pred_cont pred_coef3 pred_coef23 {G1,G2,G3} {G1},{G2,G3} {G2},{G3,G1} {G3},{G1,G2} {G1},{G2},{G3}
DEG DEG DEG DEG 0 0 2 1 621
nonDEG DEG DEG DEG 15 0 4 688 174
DEG nonDEG DEG DEG 16 664 5 3 170
nonDEG nonDEG DEG DEG 57 94 1 58 35
DEG DEG nonDEG DEG 17 1 649 2 0
nonDEG DEG nonDEG DEG 63 0 75 77 0
DEG nonDEG nonDEG DEG 56 73 84 0 0
nonDEG nonDEG nonDEG DEG 5 3 2 2 0
DEG DEG DEG nonDEG 0 0 0 0 0
nonDEG DEG DEG nonDEG 0 0 0 0 0
DEG nonDEG DEG nonDEG 0 0 0 0 0
nonDEG nonDEG DEG nonDEG 57 13 0 14 0
DEG DEG nonDEG nonDEG 0 0 0 0 0
nonDEG DEG nonDEG nonDEG 49 0 12 23 0
DEG nonDEG nonDEG nonDEG 42 15 21 0 0
nonDEG nonDEG nonDEG nonDEG 15623 137 145 132 0

参考