python版本的功能富集分析：GSEApy

关于 R 版本的 富集分析笔记太多了，R 生态的生信分析笔记超全。但是随着数据量的日益壮大，我们有必要开始学习python了。生信技能树从今年开始会大力推行 python 版本的生信生态，写超多关于 python 版本的生信分析教程。敬请关注~新专辑《python生信笔记2025》。

今天是第一篇，我们生信技能树前面推出了一个python单细胞课程：《掌握Python，解锁单细胞数据的无限可能》。群里学员问的最多的一个问题就是: python版本的功能富集分析如何做？一起来看看 python包：GSEApy。

学习笔记最好的网页当然是官网，当然也可以去看其他人写的各种笔记，也有官网没有提到的各种小细节。

GSEApy包官网：https://gseapy.readthedocs.io/en/latest/

GSEApy 包简介

这个包的名字全称为 GSEAPY: Gene Set Enrichment Analysis in Python，就是在python中的基因集富集分析。

他主要是 GSEA（https://www.gsea-msigdb.org/gsea/index.jsp） 和 Enrichr（http://amp.pharm.mssm.edu/Enrichr） 两种富集分析方法的封装。

适用数据：RNA-seq, ChIP-seq, Microarry data。

这个方法于2022年发表在生信的经典老牌杂志 Bioinformatics 上：

Zhuoqing Fang, Xinyuan Liu, Gary Peltz, GSEApy: a comprehensive package for performing gene set enrichment analysis in Python, Bioinformatics, 2022;, btac757, https://doi.org/10.1093/bioinformatics/btac757

软件安装

我们前面创建了一个python=3.9的conda小环境叫sc，然后在这个环境里面进行安装：

# bash终端
conda activate sc
# 安装 gseapy
pip install gseapy -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple

安装完成可以检查一下是否可以成功导入模块，我这里运行 python 的环境为 vscode+jupter 插件连接服务器，超级方便。

# 导入模块
import pandas as pd
import gseapy as gp
import matplotlib.pyplot as plt

# 查看 gseapy 的版本
gp.__version__

gseapy 版本为1.14版本：

差异基因列表

使用前面做过的一个芯片数据的差异结果吧：2万个基因少一半也不影响最后的差异分析富集结果啊？

数据为：https://ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE17351，拿到之后进行芯片预处理并做差异表达分析得到一个差异结果，或者 微信找我发你：Biotree123。

或者百度云盘：链接: https://pan.baidu.com/s/1sAXzlxs4jU24ZcAW4LThqQ?pwd=uavp 提取码: uavp 。

读取差异基因结果：

# 读取差异分析结果
# 记住文件路径改成自己的
# index_col=0：第一列读取为行名
deg = pd.read_csv("./GSE17351/DEG.csv", index_col=0)
deg.head()

总共有14080个基因，9列：

# 查看有多少行，有多少列
deg.shape
# (14080, 9)

获取显著差异表达基因列表：共有1895个差异基因。

gene_list = deg.loc[deg.g!="stable","name"]
gene_list

获取基因集

基因集这里有好几种获取方式，我们这里看两个来源的吧。

来源一：下载 Msigdb API

使用 get_gmt函数获取基因集：这里选择了 Msigdb数据库 2024.1.Hs 人类版本的 h.all 即hallmark通路基因集，共50个通路。

# 基因集：使用来自 Msigdb 数据库的
# gp为前面加载的gseapy模块的缩写 import gseapy as gp
msig = gp.Msigdb()
gmt = msig.get_gmt(category='h.all', dbver="2024.1.Hs")

简单看一下gmt的属性：

msig
gmt
type(gmt)
gmt.keys()
gmt.values()
print(gmt['HALLMARK_ADIPOGENESIS'])
print(gmt['HALLMARK_WNT_BETA_CATENIN_SIGNALING'])

get_gmt()函数有两个参数：

• dbver="2024.1.Hs"  这个参数的选择可以使用 list_dbver()查看，选择数据库的版本，最下面的是最新的
• category='h.all' 这个参数的选择可以使用 list_category查看，看数据库下面的基因集合名称

# list msigdb version you wanna query
msig.list_dbver()

# list categories given dbver. 列出该数据库版本下都有哪些基因集合
msig.list_category(dbver="2024.1.Hs") # human

