Identifying biological processes of DRVI factors

Identifying biological processes of DRVI factors#

Factors that do not map cleanly to a single cell type often capture biological processes (interferon response, cell cycle, stress, …). This notebook annotates them with three complementary, non-LLM enrichment tools that all consume DRVI’s per-gene interpretability scores:

Tool

Method

Input

Strengths

Enrichr (via GSEApy)

Over-representation (ORA)

Top-N gene list

Fast; huge library collection

g:Profiler

Over-representation (ORA)

Ordered gene query

g:SCS correction for hierarchical GO terms

decoupler

Activity inference (ULM/MLM)

Gene-score matrix + prior network

Reports regulators (TFs), not just terms

It ends with an appendix of databases you can swap in. A separate factor-curation notebook then pulls every tool’s stored output (including the cell-type and LLM notebooks, if you ran them) into one per-factor view for final labelling.

Note: Enrichment usually returns broad, high-level processes (e.g. “defense response to virus”, “cell cycle”). For a finer reading — the specific pathway branch or the program driving a factor — align the factor’s top genes with the primary literature, or use the LLM-based tools to suggest candidate processes and then verify those suggestions against the literature. LLM output is a hypothesis, not evidence.

This is one of four companion notebooks. See cell types from annotations, LLM-based tools, and factor curation. They share embed.h5ad, so run the cell-type notebook first so its labels feed the curation step.

A note on pre-ranked GSEA. GSEA-style methods expect signed, two-sided rankings (top = up, bottom = down). DRVI interpretability scores are non-negative by construction (each factor-direction has its own positive ranking), so there is no meaningful “bottom” for GSEA to exploit. In practice pre-ranked GSEA returned very few, no-more-informative hits here, so we restrict this notebook to ORA (Enrichr, g:Profiler) plus TF-activity inference (decoupler).

All results are guiding — interpret in context and validate against known biology.

Prerequisites#

Assumes a trained DRVI model with interpretability scores (see the general pipeline).

Adapting to your own model — change io_dir (Section 0) and, per tool, the database/organism config: gseapy_db (Section 1), organism + gp_source (Section 2), and dc_geneset / dc_organism (Section 3). See the appendix for options.

Contact#

Questions: scverse discourse. Bugs: issue tracker.

Install#

This notebook uses the tutorials-biological-processes extra of drvi-py, which adds GSEApy, g:Profiler, decoupler, and statsmodels on top of DRVI. Install it once in your environment with:

pip install "drvi-py[tutorials-biological-processes]"

On Colab, the next cell does this for you. Remove it if your environment is already set up.

import sys
import subprocess

branch = "latest"
IN_COLAB = "google.colab" in sys.modules

if IN_COLAB and branch == "stable":
    subprocess.check_call([sys.executable, "-m", "pip", "install", "drvi-py[tutorials-biological-processes]"])
elif IN_COLAB and branch != "stable":
    subprocess.check_call([sys.executable, "-m", "pip", "install",
                           "git+https://github.com/theislab/drvi.git#egg=drvi-py[tutorials-biological-processes]"])

Imports#

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
import anndata as ad
import scanpy as sc
import matplotlib.pyplot as plt
from pathlib import Path

import scvi
import drvi
from drvi.model import DRVI

print("Last run with scvi-tools version:", scvi.__version__)
print("Last run with DRVI version:", drvi.__version__)

Last run with scvi-tools version: 1.4.3
Last run with DRVI version: 0.2.6

sc.set_figure_params(dpi=100, frameon=False, figsize=(3, 3))
plt.rcParams["figure.dpi"] = 100

0. Setup#

Config#

# Input/output directory holding the trained model and embeddings. Update accordingly.
io_dir = Path("./tmp_io/drvi_immune_128/").resolve()
print(f"Using directory: {io_dir}")

# DRVI provides two per-gene scoring approaches (both precomputed by the general pipeline):
#   "OOD_combined"            — decoder-based, favors specific genes (best for identity/top genes)
#   "IND_linear_weighted_mean"— per-cell mean effect, keeps broadly shared genes
# This is the shared scoring choice for every tool below.
score_key = "OOD_combined"

Using directory: /lustre/groups/ml01/code/amirali.moinfar/projects/drvi_tutorials/tmp_io/drvi_immune_128

Load model and embeddings#

adata = sc.read_h5ad(io_dir / "adata_preprocesses.h5ad")
model = DRVI.load(io_dir / "drvi_model", adata)

embed_path = io_dir / "embed.h5ad"
embed = sc.read_h5ad(embed_path)

# Gene x factor-direction score matrix used by all three tools.
scores_df = model.get_interpretability_scores(embed, adata, key=score_key)
gene_background = adata.var_names.tolist()

scores_df.iloc[:3, :3]

