ex403.v2 - nibb-gitc/gitc2025mar-rnaseq GitHub Wiki

ex403.v2-1: Differential expression analysis with edgeR (GLM)

arab2データの遺伝子発現の2群間比較を、edgeRで行う。ex402.v2と同じデータセットを使い、ここでは、一般化線形モデル(GLM)で、感染サンプルと対照区(mock)の間でdifferential expressionを示す遺伝子(DEG)を同定する。

以下、すべてR環境。

Import library

> library(edgeR)

Import data

arab2.txt データを読み込む。arab2.txtは data/SS/arab2.txt にある。

> x <- read.delim("arab2.txt", row.names=1)

# ... dat中身の確認作業 ...

2グループ、各3繰り返し実験、という実験デザインを定義する。m1, m2, m3 が mock (control), h1, h2, h3 が　Pseudomonas感染。

> treat <- factor(c("M", "M", "M", "H", "H", "H"))
> treat <- relevel(treat, ref="M")

edgeRのDGEList関数でカウントデータを読み込む。

> y <- DGEList(x, group=treat)

Filtering genes by expression level

発現量が全体的にゼロに近い遺伝子を除去する。ここではCPM(count-per-million)が0.5以上のサンプルが３つ以上、をフィルタリングの条件とする。

> keep <- rowSums(cpm(y) > 0.5) >= 3
> table(keep)
keep
FALSE  TRUE 
 6663 19558

y <- y[keep, , keep.lib.sizes=FALSE]

参考）edgeRには発現量が低い遺伝子を判定するコマンド、filterByExprがあり、これを使っても良い。

参考）data.frame(y, keep)　を実行すると、どの遺伝子がkeepされたかわかる。

Normalization

TMM法で、ノーマライズする。normLibSizesを使う。(注：edgeRのバージョンが３系の場合は、calcNormFactors)

> y <- normLibSizes(y, method="TMM")

計算結果の確認

> y$samples
   group lib.size norm.factors
m1     M  1901107    1.0247736
m2     M  1917419    1.0537590
m3     M  3257400    0.8958933
h1     H  2126850    1.0297446
h2     H  1280838    1.1293858
h3     H  3524342    0.8887973

Design matrix for GLM

GLMにフィットさせる前に、デザインマトリックスを定義する。

> design <- model.matrix(~treat)

> rownames(design) <- colnames(y)
> design
   (Intercept) treatH
m1           1      0
m2           1      0
m3           1      0
h1           1      1
h2           1      1
h3           1      1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$treat
[1] "contr.treatment"

DE testing

Estimate dispersion

負の二項分布に基づいた確率分布モデルのdispersionを推定する。

> y <- estimateDisp(y, design, robust=TRUE)

Quasi-likelihood法を用いて一般化線型モデルへの当てはめを行う。glmQLFit関数を使う。

> fit <- glmQLFit(y, design)

結果を確認しよう。

> head(fit$coefficients)
          (Intercept)      treatH
AT1G01010   -10.67114  0.18271042
AT1G01020   -10.90379  0.11550050
AT1G01030   -11.54218  0.37257375
AT1G01040   -10.53037 -0.04228772
AT1G01050   -10.36275 -0.16222071
AT1G01060   -12.81549  0.71390043

Rでの線形モデルあてはめ結果を解釈する際の注意点: (Intercept)はこの場合、treatMの値であり、treatHはtreatMとの「違い」を示している。このような結果の表示は"contrast"をRが設定して線形モデルを解くことに由来する。上でtreat <- relevel(treat, ref="M")というように説明変数に「水準（levels）]を設定している。

それでは、いよいよDEGのstatistical testである。glmQLFTest関数を使う。

> qlf <- glmQLFTest(fit, coef=2) 
> topTags(qlf)
Coefficient:  treatH 
             logFC    logCPM        F       PValue         FDR
AT2G19190 4.538369  7.388102 154.4349 1.819165e-07 0.001480407
AT2G44370 5.428778  5.211831 152.1667 1.952620e-07 0.001480407
AT3G55150 5.777979  4.916746 147.4347 2.270794e-07 0.001480407
AT1G51820 4.376266  6.385171 133.7755 3.607117e-07 0.001763700
AT2G39380 4.943932  5.785226 122.1599 5.543690e-07 0.002025755
AT2G39530 4.346822  6.721769 115.4041 7.246031e-07 0.002025755
AT3G46280 4.911814  8.119363 112.2486 8.252758e-07 0.002025755
AT5G64120 3.746759  9.703124 112.1520 8.286144e-07 0.002025755
AT4G12500 4.389182 10.442749 104.6792 1.143799e-06 0.002251886
AT1G51800 3.980086  7.722830 104.5308 1.151389e-06 0.002251886

