Pancreas MultiScorer tutorial¶

This notebook demonstrates a three-condition ordered design. In `DEMO_MODE`, a
graded controlled transition tilt toward Beta validates omnibus, post-hoc, and
planned-contrast workflows. Replace it with genuine biological conditions and
replicates for scientific inference.

Execution provenance. The displayed outputs were generated with scCS 0.8.0.dev33. Version 0.8.0.dev34 changes documentation and tutorial explanation only; the scientific calculations are unchanged.

1. Installation, versions, and analysis settings¶

[1]:

from __future__ import annotations
import warnings

from pathlib import Path
import importlib.metadata as importlib_metadata
import platform
import sys

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scanpy as sc
import scvelo as scv
from scipy.spatial.distance import jensenshannon

import scCS

SEED = 20260714
N_JOBS = 1
sc.settings.verbosity = 2
scv.settings.verbosity = 2
sc.set_figure_params(dpi=100, facecolor="white")

print("Python", sys.version.split()[0])
print("Platform", platform.platform())
print("scCS", scCS.__version__)
print("Scanpy", importlib_metadata.version("scanpy"))
print("scVelo", scv.__version__)

OUTPUT_DIR = Path("tutorial_outputs/pancreas_multi")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

# Hide known warnings emitted by optional upstream dependencies. These do not
# change the scCS calculation and would otherwise distract from the tutorial.
warnings.filterwarnings(
    "ignore",
    message=r"This process .* is multi-threaded, use of fork\(\) may lead to deadlocks.*",
    category=DeprecationWarning,
)
warnings.filterwarnings(
    "ignore",
    category=DeprecationWarning,
    module=r"cellrank\..*",
)

Python 3.12.13
Platform Linux-6.6.114.1-microsoft-standard-WSL2-x86_64-with-glibc2.39
scCS 0.8.0.dev33
Scanpy 1.11.5
scVelo 0.3.4

[2]:

DEMO_MODE = True
CONDITION_KEY = "condition"
REPLICATE_KEY = "sample_id"
CONDITIONS = ("control", "low", "high")
REPLICATES_PER_CONDITION = 4
CONDITION_ORDER = list(CONDITIONS)
RUN_SCOPE_SENSITIVITY = False
RUN_MIXED_MODEL_SENSITIVITY = False

2. Load the dataset and velocity graph¶

This tutorial deliberately starts from the public dataset and recomputes the velocity model. No pre-existing scCS cache is assumed. Dynamical fitting can take several minutes, but showing the complete preparation makes every analysis step reproducible for a first-time user.

[3]:

def clear_velocity_outputs(adata, *, clear_dynamics=False):
    """Remove stale velocity fields before fitting a requested model."""
    for key in ("velocity", "velocity_u", "velocity_variance"):
        adata.layers.pop(key, None)
    for key in ("velocity_graph", "velocity_graph_neg", "velocity_params"):
        adata.uns.pop(key, None)
    for key in (
        "velocity_self_transition",
        "root_cells",
        "end_points",
        "velocity_pseudotime",
        "latent_time",
    ):
        if key in adata.obs:
            del adata.obs[key]
    if clear_dynamics:
        for key in list(adata.var.columns):
            if str(key).startswith("fit_"):
                del adata.var[key]
        for key in list(adata.layers):
            if str(key).startswith("fit_"):
                del adata.layers[key]
        adata.uns.pop("recover_dynamics", None)


def ensure_pca_neighbors_moments(
    adata,
    *,
    n_pcs=30,
    n_neighbors=30,
    preserve_existing_neighbors=False,
):
    """Create the unbiased PCA neighbor graph used by RNA velocity."""
    if "X_pca" not in adata.obsm:
        n_comps = min(n_pcs, adata.n_obs - 1, adata.n_vars - 1)
        sc.pp.pca(adata, n_comps=n_comps)
    if "neighbors" not in adata.uns or not preserve_existing_neighbors:
        sc.pp.neighbors(
            adata,
            n_neighbors=min(n_neighbors, adata.n_obs - 1),
            n_pcs=min(n_pcs, adata.obsm["X_pca"].shape[1]),
            use_rep="X_pca",
            random_state=SEED,
        )
    if not {"Ms", "Mu"}.issubset(adata.layers):
        scv.pp.moments(adata, n_neighbors=None, n_pcs=None)
    if "X_umap" not in adata.obsm:
        sc.tl.umap(adata, random_state=SEED)


def resolve_present_genes(adata, candidates):
    """Resolve gene symbols case-insensitively and report missing markers."""
    lookup = {str(gene).lower(): str(gene) for gene in adata.var_names}
    resolved, missing = [], []
    for candidate in candidates:
        match = lookup.get(str(candidate).lower())
        if match is None:
            missing.append(str(candidate))
        elif match not in resolved:
            resolved.append(match)
    return resolved, missing


def row_js(left, right):
    """Row-wise Jensen-Shannon divergence with base-2 logarithms."""
    return np.asarray(
        [jensenshannon(a, b, base=2.0) ** 2 for a, b in zip(left, right)],
        dtype=float,
    )


def ordering_thirds(values):
    """Return stable early/middle/late labels for a continuous coordinate."""
    series = pd.Series(np.asarray(values, dtype=float))
    ranked = series.rank(method="first")
    return (
        pd.qcut(
            ranked,
            q=3,
            labels=["early", "middle", "late"],
            duplicates="drop",
        )
        .astype(str)
        .to_numpy()
    )

[4]:

RECOVER_DYNAMICS_MAX_ITER = 20

# Load the public scVelo pancreas dataset.
adata = scv.datasets.pancreas()
adata.var_names_make_unique()

# Fit RNA velocity in the original expression/PCA manifold.
scv.pp.filter_and_normalize(adata, min_shared_counts=20)
ensure_pca_neighbors_moments(
    adata,
    n_pcs=30,
    n_neighbors=30,
    preserve_existing_neighbors=False,
)
clear_velocity_outputs(adata, clear_dynamics=True)

