sc_tools.tl — Analysis Tools#

Analysis tools for spatial omics data following the scanpy API pattern.

import sc_tools.tl as tl

hallmark = tl.load_hallmark()
combined = tl.merge_gene_signatures(project_sigs, hallmark)
tl.score_signature(adata, combined)

Gene Signature Scoring#

Scores are stored in 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.

sc_tools.tl.score_signature(adata, signatures_nested, method='scanpy', use_raw=True, ctrl_size=50, n_bins=25, min_genes=3, copy=False, save_path=None)[source]#

Score gene signatures and store in obsm (single flat key for traceability).

Supports three scoring backends selected by method: - "scanpy" (default): scanpy score_genes (AddModuleScore-style control subtraction). - "ucell": rank-based AUC (requires pyucell; robust to normalization and sparsity). - "ssgsea": single-sample GSEA (requires gseapy; slower; publication-grade).

Writes raw scores to adata.obsm["signature_score"] and z-scored (across obs) to adata.obsm["signature_score_z"]. Column names are full path (e.g. Myeloid/Macrophage_Core). Report is stored in adata.uns["signature_score_report"]. The method used is recorded in adata.uns["scoring_method"].

Parameters:

  • adata (AnnData) – Annotated data. Must have adata.raw if use_raw=True and method="scanpy".

  • signatures_nested (dict, str, or Path) – Two-level nested dict {category: {subprocess: [genes]}} or path to JSON. Keys like _meta are skipped.

  • method (str) – Scoring backend: "scanpy" | "ucell" | "ssgsea". Default "scanpy".

  • use_raw (bool) – Use adata.raw for expression (default True). Only applies to method="scanpy".

  • ctrl_size (int) – Control genes per signature gene (default 50). Only applies to method="scanpy".

  • n_bins (int) – Expression bins for control matching (default 25). Only applies to method="scanpy".

  • min_genes (int) – Minimum present genes to score (default 3); sets below this threshold receive NaN.

  • copy (bool) – If True, copy adata before modifying (default False).

  • save_path (str, Path, or None) – If set, write adata to this path after scoring (default None).

Returns:

AnnData with obsm["signature_score"], obsm["signature_score_z"], uns["signature_score_report"], and uns["scoring_method"]. Same object as input if copy=False.

Return type:

AnnData

Gene Set Loaders and Curation#

sc_tools.tl.load_hallmark(organism='human')[source]#

Load bundled MSigDB Hallmark gene sets.

Returns the 50 Hallmark gene sets as a two-level dict:

{"Hallmark": {"TNFA_SIGNALING_VIA_NFKB": ["ABCA1", ...], ...}}

The HALLMARK_ prefix is stripped from set names so that column names produced by score_signature are concise (e.g. Hallmark/HYPOXIA instead of Hallmark/HALLMARK_HYPOXIA).

Parameters:

organism (str) – Only "human" is currently supported.

Returns:

Two-level nested dict {category: {name: [genes]}}.

Return type:

dict

sc_tools.tl.load_msigdb_json(path, category_name=None)[source]#

Load an MSigDB-format JSON file.

MSigDB JSON format: a flat dict where each key is a set name and the value is either a plain list of gene symbols or an object with a geneSymbols key (the richer export format).

Parameters:

  • path (str or Path) – Path to the MSigDB JSON file.

  • category_name (str or None) – Category name for the outer dict key. If None, inferred from the filename stem (e.g. "h.all.v2025.1.Hs" becomes "h.all").

Returns:

Two-level nested dict {category_name: {set_name: [genes]}}.

Return type:

dict

sc_tools.tl.load_gmt(path, category_name=None)[source]#

Load a GMT (Gene Matrix Transposed) file.

GMT format: tab-separated; each line is:

SET_NAME\tDESCRIPTION\tGENE1\tGENE2\t...
Parameters:

  • path (str or Path) – Path to the GMT file.

  • category_name (str or None) – Category name for the outer dict key. If None, the file stem is used.

Returns:

Two-level nested dict {category_name: {set_name: [genes]}}.

Return type:

dict

sc_tools.tl.list_gene_sets()[source]#

List names of all bundled gene set collections.

Returns:

Names of bundled collections (usable as organism hints).

Return type:

list[str]

sc_tools.tl.validate_gene_signatures(signatures, var_names=None, min_genes=3)[source]#

Validate a nested gene signature dict or JSON file.

