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("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 <- 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]) = Xβ と表すことができる。ただし、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 とすれば、各実験群と処理群を組み合わせにおいて、その発現量は行列を用いて以下のように表すことができる。
この式が一般化線形モデルのモデルとされる。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 通りに絞られると思われる。この両方の場合について見ていく。
- WT 群、C1 群と C2 群を比較した時に、発現量に差が生じる遺伝子
- 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 に置換する。
このように用意した 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$con と group$con は略してそれぞれ con と treat とだけにする。
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 のデザイン行列は以下のようになる。
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
- Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014. 15(12):550. PubMed Abstract
- Differential expression analysis for sequence count data. Genome Biol. 2010, 11(10):R106. PubMed Abstract
- 一般化線形モデル. biopapyrus バイオスタティスティックス