全部流程来自：GEO数据库挖掘—生信技能树B站视频，建议去看原文！

第一步：找到相关的GEO数据集（文献/搜索），以胃癌gastric cancer为例

可去文献中查找，用于练习

第二步：运行R包GEOquery获取数据（非常看网速，尽量下载下一点的包）

library(GEOquery)

eSet <- getGEO("GSE118916",

              destdir = '.',   #下载在当前目录

              getGPL = F)  #平台信息不要

#查看下载得到的文件大小，与GEO网页中标识的是否一致

第三步：得到表达矩阵

#使用列表取子集的方法提取eSet列表里的第一个元素：eSet[[1]]；并使用exprs函数把它转化成矩阵：exp <- exprs(eSet[[1]])

eSet[[1]]

exp <- exprs(eSet[[1]])

exp[1:4,1:4]#查看前四行

#代码逻辑，把probe_id转换为--->symbol ID/entrez ID

#分两步走：过滤probe_id，得到每个基因所对应的唯一的probe_id

得到probe_id与symbol ID这件的转换关系

ID转换的第一步必须要加载特定的R包，下载哪个包，需要根据GPL来定

eSet[[1]]@annotation #查看GPL平台信息，即芯片信息

#Google直接搜索相关R包，安装并调用

#如果有对应的R包

library(R包)

#如果找不到R包

library(Biobase)

library(GEOquery)

gpl<-getGEO('GPL15297',destdir=".")

colnames(Table(gpl))

head(Table(gpl)[,c(1,15,19)])## you need to check this , which column do you need

probe2symbol=Table(gpl)[,c(1,15)]

write.csv(Table(gpl)[,c(1,15)],"GPL15207.csv")

tail(sort(table(ids$symbol)))#看高频数的样本

table(rownames(exp)%in%ids$probe_id)#查看exp与id是否一致

exp=exp[rownames(exp)%in%ids$probe_id,]#去除不一致的项

ids=ids[match(rownames(exp),ids$probe_id),]#使exp与id完全对应

tmp=by(exp,ids$symbol,function(x)rownames(x)[which.max(rowMeans(x))])

probes=as.character(tmp)

dim(exp)

exp=exp[rownames(exp)%in%probes,]# 过滤有多个探针的基因

dim(exp)

rownames(exp)=ids[match(rownames(exp),ids$probe_id),2]

head(exp)

第四步：获取分组信息

pd <-pData(eSet[[1]])

library(stringr)

group_list=ifelse(str_detect(pd$title,"Control")==TRUE,"contorl","treat")group_list

str_detect(pd$title,"Control")

方法二：自己根据CSV里文件类型自己写代码

group_list=c(rep("control",times=3),rep("treat",times=3))

第五步：检查表达矩阵

方法一：管家基因

exp['GAPDH',]#检验用的管家基因

exp['ACTB',]#同上

boxplot(exp)

#使用ggplot2画各个样本表达量的boxplot图

# 准备画图所需数据exp_L

library(reshape2)

head(exp)

exp_L=melt(exp)

head(exp_L)

colnames(exp_L)=c('symbol','sample','value')

head(exp_L)

# 获得分组信息

library(stringr)

group_list=ifelse(str_detect(pd$title,"Control")==TRUE,"contorl","treat")

group_list

exp_L$group=rep(group_list,each=nrow(exp))

head(exp_L)

# ggplot2画图 library(ggplot2)

p=ggplot(exp_L,aes(x=sample,y=value,fill=group))+geom_boxplot()

print(p)

##boxplot图精修版p=ggplot(exp_L,aes(x=sample,y=value,fill=group))+geom_boxplot()

p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")

p=p+theme_set(theme_set(theme_bw(base_size=20)))

p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())

print(p)

#查案画出的箱型图，如果均值不在一条线上，说明样本有批次效应，需要人工校正

library(limma)

exp=normalizeBetweenArrays(exp)

#除了箱型图还可以画其他的图

p=ggplot(exp_L,aes(x=sample,y=value,fill=group))+geom_violin()

print(p)

p=ggplot(exp_L,aes(value,fill=group))+geom_histogram(bins=200)+facet_wrap(~sample,nrow=4)

print(p)

p=ggplot(exp_L,aes(value,col=group))+geom_density()+facet_wrap(~sample,nrow=4)

print(p)

p=ggplot(exp_L,aes(value,col=group))+geom_density()

print(p)

第四步：检查样本信息，检查样本分组信息，一般看PCA图，hclust图

# 更改表达矩阵列名

head(exp)

colnames(exp)=paste(group_list,1:6,sep='')

head(exp)

# 定义nodePar

nodePar<-list(lab.cex=0.6,pch=c(NA,19),cex=0.7,col="blue")

