🧬 Single-Cell Drug Response Prediction

🎮 Concept: The Perturbation Vector

Before exploring specific tools, visualize the core challenge: moving a cell from a "Control" state to a "Treated" state. Modern AI methods (like CRISP) learn this vector field in high-dimensional space, while methods like State account for resistant subpopulations.

🚀 Foundation Models & Transfer Learning

🧪 Generative Perturbation (The "What If" Engines)

🧬 Mechanism & Interpretability

🎯 Synergy, Prioritization & Screening

📸 Multi-Modal & Morphological

Figure: Multi-Modal Integration — Beyond Transcriptomics

Two paradigms for integrating imaging with drug response: Left (MolPhenix): Contrastive learning aligns Cell Painting images with molecular structures for zero-shot drug retrieval — "find molecules that produce THIS phenotype." Right (GigaTIME): Translates cheap H&E slides into virtual 21-channel multiplex immunofluorescence (IF), predicting immunotherapy response biomarkers (CD8, PD-L1) as accurately as real IF. Bottom: Why multi-modal matters — complementary information at different cost points enables cross-modal prediction as cheap proxy for expensive assays.

📚 Atlases, Databases & Benchmarks

🤖 AI Agents & Clinical Decision Support

🧫 Experimental Platforms & Translational Tools

Figure: From Prediction to Clinic — The Translational Gap

The complete pipeline from computation to clinical trial:

● Computational: scRNA-seq → AI Model (CRISP/State) → Prediction ("Drug X will affect Clone A")

● Experimental: Two validation case studies — Tumor-on-a-Chip (Nat Biotech 2025) enhances CAR-T trafficking 8× through DPP4 inhibition (vildagliptin) via CXCR3-CXCL10/11 axis; MTBBench (NeurIPS 2025) benchmarks AI clinical decision support showing tool-augmented agents reach 69.1% accuracy (+9–11.2% gain).

● Clinical: Phase I/II trials — the validation hierarchy from in silico (publication) through in vitro, ex vivo, animal models, to human trials (FDA approval).

Current frontier: Tahoe-100M (95.6M cells, 379 drugs) and Tumor-on-a-Chip pushing methods to Level 2–3 (ex vivo).

🛠️ Hands-On Practice

The walkthrough below takes a QC'd tumour AnnData object all the way to a per-cell drug-sensitivity readout, along the two routes the field actually uses. Path A is training-free: score every cell against a curated drug signature collection (the beyondcell SSc/PSc idea, run here with decoupler in Python). Path B is the bulk-to-single-cell transfer route made popular by scDEAL and SCAD: train a sensitivity model on GDSC/CCLE cell-line expression with measured IC₅₀ labels, then move it onto unlabelled single cells with an explicit domain-adaptation term. Everything is computational and retrospective — no assay is run, and every conclusion is a hypothesis to be checked against existing data.

Environment & packages

The Python side (scanpy + decoupler + torch) covers signature scoring and the transfer model. beyondcell itself is an R/Seurat package, so it is installed separately if you want its native SSc/PSc collections and BCS ranking plots; the signature tables it ships can also be exported and consumed from Python.

# conda / mamba recommended
conda create -n scdrug python=3.10 -y
conda activate scdrug

pip install scanpy anndata decoupler
pip install torch --index-url https://download.pytorch.org/whl/cpu   # or a CUDA wheel

# optional: beyondcell (R/Seurat) for native SSc + PSc signature collections
# R: install.packages(c("Seurat", "remotes"))
# R: remotes::install_github("cnio-bu/beyondcell")

Hardware. Signature scoring on ~50k cells against a few thousand drug signatures runs in minutes on a laptop (8–16 GB RAM). The transfer model in Path B is small — a two-layer encoder over a few thousand shared genes — and trains on CPU in under a minute; a GPU only matters if you scale to hundreds of drugs or use a large autoencoder pretraining stage as scDEAL does.

Data structures & formats

• AnnData (target domain) — the single-cell object: adata.X normalised/log1p expression, adata.obs per-cell metadata, adata.obsm for the cells × drugs score matrices produced below. No drug labels exist here — that is the whole problem.
• Bulk expression matrix (source domain) — GDSC or CCLE cell lines × genes, typically log₂ TPM or RMA-normalised microarray. A few hundred to ~1,000 rows only, which is why overfitting is a real risk.
• Drug response table — one row per (cell line, drug) pair with LN_IC50 and/or AUC. Binarising at the per-drug median gives the sensitive/resistant labels used for classification; keeping ln(IC₅₀) continuous turns it into regression.
• Drug signature collection — long-format table of drug, gene, weight. beyondcell's SSc signatures are built from genes that separate sensitive from resistant cell lines; its PSc signatures are the transcriptional changes caused by the drug itself (LINCS L1000 / CMap). decoupler consumes exactly this long format as its net.
• Score matrix (BCS-equivalent) — cells × drugs enrichment scores, stored in adata.obsm["mlm_estimate"]; ranking cells within a drug, or drugs within a cell, is the usual downstream move.
• Model artefacts — the shared gene list, the per-domain scaling parameters, and the encoder weights must be saved together. Re-applying a model without the exact gene order and scaler it was trained with silently produces garbage.

