Skip to content
PIASO

SCALAR: ligand-receptor interaction analysis

piaso.tl.runSCALAR asks, for every ordered pair of cell types and every known ligand–receptor pair: is the ligand specific to the sender and the receptor specific to the receiver, more than expected by chance?

The word doing the work is specific. SCALAR scores cell-type specificity, not mean expression, so a ligand every cell expresses cannot win. The null is built from genes matched to each query gene by expression, so a highly expressed pair does not win either.

import numpy as np
import pandas as pd
import piaso, cosg


piaso.settings.set_figure_params(style="cell")

1. Data and an interaction database

SEA-AD middle temporal gyrus, 20,000 nuclei, 24 annotated subclasses:

adata = piaso.data.load_dataset("sea_ad_mtg_20k")
piaso.tl.infog(adata, layer="UMIs", n_top_genes=3000)

SCALAR needs an interaction database. piaso.data fetches CellChatDB on demand and caches it in ~/.piaso/data, the same way the genome and motif references work — nothing downloads until you ask:

lr = piaso.data.load_lr_database("human")     # or "mouse"
lr.shape, [c for c in lr.columns][:4]
((2951, 28), ['interaction_name', 'pathway_name', 'ligand', 'receptor'])

The mouse table is load_lr_database("mouse") — 3,105 pairs. Any table with a ligand column and a receptor column also works, so a curated in-house list can be passed straight to lr_pairs=.

What makes CellChatDB worth using rather than a bare pair list is the annotation column, which says by what mechanism each pair interacts:

lr["annotation"].value_counts()
Secreted Signaling       1200
Non-protein Signaling     746
ECM-Receptor              515
Cell-Cell Contact         490

Those are four different biological questions. A secreted ligand can act at a distance; an ECM-receptor pair means one cell is building matrix the other adheres to; cell-cell contact requires the two cells to touch. Pooling them answers none of the three cleanly, and §6 splits them.

2. The specificity matrix

SCALAR takes a genes × cell-types matrix of specificity scores. COSG produces exactly that, which means the same scores drive annotation and interaction analysis — one definition of “specific” across the analysis, not two.

cosg.cosg(adata, key_added="cosg", groupby="Subclass",
          n_genes_user=adata.n_vars, mu=1, remove_lowly_expressed=False)


names = pd.DataFrame(adata.uns["cosg"]["names"])
scores = pd.DataFrame(adata.uns["cosg"]["scores"])
spec = pd.DataFrame(0.0, index=adata.var_names, columns=names.columns)
for c in names.columns:
    spec.loc[names[c].values, c] = scores[c].values
spec.shape
(36601, 24)

n_genes_user=adata.n_vars asks COSG for every gene, not a top-N — the matrix needs a score for any gene that might appear in the database.

3. Run it

res = piaso.tl.runSCALAR(adata, specificity_matrix=spec, lr_pairs=lr,
                         ligand_col="ligand", receptor_col="receptor",
                         annotation_col="Subclass", n_permutations=1000,
                         layer="infog", random_seed=42)
res.shape, list(res.columns)
((1691712, 8),
 ['ligand', 'receptor', 'sender', 'receiver', 'interaction_score',
  'p_value', 'p_value_fdr', 'nlog10_p_value_fdr'])

1.69 million tested interactions in 57 seconds — 2,937 usable LR pairs × 576 ordered cell-type pairs, each with a 1,000-permutation null. FDR is applied per sender–receiver pair independently, which matters: correcting globally across 576 pairs would bury every result in a single enormous multiple-testing penalty for questions that were never one family.

It reports what it dropped:

Filtered out 14 LR pairs that were not found in both the AnnData object
and the specificity matrix.

4. What comes out

Top significant interactions, and each one is checkable:

ligandreceptorsenderreceiverscoreFDR
COL1A2SDC4VLMCAstrocyte0.5680.020
COL1A2ITGA2VLMCOligodendrocyte0.5610.040
COL6A2SDC4VLMCAstrocyte0.2920.020
COL1A2ITGA9VLMCL2/3 IT0.2690.019
ANGPTL4CDH5AstrocyteEndothelial0.2410.039
ICAM2ITGALEndothelialMicroglia-PVM0.2100.019
NPYGPR83Sst ChodlL2/3 IT0.2050.026
IL33IL1RL1AstrocyteEndothelial0.1880.039
VIPSCTRVipVip0.1770.048

