多群間比較（DESeq2） |

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

R を起動して DESeq2 を呼び出す。呼び出しが完了したら DESeq2 のバージョンをチェックする。

library(DESeq2)
packageVersion("DESeq2")
##[1] ‘1.8.0’

サンプルデータ

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

遺伝子	A 群で高発現	B 群で高発現	C 群で高発現
gene_1 ～ gene_50	Yes	No	No
gene_51 ～ gene_100	No	Yes	No
gene_101 ～ gene_150	No	No	Yes
gene_151 ～ gene_200	Yes	Yes	No
gene_201 ～ gene_250	No	Yes	Yes
gene_251 ～ gene_300	Yes	No	Yes
gene_201 ～ gene_1000	No	No	No

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

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

dim(count)
## [1] 1000    9

head(count)
##         A1  A2   A3  B1  B2  B3  C1  C2  C3
## gene_1   3   2   19   0   0   0  13   0   0
## gene_2 680 807 1011 168 132 138 135 188 202
## gene_3  35  38   30  11   9  11   8   9   8
## gene_4 103 165  124  35  37  31  36  36  30
## gene_5   2   9    0   0   0   0   0   0   0
## gene_6 419 373  478 114 117 109  89 142 106

実験群の識別ラベルを作成し group 変数に代入する。DESeq2 ではデータフレーム型を要求するため、データフレーム型として作成する。

group <- data.frame(con = factor(c("A", "A", "A", "B", "B", "B", "C", "C", "C")))

Wald 検定

DESeq2 には Wald 検定と尤度比検定の 2 つの検定法が実装されている。Wald 検定は多群間比較のとき ANOVA のような解析が行えずに 2 群間ずつの比較となる。ANOVA のように群間全体で比較を行う場合は尤度比検定を利用する。

A 群と B 群の比較

results 関数を利用して結果を取り出すときに contrast 引数で比較したい群を指定する。ここでは group データフレームの中の con 列の A と B を比較したいので contrast = c("con" "A", "B") と指定する。

dds <- DESeqDataSetFromMatrix(countData = count, colData = group, design = ~ con)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)
res <- results(dds, contrast = c("con", "A", "B"))
head(res[order(res$pvalue), ])
## log2 fold change (MAP): con A vs B 
## Wald test p-value: con A vs B 
## DataFrame with 6 rows and 6 columns
##           baseMean log2FoldChange     lfcSE      stat       pvalue         padj
##          <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
## gene_74   490.4866      -1.934507 0.1637137 -11.81640 3.211618e-32 1.499406e-29
## gene_89   541.8498      -2.269198 0.1922041 -11.80619 3.626132e-32 1.499406e-29
## gene_92   269.8632      -1.989457 0.1733022 -11.47970 1.668638e-30 4.599879e-28
## gene_39   305.2514       1.940684 0.1738489  11.16305 6.183365e-29 1.278411e-26
## gene_206  240.6456      -2.046675 0.1968076 -10.39937 2.495677e-25 4.127850e-23
## gene_275  375.9168       1.881827 0.1819329  10.34352 4.477670e-25 6.171722e-23

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

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

A 群と C 群の比較

ここでは group データフレームの中の con 列の A と C を比較したいので contrast = c("con" "A", "C") と指定する。

dds <- DESeqDataSetFromMatrix(countData = count, colData = group, design = ~ con)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)
res <- results(dds, contrast = c("con", "A", "C"))
head(res[order(res$pvalue), ])
## log2 fold change (MAP): con A vs C 
## Wald test p-value: con A vs C 
## DataFrame with 6 rows and 6 columns
##           baseMean log2FoldChange     lfcSE      stat       pvalue         padj
##          <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
## gene_165 1431.8341       1.925548 0.1267020  15.19745 3.676737e-52 2.514888e-49
## gene_102  696.4436      -1.953657 0.1296664 -15.06680 2.677850e-51 9.158246e-49
## gene_172  856.5890       1.911530 0.1503935  12.71019 5.190625e-37 1.183463e-34
## gene_157  370.9141       1.986548 0.1677616  11.84150 2.381559e-32 4.072465e-30
## gene_111 3015.8000      -2.199249 0.1969138 -11.16858 5.809927e-29 7.947980e-27
## gene_181  189.8772       2.348046 0.2188444  10.72930 7.415972e-27 8.412140e-25

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

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

尤度比検定

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

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

