二群間比較(edgeR)

edgeR パッケージは RNA-seq のリードカウントデータから発現変動遺伝子を検出するときによく利用されている有名なパッケージである。性能の評価においてはよく DESeq / DESeq2 などと比較される。二群間比較のデータから発現変動遺伝子を検出方法として、edgeR には 2 種類の方法が実装されている。一つは、「観測値よりも大きい値が出現する確率」を p-value とする exact test と呼ばれている方法である。もう一つは一般化線形モデルの尤度比検定を利用する方法である。後者は二群間比較のみならず複雑なモデルにも適用できる方法となっている。いずれの方法においても以下の手順で行う。

  1. カウントデータと実験群の情報をセットアップ
  2. データの正規化(正規化係数を計算)
  3. カウントデータの分布パラメーターを推測
  4. 検定

edgeR のバージョンが異なると、解析結果も少し異なってくる。このページでは 3.10.0 を利用してる。

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

以下にこのページで使用するサンプルデータについて紹介してから、exact test による発現変動遺伝子同定と一般化線形モデルの尤度比検定による発現変動遺伝子同定について解析例を紹介する。

  1. サンプルデータ
  2. exact test による二群間比較
  3. 尤度比検定による二群間比較
  4. 二群間比較(biological replicate なしの場合)

サンプルデータ

サンプルデータは、カウントデータが負の二項分布に従うように乱数を生成して作成した。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/https://stats.biopapyrus.jp/la/33matrix.html", 
                    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"))

exact test

exact test は分布の形が定められ(推定され)てから、その分布において「観測値よりも大きい値が出現する確率」を p-value とした検出方法である。例えば、ある遺伝子 i について、その遺伝子のカウントデータから平均と分散を計算し、分布を推定する。その分布の確率質量関数を Pr とおく。次に、観測値の合計値 N を、推定された分布(分散は合計値に合わせるように修正される)に当てはめ、その観測値の合計値よりも大きい値が観測される確率 Pr(X > N) を p-value としている。以下に exact test を利用した方法を示す。

はじめにカウントデータと実験群の識別ラベルを一つの変数に代入する。

d <- DGEList(counts = count, group = group)

次に、正規化係数の計算(calcNormFactors)および分散(estimateCommonDisp および estimateTagwiseDsip)の推測を行う。edgeR が採用している負の二項分布については平均を μ、分散を μ(1 + μφ) と定義している。common dispersion はすべてのデータを利用して推測された φ のことを指す、逆に言えばすべての遺伝子は 1 つだけの φ を共用している。しかし、現実にはすべての遺伝子が同じ φ を持つと考えられない。そこで各遺伝子のデータと commonDisp の関係を考慮し、各遺伝子 1 つずつについて φ を計算する。このように遺伝子ごとに計算された φ のことを tagwiseDsip という。

d <- calcNormFactors(d)
d <- estimateCommonDisp(d)
d <- estimateTagwiseDisp(d)

最後に exact test による検定を行う。検定結果を result 変数に代入する。

result <- exactTest(d)

検定結果をテーブル形式で表示させるには topTags 関数を利用する。

topTags(result)
## Comparison of groups:  B-A 
##              logFC   logCPM       PValue          FDR
## gene_188  2.445160 11.51743 4.729100e-10 4.729100e-07
## gene_189  2.394442 10.91102 3.635673e-09 1.817836e-06
## gene_136 -2.207481 11.53992 5.738655e-09 1.912885e-06
## gene_2   -2.536679 11.76044 9.836109e-09 2.459027e-06
## gene_196  2.536024 11.58500 1.329859e-08 2.659717e-06
## gene_200  2.195887 12.64491 1.690998e-08 2.818330e-06
## gene_192  2.307607 10.62127 2.283519e-08 2.911816e-06
## gene_165  2.044242 13.34506 2.351758e-08 2.911816e-06
## gene_77  -2.415902 10.90162 2.620634e-08 2.911816e-06
## gene_127 -2.280278 11.15094 3.400317e-08 3.400317e-06

サンプルデータの生成条件として、発現変動遺伝子を gene_1 ~ gene_200 としている。解析結果を見ると上位にランキングされている遺伝子はすべてこの範囲に入っていることが確認できる。

