【r<-芯片处理】基因芯片基因差异表达分析流程示例与讨论

交流地址：https://github.com/ShixiangWang/MessageBoard/issues/13

基因芯片的差异表达分析主要有 构建基因表达矩阵、构建实验设计矩阵、构建对比模型（对比矩阵）、线性模型拟合、贝叶斯检验和生成结果报表 六个关键步骤。

下面是模拟的一个示例：

# Simulate gene expression data for 100 probes and 6 microarraexprSets
# MicroarraexprSet are in two groups
# First two probes are differentiallexprSet expressed in second group
# Std deviations varexprSet between genes with prior df=4

# 构建模拟的表达矩阵，实际处理时换成自己的表达矩阵即可
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)
# 构建实验设计矩阵
group_list = c(rep("con",3), rep("G3",3))
# 这里根据实际的情况设置（表型）分组，对应表达矩阵的列：样本

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

colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
# 实验设计矩阵的每一行对应一个样品的编码，
# 每一列对应样品的一个特征。这里只考虑了一个因素两个水平，
# 如果是多因素和多水平的实验设计，会产生更多的特征，需要参考文档设计。

# 构建对比模型，比较两个实验条件下表达数据
contrast.matrix<-makeContrasts(G3-con,levels = design)
#contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix ##这个矩阵声明，我们要把G3组跟con组进行差异分析比较


##### 差异分析
##step1 线性模型拟合
fit <- lmFit(exprSet,design)
##step2 根据对比模型进行差值计算 
fit2 <- contrasts.fit(fit, contrast.matrix) 
##step3 贝叶斯检验
fit2 <- eBayes(fit2) 
##step4 生成所有基因的检验结果报告
tempOutput = topTable(fit2, coef=1, n=Inf)
##step5 用P.Value进行筛选，得到全部差异表达基因
dif <- tempOutput[tempOutput[, "P.Value"]<0.01,]
# 显示一部分报告结果
head(dif)

参考：

• 用limma对芯片数据做差异分析
• Bioconductor分析基因芯片数据
• limFit函数文档 ?limFit()

更新资料：

• 差异分析是否需要比较矩阵

相关问题请下面留言讨论。