Integrating a patient cohort with scDEF¶
Five patients, one shared biology, five different batches. The task is the one every cohort study faces: pull out the cell states that are shared across patients without mistaking a batch effect for biology, or erasing real patient-specific signal while correcting for one.
This runs the real scDEF (Ferreira et al.) — a hierarchical
Bayesian factor model — in this notebook's own kernel (Python (iscc-scdef)), so its own
plotting works and the model object is right here to inspect.
It loads a pre-generated cohort; see The analysis dataset and its ground truth for how it was made and what is held back:
| file | what it is | who sees it |
|---|---|---|
counts.csv.gz |
cells × genes, five patients pooled | the tool |
obs.csv |
each cell's patient and batch | the tool (batch only) |
truth_loading.csv.gz |
the shared program dictionary iscc planted | not read here |
truth_activity.csv.gz |
each cell's true program activities | not read here |
What is being tested. Not just that shared biology survives integration, but that both kinds of state show up: the ones every patient uses, and the ones a patient has to itself. Every patient draws from the same program dictionary while its tumour evolves privately — its own clones, its own breakout, its own copy-number landscape — so which programs end up shared and which end up patient-specific is an outcome of that evolution, not something planted.
This notebook shows what scDEF produced. The answer key is not opened here: the planted dictionary, the genes each program spans and every cell's true activities are laid out in The analysis dataset and its ground truth, and the comparison is one a reader makes between the two.
import os, json
import numpy as np, pandas as pd
import matplotlib.pyplot as plt
import anndata as ad
import scdef
from scdef import plotting as sdp
from scdef import tools as sdt
DATA = os.path.join("..", "analysis_data", "cohort")
assert os.path.isdir(DATA), (
f"dataset not found at {DATA}\n"
" generate it first: python validation/make_analysis_data.py --only cohort")
X = pd.read_csv(os.path.join(DATA, "counts.csv.gz"), index_col=0)
obs = pd.read_csv(os.path.join(DATA, "obs.csv"), index_col=0)
meta = json.load(open(os.path.join(DATA, "meta.json")))
adata = ad.AnnData(X.values.astype("float32"))
adata.obs_names = list(X.index); adata.var_names = list(X.columns)
adata.obs["patient"] = pd.Categorical(obs.loc[X.index, "patient"])
adata.obs["batch"] = pd.Categorical(obs.loc[X.index, "batch"])
adata.layers["counts"] = adata.X.copy()
print(f"cohort: {meta['n_patients']} patients, {meta['n_cells']:,} cells, "
f"{meta['n_genes']} genes, {meta['n_programs']} planted programs")
print("cells per patient:", adata.obs['patient'].value_counts().to_dict())
print(f"scDEF {scdef.__version__} running in-process (kernel: Python (iscc-scdef))")
cohort: 5 patients, 1,500 cells, 6000 genes, 6 planted programs
cells per patient: {'P0': 300, 'P1': 300, 'P2': 300, 'P3': 300, 'P4': 300}
scDEF 0.6.1 running in-process (kernel: Python (iscc-scdef))
Fit¶
scDEF is given the batch label and does its own correction — the honest way to test it, rather than pre-correcting the counts by hand and handing it something already cleaned.
# n_factors above the 6 planted programs, so scDEF has room to invent a spurious one if the data
# invites it — giving it exactly as many factors as there are programs would hide that.
sd = scdef.scDEF(adata, counts_layer="counts", n_factors=meta["n_programs"] + 2,
batch_key="batch", seed=1)
sd.fit(n_epoch=300)
print("layer sizes (finest -> coarsest):", sd.layer_sizes)
INFO:scDEF:Found 5 values for `batch` in data: ['batch0' 'batch1' 'batch2' 'batch3' 'batch4']
INFO:scDEF:Each epoch contains 6 batches of size 256
INFO:scDEF:Stopping learning: reached max epochs (n_epoch=300).
INFO:scDEF:Updated adata.var: `gene_scale_0` for batch 0.
INFO:scDEF:Updated adata.var: `gene_scale_1` for batch 1.
INFO:scDEF:Updated adata.var: `gene_scale_2` for batch 2.
INFO:scDEF:Updated adata.var: `gene_scale_3` for batch 3.
INFO:scDEF:Updated adata.var: `gene_scale_4` for batch 4.
INFO:scDEF:Updated adata.obs with layer 0: `L0` and `L0_score` for all factors in layer 0
INFO:scDEF:Updated adata.obsm with layer 0: `X_L0`
INFO:scDEF:Updated adata.uns with layer 0 signatures: `L0_signatures`.
INFO:scDEF:Updated adata.obs with layer 1: `L1` and `L1_score` for all factors in layer 1
INFO:scDEF:Updated adata.obsm with layer 1: `X_L1`
INFO:scDEF:Updated adata.uns with layer 1 signatures: `L1_signatures`.
INFO:scDEF:Updated adata.obs with layer 2: `L2` and `L2_score` for all factors in layer 2
INFO:scDEF:Updated adata.obsm with layer 2: `X_L2`
INFO:scDEF:Updated adata.uns with layer 2 signatures: `L2_signatures`.
INFO:scDEF:Updated adata.obs with layer 3: `L3` and `L3_score` for all factors in layer 3
INFO:scDEF:Updated adata.obsm with layer 3: `X_L3`
INFO:scDEF:Updated adata.uns with layer 3 signatures: `L3_signatures`.
INFO:scDEF:Updated adata.obs with layer 4: `L4` and `L4_score` for all factors in layer 4
INFO:scDEF:Updated adata.obsm with layer 4: `X_L4`
INFO:scDEF:Updated adata.uns with layer 4 signatures: `L4_signatures`.
INFO:scDEF:Updated adata.obs with layer 5: `L5` and `L5_score` for all factors in layer 5
INFO:scDEF:Updated adata.obsm with layer 5: `X_L5`
INFO:scDEF:Updated adata.uns with layer 5 signatures: `L5_signatures`.
layer sizes (finest -> coarsest): [8, 5, 3, 2, 1, 1]
Did it fit? — scdef.pl.qc¶
scDEF ships its own diagnostic panel: the ELBO trace, the biological relevance determination that decides how many factors survive, and the learned-vs-observed cell and gene scales. This is the figure to read before trusting anything downstream.
sdp.qc(sd)
The factor hierarchy — scdef.pl.make_graph¶
scDEF does not infer a flat set of factors: it infers a hierarchy, broad factors at the top splitting into finer ones below. That graph is its signature output and the thing it offers over a flat factor model, so it is what a user would look at first.
fit() builds the graph, but the drawing is done from the model's confident signatures, which are
a separate post-fit step (scdef.tl.set_confident_signatures) — calling make_graph without it
raises a KeyError telling you exactly that.
sdt.set_confident_signatures(sd)
sdp.make_graph(sd, show_all=False)
INFO:scDEF:Updated scDEF graph
Do the patients mix? — scdef.pl.umap¶
The failure mode of any cohort integration is one factor per patient: batch masquerading as biology. scDEF embeds each layer of its hierarchy separately, so this is checkable layer by layer. Coloured by patient and by batch, cells drawn from a shared program should sit together regardless of which patient they came from.
sdt.umap(sd)
sdp.umap(sd, color=["patient", "batch"], figsize=(16, 7))
Which factors does each patient use? — scdef.pl.obs_scores¶
The same question, quantified. obs_scores plots the association between an obs annotation and the
factors of every layer. A shared program is used by all five patients (a full column); a factor
confined to a single row is patient-specific — or that patient's batch.
sdp.obs_scores(sd, ["patient"], mode="fracs", figsize=(11, 3.2))
What to take from it¶
A cohort integration is only useful if the shared biology survives it. Two things worth checking on your own data, both visible above: whether patients mix in the factor space (they should, for a shared program) and whether the factor usage per patient shows structure that is neither uniform nor one-factor-per-patient — the first would mean nothing patient-specific was found, the second that batch was mistaken for biology.
Hold the hierarchy and the per-patient factor usage against
The analysis dataset and its ground truth, which shows the planted
dictionary, the genes in each program, and how each patient actually used them. iscc makes that
comparison possible at all, because the dictionary every patient drew from is recorded.