対応あり二群間比較（edgeR） |

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

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

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

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

サンプルデータ

サンプルデータは、カウントデータが負の二項分布に従うように乱数を生成して作成した。例えるならば、4 人の患者A、B、C、D から正常細胞とがん細胞のサンプルを採取し、正常細胞とがん細胞の比較を行う。サンプルデータをこのような構成で生成した。全遺伝子数は 1,000 個であり、その発現パターンは以下のように設定した。

遺伝子	正常細胞	がん細胞
gene_1 ～ gene_50	がん細胞に比べて高発現
gene_51 ～ gene_100		正常細胞に比べて高発現
gene_101 ～ gene_1000	発現量に差がない

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

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

dim(count)
## [1] 1000    8

head(count)
##        N_A N_B N_C N_D T_A T_B T_C T_D
## gene_1   1   1   0   2   4   1  54   0
## gene_2 170 212 187 195 745 630 828 807
## gene_3   8  10  11  10  39  41  45  37
## gene_4  44  29  40  49 196 208 129 158
## gene_5   0   2   0   0   0   0   1   0
## gene_6  88  95 141 123 555 270 340 561

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

cell = factor(c("N", "N", "N", "N", "T", "T", "T", "T"))
patient = factor(c("A", "B", "C", "D", "A", "B", "C", "D"))

一般化線形モデルによる二群間二因子比較の解析

対応ありの二群間比較においては、患者間の比較に興味がなく、正常細胞とがん細胞のような比較が目的である。一般化線形モデルによりこれを解析を行う方法について説明する。

尤度比検定では 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 について考える。患者 A の正常細胞の発現量の期待値を β₁ とする。患者 A と患者 B、C、D の発現量の期待値の差をそれぞれ β₂、β₃、β₄ とおく。また、正常細胞とがん細胞の発現量の期待値の差を β₅ とおく。これにより E[A-Normal] = β₁、E[B-Normal] = β₁ + β₂、E[A-Tumor] = β₁ + β₅ などである。これらを行列の形で表すと以下のようになる。（N: Normal; T: Tumor）

\[ \begin{pmatrix} E[\mbox{NA}] \\ E[\mbox{NB}] \\ E[\mbox{NC}] \\ E[\mbox{ND}] \\ E[\mbox{TA}] \\ E[\mbox{TB}] \\ E[\mbox{TC}] \\ E[\mbox{TD}] \end{pmatrix} ^{T} = \begin{pmatrix}1 & 0 & 0 & 0 & 0\\ 1 & 1 & 0 & 0&0\\ 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 &1 \\ 1 & 1 & 0 & 0 & 1\\ 1 & 0 & 1 & 0 & 1\\ 1 & 0 &0 &1 & 1 \end{pmatrix}\begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3} \\ \beta_{4} \\ \beta_{5}\end{pmatrix} \]

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

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

design <- model.matrix(~ patient + cell)
design
##   (Intercept) patientB patientC patientD cellT
## 1           1        0        0        0     0
## 2           1        1        0        0     0
## 3           1        0        1        0     0
## 4           1        0        0        1     0
## 5           1        0        0        0     1
## 6           1        1        0        0     1
## 7           1        0        1        0     1
## 8           1        0        0        1     1
## attr(,"assign")
## [1] 0 1 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$patient
## [1] "contr.treatment"
## 
## attr(,"contrasts")$cell
## [1] "contr.treatment"

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

d <- DGEList(counts = count, group = cell)
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 を必要とする。

まず、目的は正常細胞とがん細胞の比較で差異のある遺伝子を検出したい。そこで帰無仮説としては正常細胞とがん細胞の発現量は同じで E[N] = E[T]、すなわち β₅ = 0 とすることができる。β₅ = 0 となるモデルを reduced model として、full model と比較し尤度比検定を行う。そこで、もし、両者の尤度に差異が認められるならば、full model と reduced model は相違異なるモデルと考えられる。full model と reduced model の差は β₅ によってもたらしているから、両者に差異があるということは β₅ ≠ 0 を意味する。しかし、これは帰無仮説である β₅ = 0 に矛盾する。よって、帰無仮説は間違いであり、棄却するべきである。逆に、両者の尤度に差異がなければ帰無仮説は保留される。これが尤度比検定である。

reduced model は以下のように作成できる。

\[ \begin{pmatrix} E[\mbox{NA}] \\ E[\mbox{NB}] \\ E[\mbox{NC}] \\ E[\mbox{ND}] \\ E[\mbox{TA}] \\ E[\mbox{TB}] \\ E[\mbox{TC}] \\ E[\mbox{TD}] \end{pmatrix} ^{T} = \begin{pmatrix}1 & 0 & 0 & 0 & 0\\ 1 & 1 & 0 & 0&0\\ 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 &0 \\ 1 & 1 & 0 & 0 & 0\\ 1 & 0 & 1 & 0 & 0\\ 1 & 0 &0 &1 & 0 \end{pmatrix}\begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3} \\ \beta_{4} \\ \beta_{5}\end{pmatrix} \]

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

lrt <- glmLRT(fit, coef = 5)

結果をテーブル形式で表示するには topTags 関数を利用する。

topTags(lrt)
## Coefficient:  cellT
##             logFC    logCPM       LR       PValue          FDR
## gene_32  2.740074 13.567081 60.27103 8.265514e-15 4.327905e-12
## gene_62 -2.277109 14.002016 60.18021 8.655810e-15 4.327905e-12
## gene_95 -8.271141  7.723015 54.78282 1.346121e-13 4.487069e-11
## gene_16  2.342469 13.073516 43.08342 5.245480e-11 1.311370e-08
## gene_43  1.854217 13.723760 41.57287 1.135582e-10 2.271163e-08
## gene_83 -1.953828 13.281703 40.46826 1.998360e-10 3.330600e-08
## gene_74 -1.896861 12.715743 34.10152 5.231051e-09 7.472929e-07
## gene_9   2.114943 13.924096 33.47930 7.202644e-09 9.003305e-07
## gene_48  2.432970  9.978364 33.15151 8.525012e-09 9.472235e-07
## gene_46  4.386886  8.950959 32.92403 9.583118e-09 9.583118e-07

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

デフォルトでは、topTags は最初の 10 遺伝子の結果をしか表示しない。そこで、すべての結果を表示させるためには引数 n を利用する。以下は、すべての結果を 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 バイオスタティスティックス