解析結果を result.txt ファイルに保存する例。topTags 関数はデフォルトでは結果の最初の 10 件だけを表示する。すべての結果を表示させるためには引数 n に入力データの遺伝子数を与えればよい。

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

解析結果を図示する場合は以下のようにする。FDR < 0.01 の遺伝子を発現変動遺伝子として赤色でプロットする例。

is.DEG <- as.logical(table$FDR < 0.01)
DEG.names <- rownames(table)[is.DEG]
plotSmear(result, de.tags = DEG.names)

尤度比検定

尤度比検定は一般化線形モデルの検定でよく使われる検定の一つである。検定を行うには 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 とする。これを行列式で記述すると以下のようになる。

\[ \begin{pmatrix}E[Y_{iA}] & E[Y_{iB}]\end{pmatrix} = \begin{pmatrix}1 & 0 \\ 1 & 1 \end{pmatrix} \begin{pmatrix}\beta_{i1} \\ \beta_{i2}\end{pmatrix} = \mathbf{X}_{F}\mathbf{\beta_{i}} \]

この仮定で検定について考える。検定を考えるとはまず帰無仮説を立てなければならない。ここでは A 群の発現量と B 群の発現量に差異があるかどうかを検定するのが目的である。このとき帰無仮説を「A 群の発現量と B 群の発現量に差異がない」とすることができる。すると、これは「βi2 = 0」に等しい。そこで、βi2 = 0 のときの各群の期待値は以下のように記述できる。

\[ \begin{pmatrix}E[Y_{iA}] & E[Y_{iB}]\end{pmatrix} = \begin{pmatrix}1 & 0 \\ 1 & 0 \end{pmatrix} \begin{pmatrix}\beta_{i1} \\ \beta_{i2}\end{pmatrix} = \mathbf{X}_{R}\mathbf{\beta_{i}} \]

ここまで行列で記述した E[Y] = の形をした数式をモデルという。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] = と記述したが、正しくは g(E[Y]) = と記述すべきである。関数 g はリンク関数と呼ばれている。正規分布では g(x) = x であり、リードカウントデータが従う負の二項分布では g(x) = log(x) が使われている。データの従う分布が決まれば、リンク関数の形もほぼ決められているため、R で解析を行う際にリンク関数のことを考慮せずに、右辺の だけに着目すればよい。このページおいては、説明の簡潔のためにリンク関数を省略して書く。

では、実際にサンプルデータをモデル化して edgeR で解析していく。ある遺伝子 1 つだけに着目したとき、A 群の発現量を β1 とおく。次に、A 群の発現量と B 群の発現量の差を β2 とおく。このとき、E[A] = β1、E[B] = β1 + β2 である。サンプルデータは 3 × 3 の biological replicate からなるから、モデルは以下のようになる。

\[ \begin{pmatrix} E[A] & E[A] & E[A] & E[B] & E[B] & E[B]\end{pmatrix} = \begin{pmatrix}1 & 0 \\ 1 & 0 \\ 1 & 0 \\ 1 & 1 \\ 1 & 1 \\ 1 & 1\end{pmatrix}\begin{pmatrix}\beta_{1} \\ \beta_{2}\end{pmatrix} \]

遺伝子が n 個があればこのような式は n 個できる。遺伝子が n 個あるときプログラムが逐次的に処理を行ってくれるため、遺伝子が 1 個だけの時を考えればよい。

A 群と B 群の発現量を比較したいので、帰無仮説として A 群と B 群の発現量は同じであると考えればよい。そこで帰無仮説を検定するために、1 つ目のモデルと比較するためのモデルを以下のように作成できる。

\[ \begin{pmatrix} E[A] & E[A] & E[A] & E[B] & E[B] & E[B]\end{pmatrix} = \begin{pmatrix}1 & 0 \\ 1 & 0 \\ 1 & 0 \\ 1 & 0 \\ 1 & 0 \\ 1 & 0\end{pmatrix}\begin{pmatrix}\beta_{1} \\ \beta_{2}\end{pmatrix} \]

