limma差异分析

多组比较

setwd("C:\\Users\\HUAWEI\\Desktop\\proteomic_correct\\bacteria\\limma差异分析")
library(limma)
library(edgeR)
library(dplyr)

# 加载数据
rm(list = ls())
exprSet <- read.csv("bacteria_filled_scaled.csv",row.names = 1)
group = rep(c("CR_hv","CS_hv","CR_c","CS_c"),c(10,10,5,5))

# -------------------------------------------------------------------------

# 添加分组信息 
group <- factor(group,levels=c("CR_hv","CS_hv","CR_c","CS_c") )
design <- model.matrix(~0+group)    
colnames(design) = c("CR_hv","CS_hv","CR_c","CS_c")

# -------------------------------------------------------------------------

# 多组差异分析，多出makeContrasts(),contrasts.fit()

fit = lmFit (exprSet, design ) 

contrast.matrix = makeContrasts(CR_hv-CS_hv,CR_hv-CR_c,CS_hv-CS_c,CR_c-CS_c,
                                levels=design)

fit2 = contrasts.fit(fit, contrast.matrix)

fit2 = eBayes(fit2) 
plotSA(fit2)

allDiff2 = topTable(fit2, 
                    adjust = "BH",number=Inf)
head(allDiff)

results <- decideTests(fit2,lfc=log2(1.5))
summary(results)
CRhv_CShv<-topTreat(fit2,coef=1,n=Inf)
CRhv_CRc<-topTreat(fit2,coef=2,n=Inf)
CShv_CSc<-topTreat(fit2,coef=3,n=Inf)
CRc_CSc<-topTreat(fit2,coef=4,n=Inf)
vennDiagram(results[,c(1,4)])

data <- CRc_CSc
select.log2FC <- abs(data$logFC)>log2(1.5)
select.Pvalue <- data$adj.P.Val<0.05
select.vec <- select.log2FC & select.Pvalue
table(select.vec)

data$change <- "Normal"
data$change[data$logFC>=log2(1.5) & data$adj.P.Val<0.05]="Up"
data$change[data$logFC<=-log2(1.5) & data$adj.P.Val<0.05]="Down"
write.csv(data,"CRc_CSc_DEG.csv")