# Dynamical velocity is used because the pancreas dataset contains several
# kinetic regimes and non-monotonic terminal trajectories.
scv.tl.recover_dynamics(
    adata,
    max_iter=RECOVER_DYNAMICS_MAX_ITER,
    n_jobs=N_JOBS,
)
scv.tl.velocity(adata, mode="dynamical")
scv.tl.velocity_graph(adata, n_jobs=N_JOBS)
scv.tl.terminal_states(adata)
scv.tl.latent_time(
    adata,
    vkey="velocity",
    root_key="root_cells",
    end_key="end_points",
)
scv.tl.velocity_pseudotime(adata, vkey="velocity")

print(adata)
display(adata.obs["clusters"].astype(str).value_counts().to_frame("n_cells"))

Filtered out 20801 genes that are detected 20 counts (shared).
Normalized count data: X, spliced, unspliced.
computing neighbors
    finished (0:00:10)
computing moments based on connectivities
    finished (0:00:02)
recovering dynamics (using 1/24 cores)

/home/emil/miniforge3/envs/lab-py312/lib/python3.12/multiprocessing/popen_fork.py:66: DeprecationWarning: This process (pid=92183) is multi-threaded, use of fork() may lead to deadlocks in the child.
  self.pid = os.fork()

    finished (0:26:48)
computing velocities
    finished (0:00:08)
computing velocity graph (using 1/24 cores)

/home/emil/miniforge3/envs/lab-py312/lib/python3.12/multiprocessing/popen_fork.py:66: DeprecationWarning: This process (pid=92183) is multi-threaded, use of fork() may lead to deadlocks in the child.
  self.pid = os.fork()

    finished (0:00:15)
computing terminal states
    identified 2 regions of root cells and 1 region of end points .
    finished (0:00:00)
computing latent time using root_cells, end_points as prior
    finished (0:00:03)
AnnData object with n_obs × n_vars = 3696 × 7197
    obs: 'clusters_coarse', 'clusters', 'S_score', 'G2M_score', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size', 'n_counts', 'velocity_self_transition', 'root_cells', 'end_points', 'velocity_pseudotime', 'latent_time'
    var: 'highly_variable_genes', 'gene_count_corr', 'fit_r2', 'fit_alpha', 'fit_beta', 'fit_gamma', 'fit_t_', 'fit_scaling', 'fit_std_u', 'fit_std_s', 'fit_likelihood', 'fit_u0', 'fit_s0', 'fit_pval_steady', 'fit_steady_u', 'fit_steady_s', 'fit_variance', 'fit_alignment_scaling', 'velocity_genes'
    uns: 'clusters_coarse_colors', 'clusters_colors', 'day_colors', 'neighbors', 'pca', 'recover_dynamics', 'velocity_params', 'velocity_graph', 'velocity_graph_neg'
    obsm: 'X_pca', 'X_umap'
    varm: 'loss'
    layers: 'spliced', 'unspliced', 'Ms', 'Mu', 'fit_t', 'fit_tau', 'fit_tau_', 'velocity', 'velocity_u'
    obsp: 'distances', 'connectivities'

	n_cells
clusters
Ductal	916
Ngn3 high EP	642
Pre-endocrine	592
Beta	591
Alpha	481
Ngn3 low EP	262
Epsilon	142
Delta	70

3. Inspect the native velocity field¶

[5]:

native_color = "clusters"
scv.pl.velocity_embedding_stream(
    adata,
    basis="umap",
    color=native_color,
    legend_loc="right margin",
    title="Native RNA-velocity field",
)

computing velocity embedding
    finished (0:00:00)

../_images/tutorials_05_pancreas_multi_scorer_9_1.png

4. Define the furcation, target fate, and future-fate settings¶

[6]:

ROOT = ("Ngn3 high EP", "Pre-endocrine")
FATES = ["Alpha", "Beta", "Delta", "Epsilon"]
OBS_KEY = "clusters"
ORDERING_KEY = "latent_time"
TARGET_FATE = "Beta"
FUTURE_OPTIONS = {
    "effective_horizon": 64,
    "anchor_quantile": 0.90,
    "min_anchor_cells": 10,
    "progression_scale": "rank",
}

5. Real-study metadata or graded controlled demonstration¶

[7]:

def assign_stratified_demo_design(
    adata,
    *,
    annotation_key,
    ordering_key,
    conditions,
    replicates_per_condition,
    condition_key="condition",
    replicate_key="sample_id",
    n_ordering_bins=5,
    random_state=0,
):
    """Create balanced technical pseudo-conditions within state/order strata.

    This helper is for tutorial validation only. It must not replace genuine
    biological condition and replicate metadata in a scientific analysis.
    """
    annotations = adata.obs[annotation_key].astype(str)
    ordering = pd.to_numeric(adata.obs[ordering_key], errors="coerce")
    if ordering.isna().any():
        raise ValueError(f"{ordering_key!r} contains missing or nonnumeric values.")
    strata = pd.Series(index=adata.obs_names, dtype=object)
    for annotation, indices in annotations.groupby(annotations).groups.items():
        indices = pd.Index(indices)
        values = ordering.loc[indices]
        n_bins = min(n_ordering_bins, len(indices), int(values.nunique()))
        if n_bins <= 1:
            bins = pd.Series("0", index=indices)
        else:
            bins = pd.qcut(
                values.rank(method="first"),
                q=n_bins,
                labels=False,
                duplicates="drop",
            ).astype(str)
        strata.loc[indices] = str(annotation) + "::" + bins

    conditions = tuple(map(str, conditions))
    combinations = [
        (condition, f"{condition}_R{replicate + 1}")
        for condition in conditions
        for replicate in range(int(replicates_per_condition))
    ]
    rng = np.random.default_rng(random_state)
    assigned_condition = pd.Series(index=adata.obs_names, dtype=object)
    assigned_replicate = pd.Series(index=adata.obs_names, dtype=object)
    for _, indices in strata.groupby(strata, sort=True).groups.items():
        indices = np.asarray(list(indices), dtype=object)
        rng.shuffle(indices)
        for position, index in enumerate(indices):
            condition, replicate = combinations[position % len(combinations)]
            assigned_condition.loc[index] = condition
            assigned_replicate.loc[index] = replicate

    adata.obs[condition_key] = pd.Categorical(
        assigned_condition,
        categories=list(conditions),
        ordered=True,
    )
    adata.obs[replicate_key] = pd.Categorical(assigned_replicate)
    return pd.DataFrame(
        {
            "stratum": strata.astype(str),
            condition_key: adata.obs[condition_key].astype(str),
            replicate_key: adata.obs[replicate_key].astype(str),
        },
        index=adata.obs_names,
    )


