v0.6.0

What's Changed

• Merging previous PRs by @earmingol in #27
• v0.6.0 by @earmingol in #28

Full Changelog: v0.5.0...v0.6.0

v0.6.0

New Features

Complex/multi-element support for communication scoring

• New add_complexes_to_adata function in sccellfie.preprocessing that aggregates
  multi-gene complexes (e.g., protein complex subunits like ITGA4+ITGB1) or
  multi-task metabolic elements into single per-cell variables in AnnData. Supports
  min, mean, and gmean aggregation methods. The aggregated variables
  integrate seamlessly with existing compute_communication_scores and
  compute_local_colocalization_scores functions.
• New make_complex_name helper for generating canonical complex names from
  subunit lists.
• New prepare_var_pairs convenience function that accepts var_pairs with list
  elements (e.g., [(['TASK1', 'TASK2'], ['GENE1', 'GENE2'])]), auto-generates
  complex names, adds them to adata, and returns normalized string-only pairs
  ready for scoring.

Gene ablation impact analysis (sccellfie.stats.ablation)

• compute_gene_ablation_impact(gpr_source, task_by_rxn, ...) — simulates
  single-gene ablation on a synthetic uniform-expression reference and returns
  (gene × task) impact DataFrames (rel_change, abs_change,
  fraction_zeroed). Reuses scCellFie's own GPR walker
  (compute_gpr_gene_score), so results exactly reflect what the scoring
  pipeline would produce if a gene were missing — no hand-rules.
• essential_genes_from_ablation(impact, metric, threshold, topology=None, ...)
  — thresholds the impact DataFrames into {task: [genes]}, optionally
  filtered by network topology.

Optional network-topology filter (compute_reaction_topology_essentiality)

• Given a user-provided cobra Model and per-task
  (start_metabolite, end_metabolite) pairs, builds a metabolite-node /
  reaction-edge graph per task and flags reactions whose removal disconnects
  start from end.
• Combines with the ablation impact to identify genes that affect at least one
  topologically-irreplaceable reaction.
• ignore_metabolites strips currency metabolites (ATP, H2O, ...) before the
  connectivity test.

Dataset and per-cell completeness reports (sccellfie.reports.completeness)

• Evaluates an AnnData against the database using the ablation impact as a
  weighting. Flat DataFrames carry both essential-gene and all-gene
  scopes side-by-side via _essential / _all suffixed columns (binary
  fraction_present_* and impact-weighted impact_weighted_completeness_*).
• Per-cell scoring treats both dataset-absent genes and per-cell
  zero-expression genes as missing; writes adata.obs['completeness_essential']
  / adata.obs['completeness_all'].
• New functions:
  • compute_dataset_completeness(adata, gpr_source, task_by_rxn, ...)
    → {task_completeness, reaction_completeness, overall_summary}.
  • compute_cell_completeness(adata, gpr_source, task_by_rxn, ...)
    → {per_cell, matrix_essential, matrix_all} (matrices only when
    return_matrix=True).
  • generate_completeness_report(adata, ...) — runs both in one call, sharing
    the ablation impact across sub-reports.

Pipeline opt-in

• run_sccellfie_pipeline(..., compute_ablation_impact=True) runs the ablation
  on the filtered gpr_rules after preprocessing and attaches the result to
  preprocessed_db['ablation_impact'] (default False; preserves existing
  behavior).

Dataset-wise threshold estimator (sccellfie.expression.get_sccellfie_dataset_threshold)

• Streaming, memory-bounded reimplementation of the original atlas-derived
  sccellfie_threshold workflow for a single (possibly backed) AnnData — users
  can now produce thresholds tailored to their own dataset instead of relying on
  the shipped atlas values, which is particularly useful for large spatial
  assays (Visium-HD, Stereo-seq, upcoming Atera, ...) up to ~10M cells on a
  laptop.
• Replaces the original disk-based external sort with a vectorized
  Vitter Algorithm R reservoir sampler (default 5M × float32 ≈ 20 MB) for
  global P25/P75 estimation. Peak memory is bounded by
  reservoir_size + one_chunk, independent of n_cells.
• Threshold rule is a faithful port of the original clip-with-low-expression-escape
  logic: clip(nonzero_mean, P25, P75) if max > P25 or max == 0, else the
  raw nonzero_mean — keeping dropout-heavy genes from being floored up to P25.
• Respects adata.uns['normalization']['method']=='total_counts' (skips
  re-normalization); n_counts_key auto-detects among total_counts, n_counts,
  raw_sum, nCount_RNA, or falls back to computing from the full chunk before
  gene subsetting.
• Applies CORRECT_GENES[organism] gene-name fix-ups (e.g. MT-CO1 → COX1,
  MUT → MMUT) on a local copy of var_names before intersection with the
  metabolic gene set, matching the main pipeline's behavior.