Read down that column of senders. VLMC — vascular leptomeningeal cells — dominate the collagen signalling, which is what they do: they build the extracellular matrix of the meninges and perivascular space. Astrocyte → endothelial appears twice, through ANGPTL4/CDH5 and IL33/IL1RL1, which is the neurovascular unit. ICAM2 → ITGAL is the canonical endothelium-to-immune adhesion pair, and the receiver is microglia. VIP → SCTR is autocrine on VIP interneurons.

None of that was supplied. The inputs were a count matrix, a subclass label and a public interaction table.

Significant interactions by sender and receiver

5. Read it at the level of the question

Individual pairs are noisy; the aggregate is what most analyses want:

sig = res[res["p_value_fdr"] < 0.05]
sig.groupby(["sender", "receiver"]).size().sort_values(ascending=False).head(6)
sender       receiver
L2/3 IT      Endothelial    108
Pax6         L2/3 IT        101
VLMC         L2/3 IT         98
VLMC         L5 IT           97
Pax6         Sncg            96
Endothelial  L2/3 IT         94

A caution worth stating plainly: interaction counts track cell-type abundance and specificity sharpness, not just biology. L2/3 IT is the largest subclass here (4,826 nuclei), and it appears on both sides of the busiest pairs. Compare pairs of comparable size, or normalise, before concluding that one cell type “talks more”.

6. Plot one pair properly

A ranked table does not show why an interaction scored. plotLigandReceptorInteraction does: a bar for the interaction score, and directly beneath it the two specificity scores that produced it — the ligand’s in the sender, the receptor’s in the receiver. An interaction with a tall bar and a pale square underneath is one gene carrying the pair.

It wants four derived columns, all of them one line:

sig = res[res["p_value_fdr"] < 0.05].copy()


sig["CellTypeXCellType"] = sig["sender"] + "@" + sig["receiver"]
sig["ligandXreceptor"] = sig["ligand"] + "-->" + sig["receptor"]


ann = (lr.drop_duplicates(subset=["ligand", "receptor"])
         .set_index(["ligand", "receptor"])[["annotation", "pathway_name"]])
sig = sig.join(ann, on=["ligand", "receptor"])


sig["ligand_specificity"] = [spec.at[r.ligand, r.sender] for r in sig.itertuples()]
sig["receptor_specificity"] = [spec.at[r.receptor, r.receiver] for r in sig.itertuples()]
pairs = ["VLMC@Astrocyte", "Astrocyte@Endothelial", "Endothelial@Microglia-PVM"]


piaso.pl.plotLigandReceptorInteraction(
    interactions_df=sig, specificity_df=spec, cell_type_pairs=pairs,
    ligand_receptor_sep="-->", top_n=30, y_max=0.6,
    heatmap_cmap="Purples", shared_legend=True,
    fig_width=20, fig_height_per_pair=6)
Interaction scores and their specificity components

Three panels, three different mechanisms, and the colour of each bar says which:

  • VLMC → Astrocyte is blue — ECM-receptor, top to bottom. Collagens onto syndecan-4 and CD44.
  • Astrocyte → Endothelial is green — secreted. ANGPTL4 and IL33, the neurovascular signalling axis.
  • Endothelial → Microglia-PVM is pale blue — cell-cell contact. ICAM2 onto the integrins ITGAM/ITGB2/ITGAL, which is leukocyte adhesion, plus CSF3→CSF3R.

Set y_max from the data. It defaults to 10, and these scores top out at 0.57, so leaving the default draws every bar as a sliver against an empty axis: y_max=float(np.ceil(sig["interaction_score"].max() * 10) / 10).

7. Split by mechanism

annotation is a filter, and filtering before plotting turns a mixed panel into a controlled comparison. Same sender, two receivers, one mechanism:

ecm = sig[sig["annotation"] == "ECM-Receptor"]
len(ecm)
2378
piaso.pl.plotLigandReceptorInteraction(
    interactions_df=ecm, specificity_df=spec,
    cell_type_pairs=["VLMC@Astrocyte", "VLMC@L2/3 IT"],
    ligand_receptor_sep="-->", top_n=30, y_max=0.6,
    heatmap_cmap="Purples", shared_legend=True,
    fig_width=20, fig_height_per_pair=6)
ECM-receptor interactions from VLMC to two receivers

The ligands barely change — VLMC sends COL1A2, COL6A2, COL1A1, LAMC3 to both. The receptor changes completely. Astrocytes receive them on SDC4 and CD44; L2/3 IT neurons receive the same collagens on ITGA9 and ITGA3. One matrix source, two adhesion systems, and the split only appears because the mechanism class was held constant.