注：glmQLFTest(fit, coef=2) のcoef=2 は統計テストの対象のモデルの係数。ここではdesignの第２係数＝treatH を指定。帰無仮説は係数＝０。つまり、「mockと感染で違いがない」が帰無仮説。

Dump the table to a text file

> tab <- topTags(qlf, n=nrow(qlf$table))
> write.table(tab$table, "qlf.treat.txt", sep="\t", quote=F)

Inspect DE result

edgeRに組み込まれている''decideTest''で DEGの数を集計 (古いedgeR ver.3系の場合、decideTestDGE)

> summary(decideTests(qlf))
       treatH
Down       22
NotSig  19391
Up        145

Rを終了し、UNIXコマンドライン環境に戻って、先ほど保存したDEGテーブル qlf.treat.txt の中身を確認しよう。(lessなどのコマンドを使う)。

ex403-2: Differential expression analysis with edgeR (GLM; considering batch effect)

arab2データの遺伝子発現の2群間比較を、edgeRで行う。ex402.v2と同じデータセットを使い、ここでは、一般化線形モデル(GLM)で、感染サンプルと対象区の間でdifferential expressionを示す遺伝子(DEG)を同定する。ex403-1からさらに発展させ、ここではbatch effectを考慮する。

上のex403-1では、実験群 vs 対照群のペアワイズ比較を行った。それぞれ繰り返しが３つある実験であるが、原著論文を読むと、実験1 (m1, h1), 実験2 (m2, h2), 実験3(m3, h3)は異なる時期に取得されており、バッチエフェクトが存在する可能性がある。ex403-2では、バッチエフェクトによる発現変動量を推定し、時期による影響を統計的に推定しうる範囲で除いて、感染処理によるDEGを同定する。

以下全てR環境。

> library(edgeR)
> x <- read.delim("arab2.txt", row.names=1)
> treat <- factor(c("M", "M", "M", "H", "H", "H"))
> treat <- relevel(treat, ref="M")
> Time <- factor(c(1,2,3,1,2,3))
> keep <- rowSums(cpm(y) > 0.5) >= 3
> y2 <- DGEList(x, group=treat)
> y2 <- y2[keep, , keep.lib.sizes=FALSE]
> y2 <-normLibSizes(y2, method="TMM")

ここまではex403.v2-1と同じ。（> Time <- factor(c(1,2,3,1,2,3))　以外は）

MDS plot

MDS plot はサンプル間の発現パターンの類似性の概略を把握するのに便利である。

plotMDS(y2)

dim1 （Ｘ軸）に沿って、左側にバッチ2が寄っており、バッチエフェクトがありそうだと見当がつく。

Design matrix for GLM

GLMにフィットさせる前に、デザインマトリックスを定義する。

> design2 <- model.matrix(~Time+treat)
> rownames(design2) <- colnames(y2)
> design2
design2
   (Intercept) Time2 Time3 treatH
m1           1     0     0      0
m2           1     1     0      0
m3           1     0     1      0
h1           1     0     0      1
h2           1     1     0      1
h3           1     0     1      1
attr(,"assign")
[1] 0 1 1 2
attr(,"contrasts")
attr(,"contrasts")$Time
[1] "contr.treatment"

attr(,"contrasts")$treat
[1] "contr.treatment"

DE testing

estimate dispersion

> y2 <- estimateDisp(y2, design2, robust=TRUE)
> y2$common.dispersion
[1] 0.07753743

model fitting

> fit2 <- glmQLFit(y2, design2)
> head(fit2$coefficients)
          (Intercept)      Time2      Time3      treatH
AT1G01010   -10.96162  0.7369950 -0.1554951  0.21596161
AT1G01020   -10.85267  0.1467467 -0.4366442  0.14534329
AT1G01030   -11.63777  0.5221485 -0.2430494  0.27373713
AT1G01040   -10.29519 -0.5574841 -0.2542314 -0.03036512
AT1G01050   -10.39803  0.1959288 -0.1145022 -0.16614196
AT1G01060   -17.25874  5.3363641  3.4094136  0.86928320

今回の線形モデルでは:

(Intercept): Time1のMの値が基準となっている。

coefficinetの大きさは説明変数（及び各水準）の大きさである。次のプロットで横軸にTime2による影響の大きさ、縦軸にtreatHの大きさの違いを遺伝子毎に可視化してみる。

> plot(fit2$coefficients[,2], fit2$coefficients[,4], cex=0.1, pch=16, xlab="Time2", ylab="treatH")

統計的に有意なバッチエフェクトを持つ遺伝子があるかどうかチェック。design2において、バッチエフェクトつまりTimeに関する係数が第２、第３係数に相当することに留意すると、

> qlf.time <- glmLRT(fit2, coef=2:3)
> sum(p.adjust(qlf.time$table$PValue) < 0.05)
[1] 937

