edgeRでの多群間発現変動解析について、少し実験してみました。
今回はこのような仮想的なデータを用意しました。G1、G2、G3という３つのグループを含むという設定です。今回用いるデータはかなりノイズが少ない現実のデータとは少し遠いシミュレーションデータを用いています。ただ、発現変動解析の挙動を把握するにはちょうど良いと思います。
| 発現パターン | 遺伝子数 | |
|---|---|---|
| 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 |
１行目の"G1 = G2 = G3"は非発現変動遺伝子nonDEGを表します。また、2,3行目はG1対他で発現変動、 4,5行目はG2対他で発現変動、 6,7行目はG3対他で発現変動している遺伝子を表します。最後に、8,9行目は全てのグループ間で発現変動している遺伝子を表します。これらの発現変動のパターンを図に示すと以下の様になります。また、図で表現するといろいろ面倒なので、下図中の矢印で示すようにそれぞれテキストでの表現を対応させます。この表現の仕方は、baySeqというRパッケージに倣っています。
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と判定されているものを○で示すと以下の図のようになります。
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 |
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が発現量の差が最も大きいからだと考えられます。
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 |
一番理解しやすいのは１行目の情報で、全ての比較でDEGと判定された遺伝子は{G1},{G2},{G3}を予測しているという結果を示しています。逆に最終行では全ての比較でnonDEGと表現された場合は{G1,G2,G3}を予測しているといえるでしょう。また、２行目のG1対G2ではnonDEGと予測され、その他の組み合わせではDEGと判定された遺伝子が{G3},{G1,G2}という発現パターンのDEGを予測しているという結果が示されています。３行目と５行目も同様です。残った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 |
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 |
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 |
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 |
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 |
