转录组不求人系列(六):limma包分析转录组芯片数据

转录组数据差异基因的分析是转录组的核心，针对不同的数据类型有不同的分析方法，之前我们也提到过。这次要做的是芯片数据的分析，采用limma包，方法是通用的（适用于mRNA以及小RNA的芯片数据）。

一般数据挖掘分析都是下载表达矩阵进行分析，当然也可以从原始数据开始，感兴趣的可以查看limma包的详细介绍。

读入表达矩阵，对其log转化（至于为和要log化，可以百度下，如果作者上传的数据是已经log过的，则无需处理），简单QC一下，用boxplot看看数据样本之间的一致性。

setwd("E:/生物信息学")A <- read.table("GSE125424_series_matrix.txt",header = T, row.names = 1,sep = "\t",comment.char = "!")A <- log2(A)boxplot(data.frame(A))

这个数据还可以，一致性挺好。但是通常时候，样本之间的差异还是挺大的，就需要用如下的代码进行标准化再进行差异分析。使用wateRmelon包进行标准化的处理。

#library("wateRmelon")#A<- betaqn(A)#boxplot(A)#write.table(A,file="norm.txt",sep="\t",quote=F)

limma包分析差异基因还是比较简单的，我们首先需要构建分组，并构建比较矩阵，这里只有两组，虽然默认就是两组比较，但是还是需要定义，明确是哪个比哪个。

library(limma)library(dplyr)f <- c(rep("control", 3), rep("Treat", 3)) %>% as.factor()desigN <- model.matrix(~ 0 + f)colnames(desigN) <- levels(f)fit = lmFit(A, desigN)contrast.matrix <- makeContrasts(Treat - control, levels = colnames(coef(fit)))contrast.matrix # Contrasts#Levels Treat - control# control -1#  Treat                 1

差异分析：

fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2) tempOutput <- topTable(fit2,coef=1,n=Inf,adjust="BH")head(tempOutput)
# logFC AveExpr t P.Value adj.P.Val B#1446537_at 3.848140 2.164047 22.71094 2.092108e-06 0.09435616 1.339931#1440714_at -3.426720 1.887167 -17.83508 7.254608e-06 0.09706875 1.172877#1429297_at -3.198857 1.599429 -17.80127 7.325561e-06 0.09706875 1.171274#1443604_at -3.069908 1.616981 -16.22950 1.176281e-05 0.09706875 1.086811#1421100_a_at -2.234104 1.150563 -16.19503 1.189145e-05 0.09706875 1.084715#1438524_x_at 2.943123 1.471561 15.71083 1.388776e-05 0.09706875 1.054005

最后将差异分析文件导出：

write.csv(tempOutput, file="DEGs.csv")

导出的差异分析结果就可以进行下一步的筛选了，设定不同的阈值筛选差异基因，以及基因功能富集等下游分析了。