多群間二因子比較（edgeR） |

edgeR パッケージは発現変動遺伝子を検出する際に最もよく利用されているパッケージのうちの 1 つである。一般化線形モデルによる解析をサポートしているため、二群間比較のみならず多群間や多因子のデータにも柔軟に対応でき、非常に優れている。

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

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

サンプルデータ
検定

サンプルデータ

サンプルデータは条件を設定して負の二項分布に従うように乱数生成により作成した。WT 群、C1 群、C2 群が 3 つの実験群からなり、それぞれの実験群を 2 つの処理群に分け、それぞれを A 処理および B 処理を施す。いずれの群（状態）と処理において 2 つ biological replicate がある。また、全遺伝子数は 1,000 個であり、各遺伝子の発現パターンは以下のように設定した。

遺伝子	WT, C1, C2 群の比較	A, B 処理の比較
gene_1 ～ gene_25	WT のみで高発現	発現量に差がない
gene_26 ～ gene_50	C1 のみで高発現	発現量に差がない
gene_51 ～ gene_75	C2 のみで高発現	発現量に差がない
gene_76 ～ gene_100	WT と C1 で高発現	発現量に差がない
gene_101 ～ gene_125	C1 と C2 で高発現	発現量に差がない
gene_126 ～ gene_150	WT と C2 で高発現	発現量に差がない
gene_151 ～ gene_175	発現量に差がない	A で高発現
gene_176 ～ gene_200	発現量に差がない	B で高発現
gene_201 ～ gene_1000	発現量に差がない	発現量に差がない

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

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

dim(count)
## [1] 1000    12

head(count)
##        WT1_A WT2_A WT3_B WT4_B C11_A C12_A C13_B C14_B C21_A C22_A C23_B C24_B
## gene_1     3     2    40    12     2     5     9     0     0     0     0     0
## gene_2   680   807   594   701   126   213   143   188   218    85    99   103
## gene_3    35    38    31    33     8    10     9     9     8    12    16    11
## gene_4   103   165   110   196    68    51    48    36    49    32    43    53
## gene_5     2     9     1     1     1     1     0     0     3     1     0     0
## gene_6   419   373   399   498    71    60   123   142   137    84   120   116

実験群の識別ラベルを作成し group 変数に代入する。また、処理群の識別ラベルを作成し treat 変数に代入する。

group <- factor(c("WT", "WT", "WT", "WT", "C1", "C1", "C1", "C1", "C2", "C2", "C2", "C2"))
treat <- factor(c("A", "A", "B", "B", "A", "A", "B", "B", "A", "A", "B", "B"))

尤度比検定

尤度比検定では 2 つのモデルから尤度を計算し比較することによって、2 つのモデルに差があるかどうかを検定する方法である。リードカウントデータをモデル化することによって、尤度比検定を発現変動遺伝子の検出に応用できる。尤度比検定で比較する 2 つのモデルはそれぞれ full model と reduced model と呼ばれる。full model には多くのパラメーターを含み、reduced model は full model より少ないパラメーターをしか含まない。そのため、両者を比べることで、reduced model に含まれていないパラメーターの効果を調べることができる。生物実験において、野生型と遺伝子 g のノックアウト型を比較して、遺伝子 g の機能を調べることに似ている。

一般化線形モデルを利用してモデル化を行うとき、確率変数を Y とし、確率変数が従う分布のパラメーターを β とすると、Y の期待値は g(E[Y]) = Xβ と表すことができる。ただし、g はリンク関数である。リンク関数は期待値 E[Y] と右辺が等しくなるように変換を行う単調な関数であり、しかも、その形は決まっている。例えば、データの分布は正規分布ならば g(x) = x、負の二項分布ならば g(x) = log(x) である。そこで、以下、説明を簡略するためにリンク関数を省略して書く。

まず、full model について考える。三群間二因子比較データにおいて基準を WT 群 A 処理とする。WT 群 A 処理の発現量の期待値を β₁ とおく。次に、WT 群と C1 群の発現量の差を β₂ とし、WT 群と C2 群の発現量の差を β₃ とする。また、処理 A と処理 B の差を β₄ とすれば、各実験群と処理群を組み合わせにおいて、その発現量は行列を用いて以下のように表すことができる。