Checks every leaf gene list and reports coverage, duplicates, and empty sets. Optionally checks presence against a provided gene universe (e.g. adata.var_names).

Parameters:

  • signatures (dict, str, or Path) – Two-level nested dict or path to a JSON file containing signatures.

  • var_names (list[str] or None) – Gene universe (e.g. list(adata.var_names)). If provided, reports n_present and pct_coverage.

  • min_genes (int) – Minimum number of genes required per set (flags sets below threshold).

Returns:

One row per signature with columns: signature, n_genes, n_unique, n_duplicates, n_present (if var_names given), n_missing (if var_names given), pct_coverage (if var_names given), status.

Return type:

pd.DataFrame

sc_tools.tl.merge_gene_signatures(*dicts)[source]#

Combine multiple two-level signature dicts.

Later dicts overwrite earlier ones on key collision at both the category and set-name level. The special _meta key is skipped.

Parameters:

*dicts (dict) – Two-level nested dicts to merge.

Returns:

Merged two-level nested dict.

Return type:

dict

Examples

>>> project = {"Myeloid": {"Macrophage": ["CD68", "CSF1R"]}}
>>> hallmark = load_hallmark()
>>> combined = merge_gene_signatures(project, hallmark)
sc_tools.tl.update_gene_symbols(signatures, alias_map)[source]#

Replace deprecated or alias gene symbols in a signature dict.

Does NOT fetch from the internet; the caller provides the alias map (e.g. from an HGNC download or a per-project correction list).

Parameters:

  • signatures (dict) – Two-level nested dict of gene signatures.

  • alias_map (dict[str, str]) – Mapping from old symbol to new symbol, e.g. {"FAM19A5": "TAFA5"}. Case-sensitive. Symbols not in the map are left unchanged.

Returns:

New nested dict with symbols updated. Original dict is not mutated.

Return type:

dict

sc_tools.tl.save_gene_signatures(signatures, path)[source]#

Write a gene signature dict to JSON with consistent formatting.

Adds (or updates) a _meta key at the top level with a datestamp. Keys are sorted; indent is 2.

Parameters:

  • signatures (dict) – Two-level nested dict of gene signatures.

  • path (str or Path) – Output JSON path. Parent directories are created if needed.

Return type:

None

Enrichment Testing#

Group-level enrichment. Optional dependency: pip install sc-tools[geneset].

sc_tools.tl.run_ora(adata, groupby, gene_set_dict, background=None, min_genes=3, fdr_method='fdr_bh')[source]#

Over-representation analysis (ORA) using Fisher exact test.

For each (group, gene set) pair, tests whether the overlap between the gene set and the genes expressed in that group is greater than expected by chance, relative to the background gene universe.

Multiple testing correction is Benjamini-Hochberg (BH) applied across all tests within each group.

Parameters:

  • adata (AnnData) – Annotated data matrix. Must have adata.obs[groupby].

  • groupby (str) – Column in adata.obs defining groups (e.g. "leiden").

  • gene_set_dict (dict) – Gene sets as {category: {name: [genes]}} or {name: [genes]}.

  • background (list[str] or None) – Background gene universe. Defaults to adata.var_names.

  • min_genes (int) – Minimum overlap required to include a test (default 3).

  • fdr_method (str) – Multiple testing correction method passed to statsmodels.stats.multitest.multipletests (default "fdr_bh").

Returns:

Long-format results with columns: group, category, gene_set, n_overlap, n_set, n_background, n_group_genes, odds_ratio, p_val, p_adj. Sorted by p_adj ascending within each group.

Return type:

pd.DataFrame

sc_tools.tl.run_gsea_pseudobulk(adata, groupby, gene_set_dict, layer=None, method='prerank', ranking_stat='logfc', min_genes=10, n_permutations=1000)[source]#

Pseudobulk GSEA: aggregate expression per group, rank genes, run fgsea/prerank.

Requires gseapy (optional dependency). Install with:

pip install gseapy

Steps:

  1. Aggregate expression: mean per group → pseudobulk matrix (genes x groups).

  2. Compute ranking statistic (log2FC, z-score, or t-stat) relative to the mean expression across all groups.

  3. Pass ranked gene list to gseapy.prerank for each group.

  4. Collect results; apply Benjamini-Hochberg FDR across all gene sets per group.