ここまで作成した 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 <- estimateGLMCommonDisp(d, design)
d <- estimateGLMTrendedDisp(d, design)
d <- estimateGLMTagwiseDisp(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 関数で検定結果をテーブルで表示させる。

topTags(lrt)
## Coefficient:  groupB 
##              logFC   logCPM       LR       PValue          FDR
## gene_188  2.445162 11.51540 50.14539 1.427662e-12 1.427662e-09
## gene_136 -2.207486 11.54174 44.33053 2.773571e-11 1.386786e-08
## gene_189  2.394413 10.90818 42.30694 7.801510e-11 2.600503e-08
## gene_2   -2.536700 11.76207 40.23513 2.251630e-10 5.468584e-08
## gene_196  2.535953 11.58307 39.85574 2.734292e-10 5.468584e-08
## gene_127 -2.280278 11.15333 37.37626 9.740043e-10 1.623340e-07
## gene_77  -2.415938 10.90485 36.79154 1.314595e-09 1.877993e-07
## gene_192  2.307650 10.61789 36.36782 1.633785e-09 2.028805e-07
## gene_173  2.575603 11.20475 36.15113 1.825924e-09 2.028805e-07
## gene_86  -2.055414 11.09145 35.61660 2.402321e-09 2.402321e-07

二群間比較(biological replicate なし)

biological replicate がない場合、ばらつきを予測する場所だけが異なっている。それ以外は、上述と同じように実行できる。

count <- count[, c(1, 4)]
group <- factor(c("A", "B"))
head(count)
##         A1  B1
## gene_1   3   4
## gene_2 680 165
## gene_3  35   7
## gene_4 103  38
## gene_5   2   3

exact test を利用する場合は以下のようにする。

d <- DGEList(counts = count, group = group)
d <- calcNormFactors(d)
d <- estimateGLMCommonDisp(d, method = "deviance", robust = T, subset = NULL)
result <- exactTest(d)
topTags(result)
## Comparison of groups:  B-A 
##               logFC   logCPM       PValue       FDR
## gene_266 -10.002903 9.207747 0.0009889481 0.5826503
## gene_134  -7.460980 9.882436 0.0011653006 0.5826503
## gene_81   -9.461017 8.689184 0.0019705334 0.6568445
## gene_109  -6.292145 8.762635 0.0050975778 1.0000000
## gene_174   6.001967 8.337332 0.0085507710 1.0000000
## gene_179   5.485791 8.730376 0.0098292723 1.0000000
## gene_531  -8.110561 7.440717 0.0109289430 1.0000000
## gene_712   8.072989 7.346240 0.0115181135 1.0000000
## gene_862   7.980237 7.262701 0.0128526323 1.0000000
## gene_277   7.931524 7.219072 0.0136107081 1.0000000

尤度比検定の場合は以下のようにする。

design <- model.matrix(~ group)

d <- DGEList(counts = count, group = group)
d <- calcNormFactors(d)
d <- estimateGLMCommonDisp(d, method = "deviance", robust = T, subset = NULL)
fit <- glmFit(d, design)
lrt <- glmLRT(fit, coef = 2)
topTags(lrt)
## Coefficient:  groupB 
##               logFC   logCPM        LR       PValue       FDR
## gene_266 -10.002903 9.207747 14.335467 0.0001529558 0.1419003
## gene_81   -9.461017 8.689184 12.960732 0.0003180929 0.1419003
## gene_134  -7.460980 9.882436 12.415869 0.0004257008 0.1419003
## gene_712   8.072989 7.346240  9.833202 0.0017138978 0.2597371
## gene_531  -8.110561 7.440717  9.621724 0.0019228912 0.2597371
## gene_862   7.980237 7.262701  9.608605 0.0019366766 0.2597371
## gene_277   7.931524 7.219072  9.491070 0.0020647438 0.2597371
## gene_109  -6.292145 8.762635  9.479420 0.0020778965 0.2597371
## gene_193   7.774664 7.079789  9.114679 0.0025356648 0.2817405
## gene_174   6.001967 8.337332  8.794081 0.0030220943 0.2990620

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
  • 一般化線形モデル. biopapyrus バイオスタティスティックス