\[ \begin{pmatrix} E[\mbox{WT-A}] \\ E[\mbox{WT-A}] \\ E[\mbox{WT-B} \\ E[\mbox{WT-B}] \\ E[\mbox{C1-A}] \\ E[\mbox{C1-A}] \\ E[\mbox{C1-B}] \\ E[\mbox{C1-B}] \\ E[\mbox{C2-A}] \\ E[\mbox{C2-A}] \\ E[\mbox{C2-B}] \\ E[\mbox{C2-B}] \end{pmatrix}^{T} = \begin{pmatrix} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 1\\ 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 1\\ 1 & 1 & 0 & 1\\ 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 1\\ 1 & 0 & 1 & 1\\ \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3} \\ \beta_{4} \end{pmatrix} \]

この式が一般化線形モデルのモデルとされる。0 と 1 からなる行列はデザイン行列と呼び、観測値とパラメーターを関連させる重要な働きを持つ。遺伝子が n　個があればこのような式は n 個できる。遺伝子が n　個あるときプログラムが逐次的に処理を行ってくれるため、遺伝子が 1 個だけの時を考えればよい。

full model のデザイン行列は model.matrix 関数を利用すれば簡単に作成できる。R で作成したデザイン行列は上のデザイン行列の並びが異なるが、解析の効果として同じ。（ただし、上のデザイン行列の 2 列目と 3 列目それぞれ R で作成したで 3 列目と 2 列目に対応していることを注意しておく。）

design <- model.matrix(~ group + treat)
design
##    (Intercept) groupC2 groupWT treatB
## 1            1       0       1      0
## 2            1       0       1      0
## 3            1       0       1      1
## 4            1       0       1      1
## 5            1       0       0      0
## 6            1       0       0      0
## 7            1       0       0      1
## 8            1       0       0      1
## 9            1       1       0      0
## 10           1       1       0      0
## 11           1       1       0      1
## 12           1       1       0      1
## attr(,"assign")
## [1] 0 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"
## 
## attr(,"contrasts")$treat
## [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 つのモデル reduced model を必要とする。reduced model は、どのような発現変動遺伝子を期待するかによって作り方が異なる。

実験群間の比較（WT、C1、C2）

処理 A と処理 B の影響を考慮しないで、実験群同士の比較において差を持つ遺伝子を発現変動遺伝子として検出する方法を示す。このような比較において帰無仮説としては E[WT] = E[C1] = E[C2] が考えられる。これは β₂ = β₃ = 0 に対応している。そこで、「β₂ = 0 かつ β₃ = 0」に対応しているモデルを作成する。そこで、full model と reduced model を比べた時、「β₂ = 0 かつ β₃ = 0」を比べられるように reduced model を作成すればよい。例えば、以下のように作成できる。つまり、full model の 2 列目と 3 列目の要素をすべて 0 に置換する。

\[ \begin{pmatrix} E[\mbox{WT-A}] \\ E[\mbox{WT-A}] \\ E[\mbox{WT-B} \\ E[\mbox{WT-B}] \\ E[\mbox{C1-A}] \\ E[\mbox{C1-A}] \\ E[\mbox{C1-B}] \\ E[\mbox{C1-B}] \\ E[\mbox{C2-A}] \\ E[\mbox{C2-A}] \\ E[\mbox{C2-B}] \\ E[\mbox{C2-B}] \end{pmatrix}^{T} = \begin{pmatrix} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 1\\ \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3} \\ \beta_{4} \end{pmatrix} \]

このように用意した full model と reduced model から尤度を計算し比較を行う。もし両者の尤度に差異が認められた場合、full model と reduced model は相違異なるモデルと考えられる。両者の差は　β₂ と β₃ によって持たされているため、両者に差異が認められるならば「β₂ ≠ 0 または β₃ ≠ 0」である。しかし、これは最初の仮定、すなわち帰無仮説「β₂ = 0 かつ β₃ = 0」に矛盾する。つまり帰無仮説が間違いであり、棄却するべきである。このとき、WT、C1、C2 のいずれかで発現量に有意差を持つと考えられる。

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