Minimal code walkthrough

Path A scores cells against a signature collection with no training at all. Path B trains on GDSC bulk with a held-out cell-line split, then transfers to single cells while an MMD penalty pulls the two domains' latent distributions together — the step that makes bulk-trained weights usable on dropout-heavy, shallow single-cell profiles.

import scanpy as sc
import numpy as np
import pandas as pd
import decoupler as dc
import torch
import torch.nn as nn

# =====================================================================
# PATH A - signature scoring (beyondcell-style): no training required
# =====================================================================

# 1. Load a QC'd tumour object and normalise. Signature scoring expects
#    log-normalised values, NOT raw counts.
adata = sc.read_h5ad("tumor_qc.h5ad")
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# 2. Drug signature collection as a long table: drug -> gene -> weight.
#    SSc = genes separating sensitive from resistant cell lines.
#    PSc = genes changed BY the drug (LINCS/CMap perturbation profiles).
#    Restrict to genes actually measured in this dataset, or the scores
#    are diluted by signature members that can never be detected.
sig = pd.read_csv("lincs_drug_signatures.tsv", sep="\t")   # drug, gene, weight
sig = sig[sig["gene"].isin(adata.var_names)]
sig = sig.groupby("drug").filter(lambda g: len(g) >= 15)   # drop thin signatures

# 3. Score every cell against every drug signature.
dc.run_mlm(mat=adata, net=sig, source="drug", target="gene",
           weight="weight", use_raw=False, verbose=True)
bcs = adata.obsm["mlm_estimate"]            # cells x drugs, the BCS analogue
adata.obs["bortezomib_sig"] = bcs["bortezomib"].values

# 4. Immediately ask whether a "resistant-looking" group is real biology.
#    Overlay depth and cell-cycle phase on the same embedding: low-depth
#    or cycling cells routinely masquerade as a resistant subpopulation.
sc.pl.umap(adata, color=["bortezomib_sig", "total_counts", "phase", "leiden"])

# =====================================================================
# PATH B - bulk -> single-cell transfer (the scDEAL / SCAD idea)
# =====================================================================

# 5. SOURCE domain: GDSC bulk cell lines with measured drug response.
#    IC50 is a POPULATION summary - one number per cell line, per drug.
#    There is no per-cell equivalent anywhere in this pipeline.
expr = pd.read_csv("GDSC_expression.tsv", sep="\t", index_col=0)   # lines x genes
resp = pd.read_csv("GDSC_ic50.tsv", sep="\t", index_col=0)["Bortezomib"]

lines = expr.index.intersection(resp.dropna().index)
genes = expr.columns.intersection(adata.var_names)                 # shared features
Xb = expr.loc[lines, genes].to_numpy("float32")
yb = (resp.loc[lines] < resp.loc[lines].median()).to_numpy("float32")  # 1 = sensitive

# 6. TARGET domain: the single cells, same gene order, zero labels.
Xs = adata[:, genes].X
Xs = np.asarray(Xs.todense() if hasattr(Xs, "todense") else Xs, dtype="float32")

# 7. WHY DOMAIN ADAPTATION IS NOT OPTIONAL. Bulk RNA-seq of a cell line and
#    a single-cell profile of a tumour cell do not share a distribution:
#    scRNA has heavy dropout (most genes read as exact zeros), 10-100x lower
#    effective depth, and a compressed dynamic range. A model fit on bulk and
#    applied naively is evaluated far off its training manifold. Per-domain
#    standardisation is the MINIMUM, and is not sufficient on its own.
Xb = (Xb - Xb.mean(0)) / (Xb.std(0) + 1e-6)
Xs = (Xs - Xs.mean(0)) / (Xs.std(0) + 1e-6)

# 8. Held-out split ON THE SOURCE domain, so the bulk model's own skill is
#    measurable before anything is transferred. Split by cell line - never
#    let the same line appear in both halves.
rng = np.random.default_rng(0)
perm = rng.permutation(len(Xb))
cut = int(0.8 * len(Xb))
tr, te = perm[:cut], perm[cut:]

