二群間二因子比較（DESeq2） |

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

R を起動して DESeq2 を呼び出す。

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

サンプルデータ

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

遺伝子	WT 群と KO 群の比較	A 処理と B 処理の比較
gene_1 ～ gene_50	WT 群が高発現	発現量に差がない
gene_51 ～ gene_100	KO 群で高発現	発現量に差がない
gene_101 ～ gene_150	発現量に差がない	A 処理で高発現
gene_151 ～ gene_200	発現量に差がない	B 処理で高発現
gene_201 ～ gene_1000	発現量に差がない	発現量に差がない

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

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

dim(count)
## [1] 1000    8

head(count)
##        WT1_A WT2_A WT3_B WT4_B KO1_A KO2_A KO3_B KO4_B
## gene_1     1     7     0     3    10    14     0     0
## gene_2   170   182   109   124   468   832   513   417
## gene_3     8    12     6    13    54    41    40    35
## gene_4    44    38    46    74   177   221   153   130
## gene_5     0     0     0     5     2     3     0     2
## gene_6    88   112   111   153   497   270   563   344

実験群の識別ラベルを作成し con 変数に代入する。また、処理の識別ラベルを作成し treat 変数に保存する。ただし、DESeq2 ではデータフレーム型を要求するため、これらをデータフレームにまとめる。

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

Wald 検定

WT 群と KO 群の比較

group 変数には WT 群と KO 群のデータが con 列に入っている。そこで、contrast 引数を与えて、WT 群と KO 群の比較結果を取得する。

dds <- DESeqDataSetFromMatrix(countData = count, colData = group,
                              design = ~ con + treat)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)
res <- results(dds, contrast = c("con", "WT", "KO"))
head(res[order(res$pvalue), ])
## log2 fold change (MAP): con WT vs KO 
## Wald test p-value: con WT vs KO 
## DataFrame with 6 rows and 6 columns
##          baseMean log2FoldChange     lfcSE      stat       pvalue         padj
##         <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
## gene_39  403.3604      -1.997640 0.1661023 -12.02657 2.576590e-33 2.437454e-30
## gene_62 1710.3740       2.182295 0.1824834  11.95887 5.834972e-33 2.759942e-30
## gene_92  343.9807       2.111549 0.1825877  11.56457 6.229944e-31 1.964509e-28
## gene_74  632.9787       2.188305 0.2036623  10.74477 6.271960e-27 1.483318e-24
## gene_24  327.5729      -1.951868 0.1903277 -10.25530 1.120278e-24 2.119566e-22
## gene_34  155.7111      -2.682644 0.2807417  -9.55556 1.229173e-21 1.937997e-19

サンプルデータは、gene_1 ～ gene_100 を WT 群と KO 群の比較で差が出るような遺伝子として生成している。結果を眺めると確かに gene_1 ～ gene_100 の遺伝子が上位にランキングされている。

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

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

A 処理と B 処理の比較

A 処理と B 処理のラベルは group 変数の treat 列に入っているため、contrast を以下のようにして結果を取得する。

dds <- DESeqDataSetFromMatrix(countData = count, colData = group,
                              design = ~ con + treat)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)
res <- results(dds, contrast = c("treat", "A", "B"))
head(res[order(res$pvalue), ])
## log2 fold change (MAP): treat A vs B 
## Wald test p-value: treat A vs B 
## DataFrame with 6 rows and 6 columns
##           baseMean log2FoldChange     lfcSE      stat       pvalue         padj
##          <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
## gene_108 2065.6536      -2.239121 0.1200688 -18.64866 1.294884e-77 9.672785e-75
## gene_102  892.9625      -1.966384 0.1186304 -16.57572 1.043989e-61 3.899300e-59
## gene_165 1167.0530       1.964030 0.1352186  14.52486 8.431684e-48 2.099489e-45
## gene_136  357.4408      -2.076962 0.1726086 -12.03278 2.389674e-33 4.462716e-31
## gene_200  702.0213       1.981711 0.1739028  11.39551 4.402698e-30 6.577630e-28
## gene_111 3009.9793      -1.881219 0.1755862 -10.71394 8.755860e-27 1.090105e-24