ここで二群間二因子比較データをモデル化する。まず、各実験群の発現量の差を高々 3 つのパラメーターで表すことができる。すなわち、E[A] = β₁、E[B] = β₁ + β₂、E[C] = β₁ + β₃ とする。サンプルデータは 3 つの実験群で、各実験群において 3 つの biological replicate が含まれているため、これを行列の形で表すと以下のようになる。

\[ \begin{pmatrix} E[\mbox{A}]\\ E[\mbox{A}] \\ E[\mbox{A}] \\ E[\mbox{B}] \\ E[\mbox{B}] \\ E[\mbox{B}] \\ E[\mbox{C}] \\ E[\mbox{C}] \\ E[\mbox{C}] \end{pmatrix} ^{T} = \begin{pmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \\ 1 & 0 & 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 関数を利用すれば簡単に作成できる。

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

3 群間の発現量を同時に比較を行うときの検定について考える。このとき帰無仮説は「E[A] = E[B] かつ E[A] = E[C]」とすることが考えられる。尤度比検定は 2 つのモデルを比較することによって帰無仮説を検定する。そこで、full model ともう 1 つのモデルを考えて、full model とそのモデルの差を帰無仮説と同等にすればよい。

帰無仮説の「E[A] = E[B] かつ E[A] = E[C]」に着目したとき、「β₂ = 0 かつ β₃ = 0」と同等である。そこで、もう 1 つのモデルとして以下のようなモデルを考える。

\[ \begin{pmatrix} E[\mbox{A}]\\ E[\mbox{A}] \\ E[\mbox{A}] \\ E[\mbox{B}] \\ E[\mbox{B}] \\ E[\mbox{B}] \\ E[\mbox{C}] \\ E[\mbox{C}] \\ E[\mbox{C}] \end{pmatrix} ^{T} = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} \begin{pmatrix}\beta_{1} \end{pmatrix} \]

このモデルは full model に含まれているパラメーターのうち β₂ と β₃ を取り除いてできたモデルであることから、reduced model と呼ばれる。

full model と reduced model のそれぞれから尤度を計算し比較する。もし、両者の尤度に差異が認められた場合、full model と reduced model は相違異なるモデルと考えられる。full model と reduced model の差は β₂ と β₃ によって作られていることから、「β₂ ≠ 0 または β₃ ≠ 0」となり、帰無仮説である「β₂ = 0 かつ β₃ = 0」に矛盾する。すなわち、帰無仮説は間違いであり棄却すべきである。このとき、3 つの実験群の発現量の比較でいずれかに有意差を持つ。逆に、両者の尤度に差異が認められない場合は、full model と reduced model は同じモデルと考えられ、帰無仮説の判断は保留される。

reduced model のデザイン行列につていは作成できないので、1 列目だけを 1 とするという意味で ~ 1 で代用する。

ここで、DESeq2 を利用して尤度比検定を行う。nbinomLRT 関数を利用するとき、full には full model を作成する際に使用した式、reduced には reduced model を作成する際に使用した式を代入する。DESeqDataSetFromMatrix 関数を利用するときに colData 引数に group 変数が与えられているため、nbinomLRT では group$con のように入力する代わりに、省略して con を入力する。

dds <- DESeqDataSetFromMatrix(countData = count, colData = group, design = ~ con)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomLRT(dds, full = ~ con, reduced = ~ 1)
res <- results(dds)
head(res[order(res$pvalue), ])
## log2 fold change (MLE): group C vs A 
## LRT p-value: '~ group' vs '~ 1' 
## DataFrame with 6 rows and 6 columns
##           baseMean log2FoldChange     lfcSE      stat       pvalue         padj
##          <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
## gene_102  696.4436    1.963552928 0.1303548  337.3719 5.503236e-74 4.886873e-71
## gene_165 1431.8341   -1.935416127 0.1273601  294.6380 1.047556e-64 4.651147e-62
## gene_89   541.8498   -0.027051496 0.2003886  206.3506 1.554328e-45 4.600811e-43
## gene_74   490.4866    0.037043408 0.1708722  201.0331 2.219304e-44 4.926855e-42
## gene_172  856.5890   -1.925415429 0.1515039  197.0323 1.640500e-43 2.913529e-41
## gene_92   269.8632   -0.007576415 0.1855265  196.4542 2.190297e-43 3.241639e-41

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

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

References

Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014. 15(12):550. PubMed Abstract
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010, 11(10):R106. PubMed Abstract
一般化線形モデル. biopapyrus バイオスタティスティックス