python实现富集分析（Enrichment Analysis）

浅沫～

已于 2022-10-24 20:44:48 修改

阅读量4.6k

点赞数 4

分类专栏： python 生物信息学文章标签： KEGG python 1024程序员节

于 2022-09-19 19:53:46 首次发布

本文链接：https://blog.csdn.net/qq_39579290/article/details/126937057

版权

python 同时被 2 个专栏收录

6 篇文章 6 订阅

订阅专栏

生物信息学

3 篇文章 1 订阅

订阅专栏

2.5.2 KEGG pathway通路富集

一、引言

基因组可以算突变率、转录组可以算基因表达水平、表观组可以观察甲基化水平、蛋白质组/代谢组可以看蛋白/代谢物丰度，甚至是宏基因组也可以比较菌群的丰度。最终，在分子水平的出口都是在生物学中心法则的核酸水平，确切的说是在基因水平，但是基因的种类有很多，包括蛋白编码基因（mRNA）、非编码基因（miRNA、lncRNA、snRNA等），理解这些基因所代表的生物学意义的最佳途径就是基因富集分析。

参考网站：一文读懂基因富集分析：从小白到SCI - 知乎 (zhihu.com)

二、富集分析

2.1 富集分析的概念

一个生物过程（生物学过程）通常是由一组基因共同参与，而不是由单个基因独自完成。如果一个生物学过程在已知的研究中发生异常，则共同发挥功能的基因极可能被选择出来作为一个与这一过程相关的基因集合。

富集分析（Gene Set Enrichment Analysis, GSEA）通常是分析一组基因在某个功能节点上是否相比于随机水平过于出现（over-presentation）。

富集分析原理可以由单个基因的简单注释扩展到多个基因集合的成组分析。

2.2 富集分析的原理

超几何分布时一种离散型概率分布，也许中学就学过最经典的例子：假设一个袋子有10个球，其中红色球6个，白色球4个，那么我一次从袋子中抓取4个，请问这4个球中至少3个是红球的概率是多少？

富集分析原理就是超几何分布的使用。假设人类背景基因3W个，已知共500个基因参与某一通路，由已知我们有一个自己的含有200个基因的基因群，且这200个基因都来自背景基因，恰巧这200个基因里面有6个参与此通路，那么这6个基因是偶然还是显著富集在此通路上，就需要进一步做相应的显著性检验。

参考网址：python富集分析 - CSDN

富集分析，本质上是对数据的分布检验，如果分布集中在某个区域，则认为富集。

此外，个人认为富集分析：（一）概述_生信技工的博客-CSDN博客_富集分析对富集分析描述的也十分清晰。

2.3 富集分析的作用

一组基因直接注释的结果是得到大量的功能节点，这些功能具有概念上的交叠现象，导致分析结果冗余，不利于进一步的精细分析，所以研究人员希望对得到的功能节点加以过滤和筛选，以便获得更有意义的功能信息。目前最常用的方法是基于GO和KEGG的富集分析。

也就是说，富集分析是用来解读一组基因背后所代表的生物学知识，揭示其在细胞内或细胞外扮演了什么样的角色。

2.4 富集分析的步骤

首先通过多种方法找到大量的感兴趣的基因，例如差异表达基因集、共表达基因模块、蛋白质复合物基因簇等，然后寻找这些感兴趣基因集显著富集的GO节点或者KEGG通路，这有助于进一步深入细致的实验研究。

2.5 富集分析的类型

依据富集分析过程中基因选择、注释数据库的不同，常用的富集分析可以分为以下四种类型：GO term功能富集、KEGG pathway通路富集、MSigDB基因集富集和单基因富集等等。

2.5.1 GO term功能富集

基因本体（gene ontology， GO）数据库是GO组织在2000年构建的一个结构化的标准生物学模型，涵盖了细胞组分、分子功能、生物学过程三个方面，是目前应用最广泛的基因注释体系之一。GO的注释体系是一个有向无环图，包含三个分支，注释系统中每一个节点都是基因或蛋白质的一种描述，节点之间保持严格的“父子”关系。因此，一个基因或蛋白质可以从三个层面得到注释。

本体论通常处理的问题：存在哪些本质，如何将这些本质分组，在层次结构内关联以及如何根据相似性和差异进行细分。

参考网址：一文极速读懂 Gene Ontology （GO）数据库 - 腾讯云开发者社区-腾讯云 (tencent.com)

2.5.2 KEGG pathway通路富集

