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.

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.

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

This is the expensive step

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.

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.

How the next step finds these files

Step 10 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