時系列比較（edgeR） |

edgeR パッケージは発現変動遺伝子を検出する目的で作られているパッケージである。一般化線形モデルをサポートし、時系列などのようなモデルにも柔軟に対応できる。

library(edgeR)
packageVersion("edgeR")
##[1] ‘3.10.0’

以下にこのページで使用するサンプルデータを紹介してから、一般化線形モデルを利用した発現変動遺伝子の検出法について紹介する。

サンプルデータ
検定（対称的な実験デザイン）
検定（非対称的な実験デザイン）

A 群と B 群の時系列データを比較するとき、A 群に 0h、3h、6h があるのに対して B 群も 0h、3h、6h が存在するとき、ここでは対称的な実験デザインとよぶ。しかし、A 群と B 群の 0h において両者は同じと考えられ、B 群 0h の実験を省略しデータが存在しないとき、ここでは非対称的な実験デザインとよぶ。

サンプルデータ

サンプルデータは条件を設定して負の二項分布に従うように乱数生成により作成した。A 群と B 群の 2 つの実験群からなる。A 群には薬剤 A を投与し、0 時間、3 時間、6 時間の 3 つの時刻においてサンプルを採取した。B 群には薬剤 B を投与し、0 時間、3 時間、6 時間の 3 つの時刻においてサンプルを採取した。各実験の各時刻において 2 つの biological replicate があり、実験全体で計 12 biological replicate を使用したとする。また、全遺伝子数は 1,000 であり、その発現パターンは以下のように設定した。

遺伝子	A 群の発現パターン	B 群の発現パターン
gene_1 ～ gene_40	A6h > A3h > A0h, B6h > B3h > B0h A0h = B0h, A3h = B3h, A0h = B0h
gene_41 ～ gene_80	A6h = A3h > A0h	A0h = B0h = B3h = B6h
gene_81 ～ gene_120	A0h = A3h = A6h	A0h = B6h = B3h > B0h
gene_121 ～ gene_160	A0h = A6h < A3h	A0h = B0h = B3h = B6h
gene_161 ～ gene_200	A0h = A6h < A3h	A0h = B0h = B3h < B6h
gene_201 ～ gene_1000	A0h = A3h = A6h = B0h = B3h = B6h

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

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

dim(count)
## [1] 1000    12

head(count)
##        A1_0h A2_0h A3_3h A4_3h A5_6h A6_6h B1_0h B2_0h B3_3h B4_3h B5_6h B6_6h
## gene_1     1     7    11     0    22    24     0     0     9     3    22    16
## gene_2   170   182   744   683   918   637   128   112   645   517  1248  1234
## gene_3     8    12    37    46   100    94     9     8    30    46    47    63
## gene_4    44    38   109   187   354   225    40    97   103   146   505   343
## gene_5     0     0     5     0     4    10     0     0     1     0     7     4
## gene_6    88   112   372   459  1010   611   129    80   486   360   678   767

実験群の識別ラベルとサンプルの採取時刻の識別ラベルをそれぞれ group と time 変数に代入する。

group <- factor(c("A", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "B"))
time  <- factor(c("0h", "0h", "3h", "3h", "6h", "6h", "0h", "0h", "3h", "3h", "6h", "6h"))

尤度比検定による発現変動遺伝子の検出

二群間時系列比較の解析において、次に上げるような比較等を行うことができる。

時刻差を考慮してないで A 群と B 群の比較
実験群を考慮してないで 0h、3h、6h の比較
A 群内において 0h、3h、6h の比較

尤度比検定を利用して上記の場合における発現変動遺伝子を検出する例を示す。

尤度比検定は 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 はデザイン行列あるいは計画行列と呼ばれ、実験群や処理群の情報などが行列で記述される。β はパラメーターベクトルとよばれ、モデル化する際に想定されるパラメーターからなるベクトルである。

