🧐 Limma | 三个组的差异分析怎么分析做呢！？~

rm(list = ls())
library(tidyverse)
library(limma)
library(GEOquery)

GSE194314 <- getGEO('GSE194314', destdir=".",getGPL = F)

exprSet <- exprs(GSE194314[[1]])

boxplot(log2(exprSet))

exprSet <- normalizeBetweenArrays(exprSet) %>% 
  log2(.)

boxplot(exprSet)

pdata <- data.frame(
  sample = colnames(exprSet), 
  group = rep(c("A", "B", "C"), each = 4)
)

DT::datatable(pdata)

group <- factor(pdata$group,levels = unique(pdata$group))

#===========整理分组数据为分组矩阵
design <- model.matrix(~0 + group)

rownames(design) <- colnames(exprSet)
colnames(design) <- levels(group)

# 差异比较矩阵
cont.matrix <- makeContrasts("A vs B" = A-B, 
                             "A vs C" = A-C, 
                             "B vs C" = B-C,
                             levels = colnames(design)
                             )
cont.matrix

#========== 定义阈值
logFCcutoff <- log2(0)
adjPvalueCutoff <- 0.3

#=========== 进行差异分析
fit <- lmFit(exprSet, design)

#=========== 针对给定的对比计算估计系数和标准误差
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
plotSA(fit2)

summary(decideTests(fit2, p.value = adjPvalueCutoff, lfc = logFCcutoff))

dt <- decideTests(fit2,p.value = adjPvalueCutoff, lfc = logFCcutoff)
de.common <- which(dt[,1]!=0 & dt[,2]!=0 & dt[,3]!=0)
length(de.common)
vennDiagram(dt[,1:3], circle.col=c("#3C5488E5", "#91D1C2E5", "#4DBBD5E5"))

colnames(fit2)

A_VS_B <- topTreat(fit2, coef=1, n=Inf)
A_VS_C <- topTreat(fit2, coef=2, n=Inf)
B_VS_C <- topTreat(fit2, coef=3, n=Inf)

DT::datatable(A_VS_B)