其他的简单探索，不同物种，不同基因集，随便看！

    
# mouse hallmark gene sets
gmt = msig.get_gmt(category='mh.all', dbver="2023.1.Mm")
# 改成 2024.1.Hs
gmt = msig.get_gmt(category='h.all', dbver="2024.1.Hs")
# 改成 2024.1.Hs,c5.go.bp
gmt = msig.get_gmt(category='c5.go.bp', dbver="2024.1.Hs")
print(gmt['HALLMARK_WNT_BETA_CATENIN_SIGNALING'])

报错：

---

ModuleNotFoundError                       Traceback (most recent call last) File /nas2/zhangj/biosoft/miniconda3/envs/sc/lib/python3.9/site-packages/pandas/compat/_optional.py:135, in import_optional_dependency(name, extra, errors, min_version) ...

File /nas2/zhangj/biosoft/miniconda3/envs/sc/lib/python3.9/site-packages/pandas/compat/_optional.py:138, in import_optional_dependency(name, extra, errors, min_version) 136 except ImportError: 137     if errors == "raise": --> 138         raise ImportError(msg) 139     return None 141 # Handle submodules: if we have submodule, grab parent module from sys.modules

ImportError: Missing optional dependency 'lxml'.  Use pip or conda to install lxml.

缺个模块 lxml，那就安装：

# bash终端
pip install lxml -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple

来源二：查看所有可支持的 enrichr library names

来自 https://maayanlab.cloud/Enrichr/#libraries 的基因集合名称，可以使用 gp.get_library_name() 获取

支持物种：{ ‘Human’, ‘Mouse’, ‘Yeast’, ‘Fly’, ‘Fish’, ‘Worm’ }

以下命令去试试看吧：

# default: Human
# 这些名称也可以用 gp.get_library_name() 获取
names = gp.get_library_name()
names[:10]

# yeast
yeast = gp.get_library_name(organism='Yeast')
yeast[:10]

# mouse
mouse = gp.get_library_name(organism='Mouse')
mouse[:10]

## download library or read a .gmt file
go_mf = gp.get_library(name='GO_Molecular_Function_2018', organism='Human')
print(go_mf['ATP binding (GO:0005524)'])

富集分析：Over-representation analysis

现在基因列表有了，基因集也有了，可以做ORA类的功能富集分析了，使用的函数为 enrichr()，其参数细节可以使用下面的命令查看：

gp.enrichr?

非常简单吧：

选择 hallmark通路进行分析：

msig = gp.Msigdb()
gmt = msig.get_gmt(category='h.all', dbver="2024.1.Hs")

enr_h = gp.enrichr(gene_list=gene_list, 
                   gene_sets=gmt,
                   organism='Human',  # 设置物种
                   outdir='./', # 默认输出结果到当前目录
                   cutoff=0.99, # 设置 校正后pvalue的阈值
                   verbose=True)

运行成功：

看一下结果：目录中会生成两个文件 gs_ind_0.Human.enrichr.reports.txt、gs_ind_0.Human.enrichr.reports.pdf

enr_h.results.head()

绘制富集结果

1、条形图 barplot

• figsize=(4,5)：控制绘图的宽和高
• color = 'darkblue'：设置绘制图形的颜色

# 条形图
# simple plotting function
from gseapy import barplot, dotplot
barplot(enr_h.res2d, 
  column='P-value', 
  cutoff=0.9,
  top_term=15,
  figsize=(4,5),
  color = 'darkblue'
  )

结果如下：

2、气泡图 dotplot

我们今天的次页给大家介绍了一个超棒的配色资源：https://python-graph-gallery.com/color-palette-finder/。

试试看：

# https://python-graph-gallery.com/color-palette-finder/
from pypalettes import load_cmap
#cmap = load_cmap("Althoff")
cmap = load_cmap("CeriseLimon")
dotplot(enr_h.res2d,
  size=10,
  column='P-value', 
  cutoff=0.9,
  top_term=15,
  cmap=cmap,
  marker='o',
  x='Combined Score', # set x axis, 
  title = "HALLMARK Pathway",
  figsize=(4,6)
  )

分享到此，你开始学习python了吗？

友情宣传：

生信入门&数据挖掘线上直播课2025年1月班

时隔5年，我们的生信技能树VIP学徒继续招生啦

满足你生信分析计算需求的低价解决方案