关于 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)
)