# 9. Shared encoder + sensitivity head. The adaptation term forces the encoder
#    to produce features whose bulk and single-cell distributions overlap.
#    scDEAL uses a DaNN/MMD penalty; SCAD uses an adversarial discriminator.
#    A linear-kernel MMD is shown here for readability.
class Net(nn.Module):
    def __init__(self, p):
        super().__init__()
        self.enc = nn.Sequential(
            nn.Linear(p, 512), nn.ReLU(), nn.Dropout(0.3),
            nn.Linear(512, 128), nn.ReLU(),
        )
        self.head = nn.Linear(128, 1)

    def forward(self, x):
        z = self.enc(x)
        return z, self.head(z).squeeze(-1)

def mmd(a, b):                      # linear-kernel MMD between batch means
    return ((a.mean(0) - b.mean(0)) ** 2).sum()

net = Net(len(genes))
opt = torch.optim.Adam(net.parameters(), lr=1e-3)
bce = nn.BCEWithLogitsLoss()
Xb_t, yb_t, Xs_t = torch.tensor(Xb), torch.tensor(yb), torch.tensor(Xs)

for epoch in range(300):
    net.train()
    bi = torch.randint(0, len(tr), (64,))        # labelled bulk lines
    si = torch.randint(0, len(Xs_t), (256,))     # UNLABELLED single cells
    zb, logit = net(Xb_t[tr][bi])
    zs, _     = net(Xs_t[si])
    loss = bce(logit, yb_t[tr][bi]) + 1.0 * mmd(zb, zs)   # lambda = 1.0
    opt.zero_grad()
    loss.backward()
    opt.step()

# 10. Sanity-check the SOURCE fit first; only then transfer to the cells.
net.eval()
with torch.no_grad():
    _, lo = net(Xb_t[te])
    acc_bulk = ((torch.sigmoid(lo) > 0.5).float() == yb_t[te]).float().mean()
    _, ls = net(Xs_t)
    adata.obs["pred_sensitivity"] = torch.sigmoid(ls).numpy()
print(f"held-out cell-line accuracy: {acc_bulk:.3f}")

# 11. Validation is necessarily INDIRECT - no per-cell ground truth exists.
#     (a) do the two paths agree? (b) are predicted-resistant cells enriched
#     for plausible resistance programmes? (c) does a known-sensitive vs
#     known-resistant cell line pair, pushed through the same model, come
#     out in the expected order?
print(np.corrcoef(adata.obs["pred_sensitivity"], adata.obs["bortezomib_sig"])[0, 1])
sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon")
sc.pl.umap(adata, color=["pred_sensitivity", "bortezomib_sig", "leiden"])

Common pitfalls & tips

• IC₅₀ has no single-cell ground truth. Every label in GDSC/CCLE is a population-level summary of a bulk viability curve; nothing in a scRNA-seq experiment measures the sensitivity of an individual cell. Per-cell predictions can therefore only be validated indirectly — agreement between independent methods, enrichment of known resistance programmes, or consistency with a downstream clinical/response annotation — never against a per-cell IC₅₀.
• Cell lines do not contain a tumour microenvironment. A model trained on GDSC has never seen stroma, immune infiltrate, hypoxia gradients, or drug-penetration limits. Resistance that is microenvironment-mediated is invisible to it by construction, so restrict scoring to the malignant compartment and do not read the scores of T cells or fibroblasts as drug predictions.
• Batch/domain shift is the dominant error source. GDSC bulk and patient scRNA differ in platform, depth, dropout, and dynamic range at once. Always check that the encoder's latent bulk and single-cell clouds actually overlap (UMAP the latent space, colour by domain); if they sit in separate blobs, the adaptation weight is too low and the transferred predictions are extrapolation.
• Signature collections carry their own context. SSc/PSc and LINCS signatures were derived in specific cell lines at specific doses and time points. A signature built in a haematological line at 10 µM says little about an epithelial tumour at physiological exposure — check the provenance of every collection before interpreting its scores as drug response.
• A "resistant" cluster is often a technical cluster. Low-depth cells score flat on every signature, and cycling cells score high on anything with a proliferation component. Before believing a resistant subpopulation, colour it by total_counts, n_genes_by_counts, mitochondrial fraction, and cell-cycle phase; regress or stratify if it tracks any of them.
• Always run a known-sensitive / known-resistant control pair. Push two cell lines with well-documented, opposite response to your drug through the identical pipeline. If the model cannot separate that pair, its single-cell output is not interpretable — this is the cheapest and most informative sanity check available, and it catches gene-order and scaler mismatches immediately.
• Do not compare raw scores across drugs. Signature scores and model logits are on arbitrary, per-drug scales driven by signature size and label balance. Rank cells within a drug, or z-score each drug across cells, before making any "drug X beats drug Y" statement.