gseapy

gseapy

本文比较长，长到需要一个目录

ORA

1.输入数据之基因

2.输入数据之基因集合

代码语言：javascript代码运行次数：0运行复制

2.1 可以是预设的基因集合名称

2.2 可以是来自msigdb网站的gmt

2.也可以是完全自定义的基因集合

.完成富集分析

4.可视化

GSEA

1.输入数据之预排序的基因

2.输入数据之基因集合

.完成GSEA分析

4.可视化

富集分析分为超几何分布检验（ORA）和基因集富集分析（GSEA）。R语言有clusterprofiler包可以做富集，python是用gseapy。这个包功能强大，既支持两种富集分析，还支持基因集gmt的直接获取。

ORA

1.输入数据之基因

就是要一组基因，这组基因是怎么来的都可以，比如是单细胞筛选的marker基因，芯片和转录组的差异基因等。至于格式，不只是支持列表，还支持各种其他格式:param gene_list: str, list, tuple, series, dataframe. Also support input txt file with one gene id per row. 下面的gene_list是用作示例的一组基因。来自limma差异分析的结果。如需示例数据请在【生信星球】聊天框回复【1022gseapy】

代码语言：javascript代码运行次数：0运行复制

import pandas as pd
import gseapy as gp
deg = pd.read_csv("")
deg.head()

代码语言：javascript代码运行次数：0运行复制

gene_list = deg.loc[!="OT","gene"]
gene_list

代码语言：javascript代码运行次数：0运行复制

0          USH2A
1          KC2
2          RMD2
          ASPDH
4           ESR1
          ...   
147       GAS7
14409    S100A14
14560       IL18
14889       IGHM
15096      HABP2
ame: gene, Length: 4475, dtype: object

2.输入数据之基因集合

gene_sets参数，也支持很多格式: str, list, tuple of Enrichr Library name(s).or custom defined gene_sets (dict, or gmt file).

2.1 可以是预设的基因集合名称

KEGG_2021_Human是来自 /Enrichr/#libraries 的基因集合名称之一;这些名称也可以用gp.get_library_name()获取

代码语言：javascript代码运行次数：0运行复制

gene_sets_kegg = 'KEGG_2021_Human'

gp.get_library_name()[0:5]

代码语言：javascript代码运行次数：0运行复制

['ARCHS4_Cell-lines',
 'ARCHS4_IDG_Coexp',
 'ARCHS4_Kinases_Coexp',
 'ARCHS4_TFs_Coexp',
 'ARCHS4_Tissues']

2.2 可以是来自msigdb网站的gmt

代码语言：javascript代码运行次数：0运行复制

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

列出都有哪些版本文件夹

代码语言：javascript代码运行次数：0运行复制

msig.list_dbver()

列出该文件夹下都有哪些基因集合

代码语言：javascript代码运行次数：0运行复制

msig.list_category(dbver="2024.1.Hs") 

代码语言：javascript代码运行次数：0运行复制

['c1.all',
 'c2.all',
 '',
 'biocarta',
 'kegg_legacy',
 'kegg_medicus',
 'pid',
 'reactome',
 '',
 'wikipathways',
 'c.all',
 '_legacy',
 'db',
 '',
 'gtrd',
 'tft_legacy',
 '',
 '',
 'c4.all',
 '',
 '',
 'c5.all',
 'c5.go.bp',
 'c5.go',
 'c5.',
 'c5.go',
 'c5.hpo',
 'c6.all',
 'c7.all',
 'c7.immunesigdb',
 'c7.vax',
 'c8.all',
 'h.all',
 'msigdb']

可以自由的选择啦！

2.也可以是完全自定义的基因集合

写成字典格式即可，形如：dict: gene_sets={'A':['gene1', 'gene2',...], 'B':['gene2', 'gene4',...], ...}

.完成富集分析

代码语言：javascript代码运行次数：0运行复制

enr_kegg = (
    gene_list=gene_list,
    gene_sets=gene_sets_kegg,
    organism='Human',
    outdir='./',
    cutoff=0.05
)

代码语言：javascript代码运行次数：0运行复制

enr_kegg.results.head()

代码语言：javascript代码运行次数：0运行复制