def target_destination_score(adata, result, *, annotation_key, target_fate):
    """Map baseline target-fate probability to every destination cell."""
    score = np.zeros(adata.n_obs, dtype=float)
    selected_indices = adata.obs_names.get_indexer(result.cell_ids)
    if np.any(selected_indices < 0):
        raise ValueError("Result cells could not be aligned to AnnData.")
    fate_index = result.fate_names.index(str(target_fate))
    score[selected_indices] = result.future_fate_contribution[:, fate_index]
    terminal = adata.obs[annotation_key].astype(str).eq(str(target_fate)).to_numpy()
    score[terminal] = 1.0
    return np.clip(score, 0.0, 1.0)


def reweight_root_transitions(
    transition_matrix,
    *,
    source_mask,
    condition_labels,
    replicate_labels,
    destination_score,
    log_shift_by_condition,
    replicate_shift_sd=0.08,
    random_state=0,
):
    """Inject a controlled condition effect without changing graph support.

    Outgoing root-cell transition probabilities are tilted toward destinations
    with high baseline target-fate probability and then renormalized. This is a
    software-validation perturbation, not a model for a particular experiment.
    """
    normalized = scCS.canonicalize_transition_matrix(transition_matrix)
    transition = normalized.matrix.tocsr(copy=True)
    source_mask = np.asarray(source_mask, dtype=bool)
    conditions = np.asarray(condition_labels, dtype=str)
    replicates = np.asarray(replicate_labels, dtype=str)
    destination_score = np.asarray(destination_score, dtype=float)
    if source_mask.shape != (transition.shape[0],):
        raise ValueError("source_mask must align to transition rows.")
    if destination_score.shape != (transition.shape[0],):
        raise ValueError("destination_score must align to transition columns.")

    rng = np.random.default_rng(random_state)
    replicate_effect = {
        replicate: float(rng.normal(0.0, replicate_shift_sd)) for replicate in np.unique(replicates)
    }
    affected = np.zeros(transition.shape[0], dtype=bool)
    for row in np.flatnonzero(source_mask):
        condition = conditions[row]
        if condition not in log_shift_by_condition:
            raise KeyError(f"Missing log shift for condition {condition!r}.")
        start, stop = transition.indptr[row], transition.indptr[row + 1]
        columns = transition.indices[start:stop]
        if len(columns) == 0:
            continue
        shift = float(log_shift_by_condition[condition]) + replicate_effect[replicates[row]]
        multipliers = np.exp(shift * destination_score[columns])
        if np.max(np.abs(multipliers - 1.0)) > 1e-12:
            affected[row] = True
        transition.data[start:stop] *= multipliers
        row_sum = float(transition.data[start:stop].sum())
        if row_sum > 0:
            transition.data[start:stop] /= row_sum
    audit = {
        "affected_root_fraction": float(np.mean(affected[source_mask])),
        "max_row_sum_error": float(
            np.max(np.abs(np.asarray(transition.sum(axis=1)).ravel() - 1.0))
        ),
        "replicate_shift_sd": float(replicate_shift_sd),
        "log_shift_by_condition": dict(log_shift_by_condition),
    }
    return transition, audit

[8]:

original_transition = scCS.get_scvelo_transition_matrix(adata)
if DEMO_MODE:
    assignments = assign_stratified_demo_design(
        adata,
        annotation_key=OBS_KEY,
        ordering_key=ORDERING_KEY,
        conditions=CONDITIONS,
        replicates_per_condition=REPLICATES_PER_CONDITION,
        condition_key=CONDITION_KEY,
        replicate_key=REPLICATE_KEY,
        random_state=SEED,
    )
    baseline_scorer = scCS.SingleScorer(
        adata,
        root=ROOT,
        branches=FATES,
        obs_key=OBS_KEY,
    )
    baseline_scorer.build_embedding(ordering_metric=ORDERING_KEY, verbose=False)
    baseline_scorer.fit(
        transition_matrix=original_transition,
        scoring_mode="future_fate",
        future_fate_options=FUTURE_OPTIONS,
        verbose=False,
    )
    baseline_result = baseline_scorer.score(write_to_adata=False, verbose=False)
    destination_score = target_destination_score(
        adata,
        baseline_result,
        annotation_key=OBS_KEY,
        target_fate=TARGET_FATE,
    )
    labels = adata.obs[OBS_KEY].astype(str)
    root_labels = {ROOT} if isinstance(ROOT, str) else set(ROOT)
    root_mask_full = labels.isin(root_labels).to_numpy()
    analysis_transition, perturbation_audit = reweight_root_transitions(
        original_transition,
        source_mask=root_mask_full,
        condition_labels=adata.obs[CONDITION_KEY].astype(str).to_numpy(),
        replicate_labels=adata.obs[REPLICATE_KEY].astype(str).to_numpy(),
        destination_score=destination_score,
        log_shift_by_condition={"control": 0.0, "low": 0.65, "high": 1.30},
        replicate_shift_sd=0.08,
        random_state=SEED + 1,
    )
    display(pd.Series(perturbation_audit, name="value").to_frame())