発現量の違いがバッチの違いで説明できる遺伝子が、937個もあるので、バッチエフェクトがあると判定できる。従ってバッチエフェクトを考慮してDEGを同定すべきであるといえる。

それではいよいよ、treatの有無でのDEGの同定である。design2において、treatHが第４係数であることに留意して、

> qlf.treat <- glmQLFTest(fit2, coef=4)
> topTags(qlf.treat)
Coefficient:  treatH 
             logFC    logCPM        F       PValue          FDR
AT5G48430 6.350583  6.721830 300.4095 4.144268e-67 1.077841e-62
AT3G46280 4.825892  8.134882 252.6725 9.240655e-57 1.201655e-52
AT2G19190 4.540604  7.398586 228.1022 1.988262e-51 1.723691e-47
AT4G12500 4.363585 10.445676 224.2540 1.361825e-50 8.854583e-47
AT2G39530 4.381132  6.731603 201.2569 1.351163e-45 7.028209e-42
AT1G51800 4.012131  7.732901 194.7135 3.575137e-44 1.549703e-40
AT2G17740 4.245578  6.789176 193.7168 5.888743e-44 2.187920e-40
AT2G39200 3.938831  9.023060 190.9090 2.402156e-43 7.809409e-40
AT2G39380 4.985793  5.800196 189.8356 4.111883e-43 1.188243e-39
AT1G51820 4.377731  6.393980 189.4155 5.074753e-43 1.319842e-39

Dump the table to a text file

> tab2 <- topTags(qlf.treat, n=nrow(qlf.treat$table))
> write.table(tab2$table, "qlf.treat.txt", sep="\t", quote=F)

Inspect DE result

edgeRに組み込まれている''decideTest''で DEGの数を集計 (古いedgeR 3系の場合、decideTestsDGE)

>  summary(decideTests(qlf.treat))
       treatH
Down      937
NotSig  23849
Up       1222

バッチエフェクトを考慮しないex403-1 よりも、より多くのDEGが同定された。

(発展)バッチエフェクトの影響をより深く検討する

上記の通り、バッチエフェクトを考慮しないex403-1 よりも、より多くのDEGが同定された。では、fit、fit2およびlrt.treat2の内容を掘り下げて、バッチエフェクトを考慮した影響をより深く検討してみよう。

以降、同様なプロットを複数描画するので、描画作業の関数を作成する。また、各モデル、各要因で有意に変動した遺伝子の番号を取得しておく。

plotAllPoints <- function(xlab, ylab){
  plot(
    x=x_values, y=y_values,
    xlim=c(-10, 10), ylim=c(-7, 10), 
    xlab=xlab, ylab=ylab, 
    col="grey80", cex=0.5, asp=1
  )
  abline(a=0, b=1, col="blue")
  abline(h=c(-1,1), lty=2); abline(v=c(-1,1), lty=2)
  abline(h=0); abline(v=0)
}

specifiedPoints <- function(indexes, col){
  points(
    x=x_values[indexes], y=y_values[indexes],
    col=col, pch=16, cex=0.5
  )
}

lrt_q0.05_indexes <- which(p.adjust(lrt$table$PValue, "BH") < 0.05)
lrt2_q0.05_indexes <- which(p.adjust(lrt.treat2$table$PValue, "BH") < 0.05)
lrt2_only_significant_indexes <- setdiff(lrt2_q0.05_indexes, lrt_q0.05_indexes)
time_q0.05_indexes <- which(p.adjust(lrt.time$table$PValue) < 0.05)
lrt2_time_treat_significant_indexes <- intersect(
  lrt2_q0.05_indexes,
  time_q0.05_indexes
)

まずは、バッチエフェクトを考慮しなかったモデル、考慮したモデルのtreatHの係数を比較する。

バッチエフェクトを考慮しなかったモデルのtreatHの係数（横軸）、バッチエフェクトを考慮したモデルのtreatHの係数（縦軸）を散布図として比較してみる。

fit_coefs  <- fit$coefficients
fit2_coefs <- fit2$coefficients

x_values   <- fit_coefs[, "treatH"]
y_values   <- fit2_coefs[, "treatH"]
plotAllPoints(xlab="treatH: treat", ylab="treatH: Time + treat")

両者はおおむね同様な値が各遺伝子について推定されている。

なお、破線は自然対数1, -1の値を示している。青の線は切片0、傾き1の直線である。

次に、バッチエフェクトを考慮したモデルのtreatHとTime2の係数を取り上げて分析する。

以下の図では、横軸がtreatHの係数、縦軸はtreatHの効果にTime2の効果が同時に存在する例として、treatH係数にTime2係数を足した値を示している。