That is the argument for keeping annotation rather than reducing the database to a bare pair list.

8. One pair in full detail

plotLigandReceptorLollipop puts the ligand’s specificity above the axis and the receptor’s below, so the two halves of each interaction are read at once, and adds a third variable as circle size.

Circle size comes from col_circle_size, default avg_log2FC — any per-row effect size you care about. Here it is the ligand’s fold change between high-pathology and not-AD donors, which turns a description of the cortex into a question about the disease:

adc = adata.obs["Overall AD neuropathological Change"].astype(str)
hi, lo = adc.isin(["High", "Intermediate"]), adc.isin(["Not AD", "Low"])
# CP10K means per sender cell type, high-pathology vs not-AD
sig["avg_log2FC"] = ...     # see the note below
piaso.pl.plotLigandReceptorLollipop(
    sig, cell_type_pairs=["Endothelial@Microglia-PVM"], top_n=30,
    col_cell_type_pair="CellTypeXCellType", sort_by_category=True,
    fig_height_per_pair=4.5, fig_width=16, vertical_layout=False,
    background_colors=True, logfc_range=1, base_circle_size=30,
    color_labels_by_annotation=True)
Endothelial to microglia interactions, ligand above and receptor below

sort_by_category=True groups by mechanism, so the panel reads left to right as contact → secreted → ECM → non-protein. The four leftmost stems are the tall ones: ICAM2→ITGAM, ICAM2→ITGB2, ICAM2→ITGAL, TGM2→ADGRG1, all cell-cell contact, all with substantial specificity on both ends. Further right the stems collapse onto the axis — one side of the pair is carrying the score alone, which is exactly the case where a nominally significant interaction is worth less than its p-value suggests.

Without a col_circle_size column the function warns and draws every circle at one size; the plot is still correct, it just has one variable fewer.

The fold change used above:

X = adata.layers["UMIs"].tocsr()
cp10k = X.multiply(1e4 / np.asarray(X.sum(1)).ravel()[:, None]).tocsr()
gidx = {g: i for i, g in enumerate(adata.var_names)}


lfc = {}
for ct in spec.columns:
    m = (adata.obs["Subclass"] == ct).values
    a = np.asarray(cp10k[m & hi.values].mean(0)).ravel()
    b = np.asarray(cp10k[m & lo.values].mean(0)).ravel()
    lfc[ct] = np.log2((a + 0.1) / (b + 0.1))


sig["avg_log2FC"] = [lfc[r.sender][gidx[r.ligand]] for r in sig.itertuples()]

14,600 nuclei from high or intermediate pathology against 5,400 from not-AD or low. The strongest AD-up ligand among significant interactions is endothelial SEMA3G (log2FC 3.7), signalling onto neuronal and microglial plexins and neuropilins — but note its interaction scores are small (0.02–0.10). A large fold change and a small specificity score are different statements, and the lollipop is showing you both at once rather than letting one stand in for the other.

9. Narrow the question

Scoring all 576 pairs is cheap, but if you have a hypothesis, say so — the p-values are then spent on the comparison you care about:

res = piaso.tl.runSCALAR(adata, specificity_matrix=spec, lr_pairs=lr,
                         ligand_col="ligand", receptor_col="receptor",
                         annotation_col="Subclass",
                         sender_cell_types=["Astrocyte", "VLMC", "Endothelial"],
                         receiver_cell_types=["Microglia-PVM"],
                         n_permutations=1000, layer="infog")

Parameters worth knowing

parameterwhat it changes
n_permutationssize of the null (default 1000). The smallest reachable p-value is 1/(n+1), so 1,000 bottoms out at 0.000999 — raise it if you need finer resolution at the top.
n_nearest_neighborshow control genes are matched to query genes (default 30).
sender_cell_types / receiver_cell_typesrestrict the comparison.
prefilter_fdrdrop pairs that cannot reach significance before permuting. On by default.
layerwhich matrix supplies expression. infog here.
random_seedthe null is sampled; fix it.
  • Gene set scoring (PIASOscore) — the same control-set idea, applied to gene sets rather than gene pairs.
  • LARIS is the spatial counterpart: when cells have coordinates, physical proximity constrains which interactions are possible, and LARIS uses it. Same databases, per-cell answer.
  • Datasets and genome references — the rest of what piaso.data fetches on demand.