else:
    required = {CONDITION_KEY, REPLICATE_KEY}
    missing = sorted(required - set(adata.obs.columns))
    if missing:
        raise KeyError(f"Missing real-study metadata columns: {missing}")
    analysis_transition = original_transition

design_balance = (
    adata.obs.groupby([CONDITION_KEY, REPLICATE_KEY], observed=True)
    .size()
    .rename("n_cells")
    .reset_index()
)
display(design_balance)

	value
affected_root_fraction	1.0
max_row_sum_error	0.0
replicate_shift_sd	0.08
log_shift_by_condition	{'control': 0.0, 'low': 0.65, 'high': 1.3}

	condition	sample_id	n_cells
0	control	control_R1	326
1	control	control_R2	325
2	control	control_R3	320
3	control	control_R4	316
4	low	low_R1	309
5	low	low_R2	305
6	low	low_R3	305
7	low	low_R4	305
8	high	high_R1	302
9	high	high_R2	300
10	high	high_R3	293
11	high	high_R4	290

6. Construct MultiScorer, preflight, and fit the pooled model¶

[9]:

multi = scCS.MultiScorer(
    adata,
    root=ROOT,
    branches=FATES,
    obs_key=OBS_KEY,
    condition_obs_key=CONDITION_KEY,
    replicate_obs_key=REPLICATE_KEY,
    condition_order=CONDITION_ORDER,
)
preflight = multi.preflight(ordering_metric=ORDERING_KEY, check_velocity=True)
preflight.display()
preflight.raise_for_errors()
multi.build_embedding(ordering_metric=ORDERING_KEY)
multi.fit(
    transition_matrix=analysis_transition,
    scoring_mode="future_fate",
    transition_scope="pooled",
    future_fate_options=FUTURE_OPTIONS,
)
results = multi.score_all_conditions(
    population="root",
    min_cells=20,
    min_replicates=4,
)
replicate_table = multi.replicate_table(results)
display(replicate_table)
replicate_table.to_csv(OUTPUT_DIR / "replicate_outcomes.csv", index=False)

	level	code	message	value
0	info	furcation_valid	Root and terminal annotations are valid.	2518.0
1	info	root_cells	Root population contains 1234 cells.	1234.0
2	info	terminal_Alpha	Terminal 'Alpha' contains 481 cells.	481.0
3	info	terminal_Beta	Terminal 'Beta' contains 591 cells.	591.0
4	info	terminal_Delta	Terminal 'Delta' contains 70 cells.	70.0
5	info	terminal_Epsilon	Terminal 'Epsilon' contains 142 cells.	142.0
6	info	ordering_valid	Ordering metric is finite and non-constant amo...	1.0
7	info	ordering_resolution	Root ordering has 1234 unique values across 12...	1.0
8	info	velocity_available	Velocity information is available.	NaN
9	info	condition_control_cells	Condition 'control' contains 1287 cells.	1287.0
10	info	condition_low_cells	Condition 'low' contains 1224 cells.	1224.0
11	info	condition_high_cells	Condition 'high' contains 1185 cells.	1185.0
12	info	condition_control_replicates	Condition 'control' contains 4 biological repl...	4.0
13	info	condition_low_replicates	Condition 'low' contains 4 biological replicates.	4.0
14	info	condition_high_replicates	Condition 'high' contains 4 biological replica...	4.0

[scCS] Scientific star built for 2518 cells in 4 dimensions.
       Root radial clipping: 0.050 low / 0.050 high
       Terminal scientific coordinates: fixed equal-radius simplex vertices (radius=1.000).
scCS FutureFateScoreResult
  Furcation: root -> ['Alpha', 'Beta', 'Delta', 'Epsilon']
  Effective horizon: 64 (gamma=0.984615)
  Cells: 2518 selected; 1234 root
  Root affinity coverage: 1.000
  Root mean future-fate reach: 0.636
  Root mean future-fate entropy: 0.551
  Root mean future-fate specificity: 0.449
  Root mean reach-supported specificity: 0.288
  Root mean unresolved probability: 0.364
  Root mean signed progression: 0.037
  Root future-fate composition: Alpha=0.250, Beta=0.684, Delta=0.000, Epsilon=0.065
  Solver: direct; iterations=1; residual=1.731e-15
[scCS] 'control': 420 root cells; 4 replicates.
[scCS] 'low': 420 root cells; 4 replicates.
[scCS] 'high': 394 root cells; 4 replicates.

	condition	replicate_id	replicate_label	n_cells	n_valid_projection	mean_commitment_strength	mean_directional_entropy	mean_commitment_entropy	mean_directional_specificity	mean_nearest_fate_angle_degrees	...	pairwise_log_commitment_ratio::Alpha::Delta	pairwise_log_commitment_ratio::Delta::Alpha	pairwise_log_commitment_ratio::Alpha::Epsilon	pairwise_log_commitment_ratio::Epsilon::Alpha	pairwise_log_commitment_ratio::Beta::Delta	pairwise_log_commitment_ratio::Delta::Beta	pairwise_log_commitment_ratio::Beta::Epsilon	pairwise_log_commitment_ratio::Epsilon::Beta	pairwise_log_commitment_ratio::Delta::Epsilon	pairwise_log_commitment_ratio::Epsilon::Delta