ここで二群間時系列比較データをモデル化する。まず、A 群 0h の発現量を β₁ とおく。次に、A 群と B 群の発現量の差を β₂ とする。また、0h と 3h の発現量の差を β₂、0h と 6h の発現量の差を β₃ とする。このとき、E[A-0h] = β₁、E[B-0h] = β₁ + β₂、E[A-3h] = β₁ + β₃ などのように書ける。これらを行列の形で表すと以下のようになる。

\[ \begin{pmatrix} E[\mbox{A-0h}] \\ E[\mbox{A-0h}] \\ E[\mbox{A-3h}] \\ E[\mbox{A-3h}] \\ E[\mbox{A-6h}] \\ E[\mbox{A-6h}] \\ E[\mbox{B-0h}] \\ E[\mbox{B-0h}] \\ E[\mbox{B-3h}] \\ E[\mbox{B-3h}] \\ E[\mbox{B-6h}] \\ E[\mbox{B-6h}] \end{pmatrix}^{T} = \begin{pmatrix} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 1\\ 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 1 & 0\\ 1 & 1 & 1 & 0\\ 1 & 1 & 0 & 1\\ 1 & 1 & 0 & 1\\ \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3} \\ \beta_{4} \end{pmatrix} \]

この式が一般化線形モデルのモデルとされる。0 と 1 からなる行列はデザイン行列と呼び、観測値とパラメーターを関連させる重要な働きを持つ。遺伝子が n　個があればこのような式は n 個できる。遺伝子が n　個あるときプログラムが逐次的に処理を行ってくれるため、遺伝子が 1 個だけの時を考えればよい。

このモデルを full model と呼ぶことにする。full model のデザイン行列は model.matrix 関数を利用すれば簡単に作成できる。

design <- model.matrix(~ group + time)
design
##    (Intercept) groupB time3h time6h
## 1            1      0      0      0
## 2            1      0      0      0
## 3            1      0      1      0
## 4            1      0      1      0
## 5            1      0      0      1
## 6            1      0      0      1
## 7            1      1      0      0
## 8            1      1      0      0
## 9            1      1      1      0
## 10           1      1      1      0
## 11           1      1      0      1
## 12           1      1      0      1
## attr(,"assign")
## [1] 0 1 2 2
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"
## 
## attr(,"contrasts")$time
## [1] "contr.treatment"

次に、データを edgeR にセットアップし、正規化、分布のパラメーターの推定を行う。

d <- DGEList(counts = count, group = group)
d <- calcNormFactors(d)
d <- estimateGLMCommonDisp(d, design)
d <- estimateGLMTrendedDisp(d, design)
d <- estimateGLMTagwiseDisp(d, design)

推測された分布にデータを当てはめていく。

fit <- glmFit(d, design)

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

時系列データに対して様々な比較が可能であるため、すべての解析例を示すことはできない。ここではいくつかの例と考え方だけを紹介する。

実験群間の比較（A、B）

サンプルの採取時刻を考慮しないで A 群と B 群を比較し、発現量に差をもつ遺伝子を発現変動遺伝子として検出する例を示す。A 群と B 群の発現量を比較するとき、帰無仮説として E[A] = E[B] すなわち β₂ = 0 について考えればよい。ここで、full model のデザイン行列の 2 列目の要素をすべて 0 に変換したモデルを考えると、このモデルは β₂ = 0 に対応するようになる。

\[ \begin{pmatrix} E[\mbox{A-0h}] \\ E[\mbox{A-0h}] \\ E[\mbox{A-3h}] \\ E[\mbox{A-3h}] \\ E[\mbox{A-6h}] \\ E[\mbox{A-6h}] \\ E[\mbox{B-0h}] \\ E[\mbox{B-0h}] \\ E[\mbox{B-3h}] \\ E[\mbox{B-3h}] \\ E[\mbox{B-6h}] \\ E[\mbox{B-6h}] \end{pmatrix}^{T} =\begin{pmatrix} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 1\\ \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3} \\ \beta_{4} \end{pmatrix} \]

