Real-Data Trajectory Analysis

motco.analysis runs the measurement the simulation studies validated on your own data: pooled per-feature standardization, a stage-supervised PLS-DA latent space with the double-CV rank rule, delta / angle / shape between two groups in that space with the net-displacement orientation estimator, and RRPP p-values. analyze_trajectories is the Python entry point and motco analyze the command-line one; both write the same analysis artifact.

Use this entry point rather than assembling motco plsr + motco de by hand. Those commands are a different pipeline: they skip the M-value contract, choose the rank without the parsimony tie-break, measure angle with PC1 by default, and record neither the selected rank, the null summary nor the spectrum.

Input contract

• Blocks. One or more named omic blocks, each a samples × features table with string feature names. Blocks are concatenated in the order you declare them, and every joint feature is named "{block}__{feature}". Rows must equal the metadata's sample ids in the same order; every block needs at least one feature and only finite values.
• Units. Each block declares a unit transform: identity (the default) or beta_to_mvalue (the logit of methylation beta values; values must lie in the open interval (0, 1)). MOTCO does not normalize counts, filter features or impute values. Supply blocks in the continuous units you want measured, for example log-normalized pseudobulk expression. The record states each block's declared units.
• Metadata. A sample-id column, a group column with exactly two labels and a stage column, all named by you. Every group × stage cell must be observed.

• Stage order. Required and explicit. Net displacement runs from the first to the last declared stage, and shape, the spectrum and every table follow this order. Stage order is never inferred by sorting labels, because High < Intermediate < Low alphabetically. Group order is optional and defaults to the labels sorted as strings.

• Covariates (optional). covariate_cols names metadata columns, in model order; see Covariates.

Violations raise TrajectoryAnalysisError (a ValueError) naming the block, labels or cell at fault.

What is measured and recorded

The latent space is fitted once on all samples, with the rank chosen by double cross-validation (the modal rank across outer repeats, the smaller rank winning ties). The defaults are the ones the operating studies ran: cv2_splits=4, cv1_splits=3, n_repeats=5, max_components=20, random_state=1203, each clamped to what the smallest stage can support. forced_components fixes the rank and skips cross-validation for the dimensionality sensitivity analysis; the record then says component_selection: "forced".

record (written as result.json) holds:

Key	Contents
statistics, p_values	delta, angle, shape and their RRPP p-values (shape is null with two stages)
null_summary	per statistic: retained draws, mean, sd, 50/90/95/99 % quantiles of the permutation null
config_spectrum	eigenspectrum of the latent stage-mean configuration, pooled and per group, plus the permutation eigengap summary
latent	selected_lv, component_selection, CV settings, cv_mean_auroc, feature counts, feature_order_signature
design	group and ordered stage levels, cell sizes per group × stage, contrast, analysed sample count
covariates	only when covariates are declared: names, kinds, factor levels and reference levels, design columns, the LS-mean reference row, missing policy, dropped samples, coefficients per latent component
blocks	each block's name, feature count, transform and units
measurement	orientation estimator, selected rank, permutations, seed, runtime
params, software	resolved parameters; MOTCO, numpy, scipy, scikit-learn, pandas and Python versions
attribution	the compact attribution record, or a not_requested marker

The selected rank sits beside the angle result because orientation is only interpretable relative to the space it was measured in. The null summary shows how wide the observed statistic's own null was, and the relative eigengap in the spectrum indicates whether a group's orientation is resolvable at all.

Equivalence to the simulation harness

analyze_trajectories and evaluate_semisynthetic_trajectory(..., integration_method="pls") call the same block-scaling cores and the same fit_pls_latent_space. On identical data, stage order and settings they return the same latent matrix, statistics, p-values, selected rank, null summary and spectrum. tests/test_analysis.py checks this for the default estimator, PC1 and a two-block layer selection.

Attribution

With attribution=AttributionDiagnosticSettings(enabled=True, ...) (CLI: --attribution), frozen-model orientation attribution runs on the fitted pooled estimator and the standardized joint matrix it was fitted to, with your stage order and group labels. The latent space is not refitted. Real data carry no generator truth, so the truth-recovery block is reported unavailable rather than invented. The full tables are written under attribution/.

Analysis artifact

write_analysis_outputs(result, out_dir) writes:

