DESeq2 を利用して対応ありの二群間比較データを解析する例を示す。検定方法としては Wald 検定と尤度比検定の 2 種類がある。
library(DESeq2)
packageVersion("DESeq2")
## [1] ‘1.8.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 変数に保存する。ただし、DESeq2 ではデータフレーム型を要求するため、これらをデータフレームにまとめる。
group <- data.frame(
cell = factor(c("N", "N", "N", "N", "T", "T", "T", "T")),
patient = factor(c("A", "B", "C", "D", "A", "B", "C", "D"))
)Wald 検定
まず DESeq2 を利用してリードカウントデータの分布などを計算し、Wald 統計量などを計算する。
dds <- DESeqDataSetFromMatrix(countData = count, colData = group,
design = ~ patient + cell)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)正常細胞とがん細胞を比較するのが目的で、この情報は group 変数中の cell 列に保存されている。そこで、contrast で group$cell 中の N と T を比較するように指定する。
res <- results(dds, contrast = c("cell", "N", "T"))
head(res[order(res$pvalue), ])
## log2 fold change (MAP): cell N vs T
## Wald test p-value: cell N vs T
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## gene_62 1651.1748 2.218449 0.2664193 8.326907 8.298197e-17 6.207051e-14
## gene_83 1003.6712 1.911836 0.2471800 7.734588 1.037383e-14 2.867491e-12
## gene_32 1216.4954 -2.577070 0.3337543 -7.721458 1.150063e-14 2.867491e-12
## gene_74 677.1926 1.857280 0.2494440 7.445681 9.644562e-14 1.803533e-11
## gene_43 1362.9932 -1.798362 0.2425891 -7.413203 1.232845e-13 1.844336e-11
## gene_2 461.7023 -1.890825 0.2722414 -6.945399 3.773919e-12 4.704819e-10サンプルデータは、gene_1 ～ gene_100 を正常細胞とがの比較で差が出るような遺伝子として生成している。結果を眺めると確かに gene_1 ～ gene_100 の遺伝子が上位にランキングされている。
解析結果をファイルに書き出す。
write.table(res, file = "result.txt", col.names = T, row.names = T, sep = "\t")尤度比検定
対応ありの二群間比較においては、患者間の比較に興味がなく、正常細胞とがん細胞のような比較が目的である。一般化線形モデルによりこれを解析を行う方法について説明する。
尤度比検定では 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 の正常細胞の発現量の期待値を β1 とする。患者 A と患者 B、C、D の発現量の期待値の差をそれぞれ β2、β3、β4 とおく。また、正常細胞とがん細胞の発現量の期待値の差を β5 とおく。これにより E[A-Normal] = β1、E[B-Normal] = β1 + β2、E[A-Tumor] = β1 + β5 などである。これらを行列の形で表すと以下のようになる。（N: Normal; T: Tumor）
この式が一般化線形モデルのモデルとされる。0 と 1 からなる行列はデザイン行列と呼び、観測値とパラメーターを関連させる重要な働きを持つ。遺伝子が n 個があればこのような式は n 個できる。遺伝子が n 個あるときプログラムが逐次的に処理を行ってくれるため、遺伝子が 1 個だけの時を考えればよい。
このモデルを full model と呼ぶことにする。full model のデザイン行列は model.matrix 関数を利用すれば簡単に作成できる。
model.matrix(~ group$patient + group$cell)
## (Intercept) group$patientB group$patientC group$patientD group$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")$`group$patient`
## [1] "contr.treatment"
##
## attr(,"contrasts")$`group$cell`
## [1] "contr.treatment"次に尤度比検定を行い発現変動遺伝子を検出する。尤度比検定では 2 つのモデルの尤度を比較して検定を行うために、上で作成した full model の他にもう 1 つのモデル reduced model を必要とする。
まず、目的は正常細胞とがん細胞の比較で差異のある遺伝子を検出したい。そこで帰無仮説としては正常細胞とがん細胞の発現量は同じで E[N] = E[T]、すなわち β5 = 0 とすることができる。β5 = 0 となるモデルを reduced model として、full model と比較し尤度比検定を行う。そこで、もし、両者の尤度に差異が認められるならば、full model と reduced model は相違異なるモデルと考えられる。full model と reduced model の差は β5 によってもたらしているから、両者に差異があるということは β5 ≠ 0 を意味する。しかし、これは帰無仮説である β5 = 0 に矛盾する。よって、帰無仮説は間違いであり、棄却するべきである。逆に、両者の尤度に差異がなければ帰無仮説は保留される。これが尤度比検定である。
reduced model は以下のように作成できる。
R では以下のように作成できる。full model に比べ 1 列少なくなっている。DESeq2 を利用している時、デザイン行列の列名で判断しているため、以下のデザイン行列には group$cellT が含まれていないため、それに対応している β5 が自動的考慮されなくなる。パラメーターがずれることはない。
model.matrix(~ group$patient)
## (Intercept) group$patientB group$patientC group$patientD
## 1 1 0 0 0
## 2 1 1 0 0
## 3 1 0 1 0
## 4 1 0 0 1
## 5 1 0 0 0
## 6 1 1 0 0
## 7 1 0 1 0
## 8 1 0 0 1
## attr(,"assign")
## [1] 0 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$`group$patient`
## [1] "contr.treatment"ここで、full model は ~ group$patient + group$cell、reduced model は ~ group$patient で作成できることがわかった。そこで、DESeq2 を利用してこの 2 つのモデルを比較し、両モデルに差をもたらす遺伝子を検出する。コードは以下のように実行する。nbinomLRT 関数を実行するときに full model と reduced model の式を代入する。ただし、このとき、group 変数は一番最初に colData = group と代入してあるため、group$cell、group$patient はそれぞれ略してそれぞれ cell、patient とだけにする。
dds <- DESeqDataSetFromMatrix(countData = count, colData = group,
design = ~ patient + cell)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomLRT(dds, full = ~ patient + cell, reduced = ~ patient)
res <- results(dds)
head(res[order(res$pvalue), ])
## log2 fold change (MLE): cell T vs N
## LRT p-value: '~ patient + cell' vs '~ patient'
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## gene_62 1651.1748 -2.288306 0.2749449 62.71130 2.393389e-15 1.670586e-12
## gene_83 1003.6712 -1.964825 0.2539405 55.85083 7.818385e-14 2.728616e-11
## gene_74 677.1926 -1.908092 0.2564061 51.96997 5.635600e-13 9.946906e-11
## gene_32 1216.4954 2.729419 0.3512236 51.94555 5.706134e-13 9.946906e-11
## gene_43 1362.9932 1.843879 0.2489539 51.50946 7.125291e-13 9.946906e-11
## gene_2 461.7023 1.953760 0.2814114 45.00465 1.965667e-11 2.286726e-09正常細胞とがん細胞の比較における発現変動遺伝子は gene_1 ～ gene_100 に設定してある。解析結果を見ると上位にランキングされている遺伝子はすべてこの範囲に入っている。
解析結果をファイルに保存する場合は以下のようにする。
write.table(res, file = "result.txt", col.names = T, row.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 バイオスタティスティックス