INFO     File                                                                                                      
         /lustre/groups/ml01/code/amirali.moinfar/projects/drvi_tutorials/tmp_io/drvi_immune_128/drvi_model/model.p
         t already downloaded
INFO     DRVI: The model is trained with DRVI version 0.2.6.
INFO     DRVI: Updaging data setup config ...
INFO     DRVI: Done updating data source registry. Loading in DRVI version 0.2.6.
INFO     DRVI: Loading model from DRVI version 0.2.6.
INFO     DRVI: Done updating model args. Loading in 0.2.6.
INFO     DRVI: The model has been initialized
title DR 1+ DR 1- DR 2+
index
TCL1A 0.0 0.009215 5.813187e-11
IGLL5 0.0 0.001222 0.000000e+00
PTGDS 0.0 0.001911 1.868569e-08

Shared helper#

Every tool starts from the same operation: for a factor-direction, take the top genes above a score cutoff.

def top_genes(scores_df, col, cutoff, top_n):
    """Genes for a factor-direction with score >= cutoff, at most `top_n`, ranked descending."""
    s = scores_df[col]
    return s[s >= cutoff].nlargest(top_n).index.astype(str).tolist()


def factor_first(df):
    """Return `df` with the `factor` column moved to the front (no-op if empty/absent)."""
    if df.empty or "factor" not in df.columns:
        return df
    return df[["factor", *df.columns.drop("factor")]]

1. Enrichr (via GSEApy)#

GSEApy’s Enrichr ORA tests whether a factor’s top-N gene list is enriched for terms in a chosen library, against a gene background.

import gseapy

Config#

gseapy_db = "GO_Biological_Process_2023"   # any Enrichr library (see appendix)
gseapy_cutoff = 0.1                         # ~0.1 for OOD, ~0.5 for IND
gseapy_top_n = 100
padj_threshold = 0.05                       # adjusted p-value

Run enrichment#

def run_gseapy_enrichr(scores_df, gene_sets, cutoff, top_n, padj_cutoff, background):
    rows = []
    for col in scores_df.columns:
        genes = top_genes(scores_df, col, cutoff, top_n)
        if not genes:
            continue
        try:
            enr = gseapy.enrich(gene_list=genes, gene_sets=gene_sets, background=background,
                                no_plot=True, outdir=None)
        except Exception as e:
            print(f"ORA failed for {col}: {e}")
            continue
        hits = enr.results[enr.results["Adjusted P-value"] < padj_cutoff].sort_values("Adjusted P-value")
        rows.append(hits.assign(factor=col))
    return factor_first(pd.concat(rows, ignore_index=True)) if rows else pd.DataFrame()