Parameters:

  • adata (AnnData) – Annotated data. Must have adata.obs[groupby].

  • groupby (str) – Column in adata.obs defining groups.

  • gene_set_dict (dict) – Gene sets as {category: {name: [genes]}} or {name: [genes]}.

  • layer (str or None) – Layer to use for expression. Defaults to adata.X.

  • method (str) – GSEA method: "prerank" (default) or "fgsea" (also via gseapy).

  • ranking_stat (str) – Statistic used to rank genes: "logfc" (default), "zscore", "tstat".

  • min_genes (int) – Minimum genes in a set after intersecting with ranked genes (default 10).

  • n_permutations (int) – Number of permutations for enrichment scoring (default 1000).

Returns:

Long-format results with columns: group, category, gene_set, NES, p_val, p_adj, lead_edge_genes. Sorted by p_adj ascending within each group.

Return type:

pd.DataFrame

Cell-type Deconvolution#

Backend registry: cell2location, tangram, destvi. Optional dependency: pip install sc-tools[deconvolution].

sc_tools.tl.deconvolution(spatial_adata, sc_adata=None, *, method='cell2location', celltype_key='celltype', spatial_batch_key='library_id', sc_batch_key=None, reference_profiles=None, n_signature_genes=2000, use_gpu=None, qc_labels=None, method_kwargs=None, cache_dir=None, output_dir=None, output_file=None, logger_instance=None)[source]#

Run cell-type deconvolution on spatial transcriptomics data.

Processes each library (spatial_batch_key) independently using backed AnnData loading to keep peak memory low.

Parameters:

  • spatial_adata (AnnData | str | Path) – Spatial AnnData or path to h5ad file.

  • sc_adata (AnnData | str | Path | None) – Single-cell reference AnnData or path. Required unless reference_profiles is provided (Cell2location only).

  • method (str) – Backend name: "cell2location" (default), "tangram", "destvi".

  • celltype_key (str) – Column in sc_adata.obs with cell-type labels.

  • spatial_batch_key (str) – Column in spatial obs for per-library batching.

  • sc_batch_key (str | None) – Batch key in sc_adata.obs (used for HVG selection).

  • reference_profiles (DataFrame | Path | None) – Pre-computed reference profiles (DataFrame or path to pickle). Memory optimisation for Cell2location – skips regression training.

  • n_signature_genes (int) – Number of signature genes for deconvolution.

  • use_gpu (bool | None) – GPU setting. None = auto-detect.

  • qc_labels (list[str] | None) – Cell-type labels to exclude from the reference.

  • method_kwargs (dict | None) – Extra keyword arguments forwarded to the backend run() method.

  • cache_dir (str | Path | None) – Directory for caching signature genes and reference profiles.

  • output_dir (str | Path | None) – Directory for per-library intermediate results.

  • output_file (str | Path | None) – Path to save the final AnnData with proportions.

  • logger_instance (Logger | None) – Optional logger.

Returns:

The spatial AnnData with obsm['cell_type_proportions'] (n_spots x n_celltypes) and obs['{method}_argmax'] (dominant cell type per spot).

Return type:

AnnData

sc_tools.tl.extract_reference_profiles(sc_adata, celltype_key, genes=None, qc_labels=None, cache_path=None, logger_instance=None)[source]#

Compute mean expression per cell type from scRNA-seq reference.

The resulting DataFrame (genes x cell_types) is ~100x smaller than the full reference and can be passed directly to Cell2location via its cell_state_df parameter, skipping regression model training.

Parameters:

  • sc_adata (AnnData) – Single-cell reference AnnData.

  • celltype_key (str) – Column in sc_adata.obs with cell-type labels.

  • genes (list[str] | None) – Subset of genes to include. None keeps all.

  • qc_labels (list[str] | None) – Cell-type labels to exclude (e.g. ["Doublets", "QC_Filtered"]).

  • cache_path (str | Path | None) – If given, cache the result as a pickle file.

  • logger_instance (Logger | None) – Optional logger.

Returns:

Genes (rows) x cell types (columns) mean expression matrix.

Return type:

pandas.DataFrame

sc_tools.tl.select_signature_genes(sc_adata, celltype_key, sc_batch_key, n_genes_max, skip_hvg=True, cache_dir=None, force_recompute=False, sc_data_file=None, logger_instance=None)[source]#

Select signature genes for deconvolution.

This function selects genes for deconvolution by combining: 1. Highly variable genes (HVGs) - optional 2. Top marker genes per cell type (always computed)

Results can be cached to avoid recomputation.