京都基因与基因组百科全书（Kyoto encyclopedia of genes and genomes, KEGG）是系统分析基因功能、基因组信息的数据库，整合了基因组学、生物化学及系统功能组学的信息，有助于研究者把基因及表达信息作为一个整体进行研究。目前KEGG共包含了19个子数据库，富集分析常用在KEGG Pathway通路中。

2.5.3 MSigDB基因集富集

MSigDB数据库定义了已知的基因集合，包括H和C1-C7八个系列（Collection）。H: hallmark gene sets （效应）特征基因集合，共50组；C1: positional gene sets 位置基因集合，根据染色体位置，共326个；C2: curated gene sets：（专家）共识基因集合，基于通路、文献等，包括KEGG；C3: motif gene sets：模式基因集合，主要包括microRNA和转录因子靶基因两部分；C4: computational gene sets：计算基因集合，通过挖掘癌症相关芯片数据定义的基因集合；C5: GO gene sets：Gene Ontology 基因本体论；C6: oncogenic signatures：癌症特征基因集合，大部分来源于NCBI GEO 未发表芯片数据；C7: immunologic signatures: 免疫相关基因集合。可以从中获取大量的已知基因集合从而进行富集分析。

2.5.4 单基因富集

单基因富集分析并不是说拿单个基因来进行富集分析，一个基因根本没法进行统计检验。而是基于单个基因来抓取与其相关的基因，然后用这些相关的基因来进行功能富集，有两种方法：差异法和相关法。

（1)差异法：根据给定的一个基因的表达值对样本进行分组，然后计算组间的差异表达基因，进而利用差异基因进行富集分析。

（2)相关法：计算给定的一个基因的表达值与其他基因之间的相关性，将具有显著相关的基因作为一个集合进行富集分析。

参考网站：一文读懂基因富集分析：从小白到SCI - 知乎 (zhihu.com)

三、富集分析中涉及的基本概念

3.1 细胞组分

细胞组分（cellular component）：细胞的每个部分和细胞外环境。

3.2 分子功能

分子功能（molecular function）：可以描述为分子水平的活性（activity），如催化（catalytic）或结合（binding）活性。

3.3 生物过程

生物过程（biological process）：生物学过程指由一个或多个分子功能有序组合而产生的系列事件。较大的过程或通过多种分子活动完成的“生物程序”。更加广泛的生物学过程术语的例子是DNA修复或信号转导。更具体的术语的例子是葡萄糖跨膜转运的嘧啶核苷生物合成过程。

3.4 功能节点

参考网址：Biological Pathways 生物途径（通路）_ysugarr的博客-CSDN博客_生物通路

3.5 背景基因

背景基因集 常常是一个物种的基因组注释基因总和。

一组特定基因集 常常是差异表达基因集(differentially expressed genes, DEGs)。

预先构建好基因注释数据库常用GO或KEGG数据库。

常见的统计学算法包括ORA,FCS,PT,NT四种。

四、python代码

import gseapy as gp
from gseapy.plot import barplot, dotplot
gene_sets = 'KEGG_2016'
gene_list = ['DAB2IP', 'AKT1', 'ARNT', 'FLCN', 'XRCC6', 'BTG1', 'STAT1', 'STAT1']
# 富集
enr = gp.enrichr(gene_list=gene_list,#所需查询gene_list，可以是一个列表，也可为文件（一列，每行一个基因）
                 gene_sets=gene_sets,#gene set library，多个相关的gene set 。如所有GO term组成一个gene set library.
                 organism='Human',#持(human, mouse, yeast, fly, fish, worm)， 自定义gene_set 则无影响。
                 description='kegg',#工作运行描述
                 outdir='./data/enrichr',#输出目录
                 top_term=20,
                 cutoff=0.5#pvalue阈值
                 )
dot_png = "KEGG_2016" + "_" + "dot" + ".png"
bar_png = "KEGG_2016" + "_" + "bar" + ".png"
base_path = "./static/images/"
if os.path.exists(base_path+dot_png):
    os.remove(base_path+dot_png)
#dot
# to save your figure, make sure that ``ofname`` is not None
dotplot(enr.results.loc[enr.results["Gene_set"] == "KEGG_2016",], title='KEGG Top20 Pathway', cmap='viridis_r',
        top_term=20, legend="r",
        ofname=base_path+dot_png,
        )
#bar
if os.path.exists(base_path+bar_png):
    os.remove(base_path+bar_png)
barplot(enr.res2d, title='KEGG_2013',top_term=20,ofname=base_path+bar_png,)

气泡图展示如下：