GEO数据挖掘（三）— 基因通路富集全套代码分享

杭州共生网络 2022-12-28 13:15 10172人围观数据

前几次教程“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')