lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

clusterProfiler的纯无参自定义物种注释的GO、KEGG富集分析及GSEA

已有 3600 次阅读 2020-1-8 23:22 |个人分类:R统计分析|系统分类:科研笔记| _R语言, _基因功能富集分析

RclusterProfiler的纯无参自定义物种注释的GOKEGG富集分析及GSEA

 

首先感谢水生沈梦圆姐,帮忙解决了困扰已久的问题6.gif

问题就是怎样作无参(数据库中完全没有关于该物种或其近缘种的注释信息)时的GOKEGG富集分析。

Y叔的clusterProfiler包的enricher()函数,支持自定义输入KEGGGO基因注释集,提供基因名称列表即可实现该需求。

嗯嗯,很好用的一个方法。当然师姐也批评俺不好好看说明,Y叔写了那么详细的文档给介绍了该功能。此前也常使用clusterProfiler分析有参物种,还真不知道还可以适用无参......见谅。

之后顺便整理了一下这些过程,便于以后使用(尽管网上也早已有该方法的教程,为什么先前没搜到呢,可能是我找的方式不对......)。

#Bioconductor 安装 clusterProfiler
BiocManager::install('clusterProfiler')
 
#加载 clusterProfiler
library(clusterProfiler)

下文中使用的示例数据等,可在以下链接获取:

https://pan.baidu.com/s/132SP2HHWQzIiVQQo8eAVWg

 


差异表达基因的KEGG富集分析


先以KEGG来吧,因为当时有待解决的就是差异表达基因的KEGG富集的问题。

对于无参考数据库的物种而言,因为没有现成的基因功能注释信息(这是最让人难受的),首先需要使用基因的核酸序列/或蛋白序列作注释。KEGG的注释方法通常使用KAAS等。

然后可以获得如下示例的基因功能注释表,包含了基因名称及其所在的pathway对应关系。

1.png

这样该物种的KEGG功能注释集就有了。

对于差异表达基因的获得,常见方法如limmaedgeRDESeq2等。自定义阈值筛选后,即可获得一串差异表达基因的名称列表,例如下示例。

2.png

 

之后使用这两张表,作差异基因的KEGG富集就可以了。

将得到的基因-pathway对应关系表输入enricher()作为该物种的背景集,并提供用于作富集的基因名称,即可鉴定功能富集。

#“kegg_anno.txt”为预先获得的 KEGG 注释文件,包含基因和 KEGG pathway 的对应关系
#我用蛋白序列做的 kegg 注释,后续又对应的蛋白-基因关系,不用担心同一途径中基因名称出现重复的问题,enricher() 函数可自动去重
kegg_anno <- read.delim('kegg_anno.txt', colClasses = 'character')
head(kegg_anno)
 
#目标基因列表
gene_select <- read.delim('gene_select.txt', stringsAsFactors = FALSE)$gene_id
head(gene_select)
 
#KEGG 富集分析
#默认以所有注释到 KEGG 的基因为背景集,也可通过 universe 参数指定其中的一个子集作为背景集
#默认以 p<0.05 为标准,Benjamini 方法校正 p 值,q 值阈值 0.2
#默认输出 top500 富集结果
kegg_rich <- enricher(gene = gene_select,
    TERM2GENE = kegg_anno[c('pathway_id', 'gene_id')], 
    TERM2NAME = kegg_anno[c('pathway_id', 'pathway_description')], 
    pvalueCutoff = 0.05, 
    pAdjustMethod = 'BH', 
    qvalueCutoff = 0.2, 
    maxGSSize = 500)
 
#输出默认结果,即根据上述 p 值等阈值筛选后的
write.table(kegg_rich, 'kegg_rich.significant.txt', sep = '\t', row.names = FALSE, quote = FALSE)

结果中,包含了富集的pathway id和名称,在该pathway中存在差异表达的基因数量/所有差异表达的基因数量,该pathway中的背景基因数量/背景基因集总数量,富集的p值、p调整值及q值,在该pathway中存在差异表达的基因名称等信息。

3.png

 

Y叔也在clusterProfiler中提供了一系列的可视化方案用于展示本次富集分析结果,极大的便利。