• gene_set defaults to the full metabolic gene list in the shipped database;
  accepts a list/set/Index or a .json path (supports
  All_metabolic_genes.json-style files).
• Returns a single-column pd.DataFrame with first column sccellfie_threshold,
  ready for drop-in replacement of the shipped Thresholds.csv and positional
  consumption by compute_gene_scores.

Configurable percentile bounds for the dataset-wise threshold rule

• get_sccellfie_dataset_threshold now takes lower_percentile=25 and
  upper_percentile=75 arguments that drive the clip rule:
  clip(nonzero_mean, P_lower, P_upper) if max > P_lower or max == 0,
  else raw nonzero_mean. Defaults reproduce the original atlas behavior
  bit-exactly; widen (e.g. 10, 90) for a more permissive threshold or
  tighten (e.g. 30, 70) for a stricter one. Validated against
  0 ≤ lower_percentile < upper_percentile ≤ 100.
• The reporting percentiles= tuple is now automatically merged with the
  rule's bounds, so they're always available in stats['percentiles'] when
  return_stats=True. Stat keys preserve int-vs-float types, so non-integer
  percentile choices (e.g. 12.5) are reported losslessly.

Multi-panel layout and title controls in plot_segmentation

• sccellfie.plotting.plot_segmentation now accepts a list/tuple for
  color_by (e.g. ['task_A', 'task_B', 'GENE1']) and renders a panel
  grid laid out by a new ncols=4 argument, mirroring sc.pl.spatial's
  semantics. Shared geometry/crop/axis-limits are computed once and
  applied uniformly across panels; each panel is coloured independently
  with its own legend (categorical) or colorbar (continuous).
• figsize becomes the per-panel size in multi-panel mode (default
  (4, 4), total figure scales to (figsize[0]*ncols, figsize[1]*nrows));
  single-panel behavior is unchanged ((10, 10) default).
• Return type matches the input: a single Axes for color_by=str
  (back-compat), a 2D (nrows, ncols) array of Axes for color_by=list.
  Trailing unused axes are hidden automatically.
• Passing both ax= and a list of features raises a clear ValueError.
• Titles are now drawn in both single- and multi-panel modes — by
  default each panel is titled with its feature name (in single mode the
  title falls back to celltype_key when color_by=None).
  panel_titles=False suppresses titles entirely. New title= lets
  the user override: a string for single-panel mode, a list-matching-
  color_by for multi-panel mode (length-mismatch raises). Long
  titles auto-wrap via :func:textwrap.wrap at the new
  wrapped_title_length=45 characters; title_fontsize=12 controls
  the font size — both mirror the convention already used in
  plot_spatial, create_multi_violin_plots and
  create_volcano_plot.

Bug Fixes

Scalebar padding/positioning in plot_segmentation

• The scalebar in sccellfie.plotting.plot_segmentation could land too
  close to — or even overlap with — the cells, particularly in small
  figures or when y_pad_ratio was small. The bar's vertical position
  was computed in data coordinates against the full padded view height,
  and the label was offset upward in display points regardless of axis
  inversion, so for lower_* positions the label was always pushed
  toward the data instead of into the whitespace.
• The scalebar now uses a blended transform (data X, axes-fraction Y),
  so pad_frac is a reliable inset from the axes corner regardless of
  data extent or figure size. The label is always placed on the side of
  the bar away from the data (below for lower_*, above for upper_*),
  guaranteeing no overlap with cells when y_pad_ratio > 0. New
  text_pad_pts (default 2.0) controls the gap between the bar and
  its label, exposed via scalebar_kwargs.

Details

Complexes

• Complex expression is computed at the single-cell level and stored in
  adata.X (and optionally layers), ensuring correct fraction-above-threshold
  calculations in downstream communication scoring.
• Sparse matrix format (CSR) is preserved when input is sparse.
• Complex metadata (subunit composition, aggregation method) is stored in
  adata.var.

Ablation & completeness

• Under the uniform reference (gene_score=1.0 for every gene), baseline RAL
  and MTS are both 1.0, so rel_change is numerically equal to
  1 - ablated_mts.
• rel_change and fraction_zeroed are invariant to uniform_score scaling;
  abs_change scales linearly with it.
• Only reactions whose GPR contains the ablated gene are re-walked per
  ablation; only tasks touching those reactions are re-aggregated. Runtime is
  sub-second on the shipped Recon2-2 / iMM1415 databases.
• The topology pass treats reversible reactions as bidirectional by default
  (treat_reversible_as_bidirectional=True). Reactions missing from the cobra
  Model are flagged with an aggregated warning and their topology rows stay
  False.
• essential_genes_from_ablation(..., topology=..., fallback_to_ablation_only=True)
  still returns ablation-only essentials for tasks without endpoints, so a
  partial task_endpoints dict still yields a coherent result.
