GO和KEGG富集教程

Go和KEGG富集教程 前提 *** 作步骤注释中心库

前提

假设现在你已经在R官网上下载并安装好了R，并且已经有了自己的基因数据，例如一个excel表格中存放的数据。如下面这种形式。

现在需要做GO富集和KEGG的富集并生成想要的气泡或者通路图。

*** 作步骤 安装clusterProfiler

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

#这一步先安装管理包，第一次安装会让你选择镜像源，我选了中国香港

BiocManager::install("clusterProfiler")
#安装后问要不要更新旧包，我选了a,全部更新

安装绘图及注释库

BiocManager::install("org.Hs.eg.db") #人类注释数据库
#这一步又显示安装更新包，我选择了s，部分更新，会出来d窗让你选择是不是更新，我选择更新
BiocManager::install("GOplot") #绘图包
#又更新包，选择了s

载入包

library(clusterProfiler)
library(org.Hs.eg.db) #人类注释数据库,这里需要加载自己要用的数据库

查看和更改当前工作路径

getwd()
setwd("D:/R-4.0.4/workspace") #更改当前工作路径，这个路径下存放着数据excel
#windows用户注意路径要用“/”或者“\”，否则会报错不存在路径

读取数据

install.packages("readxl") #安装读取excel的包
library(readxl) #加载包
data = read_excel("D:\R-4.0.4\workspace\WZ78.vs.WZ92.fil.xls")

#如果报错 libxls error: Unable to open file，则打开xls文件另存一下，重新使用代码读取就行了。

gene_name = data$Gene_ID  #$后面跟列名，这里需要处理第一列，名称为Gene_ID

可以打印看一下基因的名称：

注意：本教程中基因名称已经做过ID的转换，所以后续不需要进行ID转换这一个步骤

GO的富集分析

ALL <- enrichGO(gene_name, "org.Hs.eg.db", keyType = "ENTREZID",ont = 'ALL',pvalueCutoff  = 0.05,pAdjustMethod = "BH",  qvalueCutoff  = 0.1, readable=T)  #一步到位
注意这里使用的db是人类
ALL <- enrichGO(gene_name,OrgDb = org.Hs.eg.db, pvalueCutoff =0.5)#阈值放宽

3种分开进行富集

BP <-enrichGO(gene_name, "org.Hs.eg.db", keyType = "ENTREZID",ont = "BP",pvalueCutoff  = 0.05,pAdjustMethod = "BH",qvalueCutoff  = 0.1, readable=T) 
MF <- enrichGO(gene_name, "org.Hs.eg.db", keyType = "ENTREZID",ont = "MF",pvalueCutoff  = 0.05,pAdjustMethod = "BH",qvalueCutoff  = 0.1, readable=T)
CC <- enrichGO(gene_name, "org.Hs.eg.db", keyType = "ENTREZID",ont = "CC",pvalueCutoff  = 0.05,pAdjustMethod = "BH",qvalueCutoff  = 0.1, readable=T)

富集结果绘图 气泡图

dotplot(ALL,x="count",showCategory = 14,colorBy="pvalue"，title="随便一个名字") #showCategory即展示几个分类，最好都展示，ALL是前面的全部分析结果

气泡示例效果图：

柱状图

barplot(ALL, showCategory=20,title="EnrichmentGO") 

kegg富集

kegg <- enrichKEGG(gene_name, 
                   organism = 'hsa',  
                   keyType = 'kegg', 
                   pvalueCutoff = 0.05,
                   pAdjustMethod = 'BH', 
                   minGSSize = 3,
                   maxGSSize = 500,
                   qvalueCutoff = 0.2,
                   use_internal_data = FALSE)
## hsa为人的简写，bta是牛的简写 ,gene_name是前面读取的文件数据

kegg富集结果绘图 显示通路图

browseKEGG(kegg,'mmu01100') 

通路图示意效果图：

画气泡图

dotplot(kegg,font.size=8，colorBy="pvalue") 

注释中心库

注释中心库中有很多注释可用，如果你要富集的物种并不在clusterProfiler包含的19个物种中，你可以安装注释中心，在其中查找是否有你需要的注释，再将其下载下来使用。

安装注释中心

BiocManager::install("AnnotationHub")
让更新包，我选择了s

查询自己的库在不在注释中心

library(AnnotationHub) #加载注释中心
ah <- AnnotationHub() #创建数据中心对象
#会让你选择是否创建一个数据中心缓存目录，选择yes

#如果出现报错reason: Failed to connect to bioconductor.org port 443: Timed out,Consider rerunning with 'localHub=TRUE'，则按照其提示修改

ah <- AnnotationHub(localHub=TRUE)

#如果报错Invalid Cache: index file, Missing entry in cache for: annotationhub.index.rds，则再进行修改
ah <- AnnotationHub(localHub=FALSE)

query(ah, "Vibrio cholerae")  #查询Vibrio cholerae是否在注释中心，ah是上一步创建的注释中心对象

如果有查询结果，例如显示一条记录[["AH10447"]]，则可以将其下载下来使用。
查询结果示例：

ath <- ah[["AH10447"]] #下载弧菌属,如果查询有结果最后会有个数字代号

注： 如果自己的物种在注释中心也查找不到，则可以自己制作自己的物种库，具体 *** 作方法可以参看https://blog.csdn.net/u012110870/article/details/102804471/?utm_medium=distribute.pc_relevant.none-task-blog-baidujs_title-4&spm=1001.2101.3001.4242