#clusterProfiler 包里的一些默认作图方法,例如
barplot(kegg_rich)  #富集柱形图
dotplot(kegg_rich)  #富集气泡图
cnetplot(kegg_rich) #网络图展示富集功能和基因的包含关系
emapplot(kegg_rich) #网络图展示各富集功能之间共有基因关系
heatplot(kegg_rich) #热图展示富集功能和基因的包含关系

富集柱形图,纵轴为显著富集的pathway,横轴为该pathway中所含差异表达基因数量,颜色为校正后p值。

4.png

富集气泡图,纵轴为显著富集的pathway,横轴为该pathway中存在差异表达的基因数量/所有差异表达的基因数量的比值,点的大小为该pathway中所含差异表达基因数量,颜色为校正后p值。

5.png

网络图,展示富集pathway和基因的包含关系。代表pathway的节点大小为该pathway中所含差异表达基因数量。

6.png

网络图,展示各富集pathway之间共有基因关系。代表pathway的节点大小为该pathway中所含差异表达基因数量,颜色为校正后p值;节点之间的边表示pathway之间存在共有基因,边的粗细代表了共有基因数量的大小。

7.png

以及用于展示差异表达基因是否存在于特定pathway的热图。

8.png

还有其它更多效果图就不再举例了,可自行参阅clusterProfiler文档。

 

另外,以上结果均为根据p值等阈值筛选后的,其它特殊情况,有时我们可能想查看所有的结果(不管它是否显著)。

#如果想输出所有富集结果(不考虑 p 值阈值等)
#可以在上述 enricher() 函数中,将 p、q 等阈值设置为 1
 
#或者在 enrichResult 类对象中直接提取需要的结果
names(attributes(kegg_rich))  #查看对象“kegg_rich”包含元素
result <- kegg_rich@result    #其中的 result,即为所有基因的富集结果
 
write.table(result, 'kegg_rich.all.txt', sep = '\t', row.names = FALSE, quote = FALSE)

 


差异表达基因的GO富集分析


对于GO,方法同上。可首先注释基因的GO功能(GO注释方法有很多,例如Blast2GO等),然后将注释好的基因-GO对应关系表,以及差异表达基因名称列表输入至enricher()执行就可以了。

#预先获得的 GO 注释文件,基因和 GO 的对应关系
go_anno <- read.delim('go_anno.txt', header = FALSE, stringsAsFactors = FALSE)
names(go_anno) <- c('gene_id', 'ID')
 
#添加 GO 注释详情
#所有 GO 的注释描述均可在 GO 官网自定义获取:http://geneontology.org/
go_class <- read.delim('go_class.txt', header = FALSE, stringsAsFactors = FALSE)
names(go_class) <- c('ID', 'Description', 'Ontology')
go_anno <- merge(go_anno, go_class, by = 'ID', all.x = TRUE)
 
#目标基因列表
gene_select <- read.delim('gene_select.txt', stringsAsFactors = FALSE)$gene_id
 
#GO 富集分析
#默认以所有注释到 GO 的基因为背景集,也可通过 universe 参数输入背景集
#默认以 p<0.05 为标准,Benjamini 方法校正 p 值,q 值阈值 0.2
#默认输出 top500 富集结果
go_rich <- enricher(gene = gene_select,
    TERM2GENE = go_anno[c('ID', 'gene_id')], 
    TERM2NAME = go_anno[c('ID', 'Description')], 
    pvalueCutoff = 0.05, 
    pAdjustMethod = 'BH', 
    qvalueCutoff = 0.2, 
    maxGSSize = 500)
 
#输出默认结果,即根据上述 p 值等阈值筛选后的
write.table(go_rich, 'go_tmp.txt', sep = '\t', row.names = FALSE, quote = FALSE)
tmp <- read.delim('go_tmp.txt')
tmp <- merge(tmp, go_class[c('ID', 'Ontology')], by = 'ID')
tmp <- tmp[c(10, 1:9)]
tmp <- tmp[order(tmp$pvalue), ]
write.table(tmp, 'go_rich.significant.txt', sep = '\t', row.names = FALSE, quote = FALSE)
 