このように full model から β₂ を取り除いてできたモデルを reduced model という。尤度比検定では full model と reduce model の尤度を比較して検定を行う。もし両者の尤度に差がなければ帰無仮説は保留される。反対に、両者の尤度に差があれば帰無仮説は棄却され、β₂ ≠ 0 と考えられ、その遺伝子は発現変動遺伝子と考えられる。

次に glmLRT 関数を利用して尤度比検定を行う。fit 変数には full model の情報が含まれているため、reduced model の情報だけを与えればよい。ここで、full model と reduced model のデザイン行列を見比べてみると、reduced model は full model のデザイン行列の 2 列目をすべて 0 にしたものに等しいことがわかる。そこで、coef = 2 と指定して、reduced model の情報を与える。

lrt <- glmLRT(fit, coef = 2)

結果は topTags 関数を利用して確認する。

topTags(lrt)
## Coefficient:  groupB 
##              logFC    logCPM       LR       PValue          FDR
## gene_52  -2.420685  9.807322 39.60147 3.114481e-10 3.114481e-07
## gene_101  1.948105  9.510770 28.26999 1.055190e-07 5.275950e-05
## gene_55  -3.278910  7.522237 22.24126 2.404505e-06 8.015018e-04
## gene_103  1.691705  9.296790 20.43294 6.175750e-06 1.543937e-03
## gene_176 -5.048721  6.058174 19.05805 1.268020e-05 2.536040e-03
## gene_97   3.886038  7.392665 17.01748 3.703737e-05 5.876608e-03
## gene_72  -2.108316  8.785828 16.81823 4.113626e-05 5.876608e-03
## gene_66  -1.885312  9.230412 15.99195 6.361229e-05 7.951536e-03
## gene_108  1.500614 13.780379 15.47176 8.374730e-05 9.305255e-03
## gene_56  -1.821414  8.460870 15.08741 1.026458e-04 9.306621e-03

サンプルデータの作成条件として、A 群と B 群に差が生じるのは gene_41 ～ gene_200 である。検定結果をみると確かにこの範囲にある遺伝子が上位にランキングされている。

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

table <- as.data.frame(topTags(lrt, n = nrow(count)))
write.table(table, file = "result.txt", col.names = T, row.names = T, sep = "\t")

時刻間の比較（0h、3h、6h）

時刻間を比較したいとき帰無仮説として E[0h] = E[3h] = E[6h] を考えればよい。これは β₃ = 0 かつ β₄ = 0 としたモデルに対応している。ここで、reduced model を以下のように作成できる。この場合、もし帰無仮説が棄却されるならば β₃ ≠ 0 または β₄ ≠ 0 である。すなわち 0h と 3h で差がある、またはは 0h と 6h で差がある。

\[ \begin{pmatrix} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 1 & 0 & 0\\ \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3} \\ \beta_{4} \end{pmatrix} \]

reduced model のデザイン行列は full model のデザイン行列の 3、4 行目の要素をすべて 0 にした行列に等しいので、coef = c(3, 4) と指定して尤度比検定を行う。

lrt <- glmLRT(fit, coef = c(3,4))
topTags(lrt)
## Coefficient:  time3h time6h 
##          logFC.time3h logFC.time6h   logCPM       LR       PValue          FDR
## gene_22     1.9391885     2.854111 11.98952 74.88674 5.477104e-17 2.797306e-14
## gene_28     1.9022955     2.685560 11.51139 74.84429 5.594612e-17 2.797306e-14
## gene_24     1.7866966     2.935829 12.39924 74.02301 8.435445e-17 2.811815e-14
## gene_10     1.7278105     2.977852 11.79024 72.47242 1.831526e-16 4.578815e-14
## gene_181   -0.6381632     1.905344 10.31697 69.48329 8.163853e-16 1.632771e-13
## gene_39     1.8378258     2.812522 12.63787 66.39999 3.814392e-15 6.357321e-13
## gene_25     2.1455918     2.718319 11.33630 61.13536 5.304279e-14 7.577541e-12
## gene_188   -0.3602758     1.692665 11.20673 60.05157 9.119415e-14 1.014511e-11
## gene_6      1.9190497     2.720222 11.86655 60.04912 9.130598e-14 1.014511e-11
## gene_200   -0.3827228     1.765076 12.27319 59.52831 1.184652e-13 1.184652e-11