# 聚类

hc=hclust(dist(t(exp)))

par(mar=c(5,5,5,10))

# 绘图

plot(as.dendrogram(hc),nodePar=nodePar,horiz=TRUE)

#绘PCA图

library(ggfortify)

# 互换行和列，再dim一下df=as.data.frame(t(exp))

# 不要view df，列太多，软件会卡住；

dim(df)

dim(exp)

exp[1:6,1:6]

df[1:6,1:6]

df$group=group_list

autoplot(prcomp(df[,1:ncol(df)-1)]),data=df,colour='group')

#保存数据

save(exp,group_list,file="step2output.Rdata")

第五步：差异分析

#limma包需要以下三个数据，表达矩阵(exp)；分组矩阵(design)；差异比较矩阵(contrast.matrix)

#表达矩阵(exp)我们早就得到了，不用再制作了；我们也得到了存放分组信息的向量group_list，用它来制作我们的分组矩阵

rm(list=ls())## 魔幻操作，一键清空

options(stringsAsFactors=F)

load(file="step2output.Rdata")#上步得到的EXP分组信息

dim(exp)

library(limma)# 做分组矩阵 

design<-model.matrix(~0+factor(group_list))

colnames(design)=levels(factor(group_list))

rownames(design)=colnames(exp)

design#查看得到的分组矩阵

#输入数据—差异比较矩阵

contrast.matrix<-makeContrasts(paste0(c("treat","contorl"),collapse="-"),levels=design)

contrast.matrix

#limma包做差异分析，只有三个步骤：lmFit，eBayes，topTable

##step1

fit<-lmFit(exp,design)

##step2

fit2<-contrasts.fit(fit,contrast.matrix)##这一步很重要，大家可以自行看看效果fit2<-eBayes(fit2) ##default no trend!!!

##eBayes()withtrend=TRUE

##step3

tempOutput=topTable(fit2,coef=1,n=Inf)nrDEG=na.omit(tempOutput)

#可以保存数据 write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)

head(nrDEG)

save(nrDEG,file="DEGoutput.Rdata") #保存差异数据nrDEG

#检查差异数据的准确性，进行可视化

#热图

rm(list=ls())## 魔幻操作，一键清空

options(stringsAsFactors=F)

load(file="DEGoutput.Rdata")

load(file="DEGinput.Rdata")

library(pheatmap)

choose_gene=head(rownames(nrDEG),25)

choose_matrix=exp[choose_gene,]

choose_matrix=t(scale(t(choose_matrix)))

pheatmap(choose_matrix)

#火山图

rm(list=ls()) ## 魔幻操作，一键清空

options(stringsAsFactors=F)

load(file="DEGoutput.Rdata")

colnames(nrDEG)

plot(nrDEG$logFC,-log10(nrDEG$P.Value))

DEG=nrDEG

logFC_cutoff<-with(DEG,mean(abs(logFC))+2*sd(abs(logFC)))DEG$change=as.factor(ifelse(DEG$P.Value<0.05&abs(DEG$logFC)>logFC_cutoff,ifelse(DEG$logFC>logFC_cutoff,'UP','DOWN'),'NOT'))

head(DEG)

g=ggplot(data=DEG,aes(x=logFC,y=-log10(P.Value),color=change))+geom_point(alpha=0.4,size=1.75)+theme_set(theme_set(theme_bw(base_size=20)))+xlab("log2 fold change")+ylab("-log10 p-value")+ggtitle(this_tile)+theme(plot.title=element_text(size=15,hjust=0.5))+scale_colour_manual(values=c('blue','black','red'))

print(g)

第六步：富集分析

#KEGG富集分析

library(clusterProfiler)

gene_up=deg[deg$change=='up','ENTREZID']

gene_down=deg[deg$change=='down','ENTREZID']

gene_diff=c(gene_up,gene_down)

gene_all=deg[,'ENTREZID']

kk.up<-enrichKEGG(gene=gene_up,organism='hsa',universe=gene_all,pvalueCutoff=0.9,qvalueCutoff=0.9)

head(kk.up)[,1:6]

dim(kk.up)

kk.down<-enrichKEGG(gene=gene_down,organism='hsa',universe=gene_all,pvalueCutoff=0.9,qvalueCutoff=0.9)

head(kk.down)[,1:6]

dim(kk.down)

kk.diff<-enrichKEGG(gene=gene_diff,organism='hsa',pvalueCutoff=0.05)

head(kk.diff)[,1:6]

kegg_diff_dt<-kk.diff@result

down_kegg<-kk.down@result%>%

    filter(pvalue<0.05) %>%mutate(group=-1)

up_kegg<-kk.up@result%>%

    filter(pvalue<0.05) %>%mutate(group=1)

