edgeR パッケージは RNA-seq のリードカウントデータから発現変動遺伝子を検出するときによく利用されている有名なパッケージである。性能の評価においてはよく DESeq2 などと比較される。edgeR では一般化線形モデルを利用した検定法 likelihood ratio test (LRT) と quasi-likelihood F-tests (QLF) が用意されている。いずれの方法も二群間比較のみならず、複雑なデザインにも対応できるようになっている。
edgeR のバージョンが異なると、解析結果も少し異なってくる。このページでは 3.10.0 を利用してる。
library(edgeR)
packageVersion("edgeR")
## [1] ‘3.30.3’サンプルデータ
サンプルデータは、カウントデータが負の二項分布に従うように乱数を生成して作成した。A 群と B 群の 2 つの実験群からなり、それぞれの実験群には 3 つの biological replicate が含まれる。また、全遺伝子数は 1,000 個であり、その発現パターンは以下のように設定した。
| gene_1 ～ gene_160 | B 群に比べ、A 群での発現量が高い |
| gene_161 ～ gene_200 | A 群に比べ、B 群で発現量が高い |
| gene_201 ～ gene_1000 | 両群において差異なし |
サンプルデータを読み込んでから行列型に変換し count 変数に代入する。
count <- read.table("https://bi.biopapyrus.jp/data/counts_3_3.txt",
sep = "\t", header = T, row.names = 1)
count <- as.matrix(count)
dim(count)
## [1] 1000 6
head(count)
## A1 A2 A3 B1 B2 B3
## gene_1 3 2 40 4 2 0
## gene_2 680 807 594 165 91 73
## gene_3 35 38 31 7 9 11
## gene_4 103 165 110 38 31 28
## gene_5 2 9 1 3 1 0
## gene_6 419 373 399 133 89 58実験群の識別ラベルを作成し、group 変数に代入する。
group <- factor(c("A", "A", "A", "B", "B", "B"))LRT （尤度比検定）
尤度比検定は一般化線形モデルの検定でよく使われる検定の一つである。検定を行うには 2 つのモデルを用意する必要がある。それぞれモデルから尤度を計算し、比較する。尤度の比較を行うことにより、2 つのモデルが同じかどうかを検定する。例えば、モデル F は β1 と β2 の 2 つのパラメーターから構成されたものとする。また、モデル R は β1 だけから構成されたものとする。ここで、モデル F とモデル R の尤度を比較し、もし両者の尤度に有意差が見られるならば、2 つのモデルが異なるモデルであることを意味する。言い換えれば、β2 の効果により、モデル F とモデル R に差異をもたらした。すなわち β2 ≠ 0 と考えられる。逆にもし両者の尤度に有意差が見られなければ、β2 の効果が認められないこととなり、β2 = 0 と考えられる。
いま、遺伝子 i に着目して尤度比検定を説明する。遺伝子 i の A 群でのカウントデータを YiA とし、その期待値は E[YiA] と表せる。同様に、B 群でのカウントデータの YiB とし、その期待値は E[YiB] と表せる。ここでパラメーター βi1 を仮定し、E[YiA] = βi1 とする。次に、もう 1 つのパラメーター βi2 を仮定し、E[YiB] = βi1 + βi2 とする。これを行列式で記述すると以下のようになる。
この仮定で検定について考える。検定を考えるとはまず帰無仮説を立てなければならない。ここでは A 群の発現量と B 群の発現量に差異があるかどうかを検定するのが目的である。このとき帰無仮説を「A 群の発現量と B 群の発現量に差異がない」とすることができる。すると、これは「βi2 = 0」に等しい。そこで、βi2 = 0 のときの各群の期待値は以下のように記述できる。
ここまで行列で記述した E[Y] = Xβ の形をした数式をモデルという。1 つ目のモデルは想定されたすべてのパラメーターが使われていることから full model と呼ばれる。2 つ目のモデルは full model からパラメーターが 1 つ取り除かれて（βi2 = 0）いることから reduced model と呼ばれる。また、行列 X をデザイン行列または計画行列という。
次に full model と reduced model のそれぞれから尤度を計算する。両者の尤度を比較する。もし両者の尤度に有意差が認められれば、これは full model と reduced model が相違異なるモデルであることを意味する。このとき帰無仮説を棄却することができ、すなわち、βi2 ≠ 0 あるいは、E[YiA] ≠ E[YB] と判定を下すことができる。逆に、両者の尤度に有意差が認められなければ、full model と reduced model が同じと考えられて帰無仮説は保留される、すなわち βi2 = 0 すなわち E[YiA] ≠ E[YB] の判定を下せなくなる。このような作業を遺伝子 1 から遺伝子 n まで繰り返していく。
（仮説検定においてここで注意したいたのは、帰無仮説が保留される（棄却されない）とき、E[YiA] = E[YB] を意味しているわけではない、あくまで判定不可能であり、判断を保留することである。）
このように方法で検定を行う方法を尤度比検定という。上述では E[Y] = Xβ と記述したが、正しくは g(E[Y]) = Xβ と記述すべきである。関数 g はリンク関数と呼ばれている。正規分布では g(x) = x であり、リードカウントデータが従う負の二項分布では g(x) = log(x) が使われている。データの従う分布が決まれば、リンク関数の形もほぼ決められているため、R で解析を行う際にリンク関数のことを考慮せずに、右辺の Xβ だけに着目すればよい。このページおいては、説明の簡潔のためにリンク関数を省略して書く。
では、実際にサンプルデータをモデル化して edgeR で解析していく。ある遺伝子 1 つだけに着目したとき、A 群の発現量を β1 とおく。次に、A 群の発現量と B 群の発現量の差を β2 とおく。このとき、E[A] = β1、E[B] = β1 + β2 である。サンプルデータは 3 × 3 の biological replicate からなるから、モデルは以下のようになる。
遺伝子が n 個があればこのような式は n 個できる。遺伝子が n 個あるときプログラムが逐次的に処理を行ってくれるため、遺伝子が 1 個だけの時を考えればよい。
A 群と B 群の発現量を比較したいので、帰無仮説として A 群と B 群の発現量は同じであると考えればよい。そこで帰無仮説を検定するために、1 つ目のモデルと比較するためのモデルを以下のように作成できる。
ここまで作成した 2 つのモデルはそれぞれ full model と reduced model という。前者は想定されるパラメーター β1 と β2 をすべて含んでいることから full model と呼ばれている。後者は full model からパラメーター β2 を取り除いてできたことから reduced model と呼ばれる。
full model のデザイン行列 model.matrix 関数を利用して簡単に作成できる。
design <- model.matrix(~ group)
design
## (Intercept) groupB
## 1 1 0
## 2 1 0
## 3 1 0
## 4 1 1
## 5 1 1
## 6 1 1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"edgeR を利用して解析を行う。分散の推定を行うときデザイン行列を与える必要がある。
d <- DGEList(counts = count, group = group)
d <- calcNormFactors(d)
d <- estimateDisp(d, design)推測された分布に観測値を当てはめていく。
fit <- glmFit(d, design)最後に尤度比検定を行う。このとき、full model と reduced model を比べるので、この情報を関数に与える。glmFit を利用して作成した 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 関数で検定結果のうち p-value がもっとも小さい 10 遺伝子を出力する。なお、edgeR のバージョンが異なると、検定結果が大きく異なることに注意。
topTags(lrt)
## Coefficient: groupB
## logFC logCPM LR PValue FDR
## gene_200 2.195858 12.644913 57.76680 2.951089e-14 2.951089e-11
## gene_188 2.445163 11.517423 55.50753 9.310189e-14 4.655094e-11
## gene_136 -2.207490 11.539922 54.61523 1.465950e-13 4.886498e-11
## gene_165 2.044239 13.345054 51.90861 5.814504e-13 1.453626e-10
## gene_114 -2.033033 12.428906 49.67530 1.814155e-12 3.628309e-10
## gene_189 2.394394 10.911017 47.18581 6.456610e-12 9.982399e-10
## gene_192 2.307743 10.621263 47.03089 6.987679e-12 9.982399e-10
## gene_139 -2.757761 9.245284 43.37945 4.508996e-11 5.636245e-09
## gene_157 -1.954365 11.570958 42.68982 6.414512e-11 7.127235e-09
## gene_135 -2.533053 9.875935 41.79947 1.011310e-10 1.011310e-08QLF （疑似尤度 F 検定）
擬似尤度（quasi-likelihood）を利用した F 検定を使用して発現変動遺伝子を検出する場合も、考え方は基本的に LRT と同じである。また、解析コードも、一部の関数名が異なること以外は、LRT とほぼ同じである。
design <- model.matrix(~ group)
d <- DGEList(counts = count, group = group)
d <- calcNormFactors(d)
d <- estimateDisp(d, design)
fit <- glmQLFit(d, design)
qlf <- glmQLFTest(fit, coef = 2)
topTags(qlf)
## Coefficient: groupB
## logFC logCPM F PValue FDR
## gene_136 -2.207476 11.539922 155.9306 4.176689e-06 0.001660442
## gene_200 2.195897 12.644913 154.0060 4.358258e-06 0.001660442
## gene_114 -2.033018 12.428906 134.4881 6.923989e-06 0.001660442
## gene_165 2.044244 13.345054 129.9113 7.789929e-06 0.001660442
## gene_188 2.445157 11.517423 126.9919 8.415577e-06 0.001660442
## gene_157 -1.954393 11.570958 115.4519 1.162290e-05 0.001660442
## gene_192 2.307515 10.621263 110.5874 1.344179e-05 0.001660442
## gene_139 -2.757970 9.245284 109.1976 1.402749e-05 0.001660442
## gene_145 -1.813436 11.596414 100.6024 1.848150e-05 0.001660442
## gene_189 2.394487 10.911017 100.4532 1.857376e-05 0.001660442LRT と QLF で検定した結果にたいして FDR = 0.1 を閾値として発現変動遺伝子を取得し、両者のオーバーラップを計算してみると、次のようになる。少数の遺伝子において、ことなる結果となっていることがわかる。
deg_LRT <- as.data.frame(topTags(lrt, n = nrow(count)))
deg_QLF <- as.data.frame(topTags(qlf, n = nrow(count)))
setdiff(rownames(deg_LRT)[deg_LRT$FDR < 0.1], rownames(deg_QLF)[deg_QLF$FDR < 0.1])
## [1] "gene_15" "gene_193" "gene_198" "gene_19" "gene_739" "gene_78"
## [7] "gene_141" "gene_64" "gene_26" "gene_170" "gene_132" "gene_97"
## [13] "gene_51" "gene_992" "gene_121" "gene_260" "gene_128" "gene_805"
## [19] "gene_7" "gene_137" "gene_613" "gene_153" "gene_101" "gene_115"
## [25] "gene_433" "gene_65" "gene_877" "gene_169" "gene_9" "gene_100"
## [31] "gene_71" "gene_12" "gene_124" "gene_125" "gene_57"
setdiff(rownames(deg_QLF)[deg_QLF$FDR < 0.1], rownames(deg_LRT)[deg_LRT$FDR < 0.1])
## [1] "gene_644"References
- edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010, 26(1):139-40. PubMed Abstract Bioconductor
- Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012, 40(10):4288-97. PubMed Abstract
- Moderated statistical tests for assessing differences in tag abundance. Bioinformatics. 2007, 23(21):2881-7. PubMed Abstract
- Small-sample estimation of negative binomial dispersion, with applications to SAGE data. Biostatistics. 2008, 9(2):321-32. PubMed Abstract
- Robustly detecting differential expression in RNA sequencing data using observation weights. Nucleic Acids Research. 2014, 42(11):e91. PubMed Abstract
- 一般化線形モデル. biopapyrus バイオスタティスティックス