您现在的位置是:首页 > 技术教程 正文

【R语言】——基因GO/KEGG富集分析!超级简单的保姆级教程!

admin 阅读: 2024-03-18
后台-插件-广告管理-内容页头部广告(手机)

上期“干货预警——原来基因功能富集分析这么简单!”和“【R语言】——基因GO/KEGG功能富集结果可视化(保姆级教程)”介绍如何使用DAVID在线分析工具对基因进行GO/KEGG功能富集分析和使用R ggplot包对获得的基因GO/KEGG功能富集结果进行可视化。本期介绍使用R clusterProfiler包和R AnnotationHub包对基因进行GO/KEGG功能富集分析、OrgDb包制作以及结果可视化。

GO/KEGG功能富集分析中重要的是背景基因的选择,使用R clusterProfiler包对基因进行富集,需要导入目的基因(前景基因)相对应物种的参考基因组(背景基因),现阶段“bioconductor”已有十几种常见动物,如人类、小鼠等物种的OrgDb。但仍然有许多物种不在Bioconductor的OrgDb列表里,但存在参考基因组,如山羊,绵羊等,这种情况则需要用到R AnnotationHub包进行索引其对应物种的参考基因组,并制作OrgDb包使用。

1 数据准备

数据输入格式(xlsx格式):

2 R包加载、数据导入及处理

  1. #下载包#
  2. if(!requireNamespace("BiocManager", quietly = TRUE))
  3. install.packages("BiocManager")
  4. BiocManager::install("clusterProfiler")
  5. BiocManager::install("topGO")
  6. BiocManager::install("Rgraphviz")
  7. BiocManager::install("pathview")
  8. install.packages("ggplot2")
  9. BiocManager::install('stringr')
  10. install.packages("openxlsx")
  11. #加载包#
  12. library(clusterProfiler)
  13. library(topGO)
  14. library(Rgraphviz)
  15. library(pathview)
  16. library(ggplot2)
  17. library(stringr)
  18. library(openxlsx)
  19. #导入数据#
  20. remove(list = ls()) #清除 Global Environment
  21. getwd() #查看当前工作路径
  22. setwd("C:/Rdata/jc") #设置需要的工作路径
  23. list.files() #查看当前工作目录下的文件
  24. data = read.xlsx("enrich-gene.xlsx",sheet= "enrich_genes",sep=',') #导入数据
  25. head(data)
  1. #数据处理-差异基因筛选#
  2. vector = abs(data$log2FC) > 1 & data$PValue < 0.05 & data$gene_name !="" ##abs绝对值;通常logFC> 1和PValue< 0.05条件进行筛选;data$gene_name != ""表示gene_name不为空白
  3. #data$gene_name<-str_to_title(data$gene_name)#用stringr将基因名称的第一个字母大写(小鼠首字母为大写)
  4. data_sgni= data[vector,]#筛选差异基因
  5. head(data_sgni)
  6. #All_gene <- rownames(data) # 提取所有基因基因名

3 背景基因选择及GO/KEGG富集分析

3.1 在“bioconductor”中已有OrgDb的物种的富集分析

在http://bioconductor.org/packages/release/BiocViews.html#___OrgDb 中寻找需要物种的OrgDb(目前仅有下图所示物种),以人类“org.Hs.eg.db”为例:

图1 “bioconductor”中已有OrgDb的物种

  1. #已有OrgDb的常见物种#
  2. BiocManager::install("org.Hs.eg.db")
  3. library(org.Hs.eg.db)
  4. #基因ID转换#
  5. keytypes(org.Hs.eg.db) #查看所有可转化类型
  6. entrezid_all = mapIds(x = org.Hs.eg.db, #id转换的比对基因组(背景基因)的物种,以人为例
  7. keys = data_sgni$gene_name, #将输入的gene_name列进行数据转换
  8. keytype = "SYMBOL", #输入数据的类型
  9. column = "ENTREZID")#输出数据的类型
  10. entrezid_all = na.omit(entrezid_all) #na省略entrezid_all中不是一一对应的数据情况
  11. entrezid_all = data.frame(entrezid_all) #将entrezid_all变成数据框格式
  12. head(entrezid_all)
  1. ###GO富集分析###
  2. GO_enrich = enrichGO(gene = entrezid_all[,1], #表示前景基因,即待富集的基因列表;[,1]表示对entrezid_all数据集的第1列进行处理
  3. OrgDb = org.Hs.eg.db,
  4. keyType = "ENTREZID", #输入数据的类型
  5. ont = "ALL", #可以输入CC/MF/BP/ALL
  6. #universe = 背景数据集 # 表示背景基因,无参的物种选择组装出来的全部unigenes作为背景基因;有参背景基因则不需要。
  7. pvalueCutoff = 1,qvalueCutoff = 1, #表示筛选的阈值,阈值设置太严格可导致筛选不到基因。可指定 1 以输出全部
  8. readable = T) #是否将基因ID映射到基因名称。
  9. GO_enrich = data.frame(GO_enrich) #将GO_enrich导成数据框格式
  10. #数据导出#
  11. write.csv(GO_enrich,'C:/Rdata/保存文件/GO_enrich.csv')
  1. ###KEGG富集分析###
  2. KEGG_enrich = enrichKEGG(gene = entrezid_all[,1], #即待富集的基因列表
  3. keyType = "kegg",
  4. pAdjustMethod = 'fdr', #指定p值校正方法
  5. organism= "human", #hsa,可根据你自己要研究的物种更改,可在https://www.kegg.jp/brite/br08611中寻找
  6. qvalueCutoff = 1, #指定 p 值阈值(可指定 1 以输出全部)
  7. pvalueCutoff=1) #指定 q 值阈值(可指定 1 以输出全部)
  8. KEGG_enrich = data.frame(KEGG_enrich)
  9. write.csv(KEGG_enrich,'C:/Rdata/保存文件/KEGG_enrich.csv') #数据导出