0	control	control::control_R1	control_R1	105	105	0.630704	0.567961	0.845688	0.432039	NaN	...	6.665268	-6.665268	1.220628	-1.220628	7.612514	-7.612514	2.167874	-2.167874	-5.444640	5.444640
1	control	control::control_R2	control_R2	105	105	0.632746	0.568043	0.844670	0.431957	NaN	...	6.877389	-6.877389	1.293436	-1.293436	7.804059	-7.804059	2.220106	-2.220106	-5.583953	5.583953
2	control	control::control_R3	control_R3	105	105	0.635236	0.558477	0.838429	0.441523	NaN	...	6.575014	-6.575014	1.308039	-1.308039	7.569297	-7.569297	2.302322	-2.302322	-5.266975	5.266975
3	control	control::control_R4	control_R4	105	105	0.629892	0.566188	0.845001	0.433812	NaN	...	6.918626	-6.918626	1.274278	-1.274278	7.828890	-7.828890	2.184541	-2.184541	-5.644349	5.644349
4	low	low::low_R1	low_R1	105	105	0.637504	0.549795	0.833030	0.450205	NaN	...	6.744802	-6.744802	1.338701	-1.338701	7.769484	-7.769484	2.363383	-2.363383	-5.406101	5.406101
5	low	low::low_R2	low_R2	105	105	0.639071	0.537751	0.825533	0.462249	NaN	...	6.711069	-6.711069	1.425986	-1.425986	7.765343	-7.765343	2.480260	-2.480260	-5.285083	5.285083
6	low	low::low_R3	low_R3	105	105	0.638748	0.545685	0.830302	0.454315	NaN	...	6.556119	-6.556119	1.297801	-1.297801	7.594684	-7.594684	2.336365	-2.336365	-5.258319	5.258319
7	low	low::low_R4	low_R4	105	105	0.639690	0.537578	0.824504	0.462422	NaN	...	6.807831	-6.807831	1.396047	-1.396047	7.872467	-7.872467	2.460683	-2.460683	-5.411784	5.411784
8	high	high::high_R1	high_R1	102	102	0.635406	0.543753	0.831095	0.456247	NaN	...	7.181730	-7.181730	1.465726	-1.465726	8.196652	-8.196652	2.480648	-2.480648	-5.716004	5.716004
9	high	high::high_R2	high_R2	100	100	0.641262	0.541285	0.826539	0.458715	NaN	...	6.925487	-6.925487	1.378131	-1.378131	7.990710	-7.990710	2.443353	-2.443353	-5.547357	5.547357
10	high	high::high_R3	high_R3	97	97	0.637808	0.542737	0.829110	0.457263	NaN	...	7.098810	-7.098810	1.361550	-1.361550	8.127643	-8.127643	2.390384	-2.390384	-5.737260	5.737260
11	high	high::high_R4	high_R4	95	95	0.635426	0.548713	0.833513	0.451287	NaN	...	6.855809	-6.855809	1.362233	-1.362233	7.884015	-7.884015	2.390439	-2.390439	-5.493576	5.493576

12 rows × 37 columns

7. Descriptive summaries and condition order¶

[10]:

print("Condition order:", multi.conditions)
print("Condition-order source:", multi.condition_order_source)
summary_rows = []
for condition, condition_result in results.items():
    population = condition_result.population_summary
    row = {
        "condition": condition,
        "n_cells": condition_result.n_cells,
        "n_replicates": condition_result.n_replicates,
        "mean_future_fate_reach": np.nanmean(condition_result.commitment_strength),
        "mean_future_fate_entropy": np.nanmean(condition_result.directional_entropy),
        "mean_future_fate_specificity": np.nanmean(condition_result.directional_specificity),
        "mean_reach_supported_specificity": np.nanmean(condition_result.specific_commitment),
        "mean_signed_progression": np.nanmean(condition_result.progression_velocity),
        "mean_selected_path_coverage": np.nanmean(condition_result.transition_coverage),
        "population_balance_entropy": population.population_balance_entropy,
        "total_future_fate_mass": population.total_mass,
    }
    for fate_index, fate in enumerate(condition_result.fate_names):
        row[f"mean_affinity_{fate}"] = np.nanmean(
            condition_result.directional_affinity[:, fate_index]
        )
        row[f"mean_contribution_{fate}"] = np.nanmean(
            condition_result.commitment_contribution[:, fate_index]
        )
    summary_rows.append(row)
condition_summary = pd.DataFrame(summary_rows)
display(condition_summary)
condition_summary.to_csv(OUTPUT_DIR / "condition_summary.csv", index=False)

Condition order: ['control', 'low', 'high']
Condition-order source: explicit

	condition	n_cells	n_replicates	mean_future_fate_reach	mean_future_fate_entropy	mean_future_fate_specificity	mean_reach_supported_specificity	mean_signed_progression	mean_selected_path_coverage	population_balance_entropy	total_future_fate_mass	mean_affinity_Alpha	mean_contribution_Alpha	mean_affinity_Beta	mean_contribution_Beta	mean_affinity_Delta	mean_contribution_Delta	mean_affinity_Epsilon	mean_contribution_Epsilon
0	control	420	4	0.632145	0.565167	0.434833	0.276102	0.042708	0.998401	0.586336	265.500714	0.260858	0.164079	0.665966	0.421956	0.000299	0.000192	0.072877	0.045917
1	low	420	4	0.638753	0.542702	0.457298	0.294968	0.036288	0.997922	0.557871	268.276365	0.246661	0.155717	0.689875	0.443007	0.000296	0.000192	0.063168	0.039838
2	high	394	4	0.637488	0.544072	0.455928	0.293351	0.032967	0.998343	0.557118	251.170435	0.248823	0.156919	0.689170	0.441431	0.000222	0.000142	0.061784	0.038996

8. Omnibus tests across all fates and core scalar outcomes¶

[11]:

omnibus_affinity = multi.compare_omnibus(
    results,
    metric="future_fate_affinity",
    n_permutations=9999,
    random_state=SEED,
)
omnibus_contribution = multi.compare_omnibus(
    results,
    metric="future_fate_contribution",
    n_permutations=9999,
    random_state=SEED + 1,
)
scalar_omnibus = []
for offset, metric in enumerate(
    (
        "future_fate_reach",
        "future_fate_specificity",
        "reach_supported_specificity",
        "future_fate_entropy",
        "signed_progression",
        "selected_path_coverage",
    )
):
    scalar_omnibus.append(
        multi.compare_omnibus(
            results,
            metric=metric,
            n_permutations=9999,
            random_state=SEED + 10 + offset,
        )
    )