#同样地,如果想输出所有富集结果(不考虑 p 值阈值等),将 p、q 等值设置为 1 即可
#或者直接在 enrichResult 类对象中直接提取需要的结果
names(attributes(go_rich))  #查看对象“go_rich”包含元素
result <- go_rich@result    #其中的 result,即为所有基因的富集结果
 
result <- merge(result, go_class[c('ID', 'Ontology')], by = 'ID')
result <- result[c(10, 1:9)]
result <- result[order(result$pvalue), ]
write.table(result, 'go_rich.all.txt', sep = '\t', row.names = FALSE, quote = FALSE)
 
#一些常见作图方法
#上文提到的 clusterProfiler 包自带的柱形图、气泡图、网络图等默认方法略
#这些参考上文 KEGG
dotplot(go_rich)

结果格式同上,包含了富集的GO termid和名称,在该GO term中存在差异表达的基因数量/所有差异表达的基因数量,该GO term中的背景基因数量/背景基因集总数量,富集的p值、p调整值及q值,在该GO term中存在差异表达的基因名称等信息。

此外,额外使用GO官网中提供的分类关系,将富集的GO归类到3大顶层功能中(第一列,即生物学过程、细胞组分和分子功能类别)。

9.png

富集分析气泡图示例。

10.png

 

此外,对于GO的包含关系,即GO的有向无环图,在这种纯自定义富集分析的情况下,仍然可以使用topGO展示。

此时需要在作图前手动修改“enrichResult”类中的元素,主要为“result”统计表和“ontology”的GO分类,便于topGO识别。

#但需注意的是一次只能展示一类(BP、CC 或 MF)GO,如下以展示富集的“biological process”(即 BP)为例
#BiocManager::install('topGO')
#BiocManager::install('Rgraphviz')
library(topGO)
 
go_rich_BP <- go_rich
go_rich_BP@result <- go_rich_BP@result[as.vector(subset(tmp, Ontology == 'biological process')$ID), ]
go_rich_BP@ontology <- 'BP'
plotGOgraph(go_rich_BP)

图中颜色从黄色到红色,对应p值从大到小,即越来越显著。

11.png

 


基因集富集分析(GSEA


clusterProfiler还提供了自定义的基因集富集分析(GSEA)的功能,方便了无参时的GSEA

#基因表达的 log2FC
gene_express <- read.delim('gene_express.txt', stringsAsFactors = FALSE)
 
genelist <- gene_express$log2FoldChange
names(genelist) <- gene_express$gene_id
genelist <- genelist[order(genelist, decreasing = TRUE)]
 
#基因集富集分析,以 GO 富集为例
#因为示例数据集的 log2FC 是我瞎凑的,所以 p<0.05 绝对不会有结果的,为了演示数据就使用的 p<1 阈值,正常情况下视情况选择合适 p 值
#默认输出 top500 富集结果
gsea <- GSEA(genelist, 
    TERM2GENE = go_anno[c('ID', 'gene_id')], 
    TERM2NAME = go_anno[c('ID', 'Description')], 
    pvalueCutoff = 1, 
    pAdjustMethod = 'BH',
    maxGSSize = 500)
 
write.table(gsea, 'gsea.txt', sep = '\t', row.names = FALSE, quote = FALSE)
 
#查看 enrichResult 类对象中的元素
#names(attributes(gsea))
#head(gsea@result)
 
#同上文提到的,clusterProfiler 包里的一些默认作图方法,例如
dotplot(gsea)  #富集气泡图
emapplot(gsea) #网络图展示各富集功能之间共有基因关系
gseaplot(gsea, 'GO:0004386')    #基因集富集得分图

结果中,展示了基因集富集的GO term,基因集数量,富集分数、标准化富集分数(NES)、p值、校正后p值等信息。

12.png

如下展示了基因集富集得分图。

13.png

 



更多精彩,可关注个人公众号“生信小白鱼”,感谢大家支持。




http://blog.sciencenet.cn/blog-3406804-1213409.html

上一篇:R包ropls的PCA、PLS-DA和OPLS-DA
下一篇:R语言中的判别分析方法

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-8-6 03:25

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部