サンプルデータ中に、時刻の比較において発現量に差が生じる遺伝子は gene_1 ～ gene_200 の全てであるから、確かにこの範囲の遺伝子が上位にランキングされている。

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

table <- as.data.frame(topTags(lrt, n = nrow(count)))
write.table(table, file = "result.txt", col.names = T, row.names = T, sep = "\t")

A 群内における時刻間の比較

A 群内において、時刻間を比較した時に発現量に差異をもつ遺伝子を検出する方法を示す。この解析を行うためには、これまで使っていたデザイン行列では対応できない。そこで、以下のようにデザイン行列を作成する。

g <- factor(paste(group, time, sep = "_"))
design <- model.matrix(~ 0 + g)
design
##    gA_0h gA_3h gA_6h gB_0h gB_3h gB_6h
## 1      1     0     0     0     0     0
## 2      1     0     0     0     0     0
## 3      0     1     0     0     0     0
## 4      0     1     0     0     0     0
## 5      0     0     1     0     0     0
## 6      0     0     1     0     0     0
## 7      0     0     0     1     0     0
## 8      0     0     0     1     0     0
## 9      0     0     0     0     1     0
## 10     0     0     0     0     1     0
## 11     0     0     0     0     0     1
## 12     0     0     0     0     0     1
## attr(,"assign")
## [1] 1 1 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$g
## [1] "contr.treatment"

このデザイン行列について少し見ていく。このデザイン行列は 6 列であるから、パラメーター数は 6 個である。それら β₁、…、β₆ とおくと、A 群 0h の発現量の期待値は E[A-0h] = β₁ である。同様に、E[A-3h] = β₂、 E[A-6h] = β₃、E[B-0h] = β₄、E[B-3h] = β₅、E[B-6h] = β₆ である。

edgeR を利用してモデルを作成する。

d <- DGEList(counts = count, group = group)
d <- calcNormFactors(d)
d <- estimateGLMCommonDisp(d, design)
d <- estimateGLMTrendedDisp(d, design)
d <- estimateGLMTagwiseDisp(d, design)
fit <- glmFit(d, design)

まず、A 群における 0h と 3h の差は β₁ - β₂ である。β₁ および β₂ はデザイン行列の 1 列目と 2 列目に該当する。つまり、デザイン行列の 1 列目から 2 列目を引き、3 列目～6 列目は考えていないので、これをベクトルで表すと、c(1, -1, 0, 0, 0, 0) と表現できる。もちろん、c(-1, 1, 0, 0, 0, 0) も可。

このような比較で尤度比検定を行うには、coef の代わりに contrast を利用する。

lrt <- glmLRT(fit, contrast = c(1, -1, 0, 0, 0, 0))
topTags(lrt)
## Coefficient:  1*gA_0h -1*gA_3h 
##              logFC    logCPM       LR       PValue          FDR
## gene_76  -2.625642 10.553694 31.67334 1.824090e-08 0.0000182409
## gene_29  -3.937059  8.361071 24.96295 5.844279e-07 0.0002562368
## gene_68  -2.322544 10.280090 24.43471 7.687103e-07 0.0002562368
## gene_127 -2.122985 10.677586 23.61492 1.176714e-06 0.0002941786
## gene_8   -3.810994  9.133612 22.85431 1.747573e-06 0.0003495147
## gene_77  -2.791684 11.164030 22.46152 2.143954e-06 0.0003573257
## gene_136 -1.952790 10.920700 22.16095 2.507220e-06 0.0003581743
## gene_25  -2.199825 11.336295 21.88895 2.888932e-06 0.0003611164
## gene_141 -2.434866 11.295334 20.95192 4.709544e-06 0.0004869277
## gene_157 -1.879084 10.884407 20.86265 4.934233e-06 0.0004869277