青の線は切片0、傾き1の直線である。Time2の影響がない（係数がゼロ）の遺伝子はこの直線状にプロットされ、縦軸方向にどちらにどのくらい移動した値にプロットされているかで各遺伝子のTime2の影響を読み取ることができる。

x_values   <- fit2_coefs[, "treatH"]
y_values   <- fit2_coefs[, "treatH"] + fit2_coefs[, "Time2"]

plotAllPoints(xlab="treatH", ylab="treatH + Time2")
specifiedPoints(lrt_q0.05_indexes,  "yellow")
legend("topleft", legend=c("treatH: treat"), pch=16, col=c("yellow"), bty='n')

黄色でハイライトされているのが、バッチエフェクトを考慮「しなかった」モデルのtreatHがq<0.05であった遺伝子である。 434遺伝子がハイライトされているが、ほとんどが傾き１の直線周辺にある。Time2の影響が存在する遺伝子はtreatHが有意と検出しにくい。

次はプロットされるデータは同じだが、バッチエフェクトを考慮「した」モデルのtreatHがq<0.05であった遺伝子をハイライトして、プロットする。

plotAllPoints(xlab="treatH", ylab="treatH + Time2")
specifiedPoints(lrt2_q0.05_indexes, "cyan")
legend("topleft", legend=c("treatH: Time+treat"), pch=16, col=c("cyan"), bty='n')

上記2つのプロットのハイライト遺伝子を同時にプロットすると、どのような係数の遺伝子がバッチエフェクトを考慮すると有意になったのかがわかりやすい。

plotAllPoints(xlab="treatH", ylab="treatH + Time2")
specifiedPoints(lrt2_q0.05_indexes, "cyan")
specifiedPoints(lrt_q0.05_indexes,  "yellow")
legend("topleft", legend=c("treatH: Time+treat", "treatH: treat"), pch=16, col=c("cyan", "yellow"), bty='n')

バッチエフェクトを考慮した場合の方が

係数のより小さい遺伝子
傾き１の直線から離れた遺伝子

についても処理間で有意な差が検出されている。

最後に、バッチエフェクト（Time)の係数が有意だった遺伝子もハイライトしてプロットする。

plotAllPoints(xlab="treatH", ylab="treatH + Time2")
specifiedPoints(lrt2_q0.05_indexes, "cyan")
specifiedPoints(lrt_q0.05_indexes,  "yellow")
specifiedPoints(time_q0.05_indexes, "magenta")
specifiedPoints(lrt2_time_treat_significant_indexes, "green")
legend(
  "topleft", 
  legend=c(
    "treatH: Time+treat", "treatH: treat", "Time: Time+treat",
     expression(paste(paste("Time", intersect(treatH)), ": ", "Time+treat", sep=""))
  ), 
  pch=16, col=c("cyan", "yellow", "magenta", "green"), bty='n'
)

Time3で有意になった遺伝子も含めてハイライトしているため、傾き１の直線周辺にTimeが有意だった遺伝子（Time2の影響が小さい遺伝子）もハイライトされている。着目すべき点として、バッチエフェクトを考慮すると、Timeの影響が有意な遺伝子でもTimeの影響を切り分けることにより、treatHの影響が有意な遺伝子（緑でハイライト）を検出できることである。バッチエフェクトを推定して取り除き、着目している/していない要因それぞれによるデータのばらつきを考慮してDEG解析を行うことによって、このような鋭敏な解析が可能になる。

Notes

edgeR のver.4へのアップでとに伴い、関数名が変更になったものがある。本文ではv3, v4系両方に対応するように記述した。

Revision history

2024-7-30

edgeR v3, v4 系両方に対応するよう記述を修正した。

2024-2-25

likelihood ratio tests法からQuasi-likelihood法へ全面的に変更。
動作確認 R:v4.3.2, edgeR: v4.0.16 on MacOS Sonoma 14.3.1

2023-2-24

本文の軽微な修正。
動作確認 R:v4.0.3, edgeR: v3.32.1 on bias5 + X11

2022-8-27

動作確認 R:v4.0.3, edgeR: v3.32.1

2021-09-12

Added the inspection on the difference of DEGs between models with or without considering the Time batch effects.

2021-3-XX

Minor modifications.

ex403.v2 - nibb-gitc/gitc2025mar-rnaseq GitHub Wiki

ex403.v2-1: Differential expression analysis with edgeR (GLM)

Import library

Import data

Filtering genes by expression level

Normalization

Design matrix for GLM

DE testing

Estimate dispersion

Inspect DE result

ex403-2: Differential expression analysis with edgeR (GLM; considering batch effect)

MDS plot

Design matrix for GLM

DE testing

Inspect DE result

(発展)バッチエフェクトの影響をより深く検討する

Links

Notes

Revision history

2024-7-30

2024-2-25

2023-2-24

2022-8-27

2021-09-12

2021-3-XX