RegVelo Schwann MultiScorer tutorial¶

This notebook demonstrates a three-condition ordered design on the curated Schwann
furcation. In `DEMO_MODE`, a graded ChC-directed transition tilt validates the
statistical API. Gut retrograde progression remains valid biological output.

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/schwann_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¶

The Schwann tutorials use dynamical RNA velocity as the primary model. The dataset is loaded directly from RegVelo and the velocity model is fitted in the original expression/PCA manifold. No tutorial cache is assumed.

The inverse-CytoTRACE coordinate is used only as the validated ordering for this particular furcation; scCS does not prescribe it as a universal pseudotime.

[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]:

import regvelo as rgv

SCHWANN_CLUSTER_ANNOTATION = {
    "15": "Common Progenitor",
    "5": "Common Progenitor",
    "16": "Gut",
    "12": "Gut neuron",
    "8": "ChC",
    "6": "ChC",
}
VELOCITY_MODEL = "dynamical"
RECOVER_DYNAMICS_MAX_ITER = 20


def prepare_schwann_annotations(adata):
    """Reproduce the curated Schwann furcation and inverse-CytoTRACE ordering."""
    if "neighbors" not in adata.uns or "connectivities" not in adata.obsp:
        raise RuntimeError("The original RegVelo neighbor graph is required.")

    if "sccs_leiden" not in adata.obs:
        sc.tl.leiden(
            adata,
            resolution=1.0,
            random_state=0,
            key_added="sccs_leiden",
        )

    leiden = adata.obs["sccs_leiden"].astype(str)
    adata.obs["cell_type_new"] = (
        leiden.map(SCHWANN_CLUSTER_ANNOTATION).fillna("Other").astype("category")
    )

    required = ("Common Progenitor", "Gut", "Gut neuron", "ChC")
    counts = adata.obs["cell_type_new"].astype(str).value_counts().reindex(required, fill_value=0)
    if np.any(counts <= 0):
        raise ValueError(
            f"The curated furcation was not reproduced. Observed counts: {counts.to_dict()}"
        )

    cytotrace = pd.to_numeric(adata.obs["CytoTRACE"], errors="coerce").to_numpy(float)
    finite = np.isfinite(cytotrace)
    lower = float(np.nanmin(cytotrace))
    upper = float(np.nanmax(cytotrace))
    inverse = np.full(adata.n_obs, np.nan, dtype=float)
    inverse[finite] = 1.0 - (cytotrace[finite] - lower) / (upper - lower)
    adata.obs["inverse_cytotrace_pseudotime"] = inverse


# Load the public RegVelo Schwann dataset and reproduce the curated annotations.
adata = rgv.datasets.schwann()
adata.var_names_make_unique()
prepare_schwann_annotations(adata)

# Fit the dynamical model in the original expression/PCA manifold.
# This is the primary Schwann model used by the tutorials.
ensure_pca_neighbors_moments(
    adata,
    n_pcs=30,
    n_neighbors=30,
    preserve_existing_neighbors=True,
)
clear_velocity_outputs(adata, clear_dynamics=True)
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)
prepare_schwann_annotations(adata)

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

running Leiden clustering
    finished (0:00:02)
recovering dynamics (using 1/24 cores)

    finished (0:23:56)
computing velocities
    finished (0:00:10)
computing velocity graph (using 1/24 cores)

    finished (0:00:19)
AnnData object with n_obs × n_vars = 8821 × 1150
    obs: 'plates', 'devtime', 'location', 'n_genes_by_counts', 'total_counts', 'total_counts_ERCC', 'pct_counts_ERCC', 'doublet_scores', 'leiden', 'CytoTRACE', 'Gut_neuron', 'Sensory', 'Symp', 'enFib', 'ChC', 'Gut_glia', 'NCC', 'Mesenchyme', 'Melanocytes', 'SatGlia', 'SC', 'BCC', 'conflict', 'assignments', 'batch', 'initial_size_unspliced', 'initial_size_spliced', 'initial_size', 'n_counts', 'sccs_leiden', 'cell_type_new', 'inverse_cytotrace_pseudotime', 'velocity_self_transition'
    var: 'ERCC', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'n_cells', 'm', 'v', 'n_obs', 'res', 'lp', 'lpa', 'qv', 'highly_variable', 'Accession', 'Chromosome', 'End', 'Start', 'Strand', 'TF', 'means', 'dispersions', 'dispersions_norm', 'velocity_genes', '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', 'fit_r2'
    uns: 'assignments_colors', 'devtime_colors', 'hvg', 'leiden', 'leiden_colors', 'leiden_sizes', 'location_colors', 'log1p', 'neighbors', 'network', 'paga', 'regulators', 'skeleton', 'targets', 'umap', 'sccs_leiden', 'recover_dynamics', 'velocity_params', 'velocity_graph', 'velocity_graph_neg'
    obsm: 'X_diff', 'X_pca', 'X_umap'
    varm: 'loss'
    layers: 'GEX', 'Ms', 'Mu', 'ambiguous', 'matrix', 'palantir_imp', 'scaled', 'spanning', 'spliced', 'unspliced', 'fit_t', 'fit_tau', 'fit_tau_', 'velocity', 'velocity_u'
    obsp: 'connectivities', 'distances'

	n_cells
