使用limma包进行差异基因分析时，做最多的是两分类的，例如control组和disease组，但也会碰到按照序列进行的分组。这时，如果逐一使用两两比较求差异基因则略显复杂。其实开发limma包的大神们已经替我们考虑到。我自己當下limma包的PDF，仔细研读并将代码运行后分享给有需要的同道，相互学习。
library(limma)##加载包load("GSE37761exp_groupfile-00.RData")#加载匹配过的基因表达矩阵eset=exp ##将基因表达矩阵赋值给esettargets<-read.csv("group_file.csv",row.names = 1)##读入样本数据，包括两列，#第一列GSM号，第二列为样本分组#该数据集中实际样本分组中存在24h_1、24h_5、24h_21等，然而这不符合R的命名原则targets$Target=gsub("_",".",targets$Target)##将"_"替换成“.”，也可以不替换##该数据集中实际样本分组中存在24h_1、24h_5、24h_21等，然而这不符合R的命名原则，所以在没个分类前加一个“F”，具体自己定targets$Target=c(paste0("F",c(targets$Target),collapse = NULL,sep=""))colnames(targets)=c("FileName","Target")#更改列名，为了和limma包中的一致lev<-unique(targets$Target)##使用unique（）函数进行去重f <- factor(targets$Target, levels=lev) design <- model.matrix(~0+f) #样本矩阵colnames(design) <- lev #更改列名为levels名###两两之间的比较，求差异基因使用topTable函数cont.wt <- makeContrasts("F12h.1-Fcontrol",#举例levels=design) fit <- lmFit(eset, design)fit2 <- contrasts.fit(fit, cont.wt) fit2 <- eBayes(fit2) tT=topTable(fit2, adjust="BH",sort.by="logFC",n=Inf)tT = subset(tT, select=c("adj.P.Val","P.Value","logFC"))colnames(tT)=c("FDR","P.Value","logFC")#####后面可以设置fdr及logFC的阈值进行筛选，进而得到差异基因列表，此处略。##两组以上的比较，寻找差异基因，使用topTableF函数cont.wt <- makeContrasts("F12h.1-Fcontrol",##举例"F24h.1-F12h.1","F5d.1-F24h.1","F10d.1-F5d.1",levels=design) fit <- lmFit(eset, design)fit2 <- contrasts.fit(fit, cont.wt) fit2 <- eBayes(fit2) tT=topTableF(fit2, adjust="BH",sort.by="F",n=Inf)##tT = subset(tT, select=c("adj.P.Val","P.Value","logFC")) 这一条和下一条代码需要根据自己需要进行更改，topTableF函数得到的结果中包括F值，而logFC则是没两组的比较，均得出了##colnames(tT)=c("FDR","P.Value","logFC")
因为不支持PDF上传，最后附上limma包usersguide.pdf的下载地址。/packages/3.2/bioc/vignettes/limma/inst/doc/usersguide.pdf
学习过程的就是分享的过程，分享的过程也是交流的过程，交流的过程就是进步的过程。
如果觉得《python 分析两组数据的差异_R语言limma包差异基因分析（两组或两组以上）》对你有帮助，请点赞、收藏，并留下你的观点哦！
