多群間比較（edgeR） |

edgeR パッケージを利用して多群間比較データの解析方法を示す。edgeR は一般化線形モデルによる解析をサポートしているため、複雑なモデルにも柔軟に対応でき、非常に優れたパッケージである。

このページに掲載した解析結果は edgeR バージョン 3.10.0 の結果である。

library(edgeR)
packageVersion("edgeR")
##[1] ‘3.10.0’

以下にこのページで使用するサンプルデータを紹介してから、一般化線形モデルを利用した発現変動遺伝子の検出法について紹介する。

サンプルデータ
検定

サンプルデータ

サンプルデータは、カウントデータが負の二項分布に従うように乱数を生成して作成した。A 群、B 群および C 群の 3 つの実験群からなり、それぞれの実験群には 3 つの biological replicate がある。また、全遺伝子数は 1,000 個であり、その発現パターンは以下のように設定した。

遺伝子	A 群で高発現	B 群で高発現	C 群で高発現
gene_1 ～ gene_50	Yes	No	No
gene_51 ～ gene_100	No	Yes	No
gene_101 ～ gene_150	No	No	Yes
gene_151 ～ gene_200	Yes	Yes	No
gene_201 ～ gene_250	No	Yes	Yes
gene_251 ～ gene_300	Yes	No	Yes
gene_201 ～ gene_1000	No	No	No

サンプルデータを読み込んでから行列型に変換し、count 変数に代入する。

count <- read.table("https://bi.biopapyrus.jp/data/counts_3_3_3.txt", sep = "\t", header = T, row.names = 1)
count <- as.matrix(count)

dim(count)
## [1] 1000    9

head(count)
##         A1  A2   A3  B1  B2  B3  C1  C2  C3
## gene_1   3   2   19   0   0   0  13   0   0
## gene_2 680 807 1011 168 132 138 135 188 202
## gene_3  35  38   30  11   9  11   8   9   8
## gene_4 103 165  124  35  37  31  36  36  30
## gene_5   2   9    0   0   0   0   0   0   0
## gene_6 419 373  478 114 117 109  89 142 106

実験群の識別ラベルを作成し group 変数に代入する。

group <- factor(c("A", "A", "A", "B", "B", "B", "C", "C", "C"))

尤度比検定による多群間比較

尤度比検定は 2 つの統計モデルのそれぞれから計算される尤度を比較することにより検定を行う方法であり、一般化線形モデルでよく使われる検定の一つである。この考え方は発現変動遺伝子の検出にも適用できる。そこで、2 つの統計モデルをどのように作成すれば、「A 群、B 群と C 群」の発現量の差の検定に応用できるのかについて見ていく。

一般化線形モデルは、確率変数を Y とし、確率変数が従う分布のパラメーターを β としたとき、Y の期待値は g(E[Y]) = Xβ と表すことのできる統計モデルをいう。モデルに含まれる関数 g はリンク関数とよばれる、E[Y] と右辺の Xβ が等しくなるように変換を行う関数である。データが正規分布であるときリンク関数は g(x) = x、負の二項分布のときリンク関数は g(x) = log(x) のように、分布が決まればリンク関数の形もほぼ決まる。Y は観測値であり、g は観測値に依存して解析前にその形がすでにわかるため、実質的には g(E[Y]) について考えなくてよい。以下、説明を簡略するためにリンク関数を省略して書く。重要な事は右辺の Xβ をいかに設計することである。X はデザイン行列あるいは計画行列と呼ばれ、実験群や処理群の情報などが行列で記述される。β はパラメーターベクトルとよばれ、モデル化する際に想定されるパラメーターからなるベクトルである。

ここで二群間二因子比較データをモデル化する。まず、各実験群の発現量の差を高々 3 つのパラメーターで表すことができる。すなわち、E[A] = β₁、E[B] = β₁ + β₂、E[C] = β₁ + β₃ とする。サンプルデータは 3 つの実験群で、各実験群において 3 つの biological replicate が含まれているため、これを行列の形で表すと以下のようになる。

