1、下载GSE115850数据

2、读取数据且处理数据，数据格式是这样的：

1c02eb8a0eb360136b8c88dd9a3ab45.png

接下来是要把行名改成基因名，又因为行名只能是唯一的不能有重复，所以要把多个Ensembl对应的同一个基因Symbol，取平均值(一般都是取平均)，这一步代码如下：

data<-read.delim('./GSE115850.txt') ##读取数据
head(data)  ##查看数据头6行
rownames(data)<-data$Ensembl  ##把Ensembl变为行名
data<-subset(data,select = -1)  ##把Ensembl一行可以删掉了
data<-aggregate(.~Symbol,data,mean)  ##把重复的Symbol取平均值
rownames(data)<-data$Symbol  ##把Symbol变为行名
data<-subset(data,select = -1)  ##删除Symbol行
data<-log2(data+1) ##对表达矩阵log2处理

处理完数据后格式是这样的：

2c2fa06f0eaffef11f179d8c121de02.png

接下来是用limma包对这个表达矩阵做差异分析，列名就是它的分组信息

#加载R包
library(limma)
library(dplyr)
library(pheatmap)
library(ggplot2)
#创建分组信息
phe<-data.frame(sample=colnames(p),group=c(rep('ctr',3),rep('gcm',3))) 
rownames(phe)<-phe$sample #改行名
phe<-subset(phe,select = -1)
design<-model.matrix(~0+factor(phe$group))
rownames(design)<-rownames(phe)
colnames(design)<-levels(factor(phe$group))
contrast.matrix<-makeContrasts('ctr-gcm',levels=design)

常规的limma包差异分析

fit<-lmFit(p,design)
fit2<-contrasts.fit(fit,contrast.matrix)
fit2<-eBayes(fit2)
DEG<-topTable(fit2,coef=1,n=Inf) %>% na.omit()
write.table(DEG,file='./GSE115850allgene.txt',sep='\t')
a<-DEG[which(DEG$P.Value<0.05&abs(DEG$logFC)>0.58),]##筛选FC>1.5和p<0.05
write.table(a,file='./GSE115850deg.txt',sep='\t')

3、输出的差异基因列表

8f70a3d175e537c11427ff08458d7da.png

4、开始画热图

exp<-p[intersect(rownames(a),rownames(p)),]
pheatmap<- pheatmap(exp,cluster_cols=F,annotation_col=phe,show_colnames=F,scale='row',show_rownames=F)
ggsave(file='DEG_heatmap.pdf',pheatmap,width=30,height=30,units='cm')

5、结果展示，效果蛮不错

82ed70554ac9aea4295be270cdb09cb.png