scalar_omnibus = pd.concat(scalar_omnibus, ignore_index=True)
display(omnibus_affinity)
display(omnibus_contribution)
display(scalar_omnibus)
omnibus_affinity.to_csv(OUTPUT_DIR / "omnibus_affinity.csv", index=False)
omnibus_contribution.to_csv(OUTPUT_DIR / "omnibus_contribution.csv", index=False)
scalar_omnibus.to_csv(OUTPUT_DIR / "omnibus_scalar_metrics.csv", index=False)

[scCS] Multi-condition replicate omnibus test; metric='directional_affinity'.
[scCS] Multi-condition replicate omnibus test; metric='mean_commitment_contribution'.
[scCS] Multi-condition replicate omnibus test; metric='commitment_strength'.
[scCS] Multi-condition replicate omnibus test; metric='directional_specificity'.
[scCS] Multi-condition replicate omnibus test; metric='specific_commitment'.
[scCS] Multi-condition replicate omnibus test; metric='directional_entropy'.
[scCS] Multi-condition replicate omnibus test; metric='progression_velocity'.
[scCS] Multi-condition replicate omnibus test; metric='transition_coverage'.

	metric	metric_public	metric_label	fate	fate_a	fate_b	statistic	pvalue	permutation_method	n_permutations	n_conditions	condition_means	pvalue_adj
