# Advanced Features This page documents optional input files and parameters that activate advanced modes of CardamomOT. All optional files are placed in `my_project/Data/` (or set directly in `adata.obs`). ## Project directory structure (full) ``` my_project/ ├── Data/ │ ├── data.h5ad # required │ ├── gene_list.txt # optional — restrict to a gene subset │ ├── stimulus_schedule.txt # optional — custom stimulus values per timepoint │ ├── stimulus_schedule_simul.txt # optional — overrides schedule for simulation only │ ├── ref_network.csv # optional — prior interaction graph (sparsity mask) │ ├── basal_init.npy / .csv # optional — warm-start for basal parameters │ ├── inter_init.npy / .csv # optional — warm-start for interactions │ ├── basal_ref.npy / .csv # optional — regularisation target for basal │ ├── inter_ref.npy / .csv # optional — regularisation target for interactions │ ├── KO_OV_inference.txt # optional — per-sample KO/OV prior (requires dataset_id) │ ├── KO_OV_simulate.txt # optional — in-silico perturbations to simulate │ └── transition_rates.csv # optional — cell-type transition cost matrix for OT └── cardamomOT/ # generated by the pipeline ├── basal.npy # (n_samples, G, n_networks) ├── inter.npy # (G, G, n_networks) └── basal_ref_mask.npy # (n_samples, G) bool — from KO_OV_inference.txt ``` --- ## Gene subset (`Data/gene_list.txt`) By default `select_DEgenes_and_split.py` selects differentially expressed genes automatically. To fix the gene set manually, place a plain-text file with one gene name per line: ``` CHGA STMN2 PHOX2B HAND2 ... ``` ```bash cardamomot step select_DEgenes_and_split -i my_project -s full -c 0 --mean-forcing 0.5 --force-basins 1.0 --temporal-basins 1 ``` When `gene_list.txt` is present, those genes are **added** to the final selection — they are not a replacement. The full behaviour depends on the `-c` flag: - **`-c 1`**: genes variable across time, conditions, and cell types are selected automatically, then the genes in `gene_list.txt` are appended (provided they are non-constant in time). - **`-c 0`**: no automatic differential selection is performed; only the genes in `gene_list.txt` that are non-constant in time are kept. --- ## Read-depth correction (`infer_rd`) When cells have heterogeneous sequencing depths, run this optional step **before** gene selection: ```bash cardamomot step infer_rd -i my_project ``` This identifies Poisson-like genes (Chronocell method, Fang et al. 2024) and stores a per-cell read-depth factor in `adata.obs['rd']`. Subsequent steps (`infer_mixture`, `infer_network_structure`) use it automatically when the column is present. **Key parameters:** `--var_threshold 1.2` (Poisson variance threshold), `--min_mean 0.1`. > Requires at least 1 000 genes. Below this threshold the step is silently skipped. --- ## Prior interaction network (`prepare_reference_network`) To bias network inference toward known biology, first build a prior network from public databases (OmniPath, etc.): ```bash cardamomot step prepare_reference_network -i my_project -d 3 ``` This generates `Data/ref_network.csv`, a gene × gene binary (or weighted) interaction matrix. The `--prior` parameter then controls how strongly edges absent from this prior are penalised. | `--prior` value | Effect | |---|---| | `1.0` (default) | No penalisation — all edges equally possible | | `0.5` | Soft penalisation — absent edges discouraged | | `0.0` | Hard mask — only edges present in the prior are allowed | **`--prior` must be passed to both steps** — it constrains *inference* in `infer_network_structure` and the *simulation reference network* in `infer_network_simul`. Always use the same value in both. ```bash # Easiest: let the pipeline handle it (prior network with default depth 3) cardamomot pipeline -i my_project -s full --ref --prior 0.5 # With a shallower graph (faster): cardamomot pipeline -i my_project -s full --ref --ref-depth 1 --prior 0.5 # Or step by step — pass --prior to both: cardamomot step infer_network_structure -i my_project -s full --prior 0.5 cardamomot step infer_network_simul -i my_project -s full --prior 0.5 ``` --- ## Custom stimulus schedule (`Data/stimulus_schedule.txt`) By default one stimulus is `0` at the first timepoint and `1` at all subsequent ones. To define a custom schedule, create a tab/space-separated file where **rows = timepoints** and **columns = stimulus channels**: ``` # stimulus_schedule.txt 0.0 0.0 1.0 0.0 1.0 1.0 ``` - Values between 0 and 1 are allowed (partial stimulus strength). - If fewer rows than timepoints are provided, the last row is repeated. - To use a *different* schedule for simulation (e.g. a novel protocol), place `Data/stimulus_schedule_simul.txt`; it takes priority over `stimulus_schedule.txt` during `simulate_network` and `simulate_network_KOV`. The `--stimulus` parameter controls how strongly the stimulus regulates genes in the inferred network: | `--stimulus` value | Effect | |---|---| | `1.0` (default) | Full stimulus influence | | `0.0` | Stimulus has no regulatory influence | --- ## Multiple experimental samples (`dataset_id`) If your experiment contains several biological conditions that share a gene regulatory network but have **different basal transcription rates** (e.g. different cell lines, donors, or perturbation backgrounds): ```python adata.obs['dataset_id'] = ... # string or integer label per cell ``` CardamomOT then: - Builds one set of **per-sample basal parameters** `θ_basal(s)` for each sample. - Runs a single joint optimisation over shared interaction weights and all per-sample basals simultaneously. - Saves `basal.npy` with shape `(n_samples, G, n_networks)`. **Keeping per-sample basals close (`constrain_basal_uniform`):** set `model.constrain_basal_uniform = λ` (e.g. `λ = 100–1000`) to add an L2 penalty pushing each sample's basals toward their common mean. Genes with a non-zero `basal_ref` entry for a given sample (e.g. a KO prior) are excluded from the penalty for that sample. --- ## Regularisation targets (`basal_ref`, `inter_ref`) Place reference arrays in `Data/` to anchor the optimiser: | File | Shape | Role | |---|---|---| | `basal_ref.npy` / `.csv` | `(G,)` or `(n_samples, G, n_networks)` | Regularisation target for basal parameters | | `inter_ref.npy` / `.csv` | `(G, G)` or `(G, G, n_networks)` | Regularisation target for interactions | | `basal_init.npy` / `.csv` | same as `basal_ref` | Warm-start values (not regularised) | | `inter_init.npy` / `.csv` | same as `inter_ref` | Warm-start values for interactions | CSV files must have gene names matching `adata.var_names` as row and column indices. --- ## KO/OV inference with per-sample priors (`Data/KO_OV_inference.txt`) When `dataset_id` is present, you can encode prior knowledge about which genes are knocked out or overexpressed in each sample: ``` # KO_OV_inference.txt (tab-separated, header required) sample_id KO OV wt 0 0 ko_CHGA CHGA 0 ov_STMN2 0 STMN2 ko_CHGA_ov_STMN2 CHGA STMN2 ``` - `sample_id` must match `adata.obs['dataset_id']`. - `KO` genes are forced to `basal = −100` (silent during inference). - `OV` genes are forced to `basal = +100` (always active). - Multiple genes per cell: `CHGA,POSTN`. During inference, this replaces `basal_ref` for the affected genes/samples and saves `cardamomOT/basal_ref_mask.npy` for use by `simulate_network_KOV`. ```bash cardamomot step infer_network_structure -i my_project -s full # KO_OV_inference.txt is read automatically when present ``` --- ## In-silico perturbation simulation (`Data/KO_OV_simulate.txt`) After training, simulate arbitrary KO/OV combinations by defining them in: ``` # KO_OV_simulate.txt (tab-separated, header required) KO OV none none CHGA STMN2 POSTN S100B,STMN2 CHGA-80 STMN2-60 ``` Each row produces one independent simulation saved as `cardamomOT/adata_sim_KO_*_OV_*_stim*.h5ad`. **Partial KO/OV (`GENE-X` syntax):** append `-X` (0 < X < 100) for a partial perturbation of strength X%: | Mode | Creation rate factor | Effect | |---|---|---| | `KO X%` | `× (1 − X/100)` | Reduced production → lower steady-state | | `OV X%` | `× 1/(1 − X/100)` | Increased production → higher steady-state | Output files use the `pctX` suffix to distinguish partial from complete perturbations (`KO_CHGApct80_OV_STMN2pct60`). ```bash cardamomot step simulate_network_KOV -i my_project -s full ``` Visualise results with: ```python from CardamomOT import plot_results_sim_kov, compare_cell_types p, stim, prior = "my_project/", 1.0, 1.0 compare_cell_types(p, "KO_CHGA_OV_none", split="full", stim=stim, prior=prior) plot_results_sim_kov(p, "KO_CHGA_OV_none", stim=stim, prior=prior) ``` --- ## Population dynamics: proliferation, death, and transition rates ### Per-cell proliferation/death (`adata.obs`) ```python adata.obs['prolif_rate'] = ... # net proliferation rate per cell adata.obs['death_rate'] = ... # net death rate per cell ``` When both columns are present, the OT marginals are corrected so that faster-growing cells carry proportionally more weight as trajectory sources. This implements demographically corrected OT (as in Waddington OT, Schiebinger et al. 2019). ### Cell-type transition rates (`Data/transition_rates.csv`) To bias the OT cost toward biologically plausible transitions: ``` # transition_rates.csv — rows = source type at t1, cols = target type at t2 ,TypeA,TypeB,TypeC TypeA , 0.3, 0.1, 0.01 TypeB , 0.05, 0.2, 0.05 TypeC , 0.01, 0.05, 0.3 ``` Row/column names must match `adata.obs['cell_type']`. At each consecutive timepoint pair, transition probabilities are computed as `exp(rate × Δt)` and rescaled so that the mean weight per row equals 1 — transitions with weight > 1 become cheaper (preferred), weight < 1 become more expensive (penalised). --- ## Proliferation-aware simulation (`--proliferation`) CardamomOT can learn and simulate net proliferation rates (R = birth − death) directly from the inferred optimal-transport couplings, without requiring external annotations. ### How it works **Learning R from the OT coupling.** During trajectory inference (`infer_network_structure`), at each consecutive timepoint pair (t → t+1), the row marginals of the optimal coupling encode the effective mass transported out of each source cell. The optimal net growth rate per cell is recovered via the self-consistent zero-mean inversion of the WOT marginal formula: ``` m_n = (1/N) · exp(R_n · Δt/2) / C (WOT: source marginal, normalised) R_n = (2/Δt) · (log(m_n) − mean_k(log(m_k))) ``` The factor `2/Δt` comes from the WOT convention where both the source **and** the target marginals are corrected by `exp(R · Δt/2)` each — the formula recovers the true net rate R from the source marginal alone. With balanced OT (default) all `m_n = 1/N` so `R_n = 0` exactly; with unbalanced OT (`model.unbalanced_reg > 0`) they vary and encode genuine differential growth. The resulting `R_opt` array (one value per cell per timepoint transition) is always saved to `cardamomOT/data_R_opt.npy`. **Training the ProliferationMLP.** `infer_network_simul --proliferation` reads `data_R_opt.npy` and trains a lightweight two-hidden-layer MLP (64 units, Tanh activation, MSE loss): ``` R(P) : protein levels → net proliferation rate ``` The network is trained on `(prot[:, ns:], R_opt)` pairs *after* any recomputation of protein trajectories, so it captures the optimal R consistent with the final inferred network. The weights are saved to `cardamomOT/prolif_network.pt`. **Branching simulation.** `simulate_network --proliferation` (and `simulate_network_KOV --proliferation`) loads the MLP and applies a **trapezoidal forward+backward resampling** after each simulated interval — mirroring the structure of the WOT correction used during inference (which corrects both source and target marginals): ``` log_weight_n = (R(P_start_n) + R(P_end_n)) / 2 · Δt ``` `P_start` and `P_end` are the protein states of cell n at the beginning and end of the interval. The N cells are then **resampled multinomially** (with replacement, weighted by `exp(log_weight_n)`). This avoids the catastrophic variance of Poisson sampling (which would kill ~37 % of cells per step even at R = 0) and is consistent with the inference-time coupling. ### Pipeline usage ```bash # Step 4 — adapt parameters and learn R cardamomot step infer_network_simul -i my_project -s full --proliferation # Step 5 — simulate with branching cardamomot step simulate_network -i my_project -s full --proliferation # KO/OV perturbations with branching cardamomot step simulate_network_KOV -i my_project -s full --proliferation ``` Or directly: ```bash python infer_network_simul.py -i my_project -s full --proliferation python simulate_network.py -i my_project -s full --proliferation ``` ### Notes - `--proliferation` is **off by default**. `infer_network_structure` always computes and saves `R_opt` as a by-product regardless — enabling the flag later has no cost. - The flag must be passed to **both** `infer_network_simul` (trains the MLP) and `simulate_network` / `simulate_network_KOV` (uses it). Passing it only to the simulation scripts while `prolif_network.pt` is absent prints a warning and falls back to standard simulation. - Resampling is currently applied to the PDMP stochastic mode only; ODE mode ignores it silently. - For balanced OT the estimated R values are exactly zero for all cells; meaningful proliferation signals require unbalanced OT (`model.unbalanced_reg > 0`) or prior annotations via `adata.obs['prolif_rate']` / `adata.obs['death_rate']` (see section above). - The trapezoidal rule `(R_start + R_end)/2` naturally handles cells whose protein levels — and thus growth rate — change significantly within the interval.