3.2 使用“AnnotationHub”获取在线注释并创建OrgDb对象

如果你所研究的物种不在Bioconductor的OrgDb列表里,但存在参考基因组,如山羊(Capra hircus/goat/chx),绵羊(sheet/Ovis aries)等,这种情况则需要用到AnnotationHub函数进行索引其对应物种的参考基因组(背景基因),并制作OrgDb包使用。

注意:AnnotationHub包连接的Bioconductor数据库是实时更新的,所以需要用到的时候再在线查询和使用。

  1. ###制作可索引到物种的OrgDb包###
  2. #下载和加载包#
  3. BiocManager::install("AnnotationHub")
  4. BiocManager::install("AnnotationDbi")
  5. BiocManager::install("rtracklayer")
  6. library(AnnotationHub)
  7. library(AnnotationDbi)
  8. library(rtracklayer)
  9. #索引与制作OrgDb#
  10. hub <- AnnotationHub() #建立AnnotationHub对象保存到hub
  11. query(hub, 'Capra hircus') #查询包含山羊(Capra hircus)的物种信息;结果有物种的各类信息需要进一步筛选
  12. query(hub[hub$rdataclass == "OrgDb"] , "Capra hircus") #筛选我们需要OrgDb类型;也可将上一步与这一步合并成query(hub,'org.Capra hircus')进行搜索
  13. goat <- hub[['AH101444']] #制作Capra hircus的OrgDb库;AH101444是Capra hircus对应的编号。
  14. goat #查看goat
  15. #help('select')
  1. #保存、载入与查看-AnnotationDbi#
  2. saveDb(goat,file="goat.OrgDb") #把goat对象保存成goat.OrgDb文件
  3. goat = loadDb(file="goat.OrgDb") #载入goat.OrgDb文件,保存到goat
  4. length(keys(goat)) #查看包含的基因数量
  5. columns(goat) #查看goat的数据类型
  6. keys(goat, keytype = "SYMBOL") #查看SYMBOL数据集下的ID
  1. # 查看AnnotationHub内容——根据自己兴趣了解#
  2. #display(hub) #调出网页
  3. #unique(hub$species) #查看hub里包含的所有物种
  4. #unique(hub$rdataclass) #查看hub里的数据类型
  5. #hub[hub$rdataclass == "OrgDb"] #查看hub里OrgDb类型的数据
  1. #基因ID转换#
  2. keytypes(goat) #查看所有可转化类型
  3. entrezid_all = mapIds(x = goat, #id转换的比对基因组(背景基因)所属物种,这边为山羊
  4. keys = data_sgni$gene_name, #将输入的gene_name列进行数据转换
  5. keytype = "SYMBOL", #输入数据的类型
  6. column = "ENTREZID")#输出数据的类型
  7. entrezid_all = na.omit(entrezid_all) #na省略entrezid_all中不是一一对应的数据情况
  8. entrezid_all = data.frame(entrezid_all) #将entrezid_all变成数据框格式
  9. head(entrezid_all)
  10. #GO富集分析#
  11. GO_enrich = enrichGO(gene = entrezid_all[,1], #待富集的基因列表
  12. OrgDb = goat, #指定物种的基因数据库,goat直接赋值给OrgDb参数即可
  13. keyType = 'ENTREZID', #输入数据的类型
  14. ont = 'ALL', #可指定 BP\MF\CC\ALL
  15. pAdjustMethod = 'fdr', #指定 p 值校正方法
  16. pvalueCutoff = 1, #指定 p 值阈值(指定 1 以输出全部)
  17. qvalueCutoff = 1, #指定 q 值阈值(指定 1 以输出全部)
  18. readable = FALSE)
  19. GO_enrich = data.frame(GO_enrich)
  20. write.csv(GO_enrich,'C:/Rdata/保存文件/GO_enrich.csv') #数据导出#
  21. ###KEGG富集分析###
  22. KEGG_enrich = enrichKEGG(gene = entrezid_all[,1], #即待富集的基因列表
  23. keyType = "kegg",
  24. pAdjustMethod = 'fdr', #指定p值校正方法
  25. organism= "chx", #山羊,可根据你自己要研究的物种更改,可在https://www.kegg.jp/brite/br08611中寻找
  26. qvalueCutoff = 1, #指定 p 值阈值(可指定 1 以输出全部)
  27. pvalueCutoff=1) #指定 q 值阈值(可指定 1 以输出全部)
  28. KEGG_enrich = data.frame(KEGG_enrich)
  29. write.csv(KEGG_enrich,'C:/Rdata/保存文件/KEGG_enrich.csv') #数据导出

输出的GO/KEGG富集结果各列内容:

ONTOLOGY:GO的BP(生物学过程)、CC(细胞组分)或MF(分子功能)三个方面内容;

ID:富集到的GO term/KEGG term;

Description:对GO term/KEGG term的生物学功能和意义进行描述;

GeneRatio:富集到该GO term/KEGG term中的基因数目/给定基因的总数目;

BgRatio:该GO term/KEGG term中背景基因总数目/该物种所有已知GO功能基因的数目;

pvalue、p.adjust和qvalue:p值、校正后p值和q值;

geneID和Count:富集到该GO term/KEGG term中的基因名称和数目。

4 GO/KEGG富集结果可视化

  1. ###GO/KEGG富集结果可视化###
  2. #数据载入与处理#
  3. install.packages("ggplot2")
  4. library(ggplot2)
  5. go_enrich = read.xlsx("enrich-gene.xlsx",sheet= "ONTOLOGY",sep=',')
  6. go_enrich$term <- paste(go_enrich$ID, go_enrich$Description, sep = ': ') #将ID与Description合并成新的一列
  7. go_enrich$term <- factor(go_enrich$term, levels = go_enrich$term,ordered = T)
  1. #纵向柱状图#
  2. ggplot(go_enrich,
  3. aes(x=term,y=Count, fill=ONTOLOGY)) + #x、y轴定义;根据ONTOLOGY填充颜色
  4. geom_bar(stat="identity", width=0.8) + #柱状图宽度
  5. scale_fill_manual(values = c("#6666FF", "#33CC33", "#FF6666") ) + #柱状图填充颜色
  6. facet_grid(ONTOLOGY~., scale = 'free_y', space = 'free_y')+
  7. coord_flip() + #让柱状图变为纵向
  8. xlab("GO term") + #x轴标签
  9. ylab("Gene_Number") + #y轴标签
  10. labs(title = "GO Terms Enrich")+ #设置标题
  11. theme_bw()
  12. #help(theme) #查阅这个函数其他具体格式
  1. #横向柱状图#
  2. ggplot(go_enrich,
  3. aes(x=term,y=Count, fill=ONTOLOGY)) + #x、y轴定义;根据ONTOLOGY填充颜色
  4. geom_bar(stat="identity", width=0.8) + #柱状图宽度
  5. scale_fill_manual(values = c("#6666FF", "#33CC33", "#FF6666") ) + #柱状图填充颜色
  6. facet_grid(.~ONTOLOGY, scale = 'free_x', space = 'free_x')+
  7. xlab("GO term") + #x轴标签
  8. ylab("Gene_Number") + #y轴标签
  9. labs(title = "GO Terms Enrich")+ #设置标题
  10. theme_bw() +
  11. theme(axis.text.x=element_text(family="sans",face = "bold", color="gray50",angle = 70,vjust = 1, hjust = 1 )) #对字体样式、颜色、还有横坐标角度()
  1. #气泡图#
  2. ggplot(go_enrich,
  3. aes(y=term,x=Count))+
  4. geom_point(aes(size=Count,color=p.adjust))+
  5. facet_grid(ONTOLOGY~., scale = 'free_y', space = 'free_y')+
  6. scale_color_gradient(low = "red",high ="blue")+
  7. labs(color=expression(PValue,size="Count"),
  8. x="Gene Ratio",y="GO term",title="GO Enrichment")+
  9. theme_bw()

图1 为GO富集结果图

KEGG富集结果与GO富集结果可视化类似可参考上一期“【R语言】——基因GO/KEGG功能富集结果可视化(保姆级教程)”内容。

好了本次分享就到这里,下期有更精彩内容,敬请期待。

关注“在打豆豆的小潘学长”公众号,发送“富集分析2”获得完整代码包和演示数据。

标签:
声明

1.本站遵循行业规范,任何转载的稿件都会明确标注作者和来源;2.本站的原创文章,请转载时务必注明文章作者和来源,不尊重原创的行为我们将追究责任;3.作者投稿可能会经我们编辑修改或补充。

在线投稿:投稿 站长QQ:1888636

后台-插件-广告管理-内容页尾部广告(手机)
关注我们

扫一扫关注我们,了解最新精彩内容

搜索
排行榜