# 9. Ridge models and KO simulations

**Language:** Python · **Input:** `MSCs_diffmap_KO_subset.h5ad`, `MSCs_receptome.csv` ·
**Output:** `coeff_matrix_a10.csv`, `{gene}_delta_matrix.csv` ·
**Runtime:** ~20 min on a normal desktop computer — about 10 min for the model training (9.1)
and about 10 min for the knockout simulations (9.2)

The receptome tells us which receptors act on which target genes in the MSCs, but not how
strongly. In this step we learn that strength from the data, and then use it to ask what the
MSCs would look like if one of the receptors were switched off. We work with the KO condition here;
the WT half of the analysis is the same code on the WT subset.

## 9.1 Training the receptor to target gene model

We fit one ridge regression per target gene, in which the expression of the gene is explained by
the expression of the receptors that the receptome places upstream of it. The receptome takes
the role of the regulon, so we point the function at its `Receptor` and `Target_Gene` columns,
and we keep the package default of α = 10 for the ridge penalty.

```python
import scanpy as sc, pandas as pd
from intratalkerpy.perturbation.mt import calculate_coef_matrix_ridge

COND = "KO"

data      = sc.read_h5ad(PERT / f"MSCs_diffmap_{COND}_subset.h5ad")
receptome = pd.read_csv(PERT / "MSCs_receptome.csv")

outpath = PERT / f"MSCs_{COND}" / "model"
outpath.mkdir(parents=True, exist_ok=True)

data = calculate_coef_matrix_ridge(
    data, receptome,
    alpha      = 10.0,
    source_col = "Receptor",
    target_col = "Target_Gene",
)

coeff = data.uns["regression_coef_matrix"]
coeff.to_csv(outpath / "coeff_matrix_a10.csv")
```

The models come back inside the object. `data.uns["regression_coef_matrix"]` is a gene × gene
matrix that is non-zero only at the receptor → target pairs of the receptome; this is the matrix
we perturb below, which is why we write it out. Two further tables come with it:
`data.uns["regression_statistics"]` reports R², RMSE and the number of upstream receptors for
every target gene, and is where we look when we want to know how well a gene is explained by its
receptors, and `data.uns["regression_coef_statistics"]` summarises the fitted coefficients per
receptor.

Running the same call on `MSCs_diffmap_WT_subset.h5ad` gives us the WT model. The two
coefficient matrices are not the same, and that difference is what makes the Il1r1 effect
condition-specific in [step 10](10_receptor_ranking.md).

## 9.2 Simulating the receptor knockouts

We now switch off one receptor at a time. For each receptor we set its expression to zero in
every cell, let that change spread through the coefficient matrix for five propagation steps,
and translate the resulting shift in expression into a movement in the diffusion map of
[step 8](08_trajectory.md). The projection follows the scheme of Velocyto and CellOracle: the
expression shift of a cell is correlated with the shifts towards its neighbours, which turns it
into a displacement vector in the embedding — one delta matrix per receptor.

We simulate the 24 receptors that the publication reports. They are a subset of the MSC
receptome which were associated to the KO condition.

:::{admonition} This is the expensive step
:class: important

On 8 cores the 24 receptors take only about 10 minutes on a Desktop PC, but only due to the low cell count. If you are running
multiple thousand of cells in the analysis we recommend using a HPC to reduce the run time.  
:::

```python
import scanpy as sc, pandas as pd, numpy as np
import scipy.sparse as sparse
from intratalkerpy.perturbation import mt

COND = "KO"

data    = sc.read_h5ad(PERT / f"MSCs_diffmap_{COND}_subset.h5ad")
coeff   = pd.read_csv(PERT / f"MSCs_{COND}" / "model" / "coeff_matrix_a10.csv", index_col=0)
outpath = PERT / f"MSCs_{COND}" / "simulation"
outpath.mkdir(parents=True, exist_ok=True)

genes = [
    "Ripk1", "Lrp6",  "Mertk", "Tgfbr1", "Axl",   "Fgfr1",
    "Epha2", "Pdgfra", "Tgfbr2", "Pdgfrb", "Tnfrsf21", "Fzd1",
    "Insr",  "Tlr4",  "Egfr",  "Lrp5",   "Acvr1", "Itgb1",
    "Pld1",  "Igf1r", "Dip2a", "Smo",    "Fgfr2", "Il1r1",
]

df_exp = data.to_df()

for i, gene in enumerate(genes, 1):
    print(f"{i}/{len(genes)}  {gene}", flush=True)
    np.random.seed(15037)

    perturbed = df_exp.copy()
    perturbed[gene] = 0                    # knockout

    simulated = mt.simulation_of_perturbation(
        gem              = df_exp,
        simulation_input = perturbed,
        coef_matrix      = coeff,
        n_propagation    = 5,
    )
    sparse.save_npz(outpath / f"{gene}_simulated_matrix.npz",
                    sparse.csr_matrix(simulated))

    delta_embedding = mt.project_perturbation_in_embedding(
        anndata          = data,
        original_matrix  = df_exp,
        perturbed_matrix = simulated,
        reduction_name   = "X_diffmap",
        n_cpu            = 8,
    )
    pd.DataFrame(delta_embedding).to_csv(outpath / f"{gene}_delta_matrix.csv")
```

Every receptor leaves two files behind: the simulated expression matrix, and
`{gene}_delta_matrix.csv` with one displacement vector per cell. The next step reads
those delta files back in.

:::{admonition} How the next step finds these files
:class: note

[Step 10](10_receptor_ranking.md) reads the CSV files whose name ends in `_delta_matrix` and
takes the receptor name from whatever comes before it, so anything else we keep in the folder is
ignored. Receptor complexes are safe here: the whole name up to the suffix is used, which makes
`Il1r1_Il1rap_delta_matrix.csv` the receptor `Il1r1_Il1rap`. With a different naming scheme, pass
it to `differential_pseudotime_analysis` as `file_suffix`, or `file_suffix=None` to read every
CSV in the folder and take the file name itself as the receptor.
:::

Next: [Differential pseudotime and receptor ranking](10_receptor_ranking.md)