转录组数据差异基因的分析是转录组的核心,针对不同的数据类型有不同的分析方法,之前我们也提到过。这次要做的是芯片数据的分析,采用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") 导出的差异分析结果就可以进行下一步的筛选了,设定不同的阈值筛选差异基因,以及基因功能富集等下游分析了。 |
|