サンプルデータは、gene_101 ～ gene_200 を A 処理と B 処理の比較で差が出るような遺伝子として生成している。結果を眺めると確かに gene_101 ～ gene_200 の遺伝子が上位にランキングされている。

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

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

尤度比検定

二群間二因子比較の解析において、興味のある遺伝子は次の 2 通りに絞られると思われる。この両方の場合について見ていく。

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

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

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

ここで二群間二因子比較データをモデル化する。まず、各実験群と処理群の組み合わせによる観測できる発現量の期待値は高々 4 通りである。すなわち、E[WT-A]、E[WT-B]、E[KO-A]、E[KO-B] の 4 通りである。この 4 通りを表すにはパラーメーターを 3 つ必要とする。そこで、3 つのパラメーター β₁、β₂、β₃ を考える。

まず WT 群 A 処理の発現量を基準し、E[WT-A] = β₁ とおく。次に WT 群と KO 群の発現量の差を β₂、A 処理と B 処理の発現量の差を β₃ とすると、E[KO-A] = β₁ + β₂、E[WT-B] = β₁ + β₃、E[KO-B] = β₁ + β₂ + β₃ である。これを行列の形で表すと以下のようになる。

\[ \begin{pmatrix} E[\mbox{WT-A}] \\ E[\mbox{WT-A}] \\ E[\mbox{WT-B}] \\ E[\mbox{WT-B}] \\ E[\mbox{KO-A}] \\ E[\mbox{KO-A}] \\ E[\mbox{KO-B}] \\ E[\mbox{KO-B}] \end{pmatrix} ^{T} = \begin{pmatrix}1 & 0 & 0\\ 1 & 0 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \\ 1 & 1 & 0\\ 1 & 1 & 1 \\ 1 & 1 & 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 関数を利用すれば簡単に作成できる。（groupWT の列は上の行列と対称になっているが、解析にあたって効果は同じである。）

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

次に尤度比検定を行い発現変動遺伝子を検出する。尤度比検定では 2 つのモデルの尤度を比較して検定を行うために、上で作成した full model の他にもう 1 つのモデルを必要とする。そのもう 1 つの方のモデルは、どのような発現変動遺伝子を期待するかによって作り方が異なる。

WT 群と KO 群の比較

（A 処理と B 処理を考慮しないで、）WT 群と KO 群を比較した際に発現量に差をもつ遺伝子を発現変動遺伝子として検出する例を示す。WT 群と KO 群の発現量を比較するための帰無仮説を考える。帰無仮説は「WT 群と KO 群の発現量は同じ」とすることができる。そこで、full model ともう 1 つのモデルをどのように作成し比べれば帰無仮説を表せるかについて考える。WT 群と KO 群の発現量の差はパラメーター β₂ によって支配されている。帰無仮説が成り立つとき、β₂ = 0 も成り立つ。そこで、帰無仮説を言い換えれば「β₂ = 0」である。例えば、β₂ = 0 に対応できるように full model の 2 列目を取り除いてできたモデルを考える。

\[ \begin{pmatrix} E[\mbox{WT-A}] \\ E[\mbox{WT-A}] \\ E[\mbox{WT-B}] \\ E[\mbox{WT-B}] \\ E[\mbox{KO-A}] \\ E[\mbox{KO-A}] \\ E[\mbox{KO-B}] \\ E[\mbox{KO-B}] \end{pmatrix} ^{T} = \begin{pmatrix}1 & 0\\ 1 & 0 \\ 1 & 1 \\ 1 & 1 \\ 1 & 0 \\ 1 & 0\\ 1 & 1 \\ 1 & 1 \end{pmatrix}\begin{pmatrix}\beta_{1} \\ \beta_{3} \end{pmatrix} \]

このモデルは、full model からパラメーター β₂ を取り除いてできたモデルであり、full model に対して reduced model と呼ぶことにする。full model と reduced model のそれぞれから尤度を計算し比較する。もし両者の尤度に差異が認められた場合、full model と reduced モデルは相違異なるモデルと判定できる。full model と reduced model の差は β₂ だけであるから、β₂ ≠ 0 であることを意味する。しかし、これは最初に仮定した β₂ = 0 （帰無仮説）に対し矛盾である。したがって、最初の仮説（帰無仮説）が間違いであり、棄却するべきである。このとき、E[WT-A] ≠ E[KO-A] または E[WT-B] ≠ E[KO-B] が成り立つ。逆にもし両者の尤度に差異が認められない場合は β₂ = 0 となり、帰無仮説を棄却できなくなる。このとき、「E[WT-A] ≠ E[KO-A] または E[WT-B] ≠ E[KO-B]」の判定が保留される（判断できなくなる）。

reduced model は以下のように作成できる。R ではデザイン行列の列名で各パラメーターと対応付けを行っている。例えば reduced model の 2 列目は group$treatB となっているので、これは full model の 3 列目に対応している、と自動的に判断される。そのため、列数がずれることによってパラメーターがずれることはない。

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
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$`group$treat`
## [1] "contr.treatment"

ここで DESeq2 による解析に話を戻す。以上まとめると 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_39  403.3604   -0.016131061 0.1673584 138.13621 6.804579e-32 5.756674e-29
## gene_62 1710.3740   -0.234348264 0.1849182 131.38693 2.037555e-30 8.618857e-28
## gene_92  343.9807   -0.042413729 0.1843977 126.71295 2.146840e-29 6.054087e-27
## gene_74  632.9787    0.205652006 0.2069727 107.04399 4.353962e-25 9.208630e-23
## gene_24  327.5729   -0.013511205 0.1925925 100.33911 1.284156e-23 2.172791e-21
## gene_28  185.5780   -0.007838902 0.1889795  87.32559 9.204898e-21 1.297891e-18

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

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

A 処理と B 処理の比較

A 処理と B 処理の差を比較するときに考える帰無仮説とし、A 処理と B 処理に差がないとすればよい。full model から 3 列目を取り除いてできたモデルを reduced model とすればよい。

\[ \begin{pmatrix} E[\mbox{WT-A}] \\ E[\mbox{WT-A}] \\ E[\mbox{WT-B}] \\ E[\mbox{WT-B}] \\ E[\mbox{KO-A}] \\ E[\mbox{KO-A}] \\ E[\mbox{KO-B}] \\ E[\mbox{KO-B}] \end{pmatrix} ^{T} =\begin{pmatrix}1 & 0\\ 1 & 0 \\ 1 & 0 \\ 1 & 0 \\ 1 & 1 \\ 1 & 1\\ 1 & 1 \\ 1 & 1 \end{pmatrix}\begin{pmatrix}\beta_{1} \\ \beta_{2} \end{pmatrix} \]

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

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

以上まとめると full model は ~ group$con + group$treat、reduced model は ~ group$con で作成できることがわかった。そこで、DESeq2 を利用して尤度比検定を行う。コードは以下のように実行する。nbinomLRT 関数を実行するときに full model と reduced model の式を代入する。ただし、このとき、group 変数は一番最初に colData = group と代入してあるので、group$con と group$treat は略してそれぞれ con と treat とだけにする。

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_108 2065.6536       2.251938 0.1207667  321.7581 5.997197e-72 4.479906e-69
## gene_102  892.9625       1.977365 0.1193118  261.9719 6.378165e-59 2.382245e-56
## gene_165 1167.0530      -1.978333 0.1362185  199.0246 3.409543e-45 8.489761e-43
## gene_136  357.4408       2.101873 0.1747580  137.9950 7.305910e-32 1.364379e-29
## gene_200  702.0213      -2.005846 0.1760517  122.2135 2.072848e-28 3.096836e-26
## gene_157  293.9956      -1.871418 0.1764969  108.5903 1.995506e-25 2.484405e-23

解析結果を 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 バイオスタティスティックス