R语言GEO数据挖掘：步骤三：进行基因差异分析

Python016

R语言GEO数据挖掘：步骤三：进行基因差异分析,第1张

用limma包，这里注意，limma包是对基因芯片表达矩阵的分析，不能对逆转录RNAseq表达矩阵进行分析（因为数据特征不同），RNAseq需要用另一种方法

解读此表

但是上面的用法做不到随心所欲的指定任意两组进行比较，所有还有下一种方法

处理好了分组信息，再自定义比较元素

自定义函数进行比较

热土和火山图都是傻瓜式的，只要的前面得出的deg数据(也就是基因差异表达数据）是正确的

#mRNA表达矩阵与GROUP文件样式，heatmap样式见文章最后

library(limma)

 mRNA <- read.table("表达矩阵.txt",sep = "\t",header = T,comment.char = "!",encoding = "UTF-8")

#mRNA数据框行名为基因名，列命为样本名称

 group <- read.table("GROUP.txt",header=T,sep = "\t",encoding = "UTF-8")

 group_CP <- group$treat 

 m_design<- model.matrix(~0+factor(group_CP))

 colnames(m_design) = levels(factor(group_CP))

 rownames(m_design)= group$ID

 contrast.matrix<-makeContrasts("P-C",levels=m_design) #注意P-C顺序，实验组要在前面否则影响上下调结果

 m_fit <- lmFit(mRNA,m_design)

 m_fit <- contrasts.fit(m_fit, contrast.matrix)

 m_fit <- eBayes(m_fit)

 m_genlist <- topTable(m_fit, coef = 1, n=Inf)  #limma结果

#将表达矩阵与差异分析结果合并

  ID_REF <- rownames(m_genlist)

  m_genlist <- data.frame(ID_REF,m_genlist)

  ID_REF <- rownames(mRNA)

  mRNA <- data.frame(ID_REF,mRNA)

  test <-merge(mRNA,m_genlist,by = "ID_REF")

  result <- subset(test,P.Value<0.05)

  row.names(result) <- result[,1]

#绘制热图

heatmap <- result[2:(nrow(group)+1)]

annotation <- data.frame(Factor = factor(group$treat)) #标注样本的分组信息

rownames(annotation) <- colnames(heatmap)

library(pheatmap)

filename <- paste("文件名",".pdf",sep="")

pdf(filename)

pheatmap(heatmap,

        annotation=annotation,

        annotation_legend = TRUE,

        main=filename ,

        scale = "row",

        show_rownames = F,

        color = colorRampPalette(c("green","black","red"))(100))

dev.off()

#表达矩阵与GROUP文件如下所示