GEO数据挖掘（三）— 基因通路富集全套代码分享
前几次教程“GEO数据挖掘（一）简单快速下载GEO数据”、“GEO数据挖掘（一）下载SRA库原始测序数据”、 “GEO数据挖掘（二）--基因差异表达分析及可视化全套代码分享”已经分享过下载数据、基因差异分析等步骤。本次将分享如何将提取出来的差异基因进行基因富集分析，并进行可视化展示。
通过前面的教程我们已经筛选出来了差异基因，但是基因的所属分类很多，那怎么才能知道这些基因所代表的生物学意义呢，这就需要利用基因富集分析来为这些基因注释功能，分析它究竟富集到哪个基因集上，参与了哪些调控，从而影响了疾病的发生。富集分析的原理简单来说就是分析一组基因在某个功能节点上是否相比于随机水平更显著，是由单个基因的简单注释扩展到多个基因集合的成组性分析。目前最常用的富集方法就是基于GO和KEGG的富集分析，来揭示一类基因所代表的生物学背景。
GO term功能富集分析
基因本体（gene ontology， GO）数据库是目前应用最广泛的基因通路注释体系之一。其原理简单理解就是计算参与调控过程的差异基因是否显著聚集在以下三个层面中：细胞成分（Cellular component，CC） 、分子功能（Molecular function，MF）、生物学过程 （Biological process，BP）。
细胞成分CC描述基因产物执行功能的具体细胞结构位置，例如某个产物蛋白可能定位在细胞核中或者核糖体中。
分子功能MF描述基因产物在分子水平上的活动，例如催化或运输。
生物学过程BP描述基因产物所关联的某个生物功能，或者多个分子活动完成的一个大的生物活动。例如有丝分裂或嘌呤代谢。
##导入数据
rm(list = ls())
load("step4_output.Rdata")
library(clusterProfiler)
library(dplyr)
library(ggplot2)
source("kegg_plot_function.R")
#(1)输入数据
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']
#(2)GO分析，分三部分
#以下步骤耗时很长，实际运行时注意把if后面的括号里F改成T
library(org.Hs.eg.db)
if(T){
#细胞组分
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_GSE42872.Rdata")
}
load(file = "ego_GSE42872.Rdata")
#(3)可视化
#条带图
barplot(ego_CC,showCategory=20)
#气泡图
dotplot(ego_CC)
geneList = deg$logFC
names(geneList)=deg$ENTREZID
geneList = sort(geneList,decreasing = T)
#(3)展示top5通路的共同基因。
#Gene-Concept Network
cnetplot(ego_CC, categorySize="pvalue", foldChange=geneList,colorEdge = TRUE)
cnetplot(ego_CC, foldChange=geneList, circular = TRUE, colorEdge = TRUE)
#Enrichment Map
emapplot(ego_CC)
#(4)展示通路关系
goplot(ego_CC)
#(5)Heatmap-like functional classification
heatplot(ego_CC,foldChange = geneList)
pdf("heatplot.pdf",width = 14,height = 5)
heatplot(ego_CC,foldChange = geneList)
dev.off()
KEGG 富集分析
KEGG pathway通路富集（Kyoto encyclopedia of genes and genomes, KEGG）是系统分析基因功能、基因组信息的数据库，整合了基因组学、生物化学及系统功能组学的信息，有助于我们把基因及表达信息作为一个整体进行研究。KEGG的每个通路图都包含一个分子相互作用和反应网络，会将基因组中的基因与通路中的基因产物联系起来，解释细胞和生物体的新陈代谢和各种其他功能的生物学背景。
#上调、下调、差异、所有基因
#（1）输入数据
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']
#（2）对上调/下调/所有差异基因进行富集分析
if(T){
kk.up <- enrichKEGG(gene = gene_up,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff = 0.9)
kk.down <- enrichKEGG(gene = gene_down,
organism = 'hsa',
universe = gene_all,
pvalueCutoff = 0.9,
qvalueCutoff =0.9)
kk.diff <- enrichKEGG(gene = gene_diff,
organism = 'hsa',
pvalueCutoff = 0.9)
save(kk.diff,kk.down,kk.up,file = "GSE21933kegg.Rdata")
}
load("GSE21933kegg.Rdata")
#(3)从富集结果中提取出结果数据框
kegg_diff_dt <- kk.diff@result
#(4)按照pvalue筛选通路
#在enrichkegg时没有设置pvaluecutoff，在此处筛选
down_kegg <- kk.down@result %>%
filter(pvalue<0.05) %>% #筛选行
mutate(group=-1) #新增列
up_kegg <- kk.up@result %>%
filter(pvalue<0.05) %>%
mutate(group=1)
#(5)可视化
g_kegg <- kegg_plot(up_kegg,down_kegg)
#g_kegg +scale_y_continuous(labels = c(20,15,10,5,0,5))
ggsave(g_kegg,filename = 'kegg_up_down.png')
目前我们已经走完“GEO数据挖掘”的基本流程，感谢大家的持续关注，但是目前只是简单分析教程，要想掌握GEO数据挖掘的精髓，还需要了解每一步分析背后所解决的生物学问题和运用算法，后续还会继续为大家分享相关知识。
请持续关注“GEO数据挖掘”系列文章，每周一个实用干货带您上手生信分析。
