gseapy-python版的富集分析

[ixxmu]

https://mp.weixin.qq.com/s/rv_zi_de1lPNPHdrEot4FA

[ixxmu]

gseapy-python版的富集分析 by 生信星球

 今天是生信星球陪你的第1022天

---

   

公众号里的文章大多数需要编程基础，如果因为代码看不懂，而跟不上正文的节奏，可以来找我学习，相当于给自己一个新手保护期。我的课程都是循环开课，点进去咨询微信↓

生信分析直播课程(每月初开一期，春节休一个月)

生信新手保护学习小组（每月两期）

单细胞陪伴学习小组（每月两期

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

ORA

1.输入数据之基因

2.输入数据之基因集合

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

2.2 可以是来自msigdb网站的gmt

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

3.完成富集分析

4.可视化

GSEA

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

2.输入数据之基因集合

3.完成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】

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

gene_list = deg.loc[deg.change!="NOT","gene"]
gene_list

0          USH2A
1          KCNN2
2          RMDN2
3          ASPDH
4           ESR1
          ...   
14373       GAS7
14409    S100A14
14560       IL18
14889       IGHM
15096      HABP2
Name: 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是来自 https://maayanlab.cloud/Enrichr/#libraries 的基因集合名称之一;这些名称也可以用gp.get_library_name()获取

gene_sets_kegg = 'KEGG_2021_Human'

gp.get_library_name()[0:5]

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

2.2 可以是来自msigdb网站的gmt

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

列出都有哪些版本文件夹

msig.list_dbver()

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

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

['c1.all',
 'c2.all',
 'c2.cgp',
 'c2.cp.biocarta',
 'c2.cp.kegg_legacy',
 'c2.cp.kegg_medicus',
 'c2.cp.pid',
 'c2.cp.reactome',
 'c2.cp',
 'c2.cp.wikipathways',
 'c3.all',
 'c3.mir.mir_legacy',
 'c3.mir.mirdb',
 'c3.mir',
 'c3.tft.gtrd',
 'c3.tft.tft_legacy',
 'c3.tft',
 'c4.3ca',
 'c4.all',
 'c4.cgn',
 'c4.cm',
 'c5.all',
 'c5.go.bp',
 'c5.go.cc',
 'c5.go.mf',
 'c5.go',
 'c5.hpo',
 'c6.all',
 'c7.all',
 'c7.immunesigdb',
 'c7.vax',
 'c8.all',
 'h.all',
 'msigdb']

可以自由的选择啦！

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

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

3.完成富集分析

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

enr_kegg.results.head()

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

enr_h.results.head()

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

import os
os.listdir()

['.ipynb_checkpoints',
 'deg.csv',
 'enrich.ipynb',
 'gs_ind_0.Human.enrichr.reports.pdf',
 'gs_ind_0.Human.enrichr.reports.txt',
 'KEGG_2021_Human.Human.enrichr.reports.pdf',
 'KEGG_2021_Human.Human.enrichr.reports.txt']

4.可视化

条带图和气泡图

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

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

dotplot(enr_h.res2d,size =5)

<AxesSubplot: xlabel='Combined Score'>

GSEA

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

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

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

2.输入数据之基因集合

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

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

3.完成GSEA分析

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

type(pre_res)

gseapy.gsea.Prerank

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

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

50

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

['HALLMARK_MYOGENESIS',
 'HALLMARK_INTERFERON_ALPHA_RESPONSE',
 'HALLMARK_ESTROGEN_RESPONSE_EARLY',
 'HALLMARK_UNFOLDED_PROTEIN_RESPONSE',
 'HALLMARK_APICAL_JUNCTION']

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

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

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

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.可视化

term_to_graph = out_df.iloc[0].Term
term_to_graph

'HALLMARK_PEROXISOME'

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

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