0h と 6h の比較において発現量に差がある遺伝子を検出する。

lrt <- glmLRT(fit, contrast = c(1, 0, -1, 0, 0, 0))
topTags(lrt)
## Coefficient:  1*gA_0h -1*gA_6h 
##             logFC   logCPM       LR       PValue          FDR
## gene_28 -2.715344 11.51139 38.62284 5.141366e-10 3.412540e-07
## gene_18 -3.642191 11.57385 37.25205 1.038056e-09 3.412540e-07
## gene_40 -3.288071 13.37848 36.84827 1.276893e-09 3.412540e-07
## gene_24 -2.914699 12.39924 36.71816 1.365016e-09 3.412540e-07
## gene_39 -2.851247 12.63787 35.16142 3.034762e-09 5.115234e-07
## gene_10 -2.898923 11.79024 35.13949 3.069141e-09 5.115234e-07
## gene_22 -2.614589 11.98952 34.44520 4.384230e-09 6.263185e-07
## gene_6  -2.821075 11.86655 32.13168 1.440691e-08 1.800864e-06
## gene_76 -2.609176 10.55369 31.31353 2.195443e-08 2.439381e-06
## gene_11 -3.439834 11.61323 29.67760 5.102103e-08 5.102103e-06

複雑な比較

A 群と B 群、各時刻間の複雑な比較について見ていく。デザイン行列を以下のように作成する。

g <- factor(paste(group, time, sep = "_"))
design <- model.matrix(~ 0 + g)
design
##    gA_0h gA_3h gA_6h gB_0h gB_3h gB_6h
## 1      1     0     0     0     0     0
## 2      1     0     0     0     0     0
## 3      0     1     0     0     0     0
## 4      0     1     0     0     0     0
## 5      0     0     1     0     0     0
## 6      0     0     1     0     0     0
## 7      0     0     0     1     0     0
## 8      0     0     0     1     0     0
## 9      0     0     0     0     1     0
## 10     0     0     0     0     1     0
## 11     0     0     0     0     0     1
## 12     0     0     0     0     0     1
## attr(,"assign")
## [1] 1 1 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$g
## [1] "contr.treatment"

edgeR を利用してモデルを作成する。

d <- DGEList(counts = count, group = group)
d <- calcNormFactors(d)
d <- estimateGLMCommonDisp(d, design)
d <- estimateGLMTrendedDisp(d, design)
d <- estimateGLMTagwiseDisp(d, design)
fit <- glmFit(d, design)

A 群 0h と B 群 0h を比較した際に、発現量に差を持つ遺伝子を検出する。このサンプルではこの発現パターンの遺伝子が含まれていない。検出結果を見ると、FDR が非常に高いことがわかる。

lrt <- glmLRT(fit, contrast = c(1, 0, 0, -1, 0, 0))
topTags(lrt)
## Coefficient:  1*gA_0h -1*gB_0h 
##              logFC   logCPM        LR       PValue        FDR
## gene_801 -2.547985 9.356602 15.785192 7.095568e-05 0.07095568
## gene_266 -4.566140 8.826286 14.151082 1.687000e-04 0.08435002
## gene_523  7.082203 6.065970 12.119940 4.988551e-04 0.16628505
## gene_277 -3.299410 8.078138 11.218605 8.098133e-04 0.16838863
## gene_52   2.322367 9.807322 11.146413 8.419431e-04 0.16838863
## gene_600 -2.068104 9.264968  9.909299 1.644458e-03 0.27407628
## gene_64  -4.152995 8.926953  8.198662 4.192128e-03 0.45869730
## gene_472  5.328702 6.123880  7.939613 4.836418e-03 0.45869730
## gene_992 -6.881579 6.014516  7.817265 5.174944e-03 0.45869730
## gene_721 -5.683493 4.996826  7.768701 5.315928e-03 0.45869730

