用python进行GO富集分析 go富集分析步骤

Go和KEGG富集教程

• 前提
• 操作步骤
• 注释中心库

前提

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

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

操作步骤

1. 安装clusterProfiler

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

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

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

2. 安装绘图及注释库

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

3. 载入包

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

4. 查看和更改当前工作路径

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

5. 读取数据

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转换这一个步骤

6. 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)

7. 富集结果绘图

• 气泡图

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

气泡示例效果图：

• 柱状图

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

8. 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是前面读取的文件数据

9. kegg富集结果绘图

• 显示通路图

browseKEGG(kegg,'mmu01100')

通路图示意效果图：

• 画气泡图

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

注释中心库

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

1. 安装注释中心

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

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

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"]] #下载弧菌属,如果查询有结果最后会有个数字代号

注： 如果自己的物种在注释中心也查找不到，则可以自己制作自己的物种库，具体操作方法可以参看