多群間二因子比較(DESeq2)

DESeq2 は RNA-seq のリードカウントデータから発現変動遺伝子を検出するためのパッケージである。一般化線形モデルをサポートしているため複雑な解析にも柔軟に対応できる。Wald 検定と尤度比検定の両方が行える。ただし、Wald 検定は 2 群間ずつ比較するときのみ利用できる。このページでは複雑なモデルにも適用できる尤度比検定を取り扱う。

library(DESeq2)
packageVersion("DESeq2")
##[1] ‘1.8.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("http://bi.biopapyrus.net/data/count3gf.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 <- data.frame(
  con   = 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]) = と表すことができる。ただし、g はリンク関数である。リンク関数は期待値 E[Y] と右辺が等しくなるように変換を行う単調な関数であり、しかも、その形は決まっている。例えば、データの分布は正規分布ならば g(x) = x、負の二項分布ならば g(x) = log(x) である。そこで、以下、説明を簡略するためにリンク関数を省略して書く。

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

\[ \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 列目に対応していることを注意しておく。)

model.matrix(~ group$con + group$treat)
##    (Intercept) group$conC2 group$conWT group$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$con`
## [1] "contr.treatment"
## 
## attr(,"contrasts")$`group$treat`
## [1] "contr.treatment"

次に尤度比検定を行い発現変動遺伝子を検出する。尤度比検定では 2 つのモデルの尤度を比較して検定を行うために、上で作成した full model の他にもう 1 つのモデル reduced model を必要とする。reduced model は、どのような発現変動遺伝子を期待するかによって作り方が異なる。三群間二因子比較の解析において、興味のある遺伝子は次の 2 通りに絞られると思われる。この両方の場合について見ていく。

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

実験群間の比較(WT、C1、C2)

処理 A と処理 B の影響を考慮しないで、実験群同士の比較において差を持つ遺伝子を発現変動遺伝子として検出する方法を示す。このような比較において帰無仮説としては E[WT] = E[C1] = E[C2] が考えられる。これは β2 = β3 = 0 に対応している。そこで、「β2 = 0 かつ β3 = 0」に対応しているモデルを作成する。そこで、full model と reduced model を比べた時、「β2 = 0 かつ β3 = 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\\ 1 & 0\\ 1 & 1\\ 1 & 1\\ 1 & 0\\ 1 & 0\\ 1 & 1\\ 1 & 1\\ 1 & 0\\ 1 & 0\\ 1 & 1\\ 1 & 1\\ \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{4} \end{pmatrix} \]

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

reduced model のデザイン行列は以下のように作成できる。reduced model の第 2 列の名前は group$treatB と表示されている、これは full model の 第 4 列の名前に対応している。このように、列名でマッチングを行えるので、列が減ることによってパラメーターがずれたりはしない。

model.matrix(~ group$treat)
##    (Intercept) group$treatB
## 1            1            0
## 2            1            0
## 3            1            1
## 4            1            1
## 5            1            0
## 6            1            0
## 7            1            1
## 8            1            1
## 9            1            0
## 10           1            0
## 11           1            1
## 12           1            1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$`group$treat`
## [1] "contr.treatment"

以上まとめると full model は ~ group$con + group$treat、reduced model は ~ group$treat で作成できることがわかった。そこで、DESeq2 を利用して尤度比検定を行う。

コードは以下のように実行する。nbinomLRT 関数を実行するときに full model と reduced model の式を代入する。ただし、このとき、group 変数は一番最初に colData = group と代入してあるため、group$congroup$con は略してそれぞれ contreat とだけにする。

dds <- DESeqDataSetFromMatrix(countData = count, colData = group,
                              design = ~ con + treat)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomLRT(dds, full = ~ con + treat, reduced = ~ treat)
res <- results(dds)
head(res[order(res$pvalue), ])
## log2 fold change (MLE): treat B vs A 
## LRT p-value: '~ con + treat' vs '~ treat' 
## DataFrame with 6 rows and 6 columns
##           baseMean log2FoldChange      lfcSE      stat        pvalue          padj
##          <numeric>      <numeric>  <numeric> <numeric>     <numeric>     <numeric>
## gene_102 1041.3360   -0.004261709 0.06363297  699.4603 1.300566e-152 1.170510e-149
## gene_92   398.9559    0.003062482 0.09380744  341.7097  6.290444e-75  2.830700e-72
## gene_114  716.1902    0.017223690 0.09283911  296.5558  4.015399e-65  1.204620e-62
## gene_22   214.7942    0.138746165 0.13238201  289.9580  1.087459e-63  2.109909e-61
## gene_24   291.8590    0.157501164 0.11740665  289.8080  1.172171e-63  2.109909e-61
## gene_136  433.6175   -0.049194473 0.10302847  271.9067  9.040711e-60  1.356107e-57

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

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

write.table(res, file = "result.txt", row.names = T, col.names = T, sep = "\t")

処理間の比較(A、B)

処理 A と処理 B を比較するとき、帰無仮説としては β4 = 0 を考えればよい。すなわち full model から β4 を取り除いたモデルを 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 \\ 1 & 0 & 0\\ 1 & 0 & 0\\ 1 & 0 & 0\\ 1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 0 & 1\\ 1 & 0 & 1\\ 1 & 0 &1 \\ 1 & 0 & 1\\ \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3} \end{pmatrix} \]

reduced model のデザイン行列は以下のように作成できる。

model.matrix(~ group$con)
##    (Intercept) group$conC2 group$conWT
## 1            1           0           1
## 2            1           0           1
## 3            1           0           1
## 4            1           0           1
## 5            1           0           0
## 6            1           0           0
## 7            1           0           0
## 8            1           0           0
## 9            1           1           0
## 10           1           1           0
## 11           1           1           0
## 12           1           1           0
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$`group$con`
## [1] "contr.treatment"

従って、DESeq2 を利用して発現変動遺伝子を検定するには以下のようにすればよい。

dds <- DESeqDataSetFromMatrix(countData = count, colData = group,
                              design = ~ con + treat)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomLRT(dds, full = ~ con + treat, reduced = ~ con)
res <- results(dds)
head(res[order(res$pvalue), ])
## log2 fold change (MLE): treat B vs A 
## LRT p-value: '~ con + treat' vs '~ con' 
## DataFrame with 6 rows and 6 columns
##           baseMean log2FoldChange      lfcSE      stat        pvalue          padj
##          <numeric>      <numeric>  <numeric> <numeric>     <numeric>     <numeric>
## gene_165 1173.2790      -1.989978 0.07738812  633.3891 9.156571e-140 7.325257e-137
## gene_157  311.4693      -2.021377 0.14618773  183.0762  1.032277e-41  4.129109e-39
## gene_192  151.7660       1.791441 0.13451421  178.4605  1.050922e-40  2.802458e-38
## gene_200  727.4578       2.068931 0.15452852  167.1306  3.132572e-38  6.265145e-36
## gene_166  166.0653      -2.056362 0.16139928  158.5830  2.308110e-36  3.692976e-34
## gene_172  701.6266      -2.066710 0.15915193  156.8271  5.583698e-36  7.444931e-34

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

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

write.table(res, file = "result.txt", row.names = T, col.names = T, sep = "\t")

References