Pathway analysis with AUCell¶
Single-cell RNA sequencing (scRNA-seq) is a powerful tool for exploring variations in cell types between conditions, tissue types, species, and individuals. When conducting scRNA-seq analysis, the differential gene expression (DEG) analysis of the single-cell data is almost always followed by gene set enrichment analysis. The aim of this analysis is to identify gene programs and biological processes, gene ontologies, or regulatory pathways that are overrepresented in a case group compared to a control group.
Paper: AUCell: Identifying cells with active gene sets
Code: https://github.com/aertslab/AUCell
Colab_Reproducibility:https://colab.research.google.com/drive/1Rk7Zopil-Ve1WFRQLV_AZOwCHRyhsOXk?usp=sharing
Warning
There are many methods to determine the enrichment pathway between two groups, and the choice of method can significantly impact the conclusion.
This tutorial focuses on using AUCell to complete the gene set enrichment in scRNA-seq data.
Part.1 The Mathematical Principles of AUCell¶
AUCell uses the “Area Under the Curve” (AUC) to calculate whether a critical subset of the input gene set is enriched within the expressed genes for each cell. The distribution of AUC scores across all the cells allows exploring the relative expression of the signature. Since the scoring method is ranking-based, AUCell is independent of the gene expression units and the normalization procedure.
In brief, the scoring method is based on a recovery analysis where the x-axis is the ranking of all genes based on expression level (genes with the same expression value, e.g., '0', are randomly sorted); and the y-axis is the number of genes recovered from the input set. AUCell then uses the AUC to calculate whether a critical subset of the input gene set is enriched at the top of the ranking for each cell. In this way, the AUC represents the proportion of expressed genes in the signature and their relative expression values compared to the other genes within the cell. The output of this step is a matrix with the AUC score for each gene set in each cell.
Part.2 Data preprocess¶
In this part, we load a test data and perform preliminary processing of the data, such as normalization and logarithmization, in order to make the data more interpretable
import omicverse as ov
ov.style(font_path='Arial')
🔬 Starting plot initialization...
Using already downloaded Arial font from: /tmp/omicverse_arial.ttf
Registered as: Arial
🧬 Detecting GPU devices…
✅ NVIDIA CUDA GPUs detected: 1
• [CUDA 0] NVIDIA GeForce RTX 3090
Memory: 23.7 GB | Compute: 8.6
🔖 Version: 1.7.9rc1 📚 Tutorials: https://omicverse.readthedocs.io/
✅ plot_set complete.
ov.utils.download_pathway_database()
ov.utils.download_geneid_annotation_pair()
......Pathway Geneset download start: GO_Biological_Process_2021 Using Stanford mirror for GO_Biological_Process_2021 ⚠️ File ./genesets/GO_Biological_Process_2021.txt already exists ......Pathway Geneset download start: GO_Cellular_Component_2021 Using Stanford mirror for GO_Cellular_Component_2021 ⚠️ File ./genesets/GO_Cellular_Component_2021.txt already exists ......Pathway Geneset download start: GO_Molecular_Function_2021 Using Stanford mirror for GO_Molecular_Function_2021 ⚠️ File ./genesets/GO_Molecular_Function_2021.txt already exists ......Pathway Geneset download start: WikiPathway_2021_Human Using Stanford mirror for WikiPathway_2021_Human ⚠️ File ./genesets/WikiPathway_2021_Human.txt already exists ......Pathway Geneset download start: WikiPathways_2019_Mouse Using Stanford mirror for WikiPathways_2019_Mouse ⚠️ File ./genesets/WikiPathways_2019_Mouse.txt already exists ......Pathway Geneset download start: Reactome_2022 Using Stanford mirror for Reactome_2022 ⚠️ File ./genesets/Reactome_2022.txt already exists ✅ Pathway Geneset download finished! ......Other Genesets can be downloaded from https://maayanlab.cloud/Enrichr/#libraries ......Geneid Annotation Pair download start: pair_GRCm39 Using Stanford mirror for pair_GRCm39 ⚠️ File ./genesets/pair_GRCm39.tsv already exists ......Geneid Annotation Pair download start: pair_T2TCHM13 Using Stanford mirror for pair_T2TCHM13 ⚠️ File ./genesets/pair_T2TCHM13.tsv already exists ......Geneid Annotation Pair download start: pair_GRCh38 Using Stanford mirror for pair_GRCh38 ⚠️ File ./genesets/pair_GRCh38.tsv already exists ......Geneid Annotation Pair download start: pair_GRCh37 Using Stanford mirror for pair_GRCh37 ⚠️ File ./genesets/pair_GRCh37.tsv already exists ......Geneid Annotation Pair download start: pair_danRer11 Using Stanford mirror for pair_danRer11 ⚠️ File ./genesets/pair_danRer11.tsv already exists ......Geneid Annotation Pair download start: pair_danRer7 Using Stanford mirror for pair_danRer7 ⚠️ File ./genesets/pair_danRer7.tsv already exists ......Geneid Annotation Pair download start: pair_hgnc_all ⚠️ File ./genesets/pair_hgnc_all.tar.gz already exists ......Extracted pair_hgnc_all.tsv from tar.gz ✅ Geneid Annotation Pair download finished!
The dataset used here is the mouse pancreas dataset that comes with scvelo
adata = ov.datasets.pancreatic_endocrinogenesis()
adata
⚠️ File ./data/endocrinogenesis_day15.h5ad already exists Loading data from ./data/endocrinogenesis_day15.h5ad ✅ Successfully loaded: 3696 cells × 27998 genes
AnnData object with n_obs × n_vars = 3696 × 27998
obs: 'clusters_coarse', 'clusters', 'S_score', 'G2M_score'
var: 'highly_variable_genes'
uns: 'clusters_coarse_colors', 'clusters_colors', 'day_colors', 'neighbors', 'pca'
obsm: 'X_pca', 'X_umap'
layers: 'spliced', 'unspliced'
obsp: 'distances', 'connectivities'
adata.X.max()
2286.0
We found that the max value of anndata object larger than 10 and type is int. We need to normalize and log1p it
adata.var['robust']=True
adata=ov.pp.preprocess(adata, target_sum=1e4,identify_robust_genes=False)
🔍 [2026-02-19 03:51:02] Running preprocessing in 'cpu' mode... Begin robust gene identification ✅ Robust gene identification completed successfully. Begin size normalization: shiftlog and HVGs selection pearson 🔍 Count Normalization: Target sum: 10000.0 Exclude highly expressed: True Max fraction threshold: 0.2 ⚠️ Excluding 1 highly-expressed genes from normalization computation Excluded genes: ['Ghrl'] ✅ Count Normalization Completed Successfully! ✓ Processed: 3,696 cells × 27,998 genes ✓ Runtime: 0.35s 🔍 Highly Variable Genes Selection (Experimental): Method: pearson_residuals Target genes: 2,000 Theta (overdispersion): 100 ✅ Experimental HVG Selection Completed Successfully! ✓ Selected: 2,000 highly variable genes out of 27,998 total (7.1%) ✓ Results added to AnnData object: • 'highly_variable': Boolean vector (adata.var) • 'highly_variable_rank': Float vector (adata.var) • 'highly_variable_nbatches': Int vector (adata.var) • 'highly_variable_intersection': Boolean vector (adata.var) • 'means': Float vector (adata.var) • 'variances': Float vector (adata.var) • 'residual_variances': Float vector (adata.var) Time to analyze data in cpu: 3.44 seconds. ✅ Preprocessing completed successfully. Added: 'highly_variable_features', boolean vector (adata.var) 'means', float vector (adata.var) 'variances', float vector (adata.var) 'residual_variances', float vector (adata.var) 'counts', raw counts layer (adata.layers) End of size normalization: shiftlog and HVGs selection pearson ╭─ SUMMARY: preprocess ──────────────────────────────────────────────╮ │ Duration: 3.4821s │ │ Shape: 3,696 x 27,998 (Unchanged) │ │ │ │ CHANGES DETECTED │ │ ──────────────── │ │ ● VAR │ ✚ highly_variable (bool) │ │ │ ✚ highly_variable_features (bool) │ │ │ ✚ highly_variable_rank (float) │ │ │ ✚ means (float) │ │ │ ✚ residual_variances (float) │ │ │ ✚ variances (float) │ │ │ │ ● UNS │ ✚ REFERENCE_MANU │ │ │ ✚ hvg │ │ │ ✚ log1p │ │ │ ✚ status │ │ │ ✚ status_args │ │ │ ╰────────────────────────────────────────────────────────────────────╯
adata.X.max()
8.059913
Now the max value of anndata object is 8.06
Part.3 Pathway anaylsis¶
In this part, we will demonstrate how to utilize the existing gene set to conduct enrichment analysis and how to create a gene set based on our own ideas for enrichment analysis in the test dataset.
First, we need to download data sets, such as the Gene Ontology(GO) or the Kyoto Encyclopedia of Genes and Genomes(KEGG). It should be noted that here we need to select the correct species, such as 'Human'or 'Mouse'.
pathway_dict=ov.utils.geneset_prepare('genesets/GO_Biological_Process_2021.txt',organism='Mouse')
When working with existing datasets, it is possible to use the ov.single.geneset_aucell to calculate the activity of a gene set that corresponds to a particular signaling pathway within the dataset.
Additionally, we can use the sc.pl.embedding function to visualize the distribution of gene set activity. By doing so, we can gain insights into the behavior of the gene set within the dataset and how it relates to the signaling pathway of interest.
##Assest one geneset
geneset_name='response to vitamin (GO:0033273)'
ov.single.geneset_aucell(adata,
geneset_name=geneset_name,
geneset=pathway_dict[geneset_name])
ov.pl.umap(
adata,
color=["{}_aucell".format(geneset_name)],
cmap='Reds',
)
We also can calculate the AUCell of more than one geneset
##Assest more than one geneset
geneset_names=['response to vitamin (GO:0033273)','response to vitamin D (GO:0033280)']
ov.single.pathway_aucell(
adata,
pathway_names=geneset_names,
pathways_dict=pathway_dict
)
ov.pl.umap(
adata,
color=[i+'_aucell' for i in geneset_names],
cmap='Reds',
)
In certain situations, the pathway we wish to investigate may not be available in the database. In such cases, we can manually define the gene set and its corresponding genes, calculate the AUCell activity of the gene set, and then visualize the results. By doing so, we can gain insights into the behavior of the gene set and its associated pathway, even in the absence of a pre-existing database entry.
##Assest test geneset
ov.single.geneset_aucell(
adata,
geneset_name='Sox',
geneset=['Sox17', 'Sox4', 'Sox7', 'Sox18', 'Sox5']
)
ov.pl.embedding(
adata,
basis='X_umap',
color=["Sox_aucell"],
cmap='Reds'
)
Occasionally, we may wish to examine clusters-specific signaling pathways from a comprehensive perspective. In such cases, we can compute the AUCell scores of all signaling pathways in the database and store the results as an anndata file. Once saved, we can utilize differential expression gene calculation functions and visualization functions from Scanpy to conduct downstream analyses. By taking this approach, we can investigate signaling pathways that are specific to certain clusters and gain insights into the underlying biological mechanisms.
##Assest all pathways
adata_aucs=ov.single.pathway_aucell_enrichment(
adata,
pathways_dict=pathway_dict,
num_workers=8,
gene_overlap_threshold=0.5,
)
00%|██████████| 3696/3696 [00:02<00:00, 1768.88it/s]
Computing AUC scores for 6036 pathways using 8 workers... Splitting 6036 pathways into 8 chunks of ~755 pathways each... Starting parallel pathway processing...
Only 36.4% of the genes in cytidine to uridine editing (GO:0016554) are present in the expression matrix (threshold: 50.0%). Only 30.8% of the genes in negative regulation of gene silencing (GO:0060969) are present in the expression matrix (threshold: 50.0%). Only 40.0% of the genes in negative regulation of histone H3-K27 methylation (GO:0061086) are present in the expression matrix (threshold: 50.0%). Only 4.0% of the genes in detection of chemical stimulus involved in sensory perception (GO:0050907) are present in the expression matrix (threshold: 50.0%). Only 25.0% of the genes in detection of chemical stimulus involved in sensory perception of bitter taste (GO:0001580) are present in the expression matrix (threshold: 50.0%). Only 0.0% of the genes in detection of chemical stimulus involved in sensory perception of smell (GO:0050911) are present in the expression matrix (threshold: 50.0%). Only 46.2% of the genes in positive regulation of T cell mediated cytotoxicity (GO:0001916) are present in the expression matrix (threshold: 50.0%). Only 31.8% of the genes in detection of chemical stimulus involved in sensory perception of taste (GO:0050912) are present in the expression matrix (threshold: 50.0%). Only 45.5% of the genes in regulation of execution phase of apoptosis (GO:1900117) are present in the expression matrix (threshold: 50.0%). Only 40.0% of the genes in alkaloid metabolic process (GO:0009820) are present in the expression matrix (threshold: 50.0%). Only 20.0% of the genes in leukotriene B4 metabolic process (GO:0036102) are present in the expression matrix (threshold: 50.0%). Only 22.2% of the genes in DNA cytosine deamination (GO:0070383) are present in the expression matrix (threshold: 50.0%). Only 28.6% of the genes in lipid hydroxylation (GO:0002933) are present in the expression matrix (threshold: 50.0%). Only 38.5% of the genes in DNA deamination (GO:0045006) are present in the expression matrix (threshold: 50.0%). Only 40.0% of the genes in rescue of stalled ribosome (GO:0072344) are present in the expression matrix (threshold: 50.0%). Only 14.3% of the genes in antigen processing and presentation of endogenous peptide antigen via MHC class I via ER pathway (GO:0002484) are present in the expression matrix (threshold: 50.0%). Only 14.3% of the genes in antigen processing and presentation of endogenous peptide antigen via MHC class I via ER pathway, TAP-independent (GO:0002486) are present in the expression matrix (threshold: 50.0%). Only 25.0% of the genes in antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP-independent (GO:0002480) are present in the expression matrix (threshold: 50.0%). Only 0.0% of the genes in antigen processing and presentation of lipid antigen via MHC class Ib (GO:0048003) are present in the expression matrix (threshold: 50.0%). Only 0.0% of the genes in antigen processing and presentation, endogenous lipid antigen via MHC class Ib (GO:0048006) are present in the expression matrix (threshold: 50.0%). Only 0.0% of the genes in antigen processing and presentation, exogenous lipid antigen via MHC class Ib (GO:0048007) are present in the expression matrix (threshold: 50.0%). Only 43.8% of the genes in negative regulation of natural killer cell mediated cytotoxicity (GO:0045953) are present in the expression matrix (threshold: 50.0%). Only 45.5% of the genes in negative regulation of natural killer cell mediated immunity (GO:0002716) are present in the expression matrix (threshold: 50.0%). Only 39.1% of the genes in drug catabolic process (GO:0042737) are present in the expression matrix (threshold: 50.0%). Only 10.5% of the genes in ATP synthesis coupled proton transport (GO:0015986) are present in the expression matrix (threshold: 50.0%). Only 16.7% of the genes in protection from natural killer cell mediated cytotoxicity (GO:0042270) are present in the expression matrix (threshold: 50.0%). Only 8.3% of the genes in energy coupled proton transport, down electrochemical gradient (GO:0015985) are present in the expression matrix (threshold: 50.0%). Only 20.0% of the genes in menaquinone metabolic process (GO:0009233) are present in the expression matrix (threshold: 50.0%). Only 31.6% of the genes in epoxygenase P450 pathway (GO:0019373) are present in the expression matrix (threshold: 50.0%). Only 5.9% of the genes in mitochondrial ATP synthesis coupled proton transport (GO:0042776) are present in the expression matrix (threshold: 50.0%). Only 33.3% of the genes in negative regulation of single stranded viral RNA replication via double stranded DNA intermediate (GO:0045869) are present in the expression matrix (threshold: 50.0%). Only 40.9% of the genes in exogenous drug catabolic process (GO:0042738) are present in the expression matrix (threshold: 50.0%). Only 40.0% of the genes in negative regulation of transcription by RNA polymerase I (GO:0016479) are present in the expression matrix (threshold: 50.0%). Only 26.8% of the genes in sensory perception of bitter taste (GO:0050913) are present in the expression matrix (threshold: 50.0%). Only 17.2% of the genes in sensory perception of chemical stimulus (GO:0007606) are present in the expression matrix (threshold: 50.0%). Only 9.9% of the genes in sensory perception of smell (GO:0007608) are present in the expression matrix (threshold: 50.0%). Only 33.3% of the genes in monoterpenoid metabolic process (GO:0016098) are present in the expression matrix (threshold: 50.0%). Only 28.6% of the genes in regulation of natural killer cell cytokine production (GO:0002727) are present in the expression matrix (threshold: 50.0%). Only 47.4% of the genes in cellular glucuronidation (GO:0052695) are present in the expression matrix (threshold: 50.0%). Only 46.7% of the genes in glucuronate metabolic process (GO:0019585) are present in the expression matrix (threshold: 50.0%). Only 20.0% of the genes in quinone catabolic process (GO:1901662) are present in the expression matrix (threshold: 50.0%). Only 44.4% of the genes in N-acylphosphatidylethanolamine metabolic process (GO:0070292) are present in the expression matrix (threshold: 50.0%). Only 47.4% of the genes in cellular response to copper ion (GO:0071280) are present in the expression matrix (threshold: 50.0%). Only 40.0% of the genes in cellular response to erythropoietin (GO:0036018) are present in the expression matrix (threshold: 50.0%). Only 40.0% of the genes in nucleosome positioning (GO:0016584) are present in the expression matrix (threshold: 50.0%). Only 42.9% of the genes in terpenoid metabolic process (GO:0006721) are present in the expression matrix (threshold: 50.0%). Only 28.6% of the genes in positive regulation of natural killer cell cytokine production (GO:0002729) are present in the expression matrix (threshold: 50.0%). Only 20.0% of the genes in regulation of barbed-end actin filament capping (GO:2000812) are present in the expression matrix (threshold: 50.0%). Only 40.0% of the genes in omega-hydroxylase P450 pathway (GO:0097267) are present in the expression matrix (threshold: 50.0%). Only 40.0% of the genes in positive regulation of neutrophil extravasation (GO:2000391) are present in the expression matrix (threshold: 50.0%). Only 38.9% of the genes in cellular response to zinc ion (GO:0071294) are present in the expression matrix (threshold: 50.0%). Only 33.3% of the genes in peptide antigen assembly with MHC protein complex (GO:0002501) are present in the expression matrix (threshold: 50.0%). Only 30.8% of the genes in negative regulation of chromatin organization (GO:1905268) are present in the expression matrix (threshold: 50.0%). Only 30.8% of the genes in negative regulation of chromatin silencing (GO:0031936) are present in the expression matrix (threshold: 50.0%). Only 37.5% of the genes in immune response-inhibiting cell surface receptor signaling pathway (GO:0002767) are present in the expression matrix (threshold: 50.0%). Only 33.3% of the genes in regulation of serine-type endopeptidase activity (GO:1900003) are present in the expression matrix (threshold: 50.0%). Only 38.9% of the genes in regulation of single stranded viral RNA replication via double stranded DNA intermediate (GO:0045091) are present in the expression matrix (threshold: 50.0%). Only 47.9% of the genes in chromatin assembly (GO:0031497) are present in the expression matrix (threshold: 50.0%). Only 41.2% of the genes in negative regulation of execution phase of apoptosis (GO:1900118) are present in the expression matrix (threshold: 50.0%). Only 42.9% of the genes in regulation of T cell tolerance induction (GO:0002664) are present in the expression matrix (threshold: 50.0%). Only 45.5% of the genes in xenobiotic glucuronidation (GO:0052697) are present in the expression matrix (threshold: 50.0%). Only 43.8% of the genes in regulation of chromatin silencing (GO:0031935) are present in the expression matrix (threshold: 50.0%). Only 16.7% of the genes in coumarin metabolic process (GO:0009804) are present in the expression matrix (threshold: 50.0%). Only 40.0% of the genes in zymogen inhibition (GO:0097341) are present in the expression matrix (threshold: 50.0%). Only 39.3% of the genes in cristae formation (GO:0042407) are present in the expression matrix (threshold: 50.0%).
Parallel processing completed! AUC calculation completed! Generated scores for 6036 pathways across 3696 cells.
We can use the anndata class to store the resulting gene set pathway and visualize it using the functions provided by the anndata class.
adata_aucs.obs=adata[adata_aucs.obs.index].obs.copy()
adata_aucs.obsm=adata[adata_aucs.obs.index].obsm.copy()
adata_aucs.obsp=adata[adata_aucs.obs.index].obsp.copy()
adata_aucs
AnnData object with n_obs × n_vars = 3696 × 6036
obs: 'clusters_coarse', 'clusters', 'S_score', 'G2M_score', 'response to vitamin (GO:0033273)_aucell'
obsm: 'X_pca', 'X_umap', 'UMAP'
obsp: 'distances', 'connectivities'
adata_aucs.write_h5ad('data/pancreas_auce.h5ad',compression='gzip')
We can load the AUCell result of scRNA-seq and visualize the geneset we test above
adata_aucs=ov.read('data/pancreas_auce.h5ad')
geneset_names=['response to vitamin (GO:0033273)','response to vitamin D (GO:0033280)']
ov.pl.embedding(
adata_aucs,
basis='X_umap',
color=geneset_names,
cmap='Reds',
)
Part4. Visualize differential enrichment pathways between different cell clusters.¶
We first read the AUCell scores files of all previously calculated signal pathways。
Given that the AUCell score is stored in the anndata structure, and is roughly equivalent to the level of gene expression, we can utilize the algorithm that calculates differential expression genes across clusters from Scanpy to determine clusters-specific signaling pathways, such as sc.tl.rank_genes_groups. By employing this approach, we can identify the pathways that are most distinctive to certain clusters and gain a better understanding of their underlying biology.
#adata_aucs.uns['log1p']['base']=None\
import scanpy as sc
sc.tl.rank_genes_groups(adata_aucs, 'clusters', method='t-test',n_genes=100)
ov.pl.rank_genes_groups_dotplot(adata_aucs,groupby='clusters',
cmap='Spectral_r',
standard_scale='var',n_genes=3)
ranking genes
finished: added to `.uns['rank_genes_groups']`
'names', sorted np.recarray to be indexed by group ids
'scores', sorted np.recarray to be indexed by group ids
'logfoldchanges', sorted np.recarray to be indexed by group ids
'pvals', sorted np.recarray to be indexed by group ids
'pvals_adj', sorted np.recarray to be indexed by group ids (0:00:01)
Furthermore, we can extract specific signaling pathways that are unique to particular clusters and then visualize them. By doing so, we can explore the differences between these clusters and gain insights into the molecular mechanisms that underlie these distinctions.
degs = sc.get.rank_genes_groups_df(adata_aucs, group='Beta', key='rank_genes_groups', log2fc_min=2,
pval_cutoff=0.05)['names'].squeeze()
degs
0 insulin metabolic process (GO:1901142) 1 negative regulation of organelle organization ... 2 amylin receptor signaling pathway (GO:0097647) 3 calcitonin family receptor signaling pathway (... 4 regulation of osteoclast differentiation (GO:0... Name: names, dtype: object
import matplotlib.pyplot as plt
#fig, axes = plt.subplots(4,3,figsize=(12,9))
axes=ov.pl.embedding(adata_aucs,ncols=3,
basis='X_umap',show=False,return_fig=True,wspace=0.35,hspace=0.65,
color=['clusters']+degs.values.tolist(),cmap='Reds',
title=[ov.utils.plot_text_set(i,3,20)for i in ['clusters']+degs.values.tolist()])
axes.tight_layout()
Part.5 Enrichment of geneset in scRNA-seq¶
In addition to using AUCell to assess gene sets, we can also use the cell's marker genes to find functions specific to each cell type. This is also commonly known as enrichment analysis.
adata.uns['log1p']['base']=None
sc.tl.rank_genes_groups(adata, 'clusters', method='t-test',n_genes=100)
ranking genes
finished: added to `.uns['rank_genes_groups']`
'names', sorted np.recarray to be indexed by group ids
'scores', sorted np.recarray to be indexed by group ids
'logfoldchanges', sorted np.recarray to be indexed by group ids
'pvals', sorted np.recarray to be indexed by group ids
'pvals_adj', sorted np.recarray to be indexed by group ids (0:00:00)
Additionally, we can calculate and visualize specific signaling pathways for all clusters with the function ov.single.pathway_enrichment.
res=ov.single.pathway_enrichment(adata,pathways_dict=pathway_dict,organism='Mouse',
group_by='clusters',plot=True)
To complete our analysis, we can use heatmaps to visualize the specific signaling pathways of each clusters with the function ov.single.pathway_enrichment_plot, with the depth of color reflecting the AUCell score for each respective clusters. This approach enables us to easily visualize and compare the activity of various signaling pathways across different clusters. By examining the heatmaps, we can identify patterns and differences between the clusters, providing us with a more in-depth understanding of the underlying biological mechanisms. Overall, this approach can be a powerful tool for identifying novel therapeutic targets and guiding the development of personalized treatments.
ax=ov.single.pathway_enrichment_plot(res,plot_title='Enrichment',cmap='Reds',
xticklabels=True,cbar=False,square=True,vmax=10,
yticklabels=True,cbar_kws={'label': '-log10(qvalue)','shrink': 0.5,})