\[ \begin{pmatrix} E[\mbox{A}]\\ E[\mbox{A}] \\ E[\mbox{A}] \\ E[\mbox{B}] \\ E[\mbox{B}] \\ E[\mbox{B}] \\ E[\mbox{C}] \\ E[\mbox{C}] \\ E[\mbox{C}] \end{pmatrix} ^{T} = \begin{pmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3}\end{pmatrix} \]

このモデルは考えられるパラメーターの 3 つすべてを含んでいるため full model とよぶことにする。また、このモデルは 1 つの遺伝子のみについて着目した時のものであり、遺伝子が n 個があるときこのようなモデルは n 個できる。遺伝子が n 個あるときプログラムが逐次的に処理を行ってくれるため、遺伝子が 1 個だけの時を考えればよい。

full model のデザイン行列は model.matrix 関数を利用すれば簡単に作成できる。

design <- model.matrix(~ group)
design
##   (Intercept) groupB groupC
## 1           1      0      0
## 2           1      0      0
## 3           1      0      0
## 4           1      1      0
## 5           1      1      0
## 6           1      1      0
## 7           1      0      1
## 8           1      0      1
## 9           1      0      1
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"

edgeR を利用して、データを正規化、分布をパラメーターの推定を行う。

d <- DGEList(counts = count, group = group)
d <- calcNormFactors(d)
d <- estimateGLMCommonDisp(d, design)
d <- estimateGLMTrendedDisp(d, design)
d <- estimateGLMTagwiseDisp(d, design)

推定された分布にデータを当てはめていく。

fit <- glmFit(d, design)

次に尤度比検定を行い発現変動遺伝子を検出する。尤度比検定では 2 つのモデルの尤度を比較して検定を行うために、上で作成した full model の他にもう 1 つのモデルを必要とする。そのもう 1 つの方のモデルは、どのような発現変動遺伝子を期待するかによって作り方が異なる。例えば以下のような比較が考えられる。

A 群と B 群の比較
B 群と C 群の比較
C 群と A 群の比較
三群を同時に比較（ANOVA-like）

A 群と B 群の比較

C 群の発現量を考慮に入れずに、A 群と B 群だけを比較し発現変動遺伝子を検出する例を示す。A 群と B 群を比較するために帰無仮説を E[A] = E[B] とすればよい。すなわち、β₂ = 0 を考えればよい。ここで、β₂ = 0 となるモデルを作成するが、full model のデザイン行列の 2 列目をすべて 0 にすれば、β₂ = 0 に対応できる。

\[ \begin{pmatrix} E[\mbox{A}]\\ E[\mbox{A}] \\ E[\mbox{A}] \\ E[\mbox{B}] \\ E[\mbox{B}] \\ E[\mbox{B}] \\ E[\mbox{C}] \\ E[\mbox{C}] \\ E[\mbox{C}] \end{pmatrix} ^{T} = \begin{pmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3}\end{pmatrix} \]

このように full model からパラメーター β₂ を取り除いてできたモデルを reduced model という。尤度比検定では full model と reduce model の尤度を比較して検定を行う。もし両者の尤度に差がなければ帰無仮説は保留される。反対に、両者の尤度に差があれば帰無仮説は棄却され、β₂ ≠ 0 と考えられ、その遺伝子は発現変動遺伝子と考えられる。

次に glmLRT 関数を利用して尤度比検定を行う。fit 変数には full model の情報が含まれているため、reduced model の情報だけを与えればよい。ここで、full model と reduced model のデザイン行列を見比べてみると、reduced model は full model のデザイン行列の 2 列をすべて 0 にしたものに等しいことがわかる。そこで、coef = 2 と指定して、reduced model の情報を与える。

lrt <- glmLRT(fit, coef = 2)

結果は topTags 関数を利用して確認する。

topTags(lrt)
## Coefficient:  groupB 
##              logFC   logCPM       LR       PValue          FDR
## gene_2   -2.512281 11.67577 39.38936 3.471837e-10 3.471837e-07
## gene_212  2.404434 11.72035 36.68405 1.389110e-09 4.330560e-07
## gene_232  2.574584 10.97919 36.21190 1.769862e-09 4.330560e-07
## gene_248  2.420228 11.21263 35.86445 2.115332e-09 4.330560e-07
## gene_89   2.278600 12.19046 35.63868 2.375238e-09 4.330560e-07
## gene_86   2.621158 10.39386 35.46381 2.598336e-09 4.330560e-07
## gene_62   2.337651 13.59625 34.86642 3.531194e-09 5.044563e-07
## gene_284 -2.362287 11.90213 34.57919 4.092582e-09 5.115728e-07
## gene_253 -2.286831 10.98825 32.07388 1.484191e-08 1.649101e-06
## gene_294 -2.315419 11.30366 31.00439 2.574460e-08 2.574460e-06