gseapy_enrichr_results = run_gseapy_enrichr(
    scores_df, gseapy_db, gseapy_cutoff, gseapy_top_n, padj_threshold, gene_background
)
display(gseapy_enrichr_results.head())
factor Gene_set Term P-value Adjusted P-value Old P-value Old adjusted P-value Odds Ratio Combined Score Genes
0 DR 1- GO_Biological_Process_2023 T Cell Activation (GO:0042110) 4.001337e-07 0.000340 0 0 8.500000 125.217470 DPP4;CD2;ITK;LCK;CD7;CD28;RHOH;CCR7;CD3G;CD247...
1 DR 1- GO_Biological_Process_2023 T Cell Receptor Signaling Pathway (GO:0050852) 4.027758e-06 0.001710 0 0 8.333333 103.519172 ITK;LCK;TXK;THEMIS;CD28;CD3G;CD247;CD3E;CD3D;S...
2 DR 1- GO_Biological_Process_2023 Gamma-Delta T Cell Activation (GO:0046629) 2.836112e-05 0.008026 0 0 79.125000 828.477629 ITK;TCF7;CD3G;CD3E
3 DR 1- GO_Biological_Process_2023 Antigen Receptor-Mediated Signaling Pathway (G... 4.462457e-05 0.009472 0 0 5.337601 53.467957 ITK;PLEKHA1;LCK;TXK;CD28;THEMIS;CD3G;CD247;CD3...
4 DR 1- GO_Biological_Process_2023 Alpha-Beta T Cell Activation (GO:0046631) 1.835750e-04 0.031171 0 0 26.347222 226.662180 CD3G;CD247;CD3E;CD3D 
# Let's look at the top enriched term for each factor-direction.
display(gseapy_enrichr_results.drop_duplicates(subset=["factor"], keep="first"))
factor Gene_set Term P-value Adjusted P-value Old P-value Old adjusted P-value Odds Ratio Combined Score Genes
0 DR 1- GO_Biological_Process_2023 T Cell Activation (GO:0042110) 4.001337e-07 3.397135e-04 0 0 8.500000 125.217470 DPP4;CD2;ITK;LCK;CD7;CD28;RHOH;CCR7;CD3G;CD247...
6 DR 2- GO_Biological_Process_2023 Defense Response To Fungus (GO:0050832) 8.511883e-07 8.571467e-04 0 0 23.759857 332.082782 CLEC4D;CLEC7A;S100A12;NRG1;CLEC4E;S100A9;S100A8
14 DR 3+ GO_Biological_Process_2023 B Cell Activation (GO:0042113) 8.097060e-09 6.558619e-06 0 0 13.500000 251.528825 FCRLA;TPD52;LAT2;CD79B;CD79A;CD40;MEF2C;BANK1;...
33 DR 5+ GO_Biological_Process_2023 Regulation Of Natural Killer Cell Mediated Cyt... 2.918578e-05 2.638394e-02 0 0 15.095745 157.627184 NCR1;NCR3;KLRB1;CD160;KLRD1;KLRC1
43 DR 6+ GO_Biological_Process_2023 Phagocytosis (GO:0006909) 6.229536e-07 6.777735e-04 0 0 13.323391 190.375185 MARCO;TYROBP;CD93;LRP1;SLC11A1;CD36;CD302;CD14...
53 DR 7- GO_Biological_Process_2023 Alpha-Beta T Cell Activation (GO:0046631) 4.715882e-06 4.470656e-03 0 0 51.657609 633.558595 NKG7;CD3G;CD247;CD3E;CD3D
54 DR 8- GO_Biological_Process_2023 Alpha-Beta T Cell Activation (GO:0046631) 3.235574e-06 2.785829e-03 0 0 56.117647 709.400250 NKG7;CD3G;CD247;CD3E;CD3D
63 DR 9- GO_Biological_Process_2023 T Cell Activation (GO:0042110) 4.285637e-08 3.599935e-05 0 0 9.640507 163.555174 ITK;CRTAM;RHOH;CD3G;CD3E;CD3D;CD2;ZAP70;CD8B;C...
68 DR 12- GO_Biological_Process_2023 Regulation Of Cell Death (GO:0010941) 4.976684e-05 1.522865e-02 0 0 65.866667 652.617585 BNIP3L;HBA1;SNCA
69 DR 13+ GO_Biological_Process_2023 B Cell Activation (GO:0042113) 8.097060e-09 7.157801e-06 0 0 13.500000 251.528825 FCRLA;CD79B;CD79A;CD40;PRKCB;BANK1;CD70;GPR183...
88 DR 14+ GO_Biological_Process_2023 T Cell Activation (GO:0042110) 4.001337e-07 3.961324e-04 0 0 8.500000 125.217470 CD2;CD84;CD8B;CD8A;LCK;CRTAM;SLAMF7;CD3G;SLAMF...
103 DR 15- GO_Biological_Process_2023 Positive Regulation Of Chromosome Organization... 1.862988e-05 5.477186e-03 0 0 inf inf TAL1;CDK1;SMC4
117 DR 17+ GO_Biological_Process_2023 B Cell Activation (GO:0042113) 5.481649e-10 4.571695e-07 0 0 15.623244 333.157007 CD86;TPD52;CD40;MEF2C;PRKCB;CD70;FCRLA;CD79B;C...
139 DR 19- GO_Biological_Process_2023 Cell Cycle G1/S Phase Transition (GO:0044843) 7.946512e-06 5.530773e-03 0 0 74.692308 877.095144 CCND2;CDK6;MYC;EIF4EBP1
141 DR 22- GO_Biological_Process_2023 Microtubule Cytoskeleton Organization Involved... 1.407260e-05 1.419925e-02 0 0 33.280702 371.788077 CENPE;CCNB1;NUSAP1;CDK1;BIRC5
148 DR 24- GO_Biological_Process_2023 Regulation Of Natural Killer Cell Mediated Cyt... 8.364444e-09 6.214782e-06 0 0 37.043478 688.981877 NCR1;NCR3;KLRB1;CD160;KLRC4;CRTAM;KLRD1;KLRC1
158 DR 25+ GO_Biological_Process_2023 T Cell Activation (GO:0042110) 4.001337e-07 3.673227e-04 0 0 8.500000 125.217470 CD2;CASP8;LCK;CD28;CD3G;GATA3;CD247;CD3E;FOXP3...
163 DR 26+ GO_Biological_Process_2023 ERAD Pathway (GO:0036503) 5.494890e-06 4.242055e-03 0 0 49.947368 604.947144 ERLEC1;XBP1;DERL3;HERPUD1;UBE2J1
168 DR 27- GO_Biological_Process_2023 Immune Response-Regulating Cell Surface Recept... 5.494890e-06 6.753219e-03 0 0 49.947368 604.947144 LYN;FGR;LILRB1;FPR2;LILRB2
177 DR 28- GO_Biological_Process_2023 Antigen Processing And Presentation Of Exogeno... 1.027028e-13 1.089677e-10 0 0 68.710843 2054.930836 CLEC4A;HLA-DMA;HLA-DRB5;HLA-DMB;FCER1G;HLA-DRA...
226 DR 29+ GO_Biological_Process_2023 Positive Regulation Of Defense Response (GO:00... 7.019148e-07 7.756158e-04 0 0 7.130268 101.032006 GRN;OSM;CD1D;IL17RA;EREG;PYCARD;NLRP12;LGALS1;...
237 DR 30- GO_Biological_Process_2023 Regulation Of Myeloid Cell Differentiation (GO... 1.407260e-05 1.577538e-02 0 0 33.280702 371.788077 MEF2C;MEIS1;CDK6;FES;ZBTB16
238 DR 31+ GO_Biological_Process_2023 B Cell Differentiation (GO:0030183) 9.371622e-05 1.557190e-02 0 0 22.920821 212.596100 CD79B;CD79A;BLNK;PAX5
241 DR 32- GO_Biological_Process_2023 Antigen Processing And Presentation Of Exogeno... 2.777335e-15 3.213377e-12 0 0 94.486590 3166.933923 HLA-DRB5;FCER1G;IFI30;CTSS;CLEC4A;HLA-DMA;HLA-...
269 DR 34+ GO_Biological_Process_2023 Defense Response To Bacterium (GO:0042742) 5.730761e-09 3.971418e-06 0 0 13.070938 248.052651 PYCARD;STAB1;HMGB2;RNASE6;AZU1;LYZ;RNASE2;S100...
275 DR 35+ GO_Biological_Process_2023 Defense Response To Gram-negative Bacterium (G... 3.352359e-07 6.067770e-05 0 0 54.666667 814.994246 CTSG;AZU1;LYZ;RNASE2;ELANE
281 DR 37+ GO_Biological_Process_2023 T Cell Receptor Signaling Pathway (GO:0050852) 4.296314e-07 1.903267e-04 0 0 15.876543 232.755494 ZAP70;LCK;CD28;CTLA4;CD3G;FYN;CD3D;SKAP1
289 DR 38+ GO_Biological_Process_2023 B Cell Receptor Signaling Pathway (GO:0050853) 3.484871e-05 1.700617e-02 0 0 12.812431 131.513124 BLK;CD79B;LAT2;CD79A;CD19;CD38
290 DR 40- GO_Biological_Process_2023 Inflammatory Response (GO:0006954) 1.935338e-10 2.113389e-07 0 0 7.625791 170.555158 CCR1;CEBPB;TNFAIP6;CCL20;NFAM1;SLC11A1;C5AR1;F...
291 DR 41+ GO_Biological_Process_2023 Regulation Of Natural Killer Cell Mediated Cyt... 1.632985e-06 1.399468e-03 0 0 20.354839 271.230282 NCR1;NCR3;KLRB1;KLRC3;CD226;KLRD1;KLRC1
295 DR 43- GO_Biological_Process_2023 Defense Response To Virus (GO:0051607) 9.184993e-28 7.173480e-25 0 0 37.678106 2345.643372 IFITM3;IFITM1;IFI6;IFIT1;DDX60;IFIT3;IFI44L;IF...
329 DR 45+ GO_Biological_Process_2023 Microtubule Cytoskeleton Organization Involved... 1.163970e-07 3.087630e-05 0 0 95.980392 1532.447821 CENPE;STMN1;CDK1;NUSAP1;BIRC5
378 DR 46+ GO_Biological_Process_2023 Antigen Processing And Presentation Of Exogeno... 2.777335e-15 3.010631e-12 0 0 94.486590 3166.933923 HLA-DRB5;FCER1G;IFI30;CLEC4A;HLA-DMA;HLA-DMB;H...
402 DR 47- GO_Biological_Process_2023 MHC Class II Protein Complex Assembly (GO:0002... 6.271388e-14 2.429118e-11 0 0 inf inf GNAO1;HLA-DRB5;HLA-DMA;HLA-DMB;HLA-DRA;HLA-DOB...
417 DR 48- GO_Biological_Process_2023 Microtubule Cytoskeleton Organization Involved... 2.334885e-11 5.930608e-09 0 0 370.500000 9070.015400 CENPE;CCNB1;STMN1;CDK1;NUSAP1;BIRC5
508 DR 49+ GO_Biological_Process_2023 Negative Regulation Of Cell Population Prolife... 2.431181e-06 2.923566e-03 0 0 5.320154 68.774345 APP;IGFBP3;ITGA1;PTPRK;IL6;ADAMTS1;CYP1B1;ZNF5...
513 DR 50- GO_Biological_Process_2023 MHC Class II Protein Complex Assembly (GO:0002... 2.125042e-07 3.520487e-05 0 0 70.071429 1076.598731 HLA-DMB;HLA-DRA;HLA-DQA1;HLA-DRB1;HLA-DPA1
551 DR 51- GO_Biological_Process_2023 Antigen Processing And Presentation Of Exogeno... 5.008729e-06 2.188732e-03 0 0 15.814815 193.009192 CLEC4A;HLA-DMA;HLA-DRB5;HLA-DMB;FCER1G;HLA-DRA...
568 DR 53+ GO_Biological_Process_2023 B Cell Activation (GO:0042113) 8.951769e-05 4.060559e-02 0 0 7.096408 66.146149 ITGB1;TPD52;CD79A;CD40;MEF2C;CD70;BTK;LAX1
572 DR 55- GO_Biological_Process_2023 Response To Zinc Ion (GO:0010043) 2.284626e-06 2.245787e-03 0 0 60.569620 786.757469 MT2A;MT1F;MT1X;CRIP1;MT1E

Store results#

embed.uns["gseapy_enrichr_results"] = gseapy_enrichr_results.convert_dtypes(
    convert_integer=False, convert_floating=False
)
n_sig = gseapy_enrichr_results["factor"].nunique() if not gseapy_enrichr_results.empty else 0
print(f"Enrichr significant factor-directions: {n_sig} / {scores_df.shape[1]}")

Enrichr significant factor-directions: 40 / 111

How to read this. ORA works well when the factor’s biology is well represented in the chosen library: hits for a lineage-aligned factor are usually specific and internally consistent (e.g. a B-cell factor returning Ig-mediated immune response, BCR signaling). Quality is bounded by the database, so factors whose biology is under-represented return loosely related terms. Interpret alongside g:Profiler below — convergent terms across the two ORA tools are better supported.

2. g:Profiler#

g:Profiler runs ORA with g:SCS multiple-testing correction, designed for hierarchical GO terms. In ordered-query mode it walks the ranked gene list to find the best-enriched prefix, which suits continuous factor scores better than a fixed top-N.

from gprofiler import GProfiler

Config#

organism = "hsapiens"        # e.g. "mmusculus", "drerio"
gp_source = ["GO:BP"]        # e.g. ["GO:MF"], ["REAC"], ["KEGG"], ["HP"]
gp_cutoff = 0.1              # ~0.1 for OOD, ~0.5 for IND
gp_top_n = 100               # max top genes per factor-direction
pval_threshold = 0.05

Run enrichment#

def run_gprofiler(scores_df, background, organism, sources, pval_threshold, cutoff, top_n):
    gp = GProfiler(return_dataframe=True)
    rows = []
    for col in scores_df.columns:
        genes = top_genes(scores_df, col, cutoff, top_n)
        if not genes:
            continue
        res = gp.profile(organism=organism, query=genes, sources=sources, ordered=True,
                         user_threshold=pval_threshold, background=background)
        if res is None or res.empty:
            continue
        rows.append(res.assign(factor=col))
    return factor_first(pd.concat(rows, ignore_index=True)) if rows else pd.DataFrame()


gprofiler_results = run_gprofiler(
    scores_df, gene_background, organism, gp_source, pval_threshold, gp_cutoff, gp_top_n
)
if not gprofiler_results.empty:
    gprofiler_results["parents"] = gprofiler_results["parents"].astype(str)

n_sig = gprofiler_results["factor"].nunique() if not gprofiler_results.empty else 0
print(f"g:Profiler significant factor-directions: {n_sig} / {scores_df.shape[1]}")
gprofiler_results.sort_values(["factor", "p_value"]).head() if not gprofiler_results.empty else gprofiler_results

g:Profiler significant factor-directions: 14 / 111
factor source native name p_value significant description term_size query_size intersection_size effective_domain_size precision recall query parents
2 DR 13+ GO:BP GO:0050853 B cell receptor signaling pathway 0.001592 True "The series of molecular signals initiated by ... 38 30 10 1977 0.333333 0.263158 query_1 ['GO:0050851']
3 DR 17+ GO:BP GO:0002399 MHC class II protein complex assembly 0.010536 True "The aggregation, arrangement and bonding toge... 10 81 7 1977 0.086420 0.700000 query_1 ['GO:0002396']
4 DR 17+ GO:BP GO:0002503 peptide antigen assembly with MHC class II pro... 0.010536 True "The binding of a peptide to the antigen bindi... 10 81 7 1977 0.086420 0.700000 query_1 ['GO:0002399', 'GO:0002495', 'GO:0002501']
5 DR 24- GO:BP GO:0002715 regulation of natural killer cell mediated imm... 0.004125 True "Any process that modulates the frequency, rat... 24 16 7 1977 0.437500 0.291667 query_1 ['GO:0002228', 'GO:0002706', 'GO:0045088']
6 DR 24- GO:BP GO:0002228 natural killer cell mediated immunity 0.004492 True "The promotion of an immune response by natura... 36 16 8 1977 0.500000 0.222222 query_1 ['GO:0002449', 'GO:0045087'] 
# Let's look at the top enriched term for each factor-direction.
if not gprofiler_results.empty:
    display(gprofiler_results.drop_duplicates(subset=["factor"], keep="first"))
factor source native name p_value significant description term_size query_size intersection_size effective_domain_size precision recall query parents
0 DR 5+ GO:BP GO:0140507 granzyme-mediated programmed cell death signal... 1.763272e-02 True "The series of molecular signals induced by gr... 6 16 4 1977 0.250000 0.666667 query_1 ['GO:0007165', 'GO:0012501']
1 DR 8- GO:BP GO:0140507 granzyme-mediated programmed cell death signal... 2.952562e-02 True "The series of molecular signals induced by gr... 6 18 4 1977 0.222222 0.666667 query_1 ['GO:0007165', 'GO:0012501']
2 DR 13+ GO:BP GO:0050853 B cell receptor signaling pathway 1.591564e-03 True "The series of molecular signals initiated by ... 38 30 10 1977 0.333333 0.263158 query_1 ['GO:0050851']
3 DR 17+ GO:BP GO:0002399 MHC class II protein complex assembly 1.053554e-02 True "The aggregation, arrangement and bonding toge... 10 81 7 1977 0.086420 0.700000 query_1 ['GO:0002396']
5 DR 24- GO:BP GO:0002715 regulation of natural killer cell mediated imm... 4.125462e-03 True "Any process that modulates the frequency, rat... 24 16 7 1977 0.437500 0.291667 query_1 ['GO:0002228', 'GO:0002706', 'GO:0045088']
10 DR 28- GO:BP GO:0019884 antigen processing and presentation of exogeno... 2.777675e-09 True "The process in which an antigen-presenting ce... 22 26 12 1977 0.461538 0.545455 query_1 ['GO:0019882']
21 DR 32- GO:BP GO:0002504 antigen processing and presentation of peptide... 1.492052e-07 True "The process in which an antigen-presenting ce... 20 34 11 1977 0.323529 0.550000 query_1 ['GO:0019882']
33 DR 35+ GO:BP GO:0070944 neutrophil-mediated killing of bacterium 4.431323e-02 True "The directed killing of a bacterium by a neut... 4 4 3 1977 0.750000 0.750000 query_1 ['GO:0042742', 'GO:0070943']
34 DR 43- GO:BP GO:0045071 negative regulation of viral genome replication 1.718009e-07 True "Any process that stops, prevents, or reduces ... 22 32 12 1977 0.375000 0.545455 query_1 ['GO:0019079', 'GO:0045069', 'GO:0048525']
45 DR 46+ GO:BP GO:0002504 antigen processing and presentation of peptide... 3.953655e-08 True "The process in which an antigen-presenting ce... 20 16 9 1977 0.562500 0.450000 query_1 ['GO:0019882']
56 DR 47- GO:BP GO:0002399 MHC class II protein complex assembly 6.278540e-08 True "The aggregation, arrangement and bonding toge... 10 30 8 1977 0.266667 0.800000 query_1 ['GO:0002396']
67 DR 48- GO:BP GO:0098813 nuclear chromosome segregation 9.918886e-05 True "The process in which genetic material, in the... 28 11 11 1977 1.000000 0.392857 query_1 ['GO:0007059']
73 DR 50- GO:BP GO:0002504 antigen processing and presentation of peptide... 9.255973e-03 True "The process in which an antigen-presenting ce... 20 26 7 1977 0.269231 0.350000 query_1 ['GO:0019882']
77 DR 51- GO:BP GO:0002495 antigen processing and presentation of peptide... 4.041777e-02 True "The process in which an antigen-presenting ce... 18 92 10 1977 0.108696 0.555556 query_1 ['GO:0002504', 'GO:0048002']

Store results#

embed.uns["gprofiler_results"] = gprofiler_results.convert_dtypes(
    convert_integer=False, convert_floating=False
)

How to read this. Top hits for a factor often organize general-to-specific, and g:Profiler can surface more mechanistic terms than top-N ORA (e.g. UPR/ERAD components on a plasma-cell factor where Enrichr returns generic B-cell terms). Coverage is typically lower than Enrichr; read the two together and trust convergent calls more than single-tool ones.

3. decoupler#

decoupler infers the activity of regulators from gene-level scores against a prior-knowledge network, rather than testing gene-set overlap. A significant hit is a transcription factor whose targets are co-enriched in the factor’s top genes — a mechanistic layer the ORA/LLM tools cannot provide. Curated networks from OmniPath:

  • CollecTRI — comprehensive TF → target regulons (recommended for TF drivers).

  • DoRothEA — TF regulons with confidence tiers A–D (tunable stringency).

  • PROGENy — pathway footprints (Hypoxia, EGFR, TGFb, …); exploratory for signaling.

import decoupler as dc
from statsmodels.stats.multitest import multipletests

Config#

dc_geneset = "collectri"        # or "dorothea", "progeny"
dc_organism = "human"           # match the data species
dc_cutoff = 0.01                # ~0.01 for OOD, ~0.05 for IND; scores below this are zeroed
fdr_threshold = 0.05

dc_methods = ["ulm", "mlm"]
dc_min = 10                     # min genes of a regulon present in the data for a valid test
dorothea_levels = ["A", "B", "C"]
fdr_method = "fdr_bh"

Load the regulatory network#

net_dispatch = {
    "collectri": lambda: dc.op.collectri(organism=dc_organism),
    "dorothea":  lambda: dc.op.dorothea(organism=dc_organism, levels=dorothea_levels),
    "progeny":   lambda: dc.op.progeny(organism=dc_organism),
}
net = net_dispatch.get(
    dc_geneset.strip().lower(),
    lambda: dc.op.resource(name=dc_geneset, organism=dc_organism),
)()

cols = ["source", "target"] + (["weight"] if "weight" in net.columns else [])
net = net[cols].dropna().drop_duplicates().reset_index(drop=True)
print(f"Network: {len(net)} interactions, {net['source'].nunique()} regulators")

Network: 42990 interactions, 1185 regulators

Run#

def run_decouple(factors_by_genes, net, methods, tmin, fdr_method):
    mat = factors_by_genes.copy()
    mat.columns = mat.columns.astype(str).str.upper()

    net_u = net.copy()
    net_u["target"] = net_u["target"].astype(str).str.upper()

    keep = mat.columns.intersection(net_u["target"].unique())
    mat = mat[keep].replace([np.inf, -np.inf], 0.0).fillna(0.0)

    res = dc.mt.decouple(data=mat, net=net_u, methods=methods, cons=False, tmin=tmin, verbose=True)
    _, pvals = dc.mt.consensus(res)

    out = pvals.stack().reset_index(name="p_value").rename(columns={"level_0": "factor", "level_1": "term"})
    out["p_adj"] = multipletests(out["p_value"].values, method=fdr_method)[1]
    return out


input_df = scores_df.copy()
input_df[input_df < dc_cutoff] = 0
decoupler_all = run_decouple(input_df.T, net, dc_methods, dc_min, fdr_method)

# Keep the most significant regulator per factor-direction for a summary view.
decoupler_results = (
    decoupler_all[decoupler_all["p_adj"] < fdr_threshold]
    .sort_values("p_adj")
    .groupby("factor", as_index=False)
    .first()
)
print(f"Significant regulators for {decoupler_results['factor'].nunique()} / {scores_df.shape[1]} directions.")
display(decoupler_results.sort_values("p_adj"))

Significant regulators for 14 / 111 directions.
factor term p_value p_adj
12 DR 7- RORC 5.442043e-20 7.722259e-16
0 DR 18- RORC 1.077976e-15 5.098825e-12
7 DR 42+ YY1 8.776332e-16 5.098825e-12
1 DR 21- ESR2 1.078305e-13 3.825288e-10
9 DR 51- NR0B2 5.085469e-08 1.030897e-04
6 DR 38+ EBF1 4.649327e-08 1.030897e-04
3 DR 28- NR1H3 6.014425e-08 1.066809e-04
10 DR 52- BMAL1 2.103195e-07 3.316037e-04
11 DR 53+ PRDM1 6.800543e-06 9.649970e-03
4 DR 35+ RUNX2 7.498029e-06 9.672457e-03
8 DR 43- IRF9 1.374995e-05 1.393656e-02
13 DR 9- SATB1 1.891204e-05 1.789079e-02
5 DR 36- RUNX2 2.597060e-05 2.167781e-02
2 DR 22- NKX3-1 2.497380e-05 2.167781e-02

Store results#

embed.uns["decoupler_results"] = decoupler_results.convert_dtypes(
    convert_integer=False, convert_floating=False
)

How to read this. Where decoupler reports a significant TF, it often lines up with the other tools (e.g. EBF1 on a B-cell-precursor factor, PRDM1/Blimp-1 on a plasma-cell factor, RFXAP alongside “MHC class II assembly”). Coverage is limited by design — many factors are not dominated by a single TF — and some hits are not obviously related to the factor’s identity, so treat regulators as candidate drivers to cross-check, not conclusions.

4. Manual Exploration#

Not every factor is a cell type. Some capture a process shared across many different cell types. DR 43− is a good example: an antiviral / interferon response. Its top genes are canonical interferon-stimulated genes (IFIT1/2/3, ISG15, MX1, CXCL10); all three tools converge on virus/interferon terms and point to the interferon regulator IRF9.

process_factor = "DR 43-"   # a program shared across cell types (edit to explore another factor)


def show_factor_evidence(factor_label, cutoff=0.1, n_genes=12):
    genes = top_genes(scores_df, factor_label, cutoff, n_genes)
    print(f"{factor_label} — top genes:\n  {', '.join(genes)}\n")
    for name, df, col in [
        ("Enrichr", gseapy_enrichr_results, "Term"),
        ("g:Profiler", gprofiler_results, "name"),
        ("decoupler TF", decoupler_results, "term"),
    ]:
        terms = df.loc[df["factor"] == factor_label, col].head(5).tolist() if not df.empty else []
        print(f"{name}: {terms or '(no hit)'}")


show_factor_evidence(process_factor)

DR 43- — top genes:
  CXCL10, CCL8, IFIT1, IFIT3, IFIT2, ISG15, RSAD2, IFI44L, APOBEC3A, HERC5, IFI6, MX1

Enrichr: ['Defense Response To Virus (GO:0051607)', 'Defense Response To Symbiont (GO:0140546)', 'Negative Regulation Of Viral Process (GO:0048525)', 'Negative Regulation Of Viral Genome Replication (GO:0045071)', 'Regulation Of Viral Genome Replication (GO:0045069)']
g:Profiler: ['negative regulation of viral genome replication', 'defense response to virus', 'response to virus', 'negative regulation of viral process', 'viral genome replication']
decoupler TF: ['IRF9']

5. Save#

Persist the enrichment results so the curation view (and the other notebooks) can read them.

ad.settings.allow_write_nullable_strings = True
embed.write_h5ad(embed_path)
print(f"Updated embedding saved to: {embed_path}")

Updated embedding saved to: /lustre/groups/ml01/code/amirali.moinfar/projects/drvi_tutorials/tmp_io/drvi_immune_128/embed.h5ad

6. Appendix: database reference#

Swap the tool-specific config variables above (gseapy_db, gp_source, dc_geneset) to use different databases. The small tables below list common choices.

For a much larger, curated catalog, this repo bundles resources/gene_set_libraries_master_v0_8_FINAL.xlsx — a registry of gene-set libraries with, per library, its authoritative source, license, curation-quality rating, and flags for which tool can consume it (usable_in_gseapy / usable_in_gprofiler / usable_in_decoupler). It also includes recommended panels for this exact workflow (e.g. Core_Biology_DRVI, Cell_Type_Annotation, Regulatory_Activity). Load it with pd.read_excel(..., sheet_name="libraries_master") and filter to the tool you are using.

Biological process databases#

Database

Description

Enrichr library

g:Profiler

MSigDB Hallmark

50 curated, non-redundant states

MSigDB_Hallmark_2020

—

GO Biological Process

Hierarchical processes

GO_Biological_Process_2025

GO:BP

GO Cellular Component

Subcellular localization

GO_Cellular_Component_2025

GO:CC

GO Molecular Function

Molecular activities

GO_Molecular_Function_2025

GO:MF

Reactome

Reaction-based pathways

Reactome_Pathways_2024

REAC

KEGG

Metabolic/signaling maps

KEGG_2026

KEGG

WikiPathways

Community-curated maps

WikiPathways_2024_Human

WP

Cell-type marker databases (Enrichr)#

Database

Description

Enrichr library

CellMarker 2.0

Curated markers from literature

CellMarker_2024

PanglaoDB

Curated scRNA-seq markers

PanglaoDB_Augmented_2021

Tabula Sapiens

Human single-cell atlas markers

Tabula_Sapiens

Human Gene Atlas

Tissue/cell-type expression

Human_Gene_Atlas

Regulatory networks (decoupler)#

Network

Description

decoupler name

Notes

CollecTRI

TF → target regulons

collectri

Recommended for TF drivers

DoRothEA

TF → target, confidence A–D

dorothea

Tunable stringency

PROGENy

Pathway-responsive signatures

progeny

Best for signaling activity

Clinical / disease phenotypes#

Database

Description

Enrichr

g:Profiler

Human Phenotype Ontology

Genes linked to clinical signs

Human_Phenotype_Ontology

HP

OMIM

Human genes and genetic disorders

OMIM_Disease

OMIM