2. Ligand–receptor inference

Language: Python · Input: integrated .h5ad · Output: {KO,WT}_lr_ready_custom.csv

In this step, we will identify ligand–receptor interactions for the knockout (KO) and wild-type (WT) samples separately using LIANA+’s CellPhoneDB implementation. First, we load the integrated data and the consensus resource before running the analysis.

Add the Il1rn interaction manually

The consensus resource does not contain the interaction Il1rn → Il1r1_Il1rap. Because Il1rn is the antagonist ligand at the centre of this study, this interaction must be added manually. Otherwise, the knockout effect cannot be scored.

Score interactions per condition

import scanpy as sc, pandas as pd, liana as li
from liana.method import cellphonedb

adata = sc.read_h5ad(DATA / "bone_marrow_Il1rn_KOvsWT_harmony_integrated.h5ad")

CELLTYPE  = "cluster_names"
CONDITION = "Condition"

consensus_db = li.resource.select_resource("mouseconsensus")
custom_db = pd.concat(
    [consensus_db, pd.DataFrame([{"ligand": "Il1rn", "receptor": "Il1r1_Il1rap"}])],
    ignore_index=True,
)

adata.raw = adata
for cond in set(adata.obs[CONDITION]):
    lr = cellphonedb(
        adata[adata.obs[CONDITION] == cond],
        groupby=CELLTYPE,
        expr_prop=0.1,          # LR expression proportion threshold
        resource=custom_db,
        use_raw=False,
        inplace=False,
        verbose=True,
    )
    lr.to_csv(LR_OUT / f"{cond}_lr_liana_consensus.csv")

Reshape for CrossTalkeR

CrossTalkeR expects a different column layout from the one produced by LIANA+. The following code keeps only statistically significant interactions (cellphone_pvals <= 0.05), renames the required columns, and writes one CrossTalkeR-compatible file for each condition.

import os

for f in os.listdir(LR_OUT):
    if not f.endswith("lr_liana_consensus.csv"):
        continue

    ev = pd.read_csv(LR_OUT / f)
    ev = ev.loc[:, ["ligand", "receptor_complex", "source", "target",
                    "lr_means", "cellphone_pvals"]]
    ev["type_gene_A"] = "Ligand"
    ev["type_gene_B"] = "Receptor"
    ev["gene_A"]  = ev["ligand"]
    ev["gene_B"]  = ev["receptor_complex"]
    ev["MeanLR"]  = ev["lr_means"]
    ev = ev.loc[ev["cellphone_pvals"] <= 0.05, :]

    cond = f[: f.find("_lr_")]
    ev.loc[:, ["source", "target", "type_gene_A", "type_gene_B",
               "gene_A", "gene_B", "MeanLR"]].to_csv(
        LR_OUT / f"{cond}_lr_ready.csv"
    )

The output files (KO_lr_ready_custom.csv and WT_lr_ready_custom.csv) are now ready to be imported into CrossTalkeR.

In the next step, we perform TF activity predictions: Next: Regulon and TF activity