A 群と B 群の比較で発現量に差が生じる遺伝子を gene_1 ～ gene_100、gene_201 ～ gene_300 として設定してある。検定結果を眺めると、この範囲にある遺伝子が上位にランキングされていることが確認できる。

解析結果を result.txt ファイルに保存する。

table <- as.data.frame(topTags(lrt, n = nrow(count)))
write.table(table, file = "result.txt", col.names = T, row.names = T, sep = "\t")

A 群と C 群の比較

A 群と B 群の比較と同様に考えて、full model から β₃ を取り除いたモデルを reduced model とすればよい。すなわち、

\[ \begin{pmatrix} E[\mbox{A}]\\ E[\mbox{A}] \\ E[\mbox{A}] \\ E[\mbox{B}] \\ E[\mbox{B}] \\ E[\mbox{B}] \\ E[\mbox{C}] \\ E[\mbox{C}] \\ E[\mbox{C}] \end{pmatrix} ^{T} = \begin{pmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3}\end{pmatrix} \]

reduced model のデザイン行列は full model のデザイン行列の 3 列目をすべて 0 にしたときのデザイン行列に等しいので、coef = 3 としてオプションを与えて尤度比検定を行う。

lrt <- glmLRT(fit, coef = 3)
topTags(lrt)
## Coefficient:  groupC 
##              logFC   logCPM       LR       PValue          FDR
## gene_248  2.472652 11.21263 37.31415 1.005523e-09 9.787191e-07
## gene_181 -2.371428 10.68939 36.01560 1.957438e-09 9.787191e-07
## gene_111  2.236228 14.67688 34.46111 4.348533e-09 1.449511e-06
## gene_232  2.459777 10.97919 33.27681 7.993005e-09 1.928945e-06
## gene_108  2.251165 13.84255 32.91158 9.644726e-09 1.928945e-06
## gene_2   -2.231827 11.67577 31.83299 1.680147e-08 2.633894e-06
## gene_212  2.215546 11.72035 31.51011 1.984051e-08 2.633894e-06
## gene_136  2.182334 11.24055 31.39326 2.107115e-08 2.633894e-06
## gene_145  2.137884 11.32728 30.25729 3.783661e-08 4.204068e-06
## gene_10  -2.308701 10.62793 30.03748 4.237752e-08 4.237752e-06

サンプルデータを生成する際に、A 群と C 群の比較において発現量に差が生じる遺伝子を gene_1 ～ gene_50、gene_101 ～ gene_200、gene_251 ～ gene_300 として設定してある。検定結果を眺めると、この範囲にある遺伝子が上位にランキングされていることが確認できる。

B 群と C 群の比較

B 群と C 群を比較したい場合、full model の 1 列目を考慮しなければよい。しかし、coef = 1 としたとき、β₁ = 0 となってしまうため、全体に影響が生じる。これは、モデルは β₁ を基準にして作成したためである。そこで、coef のかわりに contrast 引数を利用する。B 群と C 群の比較を行うので、帰無仮説としては E[B] - E[C] = 0 すなわち、β₁ × 0 + β₂ - β₃ = 0 を考えればよい。これを指定するには contrast を以下のように 0, 1, -1 の順で代入すれば良い。

