二群間二因子比較(DESeq)

DESeq は RNA-seq のリードカウントデータから発現変動遺伝子を検出するためのパッケージである。一般化線形モデルをサポートしているため複雑な解析にも柔軟に対応できる。ここでは尤度比検定を利用して二群間二因子のデータの解説について例を示す。

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

library(DESeq)
packageVersion("DESeq")
##[1] ‘1.20.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("http://bi.biopapyrus.net/data/count2gf.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 変数に保存する。ただし、DESeq では、多因子解析においてデータフレーム型を要求するため、これらをデータフレームにまとめる。

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"))
)

尤度比検定

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

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

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

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

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

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

\[ \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 関数を利用すれば簡単に作成できる。(group$conWT の列は上の行列と対称になっているが、解析にあたって効果は同じである。)

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 群の発現量の差はパラメーター β2 によって支配されている。帰無仮説が成り立つとき、β2 = 0 も成り立つ。そこで、帰無仮説を言い換えれば「β2 = 0」である。例えば、β2 = 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 からパラメーター β2 を取り除いてできたモデルであり、full model に対して reduced model と呼ぶことにする。full model と reduced model のそれぞれから尤度を計算し比較する。もし両者の尤度に差異が認められた場合、full model と reduced モデルは相違異なるモデルと判定できる。full model と reduced model の差は β2 だけであるから、β2 ≠ 0 であることを意味する。しかし、これは最初に仮定した β2 = 0 (帰無仮説)に対し矛盾である。したがって、最初の仮説(帰無仮説)が間違いであり、棄却するべきである。このとき、E[WT-A] ≠ E[KO-A] または E[WT-B] ≠ E[KO-B] が成り立つ。逆にもし両者の尤度に差異が認められない場合は β2 = 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"

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

cds <- newCountDataSet(countData = count, conditions = group)
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
f.model <- fitNbinomGLMs(cds, count ~ con + treat)
r.model <- fitNbinomGLMs(cds, count ~ treat)
res <- nbinomGLMTest(f.model, r.model)
names(res) <- rownames(count)
head(res[order(res)], n = 20)
##      gene_89      gene_34      gene_18      gene_62      gene_74      gene_50 
## 1.645526e-10 1.674688e-10 3.622016e-10 6.050001e-10 1.567991e-09 2.098149e-09 
##      gene_16      gene_76      gene_92      gene_86      gene_68      gene_82 
## 3.955669e-09 5.475852e-09 2.152089e-08 4.426078e-08 4.498999e-08 4.918935e-08 
##      gene_39      gene_43      gene_24      gene_10      gene_22       gene_2 
## 7.513697e-08 1.398676e-07 2.043491e-07 8.287327e-07 1.128385e-06 1.396328e-06 
##      gene_14      gene_61 
## 2.114047e-06 4.163467e-06

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

解析結果を整形して、result.txt ファイルに保存する。

table <- data.frame(
  gene.name = rownames(count),
  p.value = res,
  q.value = p.adjust(res, method = "BH")
)
write.table(table, file = "result.txt", col.names = T, row.names = F, 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 で作成できることがわかった。そこで、DESeq を利用して尤度比検定を行う。コードは以下のように実行する。nbinomGLMTest 関数を実行するときに full model と reduced model の式を代入する。ただし、このとき、group 変数は一番最初に conditions = group と代入してあるので、group$congroup$treat は略してそれぞれ contreat とだけにする。

cds <- newCountDataSet(countData = count, conditions = group)
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
f.model <- fitNbinomGLMs(cds, count ~ con + treat)
r.model <- fitNbinomGLMs(cds, count ~ con)
res <- nbinomGLMTest(f.model, r.model)
names(res) <- rownames(count)
head(res[order(res)], n = 20)
##     gene_125     gene_195     gene_108     gene_189     gene_188     gene_187 
## 3.869738e-11 1.329800e-10 2.576165e-10 9.458434e-10 9.627557e-09 1.416208e-08 
##     gene_165     gene_136     gene_200     gene_169     gene_102     gene_181 
## 2.985021e-08 3.233955e-08 3.532067e-08 3.796855e-08 3.934504e-08 4.344709e-08 
##     gene_111     gene_172     gene_131     gene_145     gene_114     gene_147 
## 5.311448e-08 5.881221e-08 7.133965e-08 8.193690e-08 1.214320e-07 2.367792e-07 
##     gene_122     gene_163 
## 7.085654e-07 8.554803e-07 

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

解析結果を整形して、result.txt ファイルに保存する。

table <- data.frame(
  gene.name = rownames(count),
  p.value = res,
  q.value = p.adjust(res, method = "BH")
)
write.table(table, file = "result.txt", col.names = T, row.names = F, sep = "\t")

References