0	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward Alpha	Alpha	None	None	12.758945	0.0028	monte_carlo	9999	3	{'control': 0.26085766614671557, 'low': 0.2466...	0.0112
1	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward Beta	Beta	None	None	20.323924	0.0061	monte_carlo	9999	3	{'control': 0.665965740857547, 'low': 0.689875...	0.0183
2	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward Delta	Delta	None	None	6.083894	0.0266	monte_carlo	9999	3	{'control': 0.00029927829761174226, 'low': 0.0...	0.0266
3	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward Epsilon	Epsilon	None	None	12.783965	0.0063	monte_carlo	9999	3	{'control': 0.07287731469812567, 'low': 0.0631...	0.0183

	metric	metric_public	metric_label	fate	fate_a	fate_b	statistic	pvalue	permutation_method	n_permutations	n_conditions	condition_means	pvalue_adj
0	mean_commitment_contribution	future_fate_contribution	Future-fate contribution toward Alpha	Alpha	None	None	13.905512	0.0026	monte_carlo	9999	3	{'control': 0.16407922294142538, 'low': 0.1557...	0.0104
1	mean_commitment_contribution	future_fate_contribution	Future-fate contribution toward Beta	Beta	None	None	19.133014	0.0034	monte_carlo	9999	3	{'control': 0.42195614241208695, 'low': 0.4430...	0.0104
2	mean_commitment_contribution	future_fate_contribution	Future-fate contribution toward Delta	Delta	None	None	6.052119	0.0285	monte_carlo	9999	3	{'control': 0.00019185076637403132, 'low': 0.0...	0.0285
3	mean_commitment_contribution	future_fate_contribution	Future-fate contribution toward Epsilon	Epsilon	None	None	13.076618	0.0050	monte_carlo	9999	3	{'control': 0.04591734018222482, 'low': 0.0398...	0.0104

	metric	metric_public	metric_label	fate	fate_a	fate_b	statistic	pvalue	permutation_method	n_permutations	n_conditions	condition_means	pvalue_adj
0	commitment_strength	future_fate_reach	Discounted Fate Reach (DFR)	None	None	None	10.399351	0.0062	monte_carlo	9999	3	{'control': 0.6321445563021112, 'low': 0.63875...	0.0062
1	directional_specificity	future_fate_specificity	Future-Fate Specificity (FFS)	None	None	None	28.062297	0.0046	monte_carlo	9999	3	{'control': 0.43483282412418867, 'low': 0.4572...	0.0046
2	specific_commitment	reach_supported_specificity	Resolved Commitment (RC)	None	None	None	24.892112	0.0032	monte_carlo	9999	3	{'control': 0.27610192642686493, 'low': 0.2949...	0.0032
3	directional_entropy	future_fate_entropy	Future-fate entropy	None	None	None	28.062297	0.0038	monte_carlo	9999	3	{'control': 0.5651671758758113, 'low': 0.54270...	0.0038
4	progression_velocity	signed_progression	Signed Ordering Flux (SOF)	None	None	None	3.023416	0.0964	monte_carlo	9999	3	{'control': 0.042707898703027776, 'low': 0.036...	0.0964
5	transition_coverage	selected_path_coverage	Selected-path coverage	None	None	None	0.282596	0.7349	monte_carlo	9999	3	{'control': 0.99840108962122, 'low': 0.9979221...	0.7349

9. Target-fate post-hoc tests and prespecified linear contrast¶

[12]:

target_omnibus = omnibus_affinity.loc[omnibus_affinity["fate"] == TARGET_FATE].copy()
target_posthoc = multi.compare_posthoc(
    results,
    metric="future_fate_affinity",
    fate=TARGET_FATE,
    omnibus_results=target_omnibus,
    only_significant_omnibus=False,
    n_permutations=9999,
    random_state=SEED + 100,
)
linear_contrast = multi.compare_contrast(
    {"control": -1.0, "low": 0.0, "high": 1.0},
    results,
    metric="future_fate_affinity",
    fate=TARGET_FATE,
    n_permutations=9999,
    random_state=SEED + 200,
)
display(target_omnibus)
display(target_posthoc)
display(linear_contrast)
target_posthoc.to_csv(OUTPUT_DIR / "target_posthoc.csv", index=False)
linear_contrast.to_csv(OUTPUT_DIR / "target_linear_contrast.csv", index=False)

[scCS] Multi-condition post-hoc replicate tests; metric='directional_affinity'.

	metric	metric_public	metric_label	fate	fate_a	fate_b	statistic	pvalue	permutation_method	n_permutations	n_conditions	condition_means	pvalue_adj
1	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward Beta	Beta	None	None	20.323924	0.0061	monte_carlo	9999	3	{'control': 0.665965740857547, 'low': 0.689875...	0.0183

	metric	metric_public	metric_label	fate	fate_a	fate_b	condition_a	condition_b	effect_b_minus_a	pvalue	permutation_method	n_permutations	n_replicates_a	n_replicates_b	pvalue_adj
0	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward Beta	Beta	None	None	control	low	0.023910	0.028571	exact	70	4	4	0.085714
1	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward Beta	Beta	None	None	control	high	0.023164	0.028571	exact	70	4	4	0.085714
2	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward Beta	Beta	None	None	low	high	-0.000745	1.000000	exact	70	4	4	1.000000

	metric	metric_public	metric_label	fate	fate_a	fate_b	contrast	estimate	pvalue	n_permutations	pvalue_adj
0	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward Beta	Beta	None	None	{'control': -1.0, 'low': 0.0, 'high': 1.0}	0.023164	0.0045	9999	0.0045

10. Omnibus, post-hoc, and pairwise-effect visualizations¶

[13]:

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
multi.plot_omnibus_summary(omnibus_affinity, results=results, ax=axes[0])
multi.plot_posthoc_heatmap(target_posthoc, fate=TARGET_FATE, ax=axes[1])
fig.tight_layout()
fig.savefig(OUTPUT_DIR / "omnibus_and_posthoc.png", dpi=200, bbox_inches="tight")
plt.show()

pairwise_figure = multi.plot_pairwise_delta_grid(
    target_posthoc,
    ncols=3,
    annotate=True,
)
pairwise_figure.savefig(OUTPUT_DIR / "pairwise_delta_grid.png", dpi=200, bbox_inches="tight")
plt.show()

../_images/tutorials_05_pancreas_multi_scorer_24_0.png

../_images/tutorials_05_pancreas_multi_scorer_24_1.png

11. Replicate-first target and scalar plots¶

[14]:

fig, axes = plt.subplots(1, 3, figsize=(19, 5))
multi.plot_replicate_outcomes(
    results,
    metric="future_fate_affinity",
    fate=TARGET_FATE,
    ax=axes[0],
)
multi.plot_replicate_outcomes(
    results,
    metric="future_fate_reach",
    ax=axes[1],
)
multi.plot_replicate_outcomes(
    results,
    metric="signed_progression",
    ax=axes[2],
)
fig.tight_layout()
fig.savefig(OUTPUT_DIR / "replicate_outcome_panels.png", dpi=200, bbox_inches="tight")
plt.show()

../_images/tutorials_05_pancreas_multi_scorer_26_0.png

12. Condition-specific star grids¶

[15]:

for color_by, filename in (
    ("population", "star_population.png"),
    (f"future_fate_affinity:{TARGET_FATE}", "star_target_affinity.png"),
    ("future_fate_reach", "star_reach.png"),
    ("future_fate_specificity", "star_specificity.png"),
    ("reach_supported_specificity", "star_supported_specificity.png"),
    ("signed_progression", "star_signed_progression.png"),
):
    figure = multi.plot_star_grid(
        results,
        color_by=color_by,
        population="all",
        ncols=3,
        cmap="coolwarm" if color_by == "signed_progression" else None,
    )
    figure.savefig(OUTPUT_DIR / filename, dpi=200, bbox_inches="tight")
    plt.show()

../_images/tutorials_05_pancreas_multi_scorer_28_0.png

../_images/tutorials_05_pancreas_multi_scorer_28_1.png

../_images/tutorials_05_pancreas_multi_scorer_28_2.png

../_images/tutorials_05_pancreas_multi_scorer_28_3.png

../_images/tutorials_05_pancreas_multi_scorer_28_4.png

../_images/tutorials_05_pancreas_multi_scorer_28_5.png

13. Condition summaries, composition, and trajectory shifts¶

[16]:

fig = plt.figure(figsize=(19, 5.5))
grid = fig.add_gridspec(1, 3, width_ratios=[1.15, 1.0, 1.15])
bar_ax = fig.add_subplot(grid[0, 0])
radar_ax = fig.add_subplot(grid[0, 1], projection="polar")
progression_ax = fig.add_subplot(grid[0, 2])

multi.plot_compare_conditions_bar(
    results,
    metric="future_fate_contribution",
    ax=bar_ax,
)
multi.plot_commitment_vector_radar(
    results,
    metric="commitment_composition",
    ax=radar_ax,
)
multi.plot_trajectory_shift(results, ax=progression_ax)
fig.tight_layout()
fig.savefig(OUTPUT_DIR / "condition_summary_visuals.png", dpi=200, bbox_inches="tight")
plt.show()

../_images/tutorials_05_pancreas_multi_scorer_30_0.png

14. Heatmaps, status composition, and transition coverage¶

[17]:

fig, axes = plt.subplots(1, 3, figsize=(19, 5))
multi.plot_commitment_heatmap(
    results,
    metric="future_fate_contribution",
    level="condition",
    annotate=True,
    ax=axes[0],
)
multi.plot_status_composition(results, ax=axes[1])
multi.plot_transition_coverage(results, ax=axes[2])
fig.tight_layout()
fig.savefig(OUTPUT_DIR / "condition_heatmap_and_qc.png", dpi=200, bbox_inches="tight")
plt.show()

../_images/tutorials_05_pancreas_multi_scorer_32_0.png

15. Ordering trends and display-only rose grids¶

The grid uses scVelo’s standard embedding projection on the displayed X_sccs coordinates. It is visual quality control only. Curved, loop-like, or retrograde branches may therefore point inward or turn on the star, just as they do in the native manifold.

The rose plot includes terminal cells only. Root cells are excluded because their strong forward flow toward the furcation would otherwise be assigned to the nearest terminal display sector and could falsely inflate one branch. The bar mass combines the number of cells and their display-embedded velocity magnitude; it is not a future-fate probability.

[18]:

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
multi.plot_pseudotime_trends(
    results,
    pseudotime_key=ORDERING_KEY,
    metric="future_fate_affinity",
    fate=TARGET_FATE,
    n_bins=8,
    show_replicates=True,
    ax=axes[0],
)
multi.plot_pseudotime_trends(
    results,
    pseudotime_key=ORDERING_KEY,
    metric="signed_progression",
    n_bins=8,
    show_replicates=True,
    ax=axes[1],
)
fig.tight_layout()
fig.savefig(OUTPUT_DIR / "ordering_trends.png", dpi=200, bbox_inches="tight")
plt.show()

display_multi = scCS.MultiScorer(
    adata,
    root=ROOT,
    branches=FATES,
    obs_key=OBS_KEY,
    condition_obs_key=CONDITION_KEY,
    replicate_obs_key=REPLICATE_KEY,
    condition_order=list(CONDITIONS),
)
display_multi.build_embedding(ordering_metric=ORDERING_KEY, verbose=False)
display_multi.fit(
    transition_matrix=analysis_transition,
    scoring_mode="instantaneous",
    transition_scope="pooled",
    verbose=False,
)
display_results = display_multi.score_all_conditions(
    population="root",
    min_cells=1,
    min_replicates=1,
)
rose_figure = display_multi.plot_rose_grid(
    display_results,
    population="root",
    mode="branch",
    n_bins=36,
    ncols=3,
    title="Root fate-directed velocity profiles (instantaneous QC)",
)
rose_figure.savefig(OUTPUT_DIR / "display_velocity_rose_grid.png", dpi=200, bbox_inches="tight")
plt.show()

../_images/tutorials_05_pancreas_multi_scorer_34_0.png

[scCS] 'control': 420 root cells; 4 replicates.
[scCS] 'low': 420 root cells; 4 replicates.
[scCS] 'high': 394 root cells; 4 replicates.

../_images/tutorials_05_pancreas_multi_scorer_34_2.png

16. Gene expression across conditions¶

[19]:

gene_candidates = ["Ins1", "Ins2", "Gcg", "Sst", "Ghrl", "Neurog3"]
genes, missing_genes = resolve_present_genes(adata, gene_candidates)
print("Resolved genes:", genes)
print("Missing genes:", missing_genes)
for gene in genes[:4]:
    figure = multi.plot_gene_expression_star_grid(
        gene,
        results,
        population="all",
        shared_scale=True,
        ncols=3,
    )
    figure.savefig(
        OUTPUT_DIR / f"gene_expression_{gene}.png",
        dpi=200,
        bbox_inches="tight",
    )
    plt.show()

Resolved genes: ['Ins1', 'Ins2', 'Gcg', 'Sst', 'Ghrl', 'Neurog3']
Missing genes: []

../_images/tutorials_05_pancreas_multi_scorer_36_1.png

../_images/tutorials_05_pancreas_multi_scorer_36_2.png

../_images/tutorials_05_pancreas_multi_scorer_36_3.png

../_images/tutorials_05_pancreas_multi_scorer_36_4.png

17. Optional transition-scope sensitivity¶

[20]:

if RUN_SCOPE_SENSITIVITY:
    scope_tables = []
    for scope in ("pooled", "condition", "replicate"):
        scope_multi = scCS.MultiScorer(
            adata,
            root=ROOT,
            branches=FATES,
            obs_key=OBS_KEY,
            condition_obs_key=CONDITION_KEY,
            replicate_obs_key=REPLICATE_KEY,
            condition_order=CONDITION_ORDER,
        )
        scope_multi.build_embedding(ordering_metric=ORDERING_KEY, verbose=False)
        scope_multi.fit(
            transition_matrix=analysis_transition,
            transition_scope=scope,
            scoring_mode="future_fate",
            future_fate_options=FUTURE_OPTIONS,
            verbose=False,
        )
        table = scope_multi.transition_scope_summary(population="root")
        table["transition_scope"] = scope
        scope_tables.append(table)
    scope_summary = pd.concat(scope_tables, ignore_index=True)
    display(scope_summary)
    scope_summary.to_csv(OUTPUT_DIR / "transition_scope_sensitivity.csv", index=False)
else:
    print("Set RUN_SCOPE_SENSITIVITY=True to compare pooled/blocked graphs.")

Set RUN_SCOPE_SENSITIVITY=True to compare pooled/blocked graphs.

18. Optional fail-closed mixed-model sensitivity¶

[21]:

if RUN_MIXED_MODEL_SENSITIVITY:
    mixed = multi.fit_mixed_model(
        metric="future_fate_affinity",
        fate=TARGET_FATE,
        results=results,
        on_invalid="return",
    )
    display(mixed)
    mixed.to_csv(OUTPUT_DIR / "mixed_model_sensitivity.csv", index=False)
else:
    print("Replicate-label permutation is the primary inference.")

Replicate-label permutation is the primary inference.

19. Export analysis objects and statistical tables¶

[22]:

multi.result.write_to_adata(adata)
adata.write_h5ad(OUTPUT_DIR / "multi_scorer_analysis.h5ad")
print("Saved outputs to", OUTPUT_DIR.resolve())

Saved outputs to /home/emil/notebooks/08-tutorials/tutorial_outputs/pancreas_multi

20. Interpretation and real-study adaptation¶

The planned linear contrast tests the prespecified control → low → high trend in Beta future-fate affinity. It is not equivalent to selecting whichever post-hoc comparison appears most significant.

For a real multi-condition study:

encode the intended condition order explicitly;
define the omnibus hypothesis before inspecting pairwise tests;
run post-hoc comparisons with multiplicity correction;
use planned contrasts only when prespecified by the biological design;
report biological replicate counts and effect estimates;
keep fate identity, reach, specificity, and signed progression separate.