ex12 - nibb-gitc/gitc2018july-rnaseq GitHub Wiki
ex12-1: Differential expression analysis with edgeR (GLM)
arab2データの遺伝子発現の2群間比較を、edgeRで行う。[ex5]と同じデータセットを使い、ここでは、一般化線形モデル(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
ex12-2: Differential expression analysis with edgeR (GLM; considering batch effect)
arab2データの遺伝子発現の2群間比較を、edgeRで行う。[ex5]と同じデータセットを使い、ここでは、一般化線形モデル(GLM)で、感染サンプルと対象区の間でdifferential expressionを示す遺伝子(DEG)を同定する。ex12-1からさらに発展させ、ここではbatch effectを考慮する。
上のex12-1では、実験群 vs 対照群のペアワイズ比較を行った。それぞれ繰り返しが３つある実験であるが、原著論文を読むと、実験1 (m1, h1), 実験2 (m2, h2), 実験3(m3, h3)は異なる時期に取得されており、バッチエフェクトが存在する可能性がある。ex12-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
1 1 0 0 0
2 1 1 0 0
3 1 0 1 0
4 1 0 0 1
5 1 1 0 1
6 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")
統計的に有意なバッチエフェクトを持つ遺伝子があるかどうかチェック。
> lrt.time <- glmLRT(fit2, coef=2:3)
> sum(p.adjust(lrt.time$table$PValue) < 0.05)
[1] 943
発現量の違いがバッチの違いで説明できる遺伝子が、943個もあるので、バッチエフェクトがあると判定できる。従ってバッチエフェクトを考慮してDEGを同定すべきであるといえる。
それではいよいよ、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))
[,1]
-1 914
0 24088
1 1219
バッチエフェクトを考慮しないex12-1 よりも、より多くのDEGが同定された。