百度360必应搜狗淘宝本站头条
当前位置:网站首页 > 热门文章 > 正文

GO、GSEA富集分析一网打进

bigegpt 2024-09-08 11:36 4 浏览

富集分析是生物信息分析中快速了解目标基因或目标区域功能倾向性的最重要方法之一。其中代表性的计算方式有两种:

一是基于筛选的差异基因,采用超几何检验判断上调或下调基因在哪些GO或KEGG或其它定义的通路富集。假设背景基因数目为m,背景基因中某一通路pathway中注释的基因有n个;上调基因有k个,上调基因中落于通路pathway的数目为l。简单来讲就是比较l/k是否显著高于n/m,即上调基因中落在通路pathway的比例是否高于背景基因在这一通路的比例。(实际计算时,是算的odds ratio的差异,l/(k-l) vs (n-l)/(m-k-n+l))。这就是常说的GO富集分析或KEGG富集分析,可以做的工具很多,GOEAST是其中一个最好用的在线功能富集分析工具,数据库更新实时,操作简单,并且可以直接用之前介绍的方法绘制DotPlot。

另一种方式是不硬筛选差异基因,而是对其根据表达量或与表型的相关度排序,然后判断对应的基因集是否倾向于落在有序列表的顶部或底部,从而判断基因集合对表型差异的影响和筛选有影响的基因子集。这叫GSEA富集分析,注释信息可以是GO,KEGG,也可以是其它任何符合格式的信息。GSEA富集分析 - 界面操作详细讲述了GSEA分析的原理、可视化操作和结果解读。

微信公众号biobabble的主创Y叔写有9个Bioconductor包,其中两个DOSEclusterProfiler囊括了前面提到的两种富集分析方法。并且不只支持GO、KEGG数据库,还支持Disease Ontology、MsEH enrichment analysis、Reactome通路分析等。具体可见其公众号,或软件的文档页 https://guangchuangyu.github.io/clusterProfiler/。

这么强大的工具,学习起来的路子却不是一帆风顺,最开始的拦路虎是软件的安装,系统较老配合上软件包更新较快(作者勤快也是问题),导致经常安装的是旧版本,用起来会遇到不少坑。直到有了conda,安装再也不是问题。解决了动态库依赖后,可以在Github安装最新的开发版本。

另外一个是文档较少,在R终端,直接使用help命令查看到的使用提示信息较少,寥寥几句,看过总觉得不踏实。~在线文档页内容少、更新慢。这是最开始学习时遇到的问题,这次秉着负责的精神,又重新读了文档页,发觉不需要再写一遍了,内容挺全的,主要是这一页http://guangchuangyu.github.io/2016/01/go-analysis-using-clusterprofiler/),但有几个地方需要更新下。自己对着文档页核对了下之前写的程序,再补充几点。

GO富集分析

首先还是列一个完整的例子。输入最好是用ENTREZ ID,值比较固定,不建议使用GeneSymbol,容易匹配出问题。

entrezID_text <- "4312
8318
10874
55143
55388
991
6280
2305
9493
1062
3868
4605
9833
9133
6279
10403
8685
597
7153
23397"
entrezID <- read.table(text=entrezID_text, header=F)
head(entrezID)
V1            
1    4312            
2    8318            
3    10874            
4    55143            
5    55388            
6    991

转换为向量

entrezID <- entrezID$V1
head(entrezID)
[1] 4312 8318 10874 55143 55388 991
# 这里的ENTREZ ID是从clusterProfiler里面提取是,是人的,
# 所以用了人的注释库, org.Hs.eg.db
library(org.Hs.eg.db)

开始富集分析

# readable=T: 原文档无这个参数,使用的是setReadble函数
MF <- enrichGO(entrezID, "org.Hs.eg.db", ont = "MF", keytype = "ENTREZID",
          pAdjustMethod = "BH", pvalueCutoff  = 0.05,
          qvalueCutoff  = 0.1, readable=T)
head(summary(MF))

去除冗余度高的条目

