【R>>limma】3分组差异分析

生信技能树公众号的干货多多，今天来学习下《3个分组的表达量矩阵的两两之间差异分析》这篇推文。

示例数据

rm(list = ls())
library(limma)
library(pheatmap)
y <- matrix(rnorm(10000*9),10000,9)
fivenum(y) #四分位间距
#> [1] -3.908523705 -0.667061851  0.007111134  0.676220939  4.581587505

(group=rep(LETTERS[1:3],each=3))
#> [1] "A" "A" "A" "B" "B" "B" "C" "C" "C"
y[sample(1:10000,1000),1:3]=y[sample(1:10000,1000),1:3]+2
y[sample(1:10000,1000),4:6]=y[sample(1:10000,1000),4:6]+2
y[sample(1:10000,1000),7:9]=y[sample(1:10000,1000),7:9]+2
pheatmap(y)

image.png

三组差异分析

design <- model.matrix(~0+group)
colnames(design) <- gsub("group","",colnames(design))
contr.matrix <- makeContrasts(AVSB = A-B,
                              AVSC = A-C,
                              BVSC = B-C,
                              levels = colnames(design))
contr.matrix
#>       Contrasts
#> Levels AVSB AVSC BVSC
#>      A    1    1    0
#>      B   -1    0    1
#>      C    0   -1   -1

vfit <- lmFit(y,design)
vfit <-  contrasts.fit(vfit,contrasts = contr.matrix)
efit <- eBayes(vfit)
plotSA(efit)

image.png

summary(decideTests(efit))
#>        AVSB AVSC BVSC
#> Down    322  277  297
#> NotSig 9377 9426 9389
#> Up      301  297  314

dt <- decideTests(efit)
summary(dt)
#>        AVSB AVSC BVSC
#> Down    322  277  297
#> NotSig 9377 9426 9389
#> Up      301  297  314
de.common <- which(dt[,1]!=0&dt[,2]!=0)
length(de.common)
#> [1] 165
vennDiagram(dt[,1:3],circle.col = c("royalblue","firebrick3","orange3"))

image.png

这3次差异分析的结果都是可以独立取出的：

colnames(efit)
#> [1] "AVSB" "AVSC" "BVSC"

AVSB <- topTable(efit,coef = 1,n=Inf,adjust.method = "fdr")
AVSC <- topTable(efit,coef = 2,n=Inf,adjust.method = "fdr")
BVSC <- topTable(efit,coef = 3,n = Inf,adjust.method = "fdr")

参考链接：

3个分组的表达矩阵的两两之间差异分析