lrt <- glmLRT(fit, coef = c(2, 3))

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

topTags(lrt)
## Coefficient:  groupC2 groupWT 
##         logFC.groupC2 logFC.groupWT    logCPM       LR       PValue          FDR
## gene_22  -0.006182465    2.21250976 10.916099 75.97231 3.182902e-17 3.182902e-14
## gene_70   2.084501060   -0.05159736 12.001509 71.32798 3.245829e-16 1.622915e-13
## gene_10  -0.064842824    1.98729484 10.554947 63.46559 1.654332e-14 5.514440e-12
## gene_39  -1.970732212   -2.05117152 11.554653 62.60047 2.549654e-14 6.374134e-12
## gene_18  -0.063925704    2.52670976 10.860340 61.91905 3.584671e-14 7.169342e-12
## gene_24   0.075454872    2.01656333 11.354622 61.44204 4.550209e-14 7.583681e-12
## gene_2   -0.451131078    2.06163750 11.514351 59.90225 9.826315e-14 1.403759e-11
## gene_43  -2.453106336   -2.19487534 13.115105 59.46353 1.223655e-13 1.465662e-11
## gene_50  -2.573146092   -2.51645998  9.866134 59.31332 1.319096e-13 1.465662e-11
## gene_28  -1.983041982   -1.88524513 10.359296 59.07953 1.482662e-13 1.482662e-11

実験群間の比較における発現変動遺伝子を gene_1 ～ gene_150 として設定してある。確かに、上位にランキングされている遺伝子はこの範囲に入っていることが確認できる。

解析結果を 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、B）

処理 A と処理 B を比較するとき、帰無仮説としては β₄ = 0 を考えればよい。すなわち full model から β₄ を取り除いたモデルを reduced model として尤度比検定を行えばよい。reduced model のデザイン行列は以下のようになる。

\[ \begin{pmatrix} E[\mbox{WT-A}] \\ E[\mbox{WT-A}] \\ E[\mbox{WT-B} \\ E[\mbox{WT-B}] \\ E[\mbox{C1-A}] \\ E[\mbox{C1-A}] \\ E[\mbox{C1-B}] \\ E[\mbox{C1-B}] \\ E[\mbox{C2-A}] \\ E[\mbox{C2-A}] \\ E[\mbox{C2-B}] \\ E[\mbox{C2-B}] \end{pmatrix}^{T} = \begin{pmatrix} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0\\ \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3} \\ \beta_{4} \end{pmatrix} \]

reduced model のデザイン行列は、full model のデザイン行列の 4 列目の要素をすべて 0 したデザイン行列であるから、coef = 4 と指定すればよい。

lrt <- glmLRT(fit, coef = 4)
topTags(lrt)
## Coefficient:  treatB 
##              logFC    logCPM       LR       PValue          FDR
## gene_166 -2.087503 10.555820 62.37882 2.833545e-15 1.250237e-12
## gene_165 -2.017893 13.362927 62.13994 3.198975e-15 1.250237e-12
## gene_188  2.192741 11.457769 61.50152 4.424006e-15 1.250237e-12
## gene_157 -2.051739 11.455796 61.00638 5.689010e-15 1.250237e-12
## gene_172 -2.095978 12.622683 60.82086 6.251185e-15 1.250237e-12
## gene_200  2.038794 12.656624 59.01711 1.563067e-14 2.605111e-12
## gene_196  2.432337 11.628599 56.76341 4.915275e-14 7.021821e-12
## gene_155 -2.147115 12.549064 50.08172 1.474745e-12 1.843431e-10
## gene_169 -2.289136 10.830430 48.88212 2.718177e-12 3.020196e-10
## gene_182  2.690705  9.094169 48.56092 3.201867e-12 3.085377e-10

処理 A と処理 B の比較で発現量に差が生じる遺伝子を gene_151 ～ gene_200 として設定してある。確かに上位にランキングされていることがわかる。

解析結果を 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 バイオスタティスティックス