(参考http://guangchuangyu.github.io/2015/10/use-simplify-to-remove-redundancy-of-enriched-go-terms/)

result <- simplify(MF, cutoff=0.7, by="p.adjust", select_fun=min)
# 去除前
dim(MF)
# [1] 367 9
# 去除后
dim(result)
[1] 142 9

绘制泡泡图 (DotPlot Dotplot2)

dotplot(result, showCategory = 10)

绘制网络图 (边的宽度代表两个富集的Term共有的基因数目,点大小代表条目内基因数目多少,颜色代表P-value,值越小越红;如果想改变网络布局,参考igraph文档)

enrichMap(result, vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)

另外一种网络图由cnetplot函数获得,可以映射基因的表达量。

# geneList为一个vector,每个元素的名字为基因名,值为FoldChange,用于可视化点。
# > str(geneList)
# Named num [1:12495] 4.57 4.51 4.42 4.14 3.88 ... - attr(*, "names")= chr [1:12495] "4312" "8318" "10874" "55143" ...
# 这个geneList怎么获得的,会在后面的GSEA分析时提到
cnetplot(result, foldChange=geneList)

自己尝试了下,展示的有些乱,需要调整字体和显示的条目多少。故盗图展示如下便于解释,基因与其被注释的条目连线,点的颜色代表表达变化,圈的大小代表对应注释内基因数目多少。

如果想自己调整图的布局,还是建议把输出结果转换为Cytoscape(点击查看视频教程)可以识别的两列表格形式(如下),再赋值不同的属性就可以了。

awk 'BEGIN{OFS=FS="\t"}{if(FNR==1) print "Gene\tTerm"; else {split($8,geneL,"/"); for(i in geneL) print geneL[i],$2}}' MF | head
Gene    Term
PRDX6    cell adhesion molecule binding
MFGE8    cell adhesion molecule binding
FSCN1    cell adhesion molecule binding
ATXN2L    cell adhesion molecule binding
YWHAZ    cell adhesion molecule binding
CTNNA2    cell adhesion molecule binding
ADAM15    cell adhesion molecule binding
LDHA    cell adhesion molecule binding
PKM    cell adhesion molecule binding

KEGG富集分析

输入Gene ID的格式和类型与enrichGO一致。参考

http://guangchuangyu.github.io/2015/02/kegg-enrichment-analysis-with-latest-online-data-using-clusterprofiler/。

# clusterProfiler3.4.4版本是没有readable参数的,原文档有
kk <- enrichKEGG(entrezID, organism="hsa", pvalueCutoff=0.05, pAdjustMethod="BH", 
                 qvalueCutoff=0.1)

输出结果的格式和可视化方式与上面GO富集一致,不再赘述。

另外一个没有解决的问题是setReadable函数的使用 (用测试文档提供的数据集报出如下错误)

setReadable(kk, "org.Hs.eg.db")
Error in EXTID2NAME(OrgDb, genes, keytype) : keytype is not supported...
# 即便是如下操作也没有作用
a = bitr_kegg(names(geneList), fromType = 'ncbi-geneid', toType="kegg", organism="hsa")
# Warning message:
# In bitr_kegg(names(geneList), fromType = "ncbi-geneid", toType = "kegg",  :
#  0.77% of input gene IDs are fail to map...
kk <- enrichKEGG(a$kegg, organism="hsa", keyType="kegg", pvalueCutoff=0.05, pAdjustMethod="BH", 
                  qvalueCutoff=0.1)
setReadable(kk, org.Hs.eg.db, key)
#  Error in EXTID2NAME(OrgDb, genes, keytype) : keytype is not supported...

经过多次尝试发现,可以这么解决

setReadable(kk,  org.Hs.eg.db, keytype="ENTREZID")

为什么会有这个问题呢?setReadable中自动判断keytype的语句是

if (keytype == "auto") {
    keytype <- x@keytype
    if (keytype == "UNKNOWN") {
        stop("can't determine keytype automatically; need to set 'keytype' explicitly...")
        }
}

根据富集结果中定义的keytype,也就是enrichKEGG函数中设定的keyType的值来定的。而setReadable不支持默认的keyType=kegg

没有测试小鼠,可能需要设置不同的keytype值。

另外对拟南芥来说,分析之前需要先把Entrez ID转换为kegg再用上述命令做富集分析

entrezID <- bitr_kegg(entrezID, fromType='ncbi-geneid', toType='kegg', organism="ath")
kk <- enrichKEGG(entrezID$kegg, organism="ath", pvalueCutoff=0.05, pAdjustMethod="BH", 
                 qvalueCutoff=0.1)
# 这个setreadble是可以转换成功的
result <- setReadable(kk, "org.At.tair.db", keytype="TAIR")

GSEA分析

GSEA的解释和介绍见GSEA富集分析 - 界面操作。

注意读入的基因列表是要按照表达差异降序排列 (升序也可以,相当于样品做了对调)。这里排序方式可以是表达差异,也可以是其它方式,只要方便解释即可,即从上到下,或从前到后,基因对表型的贡献有一致的变化趋势就好。不同的排序参数和排序方式需要不同的对结果的解释。

id_with_fc = "ID;FC
4312;2
8318;3
10874;4
55143;5
55388;6
991;7"
id_with_fc <- read.table(text=id_with_fc, header=T, sep=";")
id_with_fc2 <- id_with_fc[,2]
names(id_with_fc2) <- id_with_fc[,1]
# 排序是必须的,记住排序方式
id_with_fc2 <- sort(id_with_fc2, decreasing=T)
gsecc <- gseGO(geneList=id_with_fc2, ont="CC", OrgDb=org.Hs.eg.db, verbose=F)
# 昨天测试了其它数据,参数无问题。这里没有实际运行,盗用数据,每列的解释见本段开头的文章
head(as.data.frame(gsecc))
##                    ID                              Description setSize
## GO:0031982 GO:0031982                                  vesicle    2880
## GO:0031988 GO:0031988                 membrane-bounded vesicle    2791
## GO:0005576 GO:0005576                     extracellular region    3296
## GO:0065010 GO:0065010 extracellular membrane-bounded organelle    2220
## GO:0070062 GO:0070062                    extracellular exosome    2220
## GO:0044421 GO:0044421                extracellular region part    2941
##            enrichmentScore       NES      pvalue   p.adjust    qvalues
## GO:0031982      -0.2561837 -1.222689 0.001002004 0.03721229 0.02816364
## GO:0031988      -0.2572169 -1.226003 0.001007049 0.03721229 0.02816364
## GO:0005576      -0.2746489 -1.312485 0.001009082 0.03721229 0.02816364
## GO:0065010      -0.2570342 -1.222048 0.001013171 0.03721229 0.02816364
## GO:0070062      -0.2570342 -1.222048 0.001013171 0.03721229 0.02816364
## GO:0044421      -0.2744658 -1.310299 0.001014199 0.03721229 0.02816364
# 绘制GSEA图
gseaplot(gsecc, geneSetID="GO:0000779")

自定义数据集分析

如果想用clusterProfiler的函数对自己注释的数据进行功能富集分析或GSEA分析,需要提供如下格式的注释数据。后续分析就类似了。

self_anno <- "ont;gene
KEGG_GLYCOLYSIS_GLUCONEOGENESIS;gene1
KEGG_GLYCOLYSIS_GLUCONEOGENESIS;gene2
KEGG_GLYCOLYSIS_GLUCONEOGENESIS;gene3
KEGG_GLYCOLYSIS;gene1
KEGG_GLYCOLYSIS;gene4 
KEGG_CYP;gene5"
self_anno <- read.table(text=self_anno, header=T, sep=";", quote="")
# 没具体看代码怎么写的,保险期间,设置跟示例一样的列名字
colnames(self_anno) <- c("ont", "gene")
geneL <- c("gene1", "gene2", "gene4")
# self_enrich与之前enrichGO的输出结果格式一致
self_enrich <- enricher(geneL, TERM2GENE=self_anno)
# self_gsea与之前gseGO的输出结果格式一致
self_gsea <- GSEA(geneL, TERM2GENE=self_anno, verbose=F)

Reference

  1. https://guangchuangyu.github.io/clusterProfiler/
  2. http://guangchuangyu.github.io/2016/01/go-analysis-using-clusterprofiler/
  3. http://igraph.org/c/doc/igraph-Layout.html

相关推荐

恢复软件6款汇总推荐,帮你减轻数据恢复压力!

在当今数字化生活中,数据丢失的风险如影随形。无论是误删文件、硬盘故障,还是遭遇病毒攻击,丢失的数据都可能给我们带来不小的麻烦。此时,一款优秀的数据恢复软件就成为了挽救数据的关键。今天,为大家汇总推荐...

中兴星星一号刷回官方原版recovery的教程

【搞科技教程】中兴星星一号的官方recovery也来说一下了,因为之前给大家分享过了第三方的recovery了,之前给大家分享的第三方recovery也是采用一键刷入的方式,如果细心的朋友会发现,之前...

新玩机工具箱,Uotan柚坛工具箱软件体验

以前的手机系统功能比较单调,各厂商的重视程度不一样,所以喜欢玩机的朋友会解锁手机系统的读写权限,来进行刷机或者ROOT之类的操作,让使用体验更好。随着现在的手机系统越来越保守,以及自身功能的增强,...

三星g906k刷recovery教程_三星g906k中文recovery下载

【搞科技教程】看到有一些机友在找三星g906k的第三方recovery,下面就来说一下详细的recovery的刷入方法了,因为手机只有有了第三方的recovery之后才可以刷第三方的root包和系统包...

中兴星星2号刷recovery教程_星星二号中文recovery下载

【搞科技教程】咱们的中兴星星2手机也就是中兴星星二号手机的第三方recovery已经出来了,并且是中文版的,有了这个recovery之后,咱们的手机就可以轻松的刷第三方的系统包了,如果没有第三方的re...

数据恢复软件有哪些值得推荐?这 6 款亲测好用的工具汇总请收好!

在数字生活中,数据丢失的阴霾常常突如其来。无论是误删工作文档、格式化重要磁盘,还是遭遇系统崩溃,都可能让我们陷入焦虑。关键时刻,一款得力的数据恢复软件便是那根“救命稻草”。今天,为大家精心汇总6...

中兴u956刷入recovery的教程(中兴e5900刷机)

【搞科技教程】这次主要来给大家说说中兴u956手机如何刷入第三方的recovery,因为第三方的recovery工具是咱们刷第三方rom包的基础,可是很我欠却不会刷,所以太这里来给大家整理了一下详细的...

联想A850+刷recovery教程 联想A850+第三方recovery下载

【搞科技教程】联想A850+的第三方recovery出来了,这个第三方的recovery是非常的重要的,比如咱们的手机要刷第三方的系统包的时候,都是需要用到这个第三方的recovery的,在网上也是有...

工具侠重大更新 智能机上刷机一条龙完成

工具侠是针对玩机的机油开发的一款工具,不管是发烧级别的粉丝,还是普通小白用户,都可以在工具侠上找到你喜欢的工具应用。这不,最新的工具侠2.0.16版本,更新了专门为小白准备的刷机助手工具,以及MTK超...

shift+delete删除的文件找回6种硬盘数据恢复工具

硬盘作为电脑的重要存储设备,如同一个巨大的数字仓库,承载着我们日常工作、学习和生活中的各种文件,从珍贵的照片、重要的工作文档到喜爱的视频、音乐等,都依赖硬盘来安全存放。但有时,我们可能会不小心用sh...

使用vscode+Deepseek 实现AI编程 基于Cline和continue

尊敬的诸位!我是一名专注于嵌入式开发的物联网工程师。关注我,持续分享最新物联网与AI资讯和开发实战。期望与您携手探寻物联网与AI的无尽可能。这两天deepseek3.0上线,据说编程能力比肩Cl...

详解如何使用VSCode搭建TypeScript环境(适合小白)

搭建Javascript环境因为TypeScript不能直接在浏览器上运行。它需要编译器来编译并生成JavaScript文件。所以需要首先安装好javascript环境,可以参考文章:https://...

使用VSCode来书写你的Jupyter Notebooks

现在你可以在VScode里面来书写你的notebook了,使用起来十分的方便。下面来给大家演示一下环境的搭建。首先需要安装一个jupyter的包,使用下面的命令安装:pip3install-ih...

使用VSCode模板提高Vue开发效率(vscode开发vue插件)

安装VSCode安装Vetur和VueHelper插件,安装完成后需要重启VScode。在扩展插件搜索框中找到如下Vetur和VueHelper两个插件,注意看图标。添加Vue模板打...

干货!VsCode接入DeepSeek实现AI编程的5种主流插件详解

AI大模型对编程的影响非常之大,可以说首当其冲,Cursor等对话式编程工具渐渐渗透到开发者的工作中,作为AI编程的明星产品,Cursor虽然好用,但是贵啊,所以咱们得找平替,最好免费那种。俗话说,不...