4. Gene Signature Scoring#

Pipeline phase: Phase 3.5b

What you will learn:

  • Load the bundled MSigDB Hallmark gene sets (offline)

  • Merge project-specific signatures with Hallmark

  • Score signatures using score_signature (scanpy, UCell, ssGSEA methods)

  • Read back scores from adata.obsm

  • Run over-representation analysis (ORA)

  • Plot enrichment results with plot_gsea_dotplot

Storage convention: Scores are written to adata.obsm['signature_score'] (raw) and adata.obsm['signature_score_z'] (z-scored). Column names are full paths, e.g. Hallmark/HYPOXIA, Myeloid/Macrophage_Core.

import numpy as np
import pandas as pd
import anndata as ad
import sc_tools.tl as tl
import sc_tools.pl as pl
from sc_tools.utils.signatures import get_signature_columns, get_signature_df

Synthetic data with known gene sets#

np.random.seed(1)
n_obs, n_genes = 200, 500

# Include some real gene symbols so Hallmark can match them
hallmark_sample_genes = [
    "ALDOA", "CDKN3", "ENO1", "LDHA", "PGK1",  # HYPOXIA
    "VEGFA", "HK1", "HK2", "PFKL", "TPI1",
    "MYC", "CCND1", "CDK4", "E2F1", "PCNA",    # MYC_TARGETS_V1
    "TP53", "CDKN1A", "BBC3", "PUMA", "BAX",   # P53_PATHWAY
]
filler_genes = [f"GENE_{i}" for i in range(n_genes - len(hallmark_sample_genes))]
gene_names = hallmark_sample_genes + filler_genes

X = np.random.rand(n_obs, n_genes).astype(np.float32)

adata = ad.AnnData(X=X)
adata.var_names = gene_names
adata.obs["leiden"] = pd.Categorical([str(i % 4) for i in range(n_obs)])
adata.obs["condition"] = pd.Categorical(["tumor" if i % 3 == 0 else "normal" for i in range(n_obs)])

print(adata)

Load gene sets#

Bundled MSigDB Hallmark (offline)#

hallmark = tl.load_hallmark()
print(f"Loaded {len(hallmark)} Hallmark gene sets")
print("Example (HYPOXIA):", list(hallmark["Hallmark/HYPOXIA"])[:8])

Project-specific signatures#

Project signatures are typically stored in metadata/gene_signatures.json. Here we define a minimal example inline:

project_sigs = {
    "Myeloid/Macrophage_Core": ["MRC1", "CD68", "CSF1R", "ITGAM", "CD14"],
    "Myeloid/M2_Polarization": ["MRC1", "CD163", "ARG1", "IL10", "TGFB1"],
}

# Validate and merge with Hallmark
report = tl.validate_gene_signatures(project_sigs)
print("Validation report:", report)

combined = tl.merge_gene_signatures(project_sigs, hallmark)
print(f"Combined: {len(combined)} gene sets")

Score signatures#

The default method is "scanpy" (scanpy tl.score_genes). Set method="ucell" or method="ssgsea" for rank-based alternatives (requires pip install sc-tools[geneset]).

tl.score_signature(adata, combined, method="scanpy")

print("Scoring method recorded:", adata.uns.get("scoring_method"))
print("obsm keys:", [k for k in adata.obsm if "signature" in k])

Retrieve scores#

# List all scored signatures
cols = get_signature_columns(adata)
print(f"Scored {len(cols)} signatures")
print("First 5:", cols[:5])

# Get as DataFrame (spots x signatures)
score_df = get_signature_df(adata)
print("\nScore DataFrame shape:", score_df.shape)
print(score_df.iloc[:3, :3])

Over-representation analysis (ORA)#

run_ora performs Fisher exact test with Benjamini-Hochberg FDR correction for a set of query genes against the loaded gene sets.

# Query: genes upregulated in cluster 0 (here, synthetic random subset)
query_genes = ["ALDOA", "LDHA", "PGK1", "HK1", "VEGFA", "ENO1", "TPI1"]

ora_results = tl.run_ora(
    query_genes=query_genes,
    gene_sets=hallmark,
    background=adata.var_names.tolist(),
)

print(ora_results[["Term", "pval", "fdr", "overlap"]].head(5).to_string())

Plot ORA results#

import matplotlib.pyplot as plt

fig = pl.plot_gsea_dotplot(
    ora_results,
    fdr_col="fdr",
    score_col="overlap",
    top_n=10,
    title="ORA: cluster 0 upregulated genes",
)
plt.tight_layout()
plt.show()

Save Phase 3.5b checkpoint#

adata.write_h5ad("results/adata.normalized.scored.p35.h5ad")

The checkpoint must have:

  • obsm['signature_score'] — raw scores (spots x signatures)

  • obsm['signature_score_z'] — z-scored scores

  • uns['signature_score_report'] — per-signature n_present / n_missing

  • uns['scoring_method'] — the method used

Downstream scripts read scores via get_signature_df(adata) / get_signature_columns(adata).