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

edgeR パッケージは DESeq / DESeq2 パッケージと並んで発現変動遺伝子を検出する際によく利用されている。edgeR パッケージでは一般化線形モデルによる解析をサポートしているため、適当にモデルを作成すれば、二群間二因子比較も行える。

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

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

サンプルデータ
発現変動遺伝子の検出例

サンプルデータ

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

遺伝子	WT 群と KO 群の比較	A 処理と B 処理の比較
gene_1 ～ gene_50	WT 群が高発現	発現量に差がない
gene_51 ～ gene_100	KO 群で高発現	発現量に差がない
gene_101 ～ gene_150	発現量に差がない	A 処理で高発現
gene_151 ～ gene_200	発現量に差がない	B 処理で高発現
gene_201 ～ gene_1000	発現量に差がない	発現量に差がない

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

count <- read.table("https://bi.biopapyrus.jp/[[~604]", sep = "\t", header = T, row.names = 1)
count <- as.matrix(count)

dim(count)
## [1] 1000    8

head(count)
##        WT1_A WT2_A WT3_B WT4_B KO1_A KO2_A KO3_B KO4_B
## gene_1     1     7     0     3    10    14     0     0
## gene_2   170   182   109   124   468   832   513   417
## gene_3     8    12     6    13    54    41    40    35
## gene_4    44    38    46    74   177   221   153   130
## gene_5     0     0     0     5     2     3     0     2
## gene_6    88   112   111   153   497   270   563   344

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

group = factor(c("WT", "WT", "WT", "WT", "KO", "KO", "KO", "KO"))
treat = factor(c("A", "A", "B", "B", "A", "A", "B", "B"))

尤度比検定による二群間二因子比較の解析

二群間二因子比較の解析において、興味のある遺伝子は次の 2 通りに絞られると思われる。この両方の場合について見ていく。

WT 群と KO 群を比較した時に、発現量に差が生じる遺伝子
A 処理と B 処理を比較した時に、発現量に差が生じる遺伝子

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

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

ここで二群間二因子比較データをモデル化する。まず、各実験群と処理群の組み合わせによる観測できる発現量の期待値は高々 4 通りである。すなわち、E[WT-A]、E[WT-B]、E[KO-A]、E[KO-B] の 4 通りである。この 4 通りを表すにはパラーメーターを 3 つ必要とする。そこで、3 つのパラメーター β₁、β₂、β₃ を考える。

まず WT 群 A 処理の発現量を基準し、E[WT-A] = β₁ とおく。次に WT 群と KO 群の発現量の差を β₂、A 処理と B 処理の発現量の差を β₃ とすると、E[KO-A] = β₁ + β₂、E[WT-B] = β₁ + β₃、E[KO-B] = β₁ + β₂ + β₃ である。これを行列の形で表すと以下のようになる。

