ex403 - nibb-gitc/gitc2021sep-rnaseq GitHub Wiki

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

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

以下、すべてR環境。

Import library

> library(edgeR)

Import data

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

> dat <- 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(dat, group=treat)

Normalization

TMM法で、ノーマライズする。calcNormFactorsを使う。

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

計算結果の確認

> y$samples
   group lib.size norm.factors
m1     M  1902032    1.0399197
m2     M  1934029    1.0611305
m3     M  3259705    0.8841923
h1     H  2129854    1.0266944
h2     H  1295304    1.1412144
h3     H  3526579    0.8747345

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

GLMにフィットさせ、係数を確認する。

> y <- estimateDisp(y, design)
> y$common.dispersion
[1] 0.3508238

ここまでは"classic / pairwise" 法と同じ。ここからGLM特有の処理となる。

> fit <- glmFit(y, design)

> head(fit$coefficients)
          (Intercept)      treatH
AT1G01010   -10.68003  0.18653587
AT1G01020   -10.91423  0.12030392
AT1G01030   -11.54861  0.36889215
AT1G01040   -10.53809 -0.03178354
AT1G01050   -10.36875 -0.16096270
AT1G01060   -12.81986  0.72272916

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

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

> lrt <- glmLRT(fit, coef=2) 
> topTags(lrt)
Coefficient:  treatH 
             logFC    logCPM       LR       PValue          FDR
AT2G19190 4.548601  7.389099 98.88445 2.676723e-23 7.018636e-19
AT2G44370 5.435579  5.212357 90.98708 1.446132e-21 1.895952e-17
AT3G55150 5.778102  4.915724 88.87656 4.202332e-21 3.672978e-17
AT4G12500 4.391107 10.436950 86.43435 1.444542e-20 9.302502e-17
AT3G46280 4.925035  8.124269 85.73574 2.056619e-20 9.302502e-17
AT1G51820 4.388004  6.384948 85.66768 2.128638e-20 9.302502e-17
AT2G39380 4.949054  5.787376 84.02191 4.893221e-20 1.832931e-16
AT5G64120 3.755076  9.700968 81.02047 2.233916e-19 7.321938e-16
AT2G39530 4.349719  6.722165 80.31178 3.197579e-19 9.315968e-16
AT4G12490 3.877100 10.199985 79.66823 4.428640e-19 1.161234e-15

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

Dump the table to a text file

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

Inspect DE result

edgeRに組み込まれている''decideTestDGE''で DEGの数を集計

> summary(decideTestsDGE(lrt))
       treatH
Down       96
NotSig  25787
Up        338

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

arab2データの遺伝子発現の2群間比較を、edgeRで行う。ex402と同じデータセットを使い、ここでは、一般化線形モデル(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)
> dat <- 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))
> y2 <- DGEList(dat, group=treat)
> y2 <-calcNormFactors(y2, method="TMM")

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)
> y2$common.dispersion
[1] 0.07757116

model fitting

> fit2 <- glmFit(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に関する係数が第２、第３係数に相当することに留意すると、

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

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

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

> lrt.treat2 <- glmLRT(fit2, coef=4)
> topTags(lrt.treat2)
Coefficient:  treatH 
             logFC    logCPM       LR       PValue          FDR
AT5G48430 6.351014  6.712737 272.5565 3.145880e-61 8.248811e-57
AT3G46280 4.826138  8.125943 225.0412 7.191456e-51 9.428359e-47
AT2G19190 4.540899  7.389687 204.8646 1.812742e-46 1.584397e-42
AT4G12500 4.363887 10.436555 197.3463 7.924090e-45 5.194439e-41
AT2G39530 4.381421  6.722612 182.5254 1.361593e-41 7.140464e-38
AT2G39380 4.986071  5.791255 175.5897 4.450777e-40 1.729617e-36
AT2G17740 4.245929  6.780218 175.5166 4.617412e-40 1.729617e-36
AT1G51800 4.012461  7.723948 174.1586 9.139875e-40 2.995708e-36
AT1G51820 4.377955  6.384974 172.8187 1.792941e-39 5.223634e-36
AT1G51850 5.372491  5.442450 169.2718 1.067128e-38 2.798116e-35

Dump the table to a text file

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

Inspect DE result

edgeRに組み込まれている''decideTestDGE''で DEGの数を集計

> summary(decideTestsDGE(lrt.treat2))
       treatH
Down      914
NotSig  24088
Up       1219

バッチエフェクトを考慮しない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解析を行うことによって、このような鋭敏な解析が可能になる。

Links

Revision history

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.