enr_h = (
    gene_list=gene_list,
    gene_sets=gmt,
    organism='Human',
    outdir='./',
    cutoff=0.05
)

代码语言：javascript代码运行次数：0运行复制

enr_h.results.head()

完成富集后，生成了pdf和txt文件，里面分别是条带图和富集结果文件。

代码语言：javascript代码运行次数：0运行复制

import os
os.listdir()

代码语言：javascript代码运行次数：0运行复制

['.ipynb_checkpoints',
 '',
 'enrich.ipynb',
 'gs_ind_0.reports.pdf',
 'gs_ind_0.',
 'KEGG_2021_Human.reports.pdf',
 'KEGG_2021_Human.']

4.可视化

条带图和气泡图

代码语言：javascript代码运行次数：0运行复制

from gseapy.plot import barplot, dotplot
barplot(enr_kegg.res2d)

代码语言：javascript代码运行次数：0运行复制

<AxesSubplot: xlabel='$- \\log_{10}$ (Adjusted P-value)'>

代码语言：javascript代码运行次数：0运行复制

dotplot(enr_h.res2d,size =5)

代码语言：javascript代码运行次数：0运行复制

<AxesSubplot: xlabel='Combined Score'>

GSEA

1.输入数据之预排序的基因

代码语言：javascript代码运行次数：0运行复制

ranking = deg[['gene', 'logFC']]
ranking = ranking.sort_values('logFC', ascending = False).reset_index(drop = True)
ranking

注意，是全部的基因而不是只要差异基因

2.输入数据之基因集合

同前面超几何分布检验的要求↑这里只写一个gmt

代码语言：javascript代码运行次数：0运行复制

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

.完成GSEA分析

代码语言：javascript代码运行次数：0运行复制

pre_res = gp.prerank(rnk = ranking, gene_sets = gmt, seed = 6, permutation_num = 100)

代码语言：javascript代码运行次数：0运行复制

type(pre_res)

代码语言：javascript代码运行次数：0运行复制

gseapy.gsea.Prerank

这个结果就不是那么直观了，比较复杂：

代码语言：javascript代码运行次数：0运行复制

len(list(pre_res.results.keys()))

代码语言：javascript代码运行次数：0运行复制

50

代码语言：javascript代码运行次数：0运行复制

list(pre_res.results.keys())[0:5]

代码语言：javascript代码运行次数：0运行复制

['HALLMARK_MYOGEESIS',
 'HALLMARK_ITERFERO_ALPHA_RESPOSE',
 'HALLMARK_ESTROGE_RESPOSE_EARLY',
 'HALLMARK_UFOLDED_PROTEI_RESPOSE',
 'HALLMARK_APICAL_JUCTIO']

代码语言：javascript代码运行次数：0运行复制

list(pre_res.results.values())[0].keys()

代码语言：javascript代码运行次数：0运行复制

dict_keys(['name', 'es', 'nes', 'pval', 'fdr', 'fwerp', 'tag %', 'gene %', 'lead_genes', 'matched_genes', 'hits', 'RES'])

为了方便阅读，可以把结果转换成一个数据框

代码语言：javascript代码运行次数：0运行复制

out = []

for term in list(pre_res.results):
    out.append([term,
               pre_res.results[term]['fdr'],
               pre_res.results[term]['es'],
               pre_res.results[term]['nes']])

out_df = pd.DataFrame(out, columns = ['Term','fdr', 'es', 'nes']).sort_values('fdr').reset_index(drop = True)
out_df.head()

4.可视化

代码语言：javascript代码运行次数：0运行复制

term_to_graph = out_df.iloc[0].Term
term_to_graph

代码语言：javascript代码运行次数：0运行复制

'HALLMARK_PEROXISOME'

代码语言：javascript代码运行次数：0运行复制

gp.gseaplot(rank_metric = pre_res.ranking,term = term_to_graph, **pre_res.results[term_to_graph])

代码语言：javascript代码运行次数：0运行复制

[<Axes: xlabel='Gene Rank', ylabel='Ranked metric'>,
 <Axes: >,
 <Axes: >,
 <Axes: ylabel='Enrichment Score'>]

本文参与 腾讯云自媒体同步曝光计划，分享自。原始发表：2024-11-24，如有侵权请联系 cloudcommunity@tencent 删除python集合可视化排序数据