\[ \begin{pmatrix} E[\mbox{WT-A}] \\ E[\mbox{WT-A}] \\ E[\mbox{WT-B}] \\ E[\mbox{WT-B}] \\ E[\mbox{KO-A}] \\ E[\mbox{KO-A}] \\ E[\mbox{KO-B}] \\ E[\mbox{KO-B}] \end{pmatrix} ^{T} = \begin{pmatrix}1 & 0 & 0\\ 1 & 0 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \\ 1 & 1 & 0\\ 1 & 1 & 1 \\ 1 & 1 & 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 関数を利用すれば簡単に作成できる。（groupWT の列は上の行列と対称になっているが、解析にあたって効果は同じである。）

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

edgeR を利用して、データを正規化、分布をパラメーターの推定を行う。DGEList 関数には識別ラベルとして group を与えた。処理の識別ラベルである treat については考慮しなくてもよい。

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 つの方のモデルは、どのような発現変動遺伝子を期待するかによって作り方が異なる。

WT 群と KO 群の比較

（A 処理と B 処理を考慮しないで、）WT 群と KO 群を比較した際に発現量に差をもつ遺伝子を発現変動遺伝子として検出する例を示す。WT 群と KO 群の発現量を比較するための帰無仮説を考える。帰無仮説は「WT 群と KO 群の発現量は同じ」とすることができる。そこで、full model ともう 1 つのモデルをどのように作成し比べれば帰無仮説を表せるかについて考える。WT 群と KO 群の発現量の差はパラメーター β₂ によって支配されている。帰無仮説が成り立つとき、β₂ = 0 も成り立つ。そこで、帰無仮説を言い換えれば「β₂ = 0」である。例えば、β₂ = 0 に対応できるように full model の 2 列目の要素をすべて 0 に置換したモデルを考える。

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

このモデルは、full model からパラメーター β₂ を取り除いてできたモデルであり、full model に対して reduced model と呼ぶことにする。full model と reduced model のそれぞれから尤度を計算し比較する。もし両者の尤度に差異が認められた場合、full model と reduced モデルは相違異なるモデルと判定できる。full model と reduced model の差は β₂ だけであるから、β₂ ≠ 0 であることを意味する。しかし、これは最初に仮定した β₂ = 0 （帰無仮説）に対し矛盾である。したがって、最初の仮説（帰無仮説）が間違いであり、棄却するべきである。このとき、E[WT-A] ≠ E[KO-A] または E[WT-B] ≠ E[KO-B] が成り立つ。逆にもし両者の尤度に差異が認められない場合は β₂ = 0 となり、帰無仮説を棄却できなくなる。このとき、「E[WT-A] ≠ E[KO-A] または E[WT-B] ≠ E[KO-B]」の判定が保留される（判断できなくなる）。

ここで edgeR による解析に話を戻す。次のステップとして 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:  groupWT 
##             logFC   logCPM       LR       PValue          FDR
## gene_62  2.195437 13.89505 55.44743 9.599208e-14 9.599208e-11
## gene_34 -2.779385 10.46538 53.84977 2.164209e-13 1.082104e-10
## gene_74  2.211859 12.46248 50.33813 1.294103e-12 3.241939e-10
## gene_50 -2.695960 10.07772 50.33408 1.296776e-12 3.241939e-10
## gene_89  2.388014 12.14691 44.06377 3.178494e-11 6.356989e-09
## gene_18 -2.584470 10.97108 43.61257 4.002596e-11 6.670994e-09
## gene_92  2.127916 11.58842 41.94585 9.383642e-11 1.201230e-08
## gene_16 -2.156634 12.96528 41.83333 9.939462e-11 1.201230e-08
## gene_39 -2.032617 11.82491 41.66898 1.081107e-10 1.201230e-08
## gene_82  2.417651 10.00789 40.68689 1.786819e-10 1.786819e-08

WT 群と KO 群の比較における発現変動遺伝子は gene_1 ～ gene_100 に設定してある。解析結果を見ると上位にランキングされている遺伝子はすべてこの範囲に入っている。

デフォルトでは、topTags は最初の 10 遺伝子の結果をしか表示しない。そこで、すべての結果を表示させるためには引数 n　を利用する。以下は、すべての結果を result.txt ファイルに書き出す例である。

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

A 処理と B 処理の比較

A 処理と B 処理の差を比較するときに考える帰無仮説とし、A 処理と B 処理に差がないとすればよい。すなわち、β₃ = 0 を考えればよい。これは full model のデザイン行列の 3 列目の要素をすべて 0 に変換すれば対応できる。すなわち、

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

full model に対して、β₃ を取り除いてできたモデルを reduced model という。上と同様に、検定するには coef = 3 と指定する。

lrt <- glmLRT(fit, coef = 3)
topTags(lrt)
## Coefficient:  treatB 
##              logFC    logCPM       LR       PValue          FDR
## gene_108  2.240780 14.168494 63.47805 1.621636e-15 1.621636e-12
## gene_195 -2.305379 14.929074 58.78500 1.758770e-14 8.793851e-12
## gene_125  3.996920  8.835907 52.87263 3.558950e-13 1.186317e-10
## gene_165 -1.989617 13.352630 47.17433 6.494548e-12 1.623637e-09
## gene_131  2.441912  9.869362 46.40715 9.606435e-12 1.872814e-09
## gene_102  1.966615 12.960435 45.83201 1.288411e-11 1.872814e-09
## gene_111  1.897182 14.713250 45.79801 1.310970e-11 1.872814e-09
## gene_200 -2.015206 12.620309 44.38606 2.696006e-11 3.370008e-09
## gene_189 -2.447335 11.175022 43.06084 5.306372e-11 5.862047e-09
## gene_188 -2.209175 11.541669 42.86601 5.862047e-11 5.862047e-09

A 処理と B 処理の比較における発現変動遺伝子は gene_101 ～ 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 バイオスタティスティックス