lrt <- glmLRT(fit, contrast = c(0, 1, -1))
topTags(lrt)
## Coefficient:  1*groupB -1*groupC 
##              logFC   logCPM       LR       PValue          FDR
## gene_185  5.275530 10.44068 39.86899 2.715803e-10 1.368306e-07
## gene_111 -2.422004 14.67688 39.85408 2.736613e-10 1.368306e-07
## gene_145 -2.451985 11.32728 38.73505 4.854142e-10 1.575684e-07
## gene_187  2.616635 11.62883 38.22538 6.302736e-10 1.575684e-07
## gene_89   2.292847 12.19046 36.03717 1.935891e-09 3.721380e-07
## gene_181  2.362360 10.68939 35.75913 2.232828e-09 3.721380e-07
## gene_62   2.328929 13.59625 34.62890 3.989400e-09 5.699143e-07
## gene_119 -2.288232 10.59303 33.90020 5.801281e-09 7.251601e-07
## gene_274 -2.297583 12.78458 33.60577 6.749212e-09 7.499124e-07

B 群と C 群の比較で発現量に差が生じる遺伝子を gene_51 ～ gene_200、gene_251 ～ gene_300 として設定してある。検定結果を眺めると、この範囲にある遺伝子が上位にランキングされていることが確認できる。

三群間比較（ANOVA like）

3 群間の発現量を同時に比較を行うとき、帰無仮説として「E[A] = E[B] かつ E[A] = E[C]」と考えればよい。つまり、「β₂ = 0 かつ β₃ = 0」に対応するモデルを reduced model とすればよい。すなわち、

\[ \begin{pmatrix} E[\mbox{A}]\\ E[\mbox{A}] \\ E[\mbox{A}] \\ E[\mbox{B}] \\ E[\mbox{B}] \\ E[\mbox{B}] \\ E[\mbox{C}] \\ E[\mbox{C}] \\ E[\mbox{C}] \end{pmatrix} ^{T} = \begin{pmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3}\end{pmatrix} \]

reduced model のデザイン行列は full model のデザイン行列の 2 行目と 3 行目の要素をすべて 0 に変換したものに等しい。そこで、coef = c(2, 3) と指定して尤度比検定を行う。

lrt <- glmLRT(fit, coef = c(2, 3))
topTags(lrt)
## Coefficient:  groupB groupC 
##          logFC.groupB logFC.groupC   logCPM       LR       PValue          FDR
## gene_111  -0.18577617   2.23622804 14.67688 59.01008 1.535054e-13 1.535054e-10
## gene_2    -2.51228128  -2.23182721 11.67577 56.66059 4.969449e-13 1.677943e-10
## gene_89    2.27860014  -0.01424671 12.19046 56.63485 5.033829e-13 1.677943e-10
## gene_62    2.33765089   0.00872179 13.59625 55.13243 1.066950e-12 2.667375e-10
## gene_145  -0.31410103   2.13788389 11.32728 54.62748 1.373389e-12 2.746779e-10
## gene_108   0.04749834   2.25116517 13.84255 50.88143 8.937952e-12 1.489659e-09
## gene_119  -0.22673783   2.06149421 10.59303 48.57715 2.828819e-11 4.041170e-09
## gene_136   0.08517532   2.18233372 11.24055 47.38078 5.145109e-11 6.431386e-09
## gene_10   -2.22916957  -2.30870083 10.62793 45.89531 1.081334e-10 1.174930e-08
## gene_102  -0.03968937   1.97643175 12.56623 45.72928 1.174930e-10 1.174930e-08

サンプルデータでは gene_1 ～ gene_300 を発現変動遺伝子として生成してある。検定結果をみると、この範囲にある遺伝子が上位にランキングされていることがわかる。

解析結果を result.txt ファイルに保存する。

table <- as.data.frame(topTags(lrt, n = nrow(count)))
write.table(table, file = "result.txt", col.names = T, row.names = T, sep = "\t")

References

Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010, 26(1):139-40. PubMed Abstract Bioconductor
McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012, 40(10):4288-97. PubMed Abstract
Robinson MD, Smyth GK. Moderated statistical tests for assessing differences in tag abundance. Bioinformatics. 2007, 23(21):2881-7. PubMed Abstract
Robinson MD, Smyth GK. Small-sample estimation of negative binomial dispersion, with applications to SAGE data. Biostatistics. 2008, 9(2):321-32. PubMed Abstract
Zhou X, Lindsay H, Robinson MD. Robustly detecting differential expression in RNA sequencing data using observation weights. Nucleic Acids Research. 2014, 42(11):e91. PubMed Abstract
一般化線形モデル. biopapyrus バイオスタティスティックス