• Completeness at the dataset level uses gene ∉ adata.var_names (assay
  coverage); at the per-cell level it additionally treats expression == 0 in
  that cell as missing, so dropouts and low-expression cells drop the per-cell
  score.

Usage Examples

Complex communication scoring — explicit definition with add_complexes_to_adata

import numpy as np
import sccellfie

# Combine metabolic task scores and gene expression into one AnnData
adata.metabolic_tasks.var['type'] = 'metabolic score'
bdata.var['type'] = 'gene expression'

adata_updated = sccellfie.preprocessing.transfer_variables(
    adata_target=bdata,
    adata_source=adata.metabolic_tasks,
    var_names=adata.metabolic_tasks.var_names
)

# Add complexes for multi-subunit receptors or multi-task metabolic elements
complexes = {
    'taskA&taskB': ['taskA', 'taskB'],       # multiple tasks producing same metabolite
    'GENE1&GENE2': ['GENE1', 'GENE2'],       # receptor complex subunits
}
sccellfie.preprocessing.add_complexes_to_adata(adata_updated, complexes, agg_method='min')

# Use complex names alongside regular LR pairs
# Ligands can be metabolic tasks or genes (e.g., secreted proteins)
lr_pairs = [
    ('taskA&taskB', 'GENE1&GENE2'),   # complex task ligand -> complex receptor
    ('taskC', 'GENE3'),                # single task ligand -> single receptor
    ('GENE4', 'GENE5'),               # single gene ligand -> single receptor
    ('GENE4', 'GENE1&GENE2'),         # single gene ligand -> complex receptor
]

ccc_df = sccellfie.communication.compute_communication_scores(
    adata_updated,
    var_pairs=lr_pairs,
    groupby=cell_group,
    communication_score='gmean',
    agg_func='trimean',
    ligand_threshold=np.log(2),
    receptor_threshold=0.,
)

Complex communication scoring — automated with prepare_var_pairs

import sccellfie

# Same transfer step as above
adata_updated = sccellfie.preprocessing.transfer_variables(
    adata_target=bdata,
    adata_source=adata.metabolic_tasks,
    var_names=adata.metabolic_tasks.var_names
)

# Pass list elements directly — complexes are auto-detected, named, and added
# Ligands can be metabolic tasks or genes (e.g., secreted proteins)
lr_pairs = [
    (['taskA', 'taskB'], ['GENE1', 'GENE2']),  # complex task ligand -> complex receptor
    ('taskC', 'GENE3'),                          # single task ligand -> single receptor
    ('GENE4', 'GENE5'),                          # single gene ligand -> single receptor
    ('GENE4', ['GENE1', 'GENE2']),               # single gene ligand -> complex receptor
]

normalized_pairs = sccellfie.preprocessing.prepare_var_pairs(
    adata_updated, lr_pairs, agg_method='min'
)

ccc_df = sccellfie.communication.compute_communication_scores(
    adata_updated,
    var_pairs=normalized_pairs,
    groupby=cell_group,
    communication_score='gmean',
    agg_func='trimean',
)

Ablation + completeness (no cobra Model needed)

import sccellfie

# 1. Ablation impact on the raw database.
db = sccellfie.datasets.load_sccellfie_database(organism='human')
gpr_strings = db['rxn_info'].set_index('Reaction')['GPR-symbol'].to_dict()
impact = sccellfie.stats.compute_gene_ablation_impact(gpr_strings, db['task_by_rxn'])

# 2. Essential gene set per task (zero-forcing under uniform reference).
essential = sccellfie.stats.essential_genes_from_ablation(
    impact, metric='fraction_zeroed', threshold=1.0,
)

# 3. Full completeness report against the user's AnnData (writes .obs).
report = sccellfie.reports.generate_completeness_report(
    adata, gpr_strings, db['task_by_rxn'],
    ablation_impact=impact,                    # reuse the DataFrames
    metric='fraction_zeroed', threshold=1.0,
    write_to_obs=True, return_matrix=False,
)

report['dataset']['task_completeness']       # per task: n/fraction/impact for both scopes
report['dataset']['reaction_completeness']   # same columns, per reaction
report['dataset']['overall_summary']         # single-row summary (tasks / reactions / DB coverage)
report['cell']['per_cell']                   # cells × [completeness_essential, completeness_all]

# Rank cells by completeness for quality flagging
adata.obs.sort_values('completeness_essential').head(20)

Ablation via the pipeline

results = sccellfie.run_sccellfie_pipeline(
    adata, organism='human',
    compute_ablation_impact=True,     # new kwarg
    # ... existing kwargs
)
results['ablation_impact']['fraction_zeroed']   # gene × task

Dataset-wise thresholds (no atlas needed)

import sccellfie

