Advertisement

Affymetrix芯片分析:获取差异表达基因系列二_Moderated t-test

阅读量:

library(limma)

library(tcltk)

library(affy)

filters <- matrix(c("CEL file", ".[Cc][Ee][Ll]", "All", ".*"), ncol = 2, byrow = T)

cellFiles <- tkSelectFiles(title = "Select CELs", multiple = TRUE, filters = filters, initialIndex = 1);

data.raw <- ReadAffy(filenames = cel.files)

design <- model.matrix(~-1+factor(c(1,1,2,2)))%构建实验设计矩阵(实验组2个和对照组2个)

colnames(design) <- c("control", "LPS") %实验组和对照组的列名称

eset.rma <- rma(data.raw)%调用RMA算法对数据进行预处理

eset<-exprs(data.raw)%构建样品的基因表达矩阵

fit <- lmFit(eset, design)%线性模型拟合

contrast.matrix <- makeContrasts(control-LPS, levels=design) %建立对比模型以分析两个实验条件下的基因表达差异

fit1<- contrasts.fit(fit, contrast.matrix)%根据对比模型进行差值计算

fit2 <- eBayes(fit1) %贝叶斯检验

results<-decideTests(fit2, method="global", adjust.method="BH", p.value=0.01, lfc=1.5)通过基于P-value的结果筛选获得差异表达基因

summary(results)

vennCounts(results)

vennDiagram(results)

heatDiagram(results,fit2$coef)

该代码段通过调用topTable函数对fit2进行分析,并设置相关参数以实现预期效果;具体设置包括:将系数设定为1;输出结果数量设为1万;采用Benjamini-Hochberg方法校正p值;按字母顺序对结果进行排序;并在必要时维持字母顺序

write.table(x, file="limma.xls", row.names=F, sep="\t")

y <- x[xP.Value < 0.01 & (xlogFC > 1.5 | xlogFC < -1.5& xAveExpr > 10),]

全部评论 (0)

还没有任何评论哟~