• result.json: the record above.
• latent_scores.csv: samples × pls_*, indexed by sample id.
• pls_loadings.csv: block, feature, pls_* (the PLS x_rotations_, in joint order).
• stage_means.csv: group, stage, pls_*, the fitted latent LS-mean vectors in declared stage order. Pass these to motco.viz.plot_trajectories.
• attribution/*.csv, when attribution ran.

It refuses to replace an existing result.json unless overwrite=True (CLI: --overwrite). The fitted scikit-learn estimator is not pickled, because pickles are version-fragile. The scores, loadings and recorded feature order signature are enough to audit and re-project the fit.

Covariates

Donor covariates such as age at death, PMI or RIN enter the trajectory regression in the latent space, so the stage means, the RRPP reduced model and therefore what is permuted all account for them. Declare them on the dataset:

dataset = TrajectoryDataset(
    ...,
    covariate_cols=["age_at_death", "pmi", "education_level"],
    categorical_covariates=["education_level"],   # optional; non-numeric columns are categorical anyway
    missing_covariates="complete_case",           # default "error"
)

CLI: --covariate age_at_death --covariate pmi --categorical-covariate education_level --missing-covariates complete_case.

• Coding. A numeric (non-boolean) column enters as one design column holding its values. Any other column, or one listed in categorical_covariates (for numeric-coded factors), enters as drop-first indicators name[level], with levels ordered as strings (numerically when every label is an integer); the first level is the reference.
• Model. The full model is group × stage + covariates and the reduced model group + stage + covariates. Covariates are additive: no interaction with group or stage. RRPP permutes the reduced model's residuals as before, so every permuted dataset keeps the fitted covariate effects.
• LS-means. Every group × stage mean is evaluated at the pooled mean of each covariate design column (for a factor, its observed level proportions). The covariate term then adds the same vector to every stage mean, and delta, angle and shape are translation-invariant, so the statistics do not depend on the reference. Only stage_means.csv does, and the record states the reference row.
• Missing values. missing_covariates="error" (the default) rejects any missing covariate value, naming the column and up to ten samples. "complete_case" drops those samples from the metadata and every block before standardization, the latent fit and the design, and records the dropped ids. A deletion that empties a group × stage cell is an error. Covariates are never imputed; blocks must still be complete.
• Guards. A covariate that is also the sample, group or stage column, a constant numeric covariate, a factor with one observed level, a non-finite numeric value, and a model with covariates that is not of full column rank (for example a covariate collinear with group × stage) are rejected with an error naming the covariate(s).

The latent space is covariate-blind, and the adjustment is still exact. PLS-DA is fitted on the standardized features against stage alone, exactly as without covariates, so the latent scores, loadings and selected rank do not change. OLS commutes with a fixed linear map: regressing the scores T = XW on [design, covariates] gives the coefficients of regressing X on the same design, projected by W. The stage means are therefore adjusted exactly as in feature space; only the basis W was chosen without the covariates. If a covariate effect is large and aligned with stage, PLS may spend a component on it. A sensitivity analysis that residualizes the features on the covariates before the fit (outside MOTCO) addresses that case.

Attribution is not available with covariates. Frozen-model attribution works from unadjusted cell means, so a request for both raises an error before any fitting.

Keep the covariate set parsimonious. Small cells (the SEA-AD MTG design has 9–28 donors per cell) and many covariate columns can make the model near-singular or the null unstable. The record's design_columns states the degrees of freedom used.

Type I validation. scripts/covariate_type_i_check.py measures size at the Phase 6 design with covariates confounded with group, stage and their interaction (adjusted vs unadjusted arms on a shared latent fit); see the report. Report covariate-adjusted results only under its covariates_validated decision.

Without covariates the record, the artifact and every output are identical to MOTCO 0.7.0 (except the runtime); the covariates key is absent and the schema version is unchanged (tests/test_analysis_identity.py).

Example: SEA-AD-style design

Two blocks (RNA and ATAC donor pseudobulks), sex as the group, three ADNC stages:

motco analyze \
  --block rna=rna_logcpm.csv --block atac=atac_logcpm.csv \
  --metadata donors.csv --sample-col donor_id --group-col Sex --stage-col ADNC \
  --stage-order "Not AD + Low,Intermediate,High" \
  --permutations 999 --seed 20261015 \
  --out-dir results/motco_mtg_astro

For the dimensionality sensitivity analysis, re-run with --forced-components 2, 3, … into separate output directories.

from motco.analysis import (
    OmicsBlock, TrajectoryAnalysisParams, TrajectoryDataset,
    analyze_trajectories, write_analysis_outputs,
)

dataset = TrajectoryDataset(
    blocks=[OmicsBlock("rna", rna), OmicsBlock("atac", atac)],
    metadata=donors,
    sample_col="donor_id", group_col="Sex", stage_col="ADNC",
    stage_order=["Not AD + Low", "Intermediate", "High"],
)
result = analyze_trajectories(dataset, TrajectoryAnalysisParams(permutations=999, seed=20261015))
write_analysis_outputs(result, "results/motco_mtg_astro")

API

motco.analysis

Real-data trajectory analysis: the validated MOTCO measurement on your own data.

:func:analyze_trajectories runs the production pipeline the simulation harness characterizes — pooled per-feature standardization, a stage-supervised PLS-DA latent space whose rank is chosen by double cross-validation, delta / angle / shape between two groups in that space, and RRPP p-values — on caller-supplied omic blocks and sample metadata. It calls the same cores as :func:motco.simulations.evaluate_semisynthetic_trajectory (block scaling in :mod:motco.simulations.preprocessing, :func:~motco.simulations.evaluation.fit_pls_latent_space, the trajectory design and the record helpers), so the measurement exists once.

Two inputs are load-bearing and therefore explicit:

• Stage order. Net displacement runs from the first to the last declared stage, and shape and the spectrum follow that order. It is never inferred by sorting labels.
• Units. Each block declares its unit transform (identity or beta_to_mvalue). MOTCO does not normalize counts, filter features or impute values; blocks are analysed in the units supplied.

:func:write_analysis_outputs writes the result as plain JSON/CSV (record, latent scores, PLS loadings, latent stage means and, when requested, the attribution tables) — the fitted estimator is never pickled.

TrajectoryAnalysisError

Bases: ValueError

Raised when real-data analysis inputs or parameters are invalid.

Source code in src/motco/analysis.py

class TrajectoryAnalysisError(ValueError):
    """Raised when real-data analysis inputs or parameters are invalid."""

OmicsBlock dataclass

One named omic block: a samples × features table and its unit transform.

Source code in src/motco/analysis.py

@dataclass(frozen=True)
class OmicsBlock:
    """One named omic block: a samples × features table and its unit transform."""

    name: str
    values: pd.DataFrame
    transform: UnitTransform = "identity"

TrajectoryDataset dataclass

Aligned omic blocks plus sample metadata, validated on construction.

Every block's rows must equal metadata[sample_col] in order. The group column must hold exactly two labels; stage_order must list every observed stage label exactly once (labels are compared as strings), and every group × stage cell must be observed.

covariate_cols names metadata columns that enter the trajectory model additively, in declared order. A covariate is categorical when listed in categorical_covariates or when its dtype is not numeric (booleans are categorical); otherwise numeric. missing_covariates decides what a missing covariate value does: "error" rejects it, "complete_case" drops the sample from the metadata and every block before anything else is validated or fitted (dropped_sample_ids lists them). Covariates are never imputed.

Source code in src/motco/analysis.py

@dataclass(frozen=True)
class TrajectoryDataset:
    """Aligned omic blocks plus sample metadata, validated on construction.

    Every block's rows must equal ``metadata[sample_col]`` in order. The group
    column must hold exactly two labels; ``stage_order`` must list every
    observed stage label exactly once (labels are compared as strings), and
    every group × stage cell must be observed.

    ``covariate_cols`` names metadata columns that enter the trajectory model
    additively, in declared order. A covariate is categorical when listed in
    ``categorical_covariates`` or when its dtype is not numeric (booleans are
    categorical); otherwise numeric. ``missing_covariates`` decides what a
    missing covariate value does: ``"error"`` rejects it, ``"complete_case"``
    drops the sample from the metadata and every block before anything else is
    validated or fitted (``dropped_sample_ids`` lists them). Covariates are
    never imputed.
    """

    blocks: Sequence[OmicsBlock]
    metadata: pd.DataFrame
    sample_col: str
    group_col: str
    stage_col: str
    stage_order: Sequence[str]
    group_order: Sequence[str] | None = None
    covariate_cols: Sequence[str] = ()
    categorical_covariates: Sequence[str] = ()
    missing_covariates: MissingCovariatePolicy = "error"
    _dropped_sample_ids: tuple[str, ...] = field(default=(), init=False, repr=False, compare=False)
    _covariate_design: CovariateDesign | None = field(default=None, init=False, repr=False, compare=False)

    def __post_init__(self) -> None:
        object.__setattr__(self, "blocks", tuple(self.blocks))
        object.__setattr__(self, "stage_order", tuple(str(stage) for stage in self.stage_order))
        if self.group_order is not None:
            object.__setattr__(self, "group_order", tuple(str(group) for group in self.group_order))
        object.__setattr__(self, "covariate_cols", tuple(self.covariate_cols))
        object.__setattr__(self, "categorical_covariates", tuple(self.categorical_covariates))
        self._validate_metadata()
        self._validate_covariate_declarations()
        cells_before = self._observed_cells()
        self._apply_missing_covariate_policy()
        self._validate_blocks(self.sample_ids)
        self._build_covariate_design()
        self._validate_design(cells_before)

    @property
    def sample_ids(self) -> list[str]:
        return self.metadata[self.sample_col].astype(str).tolist()

    @property
    def groups(self) -> tuple[str, ...]:
        """Group labels in analysis order (declared, else sorted as strings)."""

        if self.group_order is not None:
            return tuple(self.group_order)
        return tuple(sorted(set(self.metadata[self.group_col].astype(str))))

    @property
    def stages(self) -> tuple[str, ...]:
        return tuple(self.stage_order)

    @property
    def dropped_sample_ids(self) -> tuple[str, ...]:
        """Samples removed by complete-case deletion (empty under ``"error"``)."""

        return self._dropped_sample_ids

    @property
    def covariate_design(self) -> CovariateDesign | None:
        """The covariate design columns on the retained samples, or ``None`` without covariates."""

        return self._covariate_design

    def _observed_cells(self) -> set[tuple[str, str]]:
        return set(zip(self.metadata[self.group_col].astype(str), self.metadata[self.stage_col].astype(str)))

    def _validate_covariate_declarations(self) -> None:
        if self.missing_covariates not in MISSING_COVARIATE_POLICIES:
            raise TrajectoryAnalysisError(
                f"Unknown missing_covariates policy {self.missing_covariates!r}; allowed: "
                f"{list(MISSING_COVARIATE_POLICIES)}."
            )
        covariates = list(self.covariate_cols)
        repeated = sorted({name for name in covariates if covariates.count(name) > 1})
        if repeated:
            raise TrajectoryAnalysisError(f"Covariates must be declared once; repeated: {repeated}.")
        reserved = [name for name in covariates if name in {self.sample_col, self.group_col, self.stage_col}]
        if reserved:
            raise TrajectoryAnalysisError(
                f"Covariates must differ from the sample, group and stage columns; got {reserved}."
            )
        absent = [name for name in covariates if name not in self.metadata.columns]
        if absent:
            raise TrajectoryAnalysisError(f"metadata is missing covariate column(s): {absent}.")
        undeclared = [name for name in self.categorical_covariates if name not in covariates]
        if undeclared:
            raise TrajectoryAnalysisError(
                f"categorical_covariates must also be declared in covariate_cols; not declared: {undeclared}."
            )

    def _apply_missing_covariate_policy(self) -> None:
        if not self.covariate_cols:
            return
        missing = self.metadata[list(self.covariate_cols)].isna()
        if not missing.to_numpy().any():
            return
        sample_ids = self.metadata[self.sample_col].astype(str)
        if self.missing_covariates == "error":
            column = next(name for name in self.covariate_cols if missing[name].any())
            affected = sample_ids[missing[column].to_numpy()].tolist()
            raise TrajectoryAnalysisError(
                f"Covariate {column!r} has missing values for {len(affected)} sample(s): {affected[:10]}. "
                "Pass missing_covariates='complete_case' to drop them; covariates are never imputed."
            )
        drop = missing.any(axis=1).to_numpy()
        keep = ~drop
        object.__setattr__(self, "_dropped_sample_ids", tuple(sample_ids[drop].tolist()))
        object.__setattr__(self, "metadata", self.metadata.loc[keep].reset_index(drop=True))
        retained = []
        for block in self.blocks:
            values = block.values
            if isinstance(values, pd.DataFrame) and values.shape[0] == keep.size:
                values = values.loc[keep]
            retained.append(OmicsBlock(block.name, values, block.transform))
        object.__setattr__(self, "blocks", tuple(retained))

    def _build_covariate_design(self) -> None:
        if not self.covariate_cols:
            return
        try:
            design = build_covariate_matrix(self.metadata, self.covariate_cols, self.categorical_covariates)
        except ValueError as exc:
            raise TrajectoryAnalysisError(str(exc)) from exc
        object.__setattr__(self, "_covariate_design", design)

    def _validate_metadata(self) -> None:
        if not isinstance(self.metadata, pd.DataFrame):
            raise TrajectoryAnalysisError("metadata must be a pandas DataFrame.")
        columns = {"sample_col": self.sample_col, "group_col": self.group_col, "stage_col": self.stage_col}
        missing = [f"{role}={name!r}" for role, name in columns.items() if name not in self.metadata.columns]
        if missing:
            raise TrajectoryAnalysisError(f"metadata is missing column(s): {', '.join(missing)}.")
        if len({self.sample_col, self.group_col, self.stage_col}) != 3:
            raise TrajectoryAnalysisError("sample_col, group_col and stage_col must name three different columns.")
        if self.metadata.shape[0] == 0:
            raise TrajectoryAnalysisError("metadata must contain at least one sample.")
        for name in columns.values():
            if self.metadata[name].isna().any():
                raise TrajectoryAnalysisError(f"metadata column {name!r} contains missing values.")
        sample_ids = self.sample_ids
        duplicated = sorted({sid for sid in sample_ids if sample_ids.count(sid) > 1})
        if duplicated:
            raise TrajectoryAnalysisError(f"metadata sample ids are not unique: {duplicated[:10]}.")

    def _validate_blocks(self, sample_ids: list[str]) -> None:
        if not self.blocks:
            raise TrajectoryAnalysisError("At least one omic block is required.")
        names = [block.name for block in self.blocks]
        repeated = sorted({name for name in names if names.count(name) > 1})
        if repeated:
            raise TrajectoryAnalysisError(f"Block names must be unique; repeated: {repeated}.")
        for block in self.blocks:
            if not isinstance(block.name, str) or not block.name:
                raise TrajectoryAnalysisError(f"Block names must be non-empty strings, got {block.name!r}.")
            if block.transform not in UNIT_TRANSFORMS:
                raise TrajectoryAnalysisError(
                    f"Block {block.name!r} declares unknown transform {block.transform!r}; "
                    f"allowed: {list(UNIT_TRANSFORMS)}."
                )
            values = block.values
            if not isinstance(values, pd.DataFrame):
                raise TrajectoryAnalysisError(f"Block {block.name!r} must be a pandas DataFrame.")
            if values.shape[0] != len(sample_ids):
                raise TrajectoryAnalysisError(
                    f"Block {block.name!r} has {values.shape[0]} rows but metadata has {len(sample_ids)} samples."
                )
            if values.index.astype(str).tolist() != sample_ids:
                raise TrajectoryAnalysisError(
                    f"Block {block.name!r} rows are not aligned to metadata column {self.sample_col!r} "
                    "(same sample ids in the same order are required)."
                )
            if values.shape[1] == 0:
                raise TrajectoryAnalysisError(f"Block {block.name!r} must contain at least one feature.")
            features = values.columns.astype(str).tolist()
            if len(set(features)) != len(features):
                raise TrajectoryAnalysisError(f"Block {block.name!r} has duplicate feature names.")
            try:
                matrix = values.to_numpy(dtype=float)
            except (TypeError, ValueError) as exc:
                raise TrajectoryAnalysisError(f"Block {block.name!r} holds non-numeric values.") from exc
            if not np.all(np.isfinite(matrix)):
                raise TrajectoryAnalysisError(
                    f"Block {block.name!r} contains missing or infinite values; complete matrices are "
                    "required and no imputation is performed."
                )
            if block.transform == "beta_to_mvalue" and not np.all((matrix > 0.0) & (matrix < 1.0)):
                raise TrajectoryAnalysisError(
                    f"Block {block.name!r} declares beta_to_mvalue but holds values outside the open "
                    "interval (0, 1)."
                )

    def _validate_design(self, cells_before: set[tuple[str, str]]) -> None:
        emptied = sorted(cells_before - self._observed_cells())
        if emptied:
            raise TrajectoryAnalysisError(
                f"Complete-case deletion of {len(self._dropped_sample_ids)} sample(s) with missing covariates "
                f"emptied group/stage combination(s): {emptied}."
            )
        observed_stages = sorted(set(self.metadata[self.stage_col].astype(str)))
        declared = list(self.stage_order)
        if len(set(declared)) != len(declared) or set(declared) != set(observed_stages):
            raise TrajectoryAnalysisError(
                f"stage_order must list every observed stage label exactly once; declared {declared}, "
                f"observed {observed_stages}."
            )
        if len(declared) < 2:
            raise TrajectoryAnalysisError(f"At least two stages are required; observed {observed_stages}.")
        observed_groups = sorted(set(self.metadata[self.group_col].astype(str)))
        if len(observed_groups) != 2:
            raise TrajectoryAnalysisError(
                f"Exactly two groups are required in {self.group_col!r}; found {len(observed_groups)}: "
                f"{observed_groups}."
            )
        if self.group_order is not None:
            group_order = list(self.group_order)
            if len(set(group_order)) != len(group_order) or set(group_order) != set(observed_groups):
                raise TrajectoryAnalysisError(
                    f"group_order must list both observed groups exactly once; declared {group_order}, "
                    f"observed {observed_groups}."
                )
        present = self._observed_cells()
        missing = [(group, stage) for group in self.groups for stage in declared if (group, stage) not in present]
        if missing:
            raise TrajectoryAnalysisError(f"metadata has no sample for group/stage combination(s): {missing}.")

groups property

Group labels in analysis order (declared, else sorted as strings).

dropped_sample_ids property

Samples removed by complete-case deletion (empty under "error").

covariate_design property

The covariate design columns on the retained samples, or None without covariates.

TrajectoryAnalysisParams dataclass

Measurement settings; the defaults are the simulation harness's.

forced_components fixes the latent rank (dimensionality sensitivity) and skips cross-validation; None selects it by double CV. seed seeds RRPP exactly as the harness does.

Source code in src/motco/analysis.py

@dataclass(frozen=True)
class TrajectoryAnalysisParams:
    """Measurement settings; the defaults are the simulation harness's.

    ``forced_components`` fixes the latent rank (dimensionality sensitivity) and
    skips cross-validation; ``None`` selects it by double CV. ``seed`` seeds
    RRPP exactly as the harness does.
    """

    permutations: int = 0
    seed: int | None = None
    n_jobs: int = 1
    orientation_estimator: str = "net_displacement"
    forced_components: int | None = None
    cv1_splits: int = 3
    cv2_splits: int = 4
    n_repeats: int = 5
    max_components: int = 20
    random_state: int = 1203
    attribution: AttributionDiagnosticSettings = field(default_factory=AttributionDiagnosticSettings)
    progress: bool = False

    def __post_init__(self) -> None:
        if self.orientation_estimator not in EVALUATION_ORIENTATION_ESTIMATORS:
            raise TrajectoryAnalysisError(
                f"Unknown orientation_estimator {self.orientation_estimator!r}; allowed: "
                f"{list(EVALUATION_ORIENTATION_ESTIMATORS)}."
            )
        if self.permutations < 0:
            raise TrajectoryAnalysisError("permutations must be greater than or equal to 0.")
        if self.n_jobs == 0:
            raise TrajectoryAnalysisError("n_jobs must be -1 or a non-zero integer.")
        settings = self.attribution
        if settings.enabled:
            if settings.bootstrap_replicates < 0:
                raise TrajectoryAnalysisError("attribution.bootstrap_replicates must be greater than or equal to 0.")
            if settings.top_k < 1:
                raise TrajectoryAnalysisError("attribution.top_k must be at least 1.")

    def to_dict(self) -> dict[str, Any]:
        return asdict(self)

TrajectoryAnalysisResult dataclass

The analysis record plus the fitted representation.

record is JSON-safe. latent_scores are the per-sample PLS scores, pls_loadings the joint-feature rotations (block/feature columns first), stage_means the fitted latent LS-mean vector per group × stage in declared stage order, and attribution the full attribution result when it was requested.

Source code in src/motco/analysis.py

@dataclass(frozen=True)
class TrajectoryAnalysisResult:
    """The analysis record plus the fitted representation.

    ``record`` is JSON-safe. ``latent_scores`` are the per-sample PLS scores,
    ``pls_loadings`` the joint-feature rotations (``block``/``feature`` columns
    first), ``stage_means`` the fitted latent LS-mean vector per group × stage
    in declared stage order, and ``attribution`` the full attribution result
    when it was requested.
    """

    record: dict[str, Any]
    latent_scores: pd.DataFrame
    pls_loadings: pd.DataFrame
    stage_means: pd.DataFrame
    attribution: OrientationAttributionResult | None = None

    def to_json(self, **kwargs: Any) -> str:
        return json.dumps(self.record, **kwargs)

analyze_trajectories(dataset, params=None)

Measure two groups' trajectories across ordered stages on real data.

Standardizes every block with pooled per-feature scalers (after its unit transform), concatenates the blocks in declared order as f"{block}__{feature}", builds the stage-supervised PLS-DA latent space (double-CV rank, or params.forced_components), estimates delta, angle (net displacement by default) and shape between the groups in that space and, when params.permutations > 0, their RRPP p-values. Records the selected rank, null summary, stage-mean spectrum, cell sizes, block units and software versions beside the statistics. Optionally runs frozen-model orientation attribution on the fitted estimator.

Declared covariates (dataset.covariate_cols) are appended additively to the full (group × stage) and reduced (group + stage) trajectory models, so RRPP permutes residuals of the reduced model with covariates. Every LS-mean is evaluated at the pooled covariate mean; delta, angle and shape do not depend on that reference. The latent space is fitted without the covariates. Attribution is not available with covariates.

Source code in src/motco/analysis.py

def analyze_trajectories(
    dataset: TrajectoryDataset,
    params: TrajectoryAnalysisParams | None = None,
) -> TrajectoryAnalysisResult:
    """Measure two groups' trajectories across ordered stages on real data.

    Standardizes every block with pooled per-feature scalers (after its unit
    transform), concatenates the blocks in declared order as
    ``f"{block}__{feature}"``, builds the stage-supervised PLS-DA latent space
    (double-CV rank, or ``params.forced_components``), estimates ``delta``,
    ``angle`` (net displacement by default) and ``shape`` between the groups
    in that space and, when ``params.permutations > 0``, their RRPP p-values.
    Records the selected rank, null summary, stage-mean spectrum, cell sizes,
    block units and software versions beside the statistics. Optionally runs
    frozen-model orientation attribution on the fitted estimator.

    Declared covariates (``dataset.covariate_cols``) are appended additively to
    the full (group × stage) and reduced (group + stage) trajectory models, so
    RRPP permutes residuals of the reduced model with covariates. Every LS-mean
    is evaluated at the pooled covariate mean; ``delta``, ``angle`` and
    ``shape`` do not depend on that reference. The latent space is fitted
    without the covariates. Attribution is not available with covariates.
    """

    params = params or TrajectoryAnalysisParams()
    covariates = dataset.covariate_design
    if covariates is not None and params.attribution.enabled:
        raise TrajectoryAnalysisError(
            "Attribution is not available with covariates: it works from unadjusted cell means. "
            "Disable attribution or drop the covariates."
        )
    start = perf_counter()
    blocks = {block.name: block.values for block in dataset.blocks}
    order = [block.name for block in dataset.blocks]
    transforms = {block.name: block.transform for block in dataset.blocks}

    scalers = fit_block_scalers(blocks, transforms)
    standardized = transform_blocks(blocks, scalers, transforms)
    X = concatenate_blocks(standardized, order).reset_index(drop=True)

    groups = list(dataset.groups)
    stages = list(dataset.stages)
    group_codes = {label: str(index) for index, label in enumerate(groups)}
    stage_codes = {label: str(index) for index, label in enumerate(stages)}
    working = pd.DataFrame(
        {
            _GROUP: dataset.metadata[dataset.group_col].astype(str).map(group_codes).to_numpy(),
            _STAGE: dataset.metadata[dataset.stage_col].astype(str).map(stage_codes).to_numpy(),
        }
    )

    integration_params: dict[str, Any] = {
        "cv1_splits": params.cv1_splits,
        "cv2_splits": params.cv2_splits,
        "n_repeats": params.n_repeats,
        "max_components": params.max_components,
        "random_state": params.random_state,
        "n_jobs": params.n_jobs,
    }
    if params.forced_components is not None:
        integration_params["forced_components"] = params.forced_components
    sample_ids = dataset.sample_ids
    try:
        latent = fit_pls_latent_space(
            X,
            working[_STAGE],
            integration_params,
            stage_col=dataset.stage_col,
            sample_ids=sample_ids,
            layer_feature_counts={name: int(frame.shape[1]) for name, frame in blocks.items()},
            layers_metadata={},
            original_scales=_joint_original_scales(scalers, order, X.columns),
            methylation_units=_methylation_units(transforms),
        )
        design = build_simulation_trajectory_design(working, group_col=_GROUP, stage_col=_STAGE)
    except SimulationEvaluationError as exc:
        raise TrajectoryAnalysisError(str(exc)) from exc
    group_levels = [groups[int(code)] for code in design.group_levels]
    stage_levels = [stages[int(code)] for code in design.stage_levels]
    if stage_levels != stages or group_levels != groups:
        raise TrajectoryAnalysisError(
            f"Internal level mapping mismatch: design levels {group_levels}/{stage_levels}, "
            f"declared {groups}/{stages}."
        )
    model_full, model_reduced, ls_means = _covariate_adjusted_design(
        design.model_full, design.model_reduced, design.ls_means, covariates
    )

    observed_deltas, observed_angles, observed_shapes, observed_spectra = estimate_difference(
        Y=latent.matrix,
        model_matrix=model_full,
        LS_means=ls_means,
        contrast=design.contrast,
        return_spectra=True,
        orientation_estimator=params.orientation_estimator,
    )
    config_spectrum = _config_spectrum_block(observed_spectra, group_levels)
    shape_available = len(stage_levels) >= 3
    pair_statistics = _extract_pair_statistics(
        observed_deltas, observed_angles, observed_shapes, shape_available=shape_available
    )

    p_values: dict[str, float] = {}
    null_summary: dict[str, dict[str, float]] = {}
    if params.permutations > 0:
        dist_delta, dist_angle, dist_shape, dist_eigengap = RRPP(
            latent.matrix,
            model_full,
            model_reduced,
            ls_means,
            design.contrast,
            permutations=params.permutations,
            n_jobs=params.n_jobs,
            progress=params.progress,
            seed=params.seed,
            return_eigengaps=True,
            orientation_estimator=params.orientation_estimator,
        )
        null_values = _extract_null_distributions(
            dist_delta, dist_angle, dist_shape, shape_available=shape_available
        )
        p_values = {
            statistic: _empirical_p_value(null_values[statistic], observed)
            for statistic, observed in pair_statistics.items()
            if statistic in null_values and np.isfinite(observed)
        }
        null_summary = _summarize_null_distributions(null_values)
        config_spectrum["permutation_pooled_eigengap"] = _summarize_permutation_eigengaps(dist_eigengap)

    latent_scores = latent.matrix.copy()
    latent_scores.index.name = dataset.sample_col
    artifacts = latent.artifacts
    assert artifacts is not None  # the PLS core always returns its fitted artifacts
    pls_loadings = _loadings_frame(artifacts.model, X.columns, list(latent.matrix.columns))
    betas = np.asarray(estimate_betas(model_full, latent.matrix.to_numpy(dtype=float)), dtype=float)
    stage_means = _stage_means_frame(
        ls_means,
        betas,
        group_levels,
        stage_levels,
        list(latent.matrix.columns),
    )

    attribution_record, attribution_result = _run_attribution(dataset, params, latent, X, transforms)
    cell_sizes = {
        group: {
            stage: int(
                (
                    (dataset.metadata[dataset.group_col].astype(str) == group)
                    & (dataset.metadata[dataset.stage_col].astype(str) == stage)
                ).sum()
            )
            for stage in stage_levels
        }
        for group in group_levels
    }
    record: dict[str, Any] = {
        "schema_version": ANALYSIS_SCHEMA_VERSION,
        "statistics": pair_statistics,
        "p_values": p_values,
        "null_summary": null_summary,
        "config_spectrum": config_spectrum,
        "latent": latent.metadata,
        "design": {
            "group_col": dataset.group_col,
            "stage_col": dataset.stage_col,
            "sample_col": dataset.sample_col,
            "group_levels": group_levels,
            "stage_levels": stage_levels,
            "cell_sizes": cell_sizes,
            "contrast": design.contrast,
            "n_samples": len(sample_ids),
        },
        **(
            {"covariates": _covariate_record(dataset, covariates, betas, list(latent.matrix.columns))}
            if covariates is not None
            else {}
        ),
        "blocks": [
            {
                "name": block.name,
                "n_features": int(block.values.shape[1]),
                "transform": block.transform,
                "units": _TRANSFORM_UNITS[block.transform],
            }
            for block in dataset.blocks
        ],
        "measurement": {
            "orientation_estimator": params.orientation_estimator,
            "selected_lv": latent.metadata["selected_lv"],
            "component_selection": latent.metadata["component_selection"],
            "permutations": params.permutations,
            "seed": params.seed,
            "n_jobs": params.n_jobs,
            "shape_available": shape_available,
            "runtime_seconds": perf_counter() - start,
        },
        "params": params.to_dict(),
        "software": _software_versions(),
        "attribution": attribution_record,
    }
    return TrajectoryAnalysisResult(
        record=_json_safe(record),
        latent_scores=latent_scores,
        pls_loadings=pls_loadings,
        stage_means=stage_means,
        attribution=attribution_result,
    )

write_analysis_outputs(result, out_dir, *, overwrite=False)

Write the analysis artifact to out_dir and return the written paths.

Writes result.json, latent_scores.csv, pls_loadings.csv, stage_means.csv and, when attribution ran, its tables under attribution/. Refuses to touch a directory that already holds result.json unless overwrite is set.

Source code in src/motco/analysis.py

def write_analysis_outputs(
    result: TrajectoryAnalysisResult,
    out_dir: str | Path,
    *,
    overwrite: bool = False,
) -> dict[str, Path]:
    """Write the analysis artifact to ``out_dir`` and return the written paths.

    Writes ``result.json``, ``latent_scores.csv``, ``pls_loadings.csv``,
    ``stage_means.csv`` and, when attribution ran, its tables under
    ``attribution/``. Refuses to touch a directory that already holds
    ``result.json`` unless ``overwrite`` is set.
    """

    directory = Path(out_dir)
    result_path = directory / "result.json"
    if result_path.exists() and not overwrite:
        raise FileExistsError(f"{result_path} already exists; pass overwrite=True to replace the analysis.")
    directory.mkdir(parents=True, exist_ok=True)
    paths: dict[str, Path] = {}
    result_path.write_text(json.dumps(result.record, indent=2, allow_nan=False) + "\n", encoding="utf-8")
    paths["result"] = result_path
    paths["latent_scores"] = directory / "latent_scores.csv"
    result.latent_scores.to_csv(paths["latent_scores"])
    paths["pls_loadings"] = directory / "pls_loadings.csv"
    result.pls_loadings.to_csv(paths["pls_loadings"], index=False)
    paths["stage_means"] = directory / "stage_means.csv"
    result.stage_means.to_csv(paths["stage_means"], index=False)
    if result.attribution is not None:
        for name, path in write_attribution_outputs(result.attribution, directory / "attribution").items():
            paths[f"attribution/{name}"] = path
    return paths