kegg_plot<-function(up_kegg,down_kegg)

{

dat=rbind(up_kegg,down_kegg)colnames(dat)

dat$pvalue=-log10(dat$pvalue)

dat$pvalue=dat$pvalue*dat$group

dat=dat[order(dat$pvalue,decreasing=F),]

g_kegg<-ggplot(dat,aes(x=reorder(Description,order(pvalue,decreasing=F)),y=pvalue,fill=group))+geom_bar(stat="identity")+scale_fill_gradient(low="blue",high="red",guide=FALSE)+scale_x_discrete(name="Pathway names")+scale_y_continuous(name="log10P-value")+coord_flip()+theme_bw()+theme(plot.title=element_text(hjust=0.5))+ggtitle("Pathway Enrichment")

}

g_kegg<-kegg_plot(up_kegg,down_kegg)

g_kegg

ggsave(g_kegg,filename='kegg_up_down.png')#绘图并保存

#GSEA是另一个常用的富集分析方法，目的是看看基因全局表达量的变化是否有某些特定的基因集合的倾向性。

data(geneList,package="DOSE")

head(geneList)

length(geneList)

names(geneList)

boxplot(geneList)

boxplot(deg$logFC)

geneList=deg$logFC

names(geneList)=deg$ENTREZID

geneList=sort(geneList,decreasing=T)

kk_gse<-gseKEGG(geneList=geneList,organism='hsa',nPerm=1000,minGSSize=120,pvalueCutoff=0.9,verbose=FALSE)

head(kk_gse)[,1:6]

gseaplot(kk_gse,geneSetID=rownames(kk_gse[1,]))

down_kegg<-kk_gse[kk_gse$pvalue<0.05&kk_gse$enrichmentScore<0,];down_kegg$group=-1

up_kegg<-kk_gse[kk_gse$pvalue<0.05&kk_gse$enrichmentScore>0,];up_kegg$group=1

gse_kegg=kegg_plot(up_kegg,down_kegg)

print(gse_kegg)

ggsave(gse_kegg,filename='kegg_up_down_gsea.png')

#GO富集分析

GO富集分析原理：有一个term注释了100个差异表达基因参与了哪个过程，注释完之后（模式生物都有现成的注释包，不用我们自己注释），计算相对于背景它是否显著集中在某条通路、某一个细胞学定位、某一种生物学功能。

对GO database analysis一般从三个层面进行：

Cellular component，CC  细胞成分

Biological process， BP  生物学过程

Molecular function，MF  分子功能

这三个层面具体是指：

Cellular component解释的是基因存在在哪里，在细胞质还是在细胞核？如果存在细胞质那在哪个细胞器上？如果是在线粒体中那是存在线粒体膜上还是在线粒体的基质当中？这些信息都叫Cellular component。

Biological process是在说明该基因参与了哪些生物学过程，比如，它参与了rRNA的加工或参与了DNA的复制，这些信息都叫Biological process

Molecular function在讲该基因在分子层面的功能是什么？它是催化什么反应的？

library(clusterProfiler)

gene_up=deg[deg$change=='up','ENTREZID']

gene_down=deg[deg$change=='down','ENTREZID']

gene_diff=c(gene_up,gene_down)

head(deg)

#细胞组分

ego_CC<-enrichGO(gene=gene_diff,OrgDb=org.Hs.eg.db,ont="CC",pAdjustMethod="BH",minGSSize=1,pvalueCutoff=0.01,qvalueCutoff=0.01,readable=TRUE)

#生物过程，重要

ego_BP<-enrichGO(gene=gene_diff,OrgDb=org.Hs.eg.db,ont="BP",pAdjustMethod="BH",minGSSize=1,pvalueCutoff=0.01,qvalueCutoff=0.01,readable=TRUE)

#分子功能，很重要

ego_MF<-enrichGO(gene=gene_diff,OrgDb=org.Hs.eg.db,ont="MF",pAdjustMethod="BH",minGSSize=1,pvalueCutoff=0.01,qvalueCutoff=0.01,readable=TRUE)

save(ego_CC,ego_BP,ego_MF,file="ego_GPL6244.Rdata")

rm(list=ls())

load(file="ego_GPL6244.Rdata")

#第一种，条带图，按p从小到大排的#如果运行了没出图，就dev.new()

barplot(ego_CC,showCategory=20,title="EnrichmentGO_CC")

barplot(ego_BP,showCategory=20,title="EnrichmentGO_CC")

#第二种，点图，按富集数从大到小的dotplot(ego_CC,title="EnrichmentGO_BP_dot")

#保存

pdf(file="dotplot_GPL6244.pdf")

dotplot(ego_CC,title="EnrichmentGO_BP_dot")

dev.off()