# Load (or lazily back) the user's AnnData — can be raw counts or pre-normalized
adata = sc.read_h5ad('my_spatial_dataset.h5ad', backed='r')  # backed works

# 1. Compute thresholds from the dataset itself
thr = sccellfie.expression.get_sccellfie_dataset_threshold(
    adata,
    organism='human',           # drives CORRECT_GENES and default gene_set
    reservoir_size=2_000_000,   # tune for memory/accuracy
    chunk_size=100_000,
    random_state=0,
)
thr.head()                       # pd.DataFrame with column 'sccellfie_threshold'

# 2. Drop-in replacement for the shipped thresholds in the pipeline
db = sccellfie.datasets.load_sccellfie_database(organism='human')
db['thresholds'] = thr           # overrides the atlas-derived column
results = sccellfie.run_sccellfie_pipeline(adata, organism='human', **db)

# Optional: inspect the intermediate per-gene stats
thr, stats = sccellfie.expression.get_sccellfie_dataset_threshold(
    adata, return_stats=True, verbose=False,
)
stats['percentiles']             # {10: ..., 25: ..., 50: ..., 75: ..., 90: ..., 95: ...}
stats['nonzero_mean']            # pd.Series, per gene

Path-topology filter (requires a cobra Model)

import cobra
recon = cobra.io.load_json_model('recon2_2.json')

# User-curated {task: (start_metabolite, end_metabolite)}
task_endpoints = {
    'Glycolysis': ('glc_D_c', 'pyr_c'),
    # ...
}

topology = sccellfie.stats.compute_reaction_topology_essentiality(
    db['task_by_rxn'], recon, task_endpoints,
    ignore_metabolites={'atp_c','adp_c','h_c','h2o_c','nadh_m','nad_m','pi_c','coa_c'},
)

# A gene is network-essential iff it clears the ablation threshold AND affects
# at least one topologically-irreplaceable reaction in the task.
essential_net = sccellfie.stats.essential_genes_from_ablation(
    impact, metric='fraction_zeroed', threshold=1.0,
    topology=topology, task_by_rxn=db['task_by_rxn'], gpr_source=gpr_strings,
)

Custom percentile bounds for the dataset-wise threshold

import sccellfie

# Permissive: clip between P10 and P90 (looser low-expression gating).
thr_wide = sccellfie.expression.get_sccellfie_dataset_threshold(
    adata, organism='human',
    lower_percentile=10, upper_percentile=90,
    random_state=0,
)

# Strict: clip between P30 and P70.
thr_tight = sccellfie.expression.get_sccellfie_dataset_threshold(
    adata, organism='human',
    lower_percentile=30, upper_percentile=70,
    random_state=0,
)

# Inspect the actual percentile values used by the rule.
thr, stats = sccellfie.expression.get_sccellfie_dataset_threshold(
    adata, lower_percentile=10, upper_percentile=90,
    return_stats=True, verbose=False,
)
stats['percentiles'][10], stats['percentiles'][90]

Multi-panel plot_segmentation

import sccellfie

# Render four metabolic-task scores as cell polygons, side by side.
fig, axes = sccellfie.plotting.plot_segmentation(
    adata,
    color_by=['Glycolysis', 'OXPHOS', 'Fatty acid oxidation', 'Glutaminolysis'],
    segmentation=seg,        # {cell_id: shapely.Polygon}
    ncols=2,                  # 2x2 grid
    figsize=(4, 4),           # per-panel size; total figure = 8x8 inches
    cmap='magma',
)

# Mix categorical and continuous variables in the same figure.
fig, axes = sccellfie.plotting.plot_segmentation(
    adata,
    color_by=['cell_type', 'Glycolysis', 'OXPHOS'],
    segmentation=seg,
    ncols=3,
)
# axes[0, 0] has a categorical legend; axes[0, 1] / axes[0, 2] have colorbars.

# Suppress per-panel titles if you want to add your own.
fig, axes = sccellfie.plotting.plot_segmentation(
    adata,
    color_by=['GENE1', 'GENE2'],
    panel_titles=False,
    ncols=2,
)

# Long task names auto-wrap; title_fontsize and wrapped_title_length mirror
# the convention in plot_spatial / create_volcano_plot.
fig, axes = sccellfie.plotting.plot_segmentation(
    adata,
    color_by=['Beta-oxidation of long-chain fatty acids', 'OXPHOS'],
    title=['Long-chain FA β-oxidation', 'Oxidative phosphorylation'],
    title_fontsize=14,
    wrapped_title_length=20,
    ncols=2,
)

# Titles also work for a single panel (set to the feature name by default).
fig, ax = sccellfie.plotting.plot_segmentation(
    adata, color_by='Glycolysis', title_fontsize=16,
)
ax.get_title()   # 'Glycolysis'