A 群 3h と B 群 3h を比較した際に、発現量に差を持つ遺伝子を検出する。

lrt <- glmLRT(fit, contrast = c(0, 1, 0, 0, -1, 0))
topTags(lrt)
## Coefficient:  1*gA_3h -1*gB_3h 
##              logFC    logCPM       LR       PValue          FDR
## gene_80   2.661480 10.908759 32.28732 1.329782e-08 1.329782e-05
## gene_86  -3.359239 11.000177 30.01451 4.288265e-08 2.144132e-05
## gene_76   2.391902 10.553694 27.09941 1.932565e-07 6.243973e-05
## gene_114 -2.245603 12.155892 26.60375 2.497589e-07 6.243973e-05
## gene_74   2.193330 12.023350 25.04785 5.592520e-07 1.059103e-04
## gene_92  -2.036336 11.273326 24.80156 6.354618e-07 1.059103e-04
## gene_157  1.998217 10.884407 23.43107 1.294707e-06 1.849581e-04
## gene_896  8.493779  7.057134 23.00043 1.619653e-06 2.024566e-04
## gene_66   2.822241  9.230412 21.94140 2.811036e-06 3.123373e-04
## gene_102 -2.019170 12.616172 20.15228 7.151486e-06 6.772839e-04

非対称的な実験デザイン

上で利用したサンプルデータは A 群に 3 つの処理群、B 群に 3 つの処理群が存在する対称的な実験デザインとなっている。しかし、A 群 0h と B 群 0h は同じと考えられる場合、B 群 0h の実験が省略されて行わない場合が多い。ここでは、B 群 0h がない場合を想定して解析例を示す。

サンプルデータから B 群 0h のデータを取り除く。

count2 <- count[, -c(7, 8)]
group2 <- group[-c(7, 8)]
time2 <- time[-c(7, 8)]

実験群間の比較（A、B）

まずモデルのデザイン行列を考える。デザイン行列は上と同じような考え方で作成する。ただし、B 群 0h のデータがないことに注意する。

\[ E\left[\begin{pmatrix} \mbox{A-0h} \\ \mbox{A-0h} \\ \mbox{A-3h} \\ \mbox{A-3h} \\ \mbox{A-6h} \\ \mbox{A-6h} \\ \mbox{B-3h} \\ \mbox{B-3h} \\ \mbox{B-6h} \\ \mbox{B-6h} \end{pmatrix}^{T}\ \right] = \begin{pmatrix} 1 & 0 & 0 & 0\\ 1 & 0 & 0 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 1 & 0\\ 1 & 0 & 0 & 1\\ 1 & 0 & 0 & 1\\ 1 & 1 & 1 & 0\\ 1 & 1 & 1 & 0\\ 1 & 1 & 0 & 1\\ 1 & 1 & 0 & 1\\ \end{pmatrix} \begin{pmatrix}\beta_{1} \\ \beta_{2} \\ \beta_{3} \\ \beta_{4} \end{pmatrix} \]

R でデザイン行列を作成する。

design <- model.matrix(~ group2 + time2)
design
##    (Intercept) group2B time23h time26h
## 1            1       0       0       0
## 2            1       0       0       0
## 3            1       0       1       0
## 4            1       0       1       0
## 5            1       0       0       1
## 6            1       0       0       1
## 7            1       1       1       0
## 8            1       1       1       0
## 9            1       1       0       1
## 10           1       1       0       1
## attr(,"assign")
## [1] 0 1 2 2
## attr(,"contrasts")
## attr(,"contrasts")$group2
## [1] "contr.treatment"
## 
## attr(,"contrasts")$time2
## [1] "contr.treatment"

A 群と B 群を比較したければ、β₂ = 0 を考えればよく、これはデザイン行列の 2 列目の要素をすべて 0 にしたときに対応している。よって、edgeR で解析を行うと、以下のようになる。