cell_type_new
Other	6845
ChC	762
Common Progenitor	675
Gut neuron	294
Gut	245

3. Inspect the native velocity field¶

[5]:

native_color = "assignments" if "assignments" in adata.obs else "cell_type_new"
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:01)

../_images/tutorials_06_schwann_multi_scorer_9_1.png

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

[6]:

ROOT = "Common Progenitor"
FATES = ["Gut", "Gut neuron", "ChC"]
OBS_KEY = "cell_type_new"
ORDERING_KEY = "inverse_cytotrace_pseudotime"
TARGET_FATE = "ChC"
FUTURE_OPTIONS = {
    "effective_horizon": 64,
    "anchor_quantile": 0.90,
    "min_anchor_cells": 10,
    "progression_scale": "rank",
}

5. Real-study metadata or graded controlled demonstration¶

The controlled perturbation is intentionally strong enough to recover the predeclared ChC Conditional Fate Affinity effect with replicate-level inference. It validates the software workflow only; it is not a biological treatment model.

[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 Conditional Fate Affinity 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.conditional_fate_affinity[:, 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.90, "high": 1.80},
        replicate_shift_sd=0.03,
        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.03
log_shift_by_condition	{'control': 0.0, 'low': 0.9, 'high': 1.8}

	condition	sample_id	n_cells
0	control	control_R1	750
1	control	control_R2	740
2	control	control_R3	740
3	control	control_R4	735
4	low	low_R1	735
5	low	low_R2	735
6	low	low_R3	735
7	low	low_R4	735
8	high	high_R1	732
9	high	high_R2	730
10	high	high_R3	729
11	high	high_R4	725

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.	1976.000000
1	info	root_cells	Root population contains 675 cells.	675.000000
2	info	terminal_Gut	Terminal 'Gut' contains 245 cells.	245.000000
3	info	terminal_Gut neuron	Terminal 'Gut neuron' contains 294 cells.	294.000000
4	info	terminal_ChC	Terminal 'ChC' contains 762 cells.	762.000000
5	info	ordering_valid	Ordering metric is finite and non-constant amo...	1.000000
6	info	ordering_resolution	Root ordering has 337 unique values across 675...	0.499259
7	info	velocity_available	Velocity information is available.	NaN
8	info	condition_control_cells	Condition 'control' contains 2965 cells.	2965.000000
9	info	condition_low_cells	Condition 'low' contains 2940 cells.	2940.000000
10	info	condition_high_cells	Condition 'high' contains 2916 cells.	2916.000000
11	info	condition_control_replicates	Condition 'control' contains 4 biological repl...	4.000000
12	info	condition_low_replicates	Condition 'low' contains 4 biological replicates.	4.000000
13	info	condition_high_replicates	Condition 'high' contains 4 biological replica...	4.000000

[scCS] Scientific star built for 1976 cells in 3 dimensions.
       Root radial clipping: 0.050 low / 0.050 high
       Terminal scientific coordinates: fixed equal-radius simplex vertices (radius=1.000).
scCS FutureFateScoreResult
  Furcation: Common Progenitor -> ['Gut', 'Gut neuron', 'ChC']
  Effective horizon: 64 (gamma=0.984615)
  Cells: 1976 selected; 675 root
  Root affinity coverage: 1.000
  Root mean future-fate reach: 0.195
  Root mean future-fate entropy: 0.775
  Root mean future-fate specificity: 0.225
  Root mean reach-supported specificity: 0.048
  Root mean unresolved probability: 0.805
  Root mean signed progression: 0.014
  Root future-fate composition: Gut=0.102, Gut neuron=0.267, ChC=0.631
  Solver: direct; iterations=1; residual=2.776e-15
[scCS] 'control': 235 root cells; 4 replicates.
[scCS] 'low': 220 root cells; 4 replicates.
[scCS] 'high': 220 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	...	commitment_composition::Gut neuron	mean_commitment_contribution::ChC	directional_affinity::ChC	commitment_composition::ChC	pairwise_log_commitment_ratio::Gut::Gut neuron	pairwise_log_commitment_ratio::Gut neuron::Gut	pairwise_log_commitment_ratio::Gut::ChC	pairwise_log_commitment_ratio::ChC::Gut	pairwise_log_commitment_ratio::Gut neuron::ChC	pairwise_log_commitment_ratio::ChC::Gut neuron
0	control	control::control_R1	control_R1	60	60	0.188393	0.777648	0.988267	0.222352	NaN	...	0.270647	0.119230	0.633758	0.632879	-1.031545	1.031545	-1.881007	1.881007	-0.849462	0.849462
1	control	control::control_R2	control_R2	60	60	0.194368	0.776844	0.986148	0.223156	NaN	...	0.304324	0.116747	0.626027	0.600646	-1.163903	1.163903	-1.843816	1.843816	-0.679913	0.679913
2	control	control::control_R3	control_R3	60	60	0.201036	0.779418	0.988228	0.220582	NaN	...	0.305277	0.120502	0.619972	0.599404	-1.163991	1.163991	-1.838707	1.838707	-0.674716	0.674716
3	control	control::control_R4	control_R4	55	55	0.181898	0.787638	0.991489	0.212362	NaN	...	0.255182	0.115113	0.629614	0.632842	-0.823693	0.823693	-1.731936	1.731936	-0.908242	0.908242
4	low	low::low_R1	low_R1	55	55	0.180927	0.791686	0.992195	0.208314	NaN	...	0.256081	0.113841	0.628271	0.629212	-0.803104	0.803104	-1.702079	1.702079	-0.898976	0.898976
5	low	low::low_R2	low_R2	55	55	0.191827	0.788753	0.990717	0.211247	NaN	...	0.253664	0.118518	0.622625	0.617839	-0.680103	0.680103	-1.570319	1.570319	-0.890216	0.890216
6	low	low::low_R3	low_R3	55	55	0.202472	0.778593	0.986722	0.221407	NaN	...	0.261032	0.128375	0.635641	0.634035	-0.911329	0.911329	-1.798790	1.798790	-0.887461	0.887461
7	low	low::low_R4	low_R4	55	55	0.221607	0.763342	0.983070	0.236658	NaN	...	0.305283	0.133602	0.624781	0.602879	-1.201214	1.201214	-1.881691	1.881691	-0.680476	0.680476
8	high	high::high_R1	high_R1	55	55	0.204155	0.748747	0.986673	0.251253	NaN	...	0.248604	0.136070	0.660220	0.666500	-1.074437	1.074437	-2.060615	2.060615	-0.986178	0.986178
9	high	high::high_R2	high_R2	55	55	0.187697	0.778163	0.989006	0.221837	NaN	...	0.244471	0.123197	0.644674	0.656361	-0.902284	0.902284	-1.889898	1.889898	-0.987614	0.987614
10	high	high::high_R3	high_R3	55	55	0.186746	0.770911	0.990359	0.229089	NaN	...	0.230736	0.122345	0.642442	0.655143	-0.704007	0.704007	-1.747588	1.747588	-1.043580	1.043580
11	high	high::high_R4	high_R4	55	55	0.201754	0.762375	0.986112	0.237625	NaN	...	0.256094	0.132023	0.646621	0.654379	-1.051004	1.051004	-1.989146	1.989146	-0.938142	0.938142

12 rows × 28 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_Gut	mean_contribution_Gut	mean_affinity_Gut neuron	mean_contribution_Gut neuron	mean_affinity_ChC	mean_contribution_ChC
0	control	235	4	0.191626	0.780233	0.219767	0.046199	0.019010	0.864054	0.806243	45.032201	0.112387	0.019016	0.260319	0.054654	0.627295	0.117957
1	low	220	4	0.199208	0.780593	0.219407	0.047405	0.014196	0.911330	0.811617	43.825839	0.118852	0.021750	0.253319	0.053874	0.627830	0.123584
2	high	220	4	0.195088	0.765049	0.234951	0.050300	0.009409	0.907237	0.769744	42.919407	0.109940	0.018830	0.241570	0.047849	0.648489	0.128409

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 Gut	Gut	None	None	1.139627	0.3526	monte_carlo	9999	3	{'control': 0.11254384023332382, 'low': 0.1188...	0.3526
1	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward Gut neuron	Gut neuron	None	None	4.569582	0.0472	monte_carlo	9999	3	{'control': 0.260113217224082, 'low': 0.253318...	0.0944
2	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward ChC	ChC	None	None	13.374383	0.0054	monte_carlo	9999	3	{'control': 0.6273429425425943, 'low': 0.62782...	0.0162

	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 Gut	Gut	None	None	4.059133	0.0516	monte_carlo	9999	3	{'control': 0.019044123512835267, 'low': 0.021...	0.1548
1	mean_commitment_contribution	future_fate_contribution	Future-fate contribution toward Gut neuron	Gut neuron	None	None	1.027434	0.4123	monte_carlo	9999	3	{'control': 0.054481930552133434, 'low': 0.053...	0.4123
2	mean_commitment_contribution	future_fate_contribution	Future-fate contribution toward ChC	ChC	None	None	2.509004	0.1437	monte_carlo	9999	3	{'control': 0.1178976648097616, 'low': 0.12358...	0.2874

	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	0.403599	0.6909	monte_carlo	9999	3	{'control': 0.19142371887473028, 'low': 0.1992...	0.6909
1	directional_specificity	future_fate_specificity	Future-Fate Specificity (FFS)	None	None	None	2.740926	0.1094	monte_carlo	9999	3	{'control': 0.21961301297411065, 'low': 0.2194...	0.1094
2	specific_commitment	reach_supported_specificity	Resolved Commitment (RC)	None	None	None	0.509707	0.6025	monte_carlo	9999	3	{'control': 0.04608992245698405, 'low': 0.0474...	0.6025
3	directional_entropy	future_fate_entropy	Future-fate entropy	None	None	None	2.740926	0.1128	monte_carlo	9999	3	{'control': 0.7803869870258893, 'low': 0.78059...	0.1128
4	progression_velocity	signed_progression	Signed Ordering Flux (SOF)	None	None	None	0.477786	0.6388	monte_carlo	9999	3	{'control': 0.019021900757523465, 'low': 0.014...	0.6388
5	transition_coverage	selected_path_coverage	Selected-path coverage	None	None	None	6.264084	0.0182	monte_carlo	9999	3	{'control': 0.8642379446788946, 'low': 0.91132...	0.0182

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
2	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward ChC	ChC	None	None	13.374383	0.0054	monte_carlo	9999	3	{'control': 0.6273429425425943, 'low': 0.62782...	0.0162

	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 ChC	ChC	None	None	control	low	0.000487	0.914286	exact	70	4	4	0.914286
1	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward ChC	ChC	None	None	control	high	0.021147	0.028571	exact	70	4	4	0.085714
2	directional_affinity	future_fate_affinity	Conditional Fate Affinity (CFA) toward ChC	ChC	None	None	low	high	0.020660	0.028571	exact	70	4	4	0.085714

	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 ChC	ChC	None	None	{'control': -1.0, 'low': 0.0, 'high': 1.0}	0.021147	0.0041	9999	0.0041

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_06_schwann_multi_scorer_24_0.png

../_images/tutorials_06_schwann_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_06_schwann_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_06_schwann_multi_scorer_28_0.png

../_images/tutorials_06_schwann_multi_scorer_28_1.png

../_images/tutorials_06_schwann_multi_scorer_28_2.png

../_images/tutorials_06_schwann_multi_scorer_28_3.png

../_images/tutorials_06_schwann_multi_scorer_28_4.png

../_images/tutorials_06_schwann_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_06_schwann_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_06_schwann_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_06_schwann_multi_scorer_34_0.png

[scCS] 'control': 235 root cells; 4 replicates.
[scCS] 'low': 220 root cells; 4 replicates.
[scCS] 'high': 220 root cells; 4 replicates.

../_images/tutorials_06_schwann_multi_scorer_34_2.png

16. Gene expression across conditions¶

[19]:

gene_candidates = ["Chga", "Chgb", "Th", "Dbh", "Phox2b", "Tubb3", "Sox10"]
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: ['Chga', 'Chgb', 'Dbh', 'Phox2b', 'Tubb3']
Missing genes: ['Th', 'Sox10']

../_images/tutorials_06_schwann_multi_scorer_36_1.png

../_images/tutorials_06_schwann_multi_scorer_36_2.png

../_images/tutorials_06_schwann_multi_scorer_36_3.png

../_images/tutorials_06_schwann_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/schwann_multi

20. Interpretation and real-study adaptation¶

The planned contrast concerns increasing ChC future-fate affinity. It does not require Gut progression to become positive or imply that all Schwann branches are monotonic.

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.