抱歉，您的浏览器无法访问本站
本页面需要浏览器支持（启用）JavaScript
了解详情 >

转录组数据分析笔记（10）——初识GO/KEGG富集分析

学习笔记转录组数据分析

发布于：2022年5月16日

字数：2.2k 字

时长：9 分钟
次浏览

前面介绍了如何找到差异基因，我们通过R包DESeq2获得了差异表达基因，在此基础上做了更为直观的火山图和差异表达基因热图。但是仅仅知道差异表达基因的名字还不够，我们还要知道它到底有哪些功能和特征，就比如我看到一个很养眼的动漫角色，我就要去查查出自哪部番，是怎么样的人设和背景故事，一样的道理。这里简单记录下如何使用AnnotationHub，以及怎么进行GO\KEGG富集分析。

一个基因没有注释信息，那就只是一段核苷酸序列，有了注释信息我们才能知道这个基因在染色体上的定位，在具体的某个代谢途径上发挥什么功能等等。网上能找到很多注释信息的数据库，比如模式生物拟南芥TAIR，人类基因组hg19等等，Bioconductor有一个专门用来搜集注释信息数据库的工具包——AnnotationHub。

1. AnnotationHub注释数据库搜索工具

用bioconductor下载Annotationhub包，载入（注：为演示结果，以下命令均在Rstudio终端输入）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
> library("Annotationhub")
> hub <- AnnotationHub()
C:\Users\HUAWEI\AppData\Local/R/cache/R/AnnotationHub
  does not exist, create directory? (yes/no): yes
  |===================================================| 100%

snapshotDate(): 2021-10-20
> hub
AnnotationHub with 62386 records
# snapshotDate(): 2021-10-20
# $dataprovider: Ensembl, BroadInstitute, UCSC, ftp://ftp...
# $species: Homo sapiens, Mus musculus, Drosophila melano...
# $rdataclass: GRanges, TwoBitFile, BigWigFile, EnsDb, Rl...
# additional mcols(): taxonomyid, genome,
#   description, coordinate_1_based, maintainer,
#   rdatadateadded, preparerclass, tags, rdatapath,
#   sourceurl, sourcetype 
# retrieve records with, e.g., 'object[["AH5012"]]'

第一次使用AnnotationHub需要创建一个AnnotationHub对象。为了更直观地使用，我们将AnnotationHub对象赋值给hub变量。查看这个变量，我们可以得到如下的信息。

  • 数据库版本是2021-10-20，目前有62386条记录
  • 可以用$dataprovider 方式查看数据来源，比如数据来自于Ensembl，UCSC等等
  • 可以用$species 方式查看数据库有哪些物种，比如人类、小鼠等等
  • 可以用$rdataclass 方式查看数据类型
  • 可以通过函数mcols()查看更多信息
  • 获取数据的方式是object[["AH5012"]] object指你命名的变量名