d <- DGEList(counts = count2, group = group2)
d <- calcNormFactors(d)
d <- estimateGLMCommonDisp(d, design)
d <- estimateGLMTrendedDisp(d, design)
d <- estimateGLMTagwiseDisp(d, design)
fit <- glmFit(d, design)
lrt <- glmLRT(fit, coef = 2)
topTags(lrt)
## Coefficient:  group2B 
##              logFC    logCPM       LR       PValue          FDR
## gene_86   3.051527 11.143101 51.51443 7.107279e-13 7.107279e-10
## gene_76  -2.349481 10.568337 49.93879 1.586180e-12 7.930902e-10
## gene_92   2.036837 11.368311 47.26248 6.208941e-12 2.069647e-09
## gene_66  -2.715280  9.325026 42.94349 5.634415e-11 1.408604e-08
## gene_72  -2.922505  8.909729 40.02975 2.501244e-10 4.697773e-08
## gene_114  2.041853 12.248208 39.79638 2.818664e-10 4.697773e-08
## gene_119  1.926990 10.632820 36.76081 1.335482e-09 1.907832e-07
## gene_110  1.907673 10.532696 36.39511 1.611068e-09 2.013836e-07
## gene_70  -1.923763 11.985627 36.03545 1.937598e-09 2.152886e-07
## gene_74  -1.886521 12.084908 34.75668 3.735953e-09 3.735953e-07

サンプルデータの作成条件として、A 群と B 群に差が生じるのは gene_41 ～ gene_200 である。

時刻間の比較（0h、3h、6h）

時刻間の比較では β₃ = β₄ = 0 を考えればよい。

d <- DGEList(counts = count2, group = group2)
d <- calcNormFactors(d)
d <- estimateGLMCommonDisp(d, design)
d <- estimateGLMTrendedDisp(d, design)
d <- estimateGLMTagwiseDisp(d, design)
fit <- glmFit(d, design)
lrt <- glmLRT(fit, coef = c(3, 4))
topTags(lrt)
## Coefficient:  time23h time26h 
##          logFC.time23h logFC.time26h   logCPM       LR       PValue          FDR
## gene_181    -0.5715937      1.996192 10.39799 62.19151 3.128133e-14 3.128133e-11
## gene_188    -0.3260292      1.743912 11.27695 53.61871 2.274281e-12 1.137141e-09
## gene_200    -0.3332211      1.832648 12.34453 50.91423 8.792544e-12 2.930848e-09
## gene_183    -0.2943992      2.121198  9.18451 45.75722 1.158631e-10 2.896577e-08
## gene_187     0.3096028      2.065764 11.15479 42.86327 4.924466e-10 9.848932e-08
## gene_10      1.5238915      2.807287 11.96056 39.83214 2.241618e-09 3.268758e-07
## gene_28      1.8703839      2.673122 11.66848 39.79106 2.288130e-09 3.268758e-07
## gene_24      1.7247186      2.894191 12.55939 39.23731 3.018046e-09 3.772558e-07
## gene_22      1.6994857      2.633866 12.16310 37.12518 8.677006e-09 9.641118e-07
## gene_4       1.5056331      2.818353 10.75742 36.55078 1.156374e-08 1.156374e-06

References

Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010, 26(1):139-40. PubMed Abstract Bioconductor
McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012, 40(10):4288-97. PubMed Abstract
Robinson MD, Smyth GK. Moderated statistical tests for assessing differences in tag abundance. Bioinformatics. 2007, 23(21):2881-7. PubMed Abstract
Robinson MD, Smyth GK. Small-sample estimation of negative binomial dispersion, with applications to SAGE data. Biostatistics. 2008, 9(2):321-32. PubMed Abstract
Zhou X, Lindsay H, Robinson MD. Robustly detecting differential expression in RNA sequencing data using observation weights. Nucleic Acids Research. 2014, 42(11):e91. PubMed Abstract