本文比较长，长到需要一个目录
ORA
1.输入数据之基因
2.输入数据之基因集合
2.1 可以是预设的基因集合名称
2.2 可以是来自msigdb网站的gmt
2.3也可以是完全自定义的基因集合
3.完成富集分析
4.可视化
GSEA
1.输入数据之预排序的基因
2.输入数据之基因集合
3.完成GSEA分析
4.可视化
ORA
富集分析分为超几何分布检验（ORA）和基因集富集分析（GSEA）。R语言有clusterprofiler包可以做富集，python是用gseapy。这个包功能强大，既支持两种富集分析，还支持基因集gmt的直接获取。
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'>]
