> Status: draft. Check package versions, input layers, and state choices before use.

Input: an `.h5ad` with unique cell and gene IDs, `spliced` and `unspliced`
layers, a cluster label, and sample metadata.

| Analysis | Depends on | Main output |
|---|---|---|
| RNA velocity | Spliced and unspliced counts | Velocity field and latent-time estimates |
| CellRank fate inference | Velocity and connectivity kernels | Macrostate and fate probabilities |
| Network activity | Expression and a matched gene network | Cell-level regulator scores |
| Metabolic-task scoring | Expression, metadata, and package resources | Task-activity estimates |

::: {.callout-warning}
These are model-based estimates. Fate probabilities are not observed lineage
outcomes; regulator scores are not direct binding measurements; metabolic-task
scores are not measurements of flux.
:::

## 1. Load and check the input

```{python}
from pathlib import Path

import anndata as ad
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scanpy as sc
import scvelo as scv

adata = ad.read_h5ad("input/processed_single_cell.h5ad")
assert adata.obs_names.is_unique and adata.var_names.is_unique
assert {"spliced", "unspliced"}.issubset(adata.layers)
assert "sample_id" in adata.obs and "cluster" in adata.obs
assert adata.layers["spliced"].shape == adata.shape
assert adata.layers["unspliced"].shape == adata.shape
print(adata)
```

The velocity layers must contain the kinetic counts for the same cells and
genes as `adata`. An integrated expression matrix does not replace them.

## RNA velocity and latent time

This branch estimates transcriptional direction from spliced and unspliced
counts. It requires kinetic layers from the same cells and genes; integrated
expression values do not replace them.

```{python}
scv.settings.set_figure_params(dpi=120, frameon=False)
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
sc.tl.pca(adata, n_comps=30, random_state=1234)
sc.pp.neighbors(adata, n_neighbors=30, n_pcs=30, random_state=1234)
sc.tl.umap(adata, random_state=1234)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
scv.tl.recover_dynamics(adata, n_jobs=1)
scv.tl.velocity(adata, mode="dynamical")
scv.tl.velocity_graph(adata, n_jobs=1)
scv.tl.velocity_pseudotime(adata)
scv.tl.recover_latent_time(adata)
scv.tl.velocity_confidence(adata)
scv.pl.velocity_embedding_stream(adata, basis="umap", color="cluster")
```

Inspect the kinetic fit and velocity confidence. Pseudotime and latent time are
model-derived orderings, not elapsed time.

## CellRank states and fate probabilities

This branch combines velocity and cell-connectivity kernels. Results depend on
the kernel blend and number of states; inspect sensitivity to plausible choices.

```{python}
from cellrank.estimators import GPCCA
from cellrank.kernels import ConnectivityKernel, VelocityKernel

velocity_kernel = VelocityKernel(adata).compute_transition_matrix()
connectivity_kernel = ConnectivityKernel(adata).compute_transition_matrix()
# Example blend; check whether conclusions change with other weights.
kernel = 0.7 * velocity_kernel + 0.3 * connectivity_kernel

estimator = GPCCA(kernel)
estimator.compute_schur(method="krylov")
estimator.compute_macrostates(n_states=3, cluster_key="cluster")
estimator.predict_initial_states(allow_overlap=True)
estimator.predict_terminal_states(allow_overlap=True)
estimator.compute_fate_probabilities(solver="gmres", n_jobs=1)
fate = estimator.fate_probabilities.to_df()
fate.index.name = "cell_id"
drivers = estimator.compute_lineage_drivers(seed=1234)

estimator.plot_macrostates(basis="X_umap", which="terminal")
estimator.plot_fate_probabilities(basis="X_umap", same_plot=False)
print(fate.describe())
print(drivers.head())
```

## Network activity with decoupler

This analysis is independent of RNA velocity. The network TSV needs `source`,
`target`, and `weight` columns, and its gene identifiers must match the
expression data.

```{python}
import decoupler as dc

network = pd.read_csv("input/regulatory_network.tsv", sep="\t")
assert {"source", "target", "weight"}.issubset(network.columns)
dc.mt.ulm(data=adata, net=network, layer=None, tmin=5)
activity = dc.pp.get_obsm(adata=adata, key="score_ulm")
activity.index.name = "cell_id"
print(activity.describe().T.head())
```

ULM scores summarize expression patterns relative to a prior-knowledge
network; they do not measure transcription-factor binding or pathway flux.

## Metabolic-task scores with scCellFie (optional)

This analysis is independent of velocity and network scoring. It needs a
Scanpy neighbor graph and the metadata columns below. Check the current package
API and organism resource before enabling the example.

```{python}
#| eval: false

import sccellfie

adata = sccellfie.run_sccellfie_pipeline(
    adata,
    organism="human",
    n_counts_col="n_counts",
    process_by_group=False,
    neighbors_key="neighbors",
    n_neighbors=10,
    batch_key="sample_id",
    threshold_key="sccellfie_threshold",
    smooth_cells=True,
    alpha=0.33,
    chunk_size=5000,
    save_folder=None,
    save_filename=None,
)["adata"]

report = sccellfie.reports.generate_report_from_adata(
    adata.metabolic_tasks,
    group_by="cluster",
    tissue_col="sample_id",
    feature_name="metabolic_task",
)
```

These scores estimate metabolic-task activity; they do not measure metabolic
flux.

## Save the object and analysis tables

```{python}
output = Path("output")
output.mkdir(parents=True, exist_ok=True)
adata.write_h5ad(output / "single_cell_analysis.h5ad")
fate.reset_index().to_csv(output / "cellrank_fate.tsv", sep="\t", index=False)
drivers.to_csv(output / "cellrank_drivers.tsv", sep="\t", index=False)
activity.reset_index().to_csv(output / "decoupler_ulm.tsv", sep="\t", index=False)
```