差异分析

data=read.csv("second_magic.csv",header=T,row.names=1)

# 1.构建模拟的表达矩阵，实际处理时换成自己的表达矩阵即可

sd <- 0.3*sqrt(4/rchisq(100,df=4))

exprSet <- matrix(rnorm(100*6,sd=sd),100,6)

rownames(exprSet) <- paste("Gene",1:100)

colnames(exprSet) <- c(paste0("con-",1:3), paste0("G3-",1:3))

exprSet[1:2,4:6] <- exprSet[1:2,4:6] + 2

library(limma)

# 2.构建实验设计矩阵

group_list = c(rep("con",3), rep("G3",3))#对应sample的分组列表

# 这里根据实际的情况设置（表型）分组，对应表达矩阵的列：样本

design <- model.matrix(~0+factor(group_list))

design

colnames(design)=levels(factor(group_list))

rownames(design)=colnames(exprSet)

#给design矩阵加行名，行名为表达谱矩阵的列名，即sample

design

# 实验设计矩阵的每一行对应一个样品的编码，

# 每一列对应样品的一个特征。这里只考虑了一个因素两个水平，

# 如果是多因素和多水平的实验设计，会产生更多的特征，需要参考文档设计。

# 3.构建对比模型（对比矩阵），比较两个实验条件下表达数据

contrast.matrix<-makeContrasts(G3-con,levels = design)

#contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)

contrast.matrix ##这个矩阵声明，我们要把G3组跟con组进行差异分析比较

##### 差异分析

##4.  step1 线性模型拟合

fit <- lmFit(exprSet,design)

##    step2 根据对比模型进行差值计算

fit2 <- contrasts.fit(fit, contrast.matrix)

##5.  step3 贝叶斯检验

fit2 <- eBayes(fit2)

##6.  step4 生成所有基因的检验结果报告

tempOutput = topTable(fit2, coef=1, n=Inf)

##step5 用P.Value进行筛选，得到全部差异表达基因

dif <- tempOutput[tempOutput[, "P.Value"]<0.01,]

# 显示一部分报告结果

head(dif)