Parameters:

  • sc_adata (AnnData) – Single-cell reference data

  • celltype_key (str) – Column name in sc_adata.obs containing cell type annotations

  • sc_batch_key (str) – Column name in sc_adata.obs containing batch information

  • n_genes_max (int) – Maximum number of genes to return

  • skip_hvg (bool) – If True, skip HVG computation and use only marker genes (faster, less memory)

  • cache_dir (str, optional) – Directory to cache signature genes. If None, caching is disabled.

  • force_recompute (bool) – If True, force recomputation even if cache exists.

  • sc_data_file (str, optional) – Path to single-cell data file (for cache key generation).

  • logger_instance (Logger, optional) – Custom logger instance. If None, uses module logger.

Returns:

List of signature gene names

Return type:

list of str

Colocalization#

sc_tools.tl.colocalization.truncated_similarity(score_a, score_b)[source]#

Truncated similarity: score_a * score_b where both > 0, else 0.

Parameters:

  • score_a (array-like) – First score vector (e.g. proliferation).

  • score_b (array-like) – Second score vector (e.g. macrophage).

Returns:

Same shape as inputs; product where both > 0, else 0.

Return type:

np.ndarray

sc_tools.tl.colocalization.pearson_correlation(adata, sig_columns, min_valid_ratio=0.5)[source]#

Compute Pearson correlation matrix between signatures across spots.

Parameters:

  • adata (AnnData) – Annotated data object with signature scores

  • sig_columns (list of str) – List of signature column names

  • min_valid_ratio (float) – Minimum ratio of non-NaN values required (default: 0.5)

Returns:

Correlation matrix

Return type:

DataFrame

sc_tools.tl.colocalization.morans_i(adata, sig_column, coord_key='spatial', n_perms=1000, n_jobs=1)[source]#

Compute Moran’s I spatial autocorrelation for a signature using squidpy.

Note: squidpy’s spatial_autocorr expects genes in var_names, but signature scores are typically in obs columns. This function temporarily adds the signature as a “pseudo-gene” to compute Moran’s I.

Parameters:

  • adata (AnnData) – Annotated data object

  • sig_column (str) – Signature column name in adata.obs

  • coord_key (str) – Key in adata.obsm for spatial coordinates (default: ‘spatial’)

  • n_perms (int) – Number of permutations for p-value calculation (default: 1000)

  • n_jobs (int) – Number of jobs for parallel processing (default: 1)

Returns:

Dictionary with ‘I’ (Moran’s I) and ‘pval’ (p-value)

Return type:

dict

sc_tools.tl.colocalization.morans_i_batch(adata, sig_columns, n_perms=1000, coord_key='spatial', n_jobs=1)[source]#

Compute Moran’s I for multiple signatures using squidpy.

This function efficiently computes Moran’s I for multiple signatures by temporarily adding them all as “pseudo-genes” and computing in one batch.

Parameters:

  • adata (AnnData) – Annotated data object

  • sig_columns (list of str) – List of signature column names

  • n_perms (int) – Number of permutations (default: 1000)

  • coord_key (str) – Key in adata.obsm for spatial coordinates (default: ‘spatial’)

  • n_jobs (int) – Number of jobs for parallel processing (default: 1)

Returns:

DataFrame with columns ‘Morans_I’ and ‘p_value’

Return type:

DataFrame

sc_tools.tl.colocalization.neighborhood_enrichment(adata, sig_column, threshold_low=-1.0, threshold_high=1.0, n_perms=1000)[source]#

Compute thresholded neighborhood enrichment for a signature.

Parameters:

  • adata (AnnData) – Annotated data object

  • sig_column (str) – Signature column name

  • threshold_low (float) – Threshold for low category (default: -1.0)

  • threshold_high (float) – Threshold for high category (default: 1.0)

  • n_perms (int) – Number of permutations (default: 1000)

Returns:

Neighborhood enrichment score

Return type:

float

sc_tools.tl.colocalization.neighborhood_enrichment_batch(adata, sig_columns, threshold_low=-1.0, threshold_high=1.0, n_perms=1000)[source]#

Compute neighborhood enrichment for multiple signatures.

Parameters:

  • adata (AnnData) – Annotated data object

  • sig_columns (list of str) – List of signature column names

  • threshold_low (float) – Threshold for low category (default: -1.0)

  • threshold_high (float) – Threshold for high category (default: 1.0)

  • n_perms (int) – Number of permutations (default: 1000)

Returns:

DataFrame with ‘nhood_enrichment’ column

Return type:

DataFrame