以上就是AnnotationHub的标准用法，比如我想获得拟南芥的注释数据库，我就输入以下命令查找：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
> query(hub, "Arabidopsis thaliana")
AnnotationHub with 13 records
# snapshotDate(): 2021-10-20
# $dataprovider: UCSC, PathBank, NCBI,DBCLS, FANTOM5,DLRP...
# $species: Arabidopsis thaliana
# $rdataclass: SQLiteFile, TxDb, Tibble, list, OrgDb, Inp...
# additional mcols(): taxonomyid, genome,
#   description, coordinate_1_based, maintainer,
#   rdatadateadded, preparerclass, tags, rdatapath,
#   sourceurl, sourcetype 
# retrieve records with, e.g., 'object[["AH10456"]]' 

            title                                           
  AH10456 | hom.Arabidopsis_thaliana.inp8.sqlite            
  AH52245 | TxDb.Athaliana.BioMart.plantsmart22.sqlite      
  AH52246 | TxDb.Athaliana.BioMart.plantsmart25.sqlite      
  AH52247 | TxDb.Athaliana.BioMart.plantsmart28.sqlite      
  AH87070 | pathbank_Arabidopsis_thaliana_metabolites.rda   
  ...       ...                                             
  AH91794 | wikipathways_Arabidopsis_thaliana_metabolites...
  AH95585 | Alternative Splicing Annotation for Arabidops...
  AH95951 | org.At.tair.db.sqlite                           
  AH97723 | LRBaseDb for Arabidopsis thaliana (Thale cres...
  AH97844 | MeSHDb for Arabidopsis thaliana (Thale cress,...

query()函数查找，输入拟南芥的学名Arabidopsis thaliana我们可以看到一共找出了13个数据库。可以看到AH95951这个编号的数据库来源就是最大的拟南芥数据库TAIR（OrgDb，存储不同数据库基因ID之间对应关系，以及基因与GO等注释的对应关系，后面ID转换和GO分析要用到），我们就用这个数据库的注释资源。

1
2
3
4
> hub[["AH95951"]]
downloading 1 resources
retrieving 1 resource
  |                                                   |   0%

前面说过object[["AH5012"]]是获取数据的方式。以上，就可以挂后台自动下载了。我这里因为网速的原因不下了，从前面的基因名也能看出来，我筛选的差异基因都是AT开头的，而且之前的基因组注释文件也是TAIR下载的，我可以直接用bioconductor安装org.At.tair.db包，这里用AnnotationHub只是提供一个找注释数据库的思路。

2. GO/KEGG富集分析

2.1 基因ID转换

找到和下载注释数据库只是第一步，接下来GO/KEGG富集分析需要用到R包clusterProfiler和org.At.tair.db

先来看一下我们基因名是什么格式的：

很明显，我们的基因ID是TAIR类型（废话，我从TAIR下的），org.At.tair.db包可以转换基因ID类型

可以用keytypes(org.At.tair.db)或者columns(org.At.tair.db)查看可以转换的基因ID类型

转换基因ID代码如下：

1
2
3
4
5
6
7
8
library("org.At.tair.db")
columns(org.At.tair.db)     # 查看能转换基因的ID类型
diffgen <- nDEGs[, 1]       # 注意只需要基因名
diff_gen <- bitr(diffgen,
                 fromType = "TAIR",
                 toType = "ENTREZID",   # 基因ID类型TAIR转换为ENTREZID
                 OrgDb = "org.At.tair.db")      # 该函数是基于org.At.tair.db包的
diff_gen

这一步我的基因ID转换率只有60%左右，有将近一半的TAIR基因ID不能成功转换成ENTREZID，可能是Gene ID的版本问题，同一个基因在不同版本genecode中结果不一样，下载的注释文件原始版本我这里找不到了…暂时无法解决这个问题。只能不转换基因ID先跑一遍GO/KEGG富集分析。

看了很多教程都说clusterProfiler需要的ID类型是ENTREZID，这里我持怀疑态度，不转换后续也能得到结果，我看了函数enrichGO()默认的基因ID是ENTREZID并不代表不能改变，有可能是误传。查阅了一些资料，简单来说Entrez ID是来自于NCBI旗下Entrez gene数据库的编号系统，基因编号系统之间是可以相互转换的，这些ID可以在对应的数据库找到基因注释信息，就是说也可以在网页上手动注释。

2.2 GO/KEGG分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
library("clusterProfiler")
# GO富集分析
enrich_GO <- enrichGO(gene = diffgen,   # 基因名列表
                      OrgDb = 'org.At.tair.db', # 输入OrgDb数据库（注释对象信息）
                      keyType = 'TAIR',     # 输入的基因名ID类型
                      ont = 'ALL',      # 输出的GO分类
                      pAdjustMethod = 'fdr',
                      pvalueCutoff = 0.05,
                      qvalueCutoff = 0.2,
                      readable = TRUE)
GO_result <- enrich_GO@result
write.table(GO_result, 'GO_result.csv', sep = ',', quote = FALSE, row.names = FALSE)

# KEGG富集分析
enrich_KEGG <- enrichKEGG(gene = diffgen,
                        keyType = "kegg",
                        organism = "ath",   # 输入的物种名
                        pvalueCutoff = 0.05,
                        qvalueCutoff = 0.2)
KEGG_result <- enrich_KEGG@result
write.table(KEGG_result, 'KEGG_result.csv', sep = ',', quote = FALSE, row.names = FALSE)

clusterProfiler这个包进行GO和KEGG富集分析就这两个函数

这里我的GO只富集到两条细胞组分的内容：

说一下各列代表的意思：

  • ONTOLOGY GO分类BP（生物学过程）、CC（细胞组分）或MF（分子功能）
  • ID 富集到的GO id号
  • Description 富集到的GO描述
  • GeneRatio和BgRatio 分别为富集到该GO条目中的基因数目/给定基因的总数目，以及该条目中背景基因总数目/该物种所有已知的GO功能基因数目
  • pvalue、p.adjust和qvalue p值、校正后p值和q值信息
  • geneID和Count，富集到该GO条目中的基因名称和数目

KEGG富集分析结果表如下：

  • ID和Description 分别代表富集到KEGG的ID和描述，其他和GO富集都类似

KEGG富集分析的时候有一点需要注意，输入的organism名称需要在官网的KEGG Organisms列表中能找到，否则是不能进行分析的！点击这里进入KEGG Organisms: Complete Genomes

还发现一个很奇怪的问题，我在官网的Organisms列表能找到拟南芥Arabidopsis thaliana，但是在上面的函数中对参数赋值organism = "Arabidopsis thaliana"会显示HTTP 400错误，也就是发出的url请求有问题，但是输入organism = "ath"程序可以正常运行，以后注意写缩写吧（应该是只有缩写才行，会通过联网自动获取该物种的pathway注释信息）。

3. 可视化

clusterProfiler包还提供了GO/KEGG富集结果的可视化方案，此处代码参考CSDN，作者：Tian問

这里因为我的GO结果不好，只简单写一下流程，详细作图函数参数使用方法和效果同样可以参考上面的链接~

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
## GO富集分析可视化
#barplot
barplot(enrich_GO, showCategory = 10)
#dotplot
dotplot(enrich_GO, showCategory = 10)
#DAG有向无环图
plotGOgraph(enrich_GO)  #矩形代表富集到的top10个GO terms, 颜色从黄色过滤到红色，对应p值从大到小。
#igraph布局的DAG
goplot(enrich_GO)
#GO terms关系网络图（通过差异基因关联）
emapplot(enrich_GO, showCategory = 30)
#GO term与差异基因关系网络图
cnetplot(enrich_GO, showCategory = 5)

## KEGG富集分析可视化
#barplot
barplot(enrich_KEGG, showCategory = 10)
#dotplot
dotplot(enrich_KEGG, showCategory = 10)
#pathway关系网络图（通过差异基因关联）
emapplot(enrich_KEGG,  showCategory = 30)
#pathway与差异基因关系网络图
cnetplot(enrich_KEGG, showCategory = 5)
#pathway映射
browseKEGG(enrich_KEGG, "ath03060") #在pathway通路图上标记富集到的基因，会弹出页面链接到KEGG官网

关于GO/KEGG富集分析，还有非常多的操作和应用，我只是简单做个最基础的富集分析的学习，没有涉及到手动注释，构建orgdb等等更多操作。电脑快没电了，这篇笔记先暂时记这些，以后需要用到再补充~

本文由Phantom创作，转载请注明出处

博客内容遵循 署名-非商业性使用-相同方式共享 4.0 国际 (CC BY-NC-SA 4.0) 协议

本文永久链接是：http://www.shelven.com/2022/05/16/a.html

更新于：2022年12月4日

AnnotationHub

GO/KEGG

org.At.tair.db

欢迎小伙伴们留言评论~