Quality control from cross-modal disagreement: mis-paired cells and doublets

Quality control from cross-modal disagreement: mis-paired cells and doublets#

Open in Colab

Experimental. A QC idea tested on problems we create ourselves, so that every claim has ground truth. How well it transfers to the artifacts in your own data is an empirical question.

In a paired assay both modalities of a barcode come from the same cell, so a trained UniVI model should place a cell’s RNA and ATAC at the same point of the latent space. Barcodes whose modalities disagree more than those of well-behaved cells are candidates for problems: two cells sharing a barcode, RNA and chromatin from different cells, ambient contamination, or damaged nuclei.

This notebook, on 10x Multiome PBMCs:

  1. trains a model and scores held-out cells with four measures of disagreement or unusualness;

  2. builds ground truth by corrupting held-out cells in known ways: swapping in another cell’s ATAC (from a different or the same cell type), and merging two cells into a doublet;

  3. measures which score detects which problem (AUROC);

  4. turns a score into conformal p-values against clean calibration cells and flags barcodes with Benjamini–Hochberg, which controls the false discovery rate (Bates et al., 2023), then checks the realized false discovery rate and power on the constructed data;

  5. applies the same procedure to uncorrupted held-out cells and describes what it flags.

import sys

if "google.colab" in sys.modules:
    %pip install -q "univi[tutorials]>=1.2" "pandas==2.2.3"

import anndata as ad
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scanpy as sc
import scipy.sparse as sp
import torch
from IPython.display import display
from sklearn.metrics import roc_auc_score
from sklearn.neighbors import NearestNeighbors

import univi
import univi.datasets as uds
from univi import ModalityConfig, TrainingConfig, UniVIConfig, UniVIMultiModalVAE, UniVITrainer
from univi.evaluation import cross_modal_predict, encode_fused_adata_pair
from univi.preprocessing import ATACPreprocessor, RNAPreprocessor, split_by_label
from univi.utils.seed import set_seed
from univi.workflows import make_loader

device = "cuda" if torch.cuda.is_available() else ("mps" if torch.backends.mps.is_available() else "cpu")
set_seed(0)
dense = lambda x: x.toarray() if sp.issparse(x) else np.asarray(x)
print(f"UniVI {univi.__version__} on {device}")

UniVI 1.2.1 on cuda

C:\Users\ashfo\micromamba\envs\univi-release-050\Lib\site-packages\tqdm\auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

N_EPOCHS = 400
BATCH_SIZE = 256
N_HVG = 2000
N_LSI = 101
CORRUPT_FRACTION = 0.10   # corrupted cells of each kind, as a fraction of the held-out test cells
FDR_Q = 0.10              # target false discovery rate for flagging
QC_SCORE = "standardized disagreement"   # score used for flagging; fix it before looking at results
N_KNN = 15

1. Data, split and model#

The model is trained on 70% of cells, with 10% for early stopping (validation). Another 10% are calibration cells: clean barcodes that played no part in fitting, which define what normal disagreement looks like. The last 10% are test cells, which are scored and corrupted.

Calibration cells must not double as validation cells. Early stopping selects the checkpoint that fits the validation cells best, so they look slightly more “normal” than any new barcode. In an earlier version of this notebook that did exactly this, the realized false discovery rate was 16% at a 10% target.

data = uds.pbmc_multiome_10k()
rna, atac = data["rna"], data["atac"]
labels_all = rna.obs["cell_type"].astype(str)
idx = split_by_label(labels_all, train_fraction=0.7, val_fraction=0.1, seed=0)
rest = idx.pop("test")
halves = split_by_label(labels_all.iloc[rest], train_fraction=0.5, val_fraction=0.0, seed=1)
idx["cal"], idx["test"] = rest[halves["train"]], rest[halves["test"]]
print({k: len(v) for k, v in idx.items()})

rna_prep = RNAPreprocessor(n_hvg=N_HVG, scale=True).fit(rna[idx["train"]])
atac_prep = ATACPreprocessor(n_components=N_LSI, drop_first=True, scale=True).fit(atac[idx["train"]])
parts = {k: {"rna": rna_prep.transform(rna[i]), "atac": atac_prep.transform(atac[i])} for k, i in idx.items()}

cfg = UniVIConfig(
    latent_dim=30, beta=1.25, gamma=4.35, encoder_dropout=0.10, decoder_dropout=0.05,
    kl_anneal_start=50, kl_anneal_end=85, align_anneal_start=75, align_anneal_end=110,
    modalities=[ModalityConfig("rna", parts["train"]["rna"].n_vars, [512, 256, 128], [128, 256, 512]),
                ModalityConfig("atac", parts["train"]["atac"].n_vars, [128, 64], [64, 128])])
model = UniVIMultiModalVAE(cfg, loss_mode="v1", v1_recon="avg", normalize_v1_terms=True)
trainer = UniVITrainer(model, make_loader(parts["train"], batch_size=BATCH_SIZE, shuffle=True, drop_last=True),
                       make_loader(parts["val"], batch_size=1024),
                       TrainingConfig(n_epochs=N_EPOCHS, batch_size=BATCH_SIZE, lr=1e-3, weight_decay=1e-4,
                                      device=device, early_stopping=True, patience=50, best_epoch_warmup=110,
                                      log_every=50))
trainer.fit()
print("best epoch:", trainer.best_epoch)

{'train': 6733, 'val': 956, 'cal': 965, 'test': 977}
best epoch: 137

[2026-09-23 12:35:48,173] [UniVITrainer] [INFO] TrainingConfig:
[2026-09-23 12:35:48,174] [UniVITrainer] [INFO]   n_epochs: 400
[2026-09-23 12:35:48,175] [UniVITrainer] [INFO]   batch_size: 256
[2026-09-23 12:35:48,176] [UniVITrainer] [INFO]   lr: 0.001
[2026-09-23 12:35:48,176] [UniVITrainer] [INFO]   weight_decay: 0.0001
[2026-09-23 12:35:48,177] [UniVITrainer] [INFO]   device: 'cuda'
[2026-09-23 12:35:48,178] [UniVITrainer] [INFO]   log_every: 50
[2026-09-23 12:35:48,178] [UniVITrainer] [INFO]   grad_clip: None
[2026-09-23 12:35:48,178] [UniVITrainer] [INFO]   num_workers: 0
[2026-09-23 12:35:48,179] [UniVITrainer] [INFO]   seed: 0
[2026-09-23 12:35:48,179] [UniVITrainer] [INFO]   early_stopping: True
[2026-09-23 12:35:48,180] [UniVITrainer] [INFO]   patience: 50
[2026-09-23 12:35:48,180] [UniVITrainer] [INFO]   min_delta: 0.0
[2026-09-23 12:35:48,180] [UniVITrainer] [INFO]   best_epoch_warmup: 110
Training UniVI:   0%|          | 0/400 [00:00<?, ?it/s][2026-09-23 12:35:49,222] [UniVITrainer] [INFO] [Epoch 001] Train loss=1036.5966 (beta=0.000, gamma=0.000)
[2026-09-23 12:35:49,256] [UniVITrainer] [INFO] [Epoch 001] Val loss=2376.5859 (beta=1.250, gamma=4.350)
Training UniVI:  12%|█▏        | 49/400 [00:38<04:35,  1.27it/s, beta=0.000, gamma=0.000, train_loss=873.2736, val_loss=107760.9453][2026-09-23 12:36:27,898] [UniVITrainer] [INFO] [Epoch 050] Train loss=873.9449 (beta=0.000, gamma=0.000)
[2026-09-23 12:36:27,938] [UniVITrainer] [INFO] [Epoch 050] Val loss=98432.2500 (beta=1.250, gamma=4.350)
Training UniVI:  25%|██▍       | 99/400 [01:18<03:56,  1.27it/s, beta=1.250, gamma=2.983, train_loss=894.8223, val_loss=871.1315][2026-09-23 12:37:07,290] [UniVITrainer] [INFO] [Epoch 100] Train loss=895.6535 (beta=1.250, gamma=3.107)
[2026-09-23 12:37:07,328] [UniVITrainer] [INFO] [Epoch 100] Val loss=870.8674 (beta=1.250, gamma=4.350)
Training UniVI:  27%|██▋       | 108/400 [01:26<03:48,  1.28it/s, beta=1.250, gamma=4.226, train_loss=895.4891, val_loss=870.4218][2026-09-23 12:37:14,521] [UniVITrainer] [INFO] [Epoch 109] Best tracking warmup ends next epoch (best_epoch_warmup=110).
Training UniVI:  27%|██▋       | 109/400 [01:27<03:48,  1.27it/s, beta=1.250, gamma=4.350, train_loss=895.4049, val_loss=870.1731][2026-09-23 12:37:15,548] [UniVITrainer] [INFO] [Epoch 110] New best val loss: 870.1731
Training UniVI:  28%|██▊       | 110/400 [01:28<04:08,  1.17it/s, beta=1.250, gamma=4.350, train_loss=893.6312, val_loss=869.9840][2026-09-23 12:37:16,333] [UniVITrainer] [INFO] [Epoch 111] New best val loss: 869.9840
Training UniVI:  28%|██▊       | 112/400 [01:29<03:57,  1.21it/s, beta=1.250, gamma=4.350, train_loss=893.8864, val_loss=869.1840][2026-09-23 12:37:17,923] [UniVITrainer] [INFO] [Epoch 113] New best val loss: 869.1840
Training UniVI:  29%|██▉       | 116/400 [01:32<03:47,  1.25it/s, beta=1.250, gamma=4.350, train_loss=894.3108, val_loss=868.5476][2026-09-23 12:37:21,127] [UniVITrainer] [INFO] [Epoch 117] New best val loss: 868.5476
Training UniVI:  34%|███▍      | 136/400 [01:48<03:27,  1.27it/s, beta=1.250, gamma=4.350, train_loss=887.3060, val_loss=868.5013][2026-09-23 12:37:37,073] [UniVITrainer] [INFO] [Epoch 137] New best val loss: 868.5013
Training UniVI:  37%|███▋      | 149/400 [01:58<03:18,  1.27it/s, beta=1.250, gamma=4.350, train_loss=886.3153, val_loss=869.6022][2026-09-23 12:37:47,668] [UniVITrainer] [INFO] [Epoch 150] Train loss=885.0007 (beta=1.250, gamma=4.350)
[2026-09-23 12:37:47,701] [UniVITrainer] [INFO] [Epoch 150] Val loss=869.9655 (beta=1.250, gamma=4.350)
Training UniVI:  46%|████▋     | 186/400 [02:29<02:49,  1.26it/s, beta=1.250, gamma=4.350, train_loss=877.3926, val_loss=870.0141][2026-09-23 12:38:17,287] [UniVITrainer] [INFO] Early stopping at epoch 187 (best val loss=868.5013, best epoch=137)
Training UniVI:  46%|████▋     | 186/400 [02:29<02:51,  1.25it/s, beta=1.250, gamma=4.350, train_loss=877.3926, val_loss=870.0141]
[2026-09-23 12:38:17,295] [UniVITrainer] [INFO] Restored best model from epoch 137 (val loss=868.5013)

2. Corrupt held-out cells in known ways#

All corruptions are made from raw counts of test cells and then passed through the fitted transforms, exactly like real barcodes:

  • swap, other type: RNA of one cell with the ATAC of a cell of a different type (a mis-paired barcode between populations);

  • swap, same type: RNA of one cell with the ATAC of another cell of the same type (the hard case);

  • doublet: RNA counts of two cells of different types added together, and likewise their ATAC counts (a heterotypic doublet, whose modalities still agree with each other).

rng = np.random.default_rng(0)
test_rna_raw, test_atac_raw = rna[idx["test"]], atac[idx["test"]]
test_labels = labels_all.iloc[idx["test"]].to_numpy()
n_test = len(test_labels)
K = max(1, int(CORRUPT_FRACTION * n_test))


def pick_partner(i, same):
    pool = np.flatnonzero((test_labels == test_labels[i]) if same else (test_labels != test_labels[i]))
    pool = pool[pool != i]
    return rng.choice(pool) if len(pool) else None


def counts_rows(adata, rows):
    return sp.csr_matrix(adata.layers["counts"][rows], dtype=np.float32)


def as_counts(x, obs_names, var):
    out = ad.AnnData(sp.csr_matrix(x), obs=pd.DataFrame(index=obs_names), var=var.copy())
    out.layers["counts"] = out.X.copy()
    return out


pieces_rna, pieces_atac, kinds, cell_labels = [], [], [], []
# clean test cells
pieces_rna.append(counts_rows(test_rna_raw, np.arange(n_test)))
pieces_atac.append(counts_rows(test_atac_raw, np.arange(n_test)))
kinds += ["clean"] * n_test
cell_labels += list(test_labels)
for kind in ("swap, other type", "swap, same type", "doublet"):
    same = kind == "swap, same type"
    recipients = rng.choice(n_test, size=K, replace=False)
    pairs = [(i, pick_partner(i, same)) for i in recipients]
    pairs = [(i, j) for i, j in pairs if j is not None]
    a, b = np.array([p[0] for p in pairs]), np.array([p[1] for p in pairs])
    if kind == "doublet":
        pieces_rna.append(counts_rows(test_rna_raw, a) + counts_rows(test_rna_raw, b))
        pieces_atac.append(counts_rows(test_atac_raw, a) + counts_rows(test_atac_raw, b))
    else:
        pieces_rna.append(counts_rows(test_rna_raw, a))
        pieces_atac.append(counts_rows(test_atac_raw, b))
    kinds += [kind] * len(pairs)
    cell_labels += list(test_labels[a])

names = [f"q{i}" for i in range(len(kinds))]
query = {"rna": rna_prep.transform(as_counts(sp.vstack(pieces_rna), names, rna.var)),
         "atac": atac_prep.transform(as_counts(sp.vstack(pieces_atac), names, atac.var))}
kinds = np.array(kinds)
print(pd.Series(kinds).value_counts().to_dict())

{'clean': 977, 'swap, other type': 97, 'swap, same type': 97, 'doublet': 97}

3. Four scores#

  • latent disagreement: distance between the posterior means of the RNA and the ATAC encoder;

  • standardized disagreement: the same difference scaled by the two posteriors’ variances, Σ_d (μ_rna − μ_atac)² / (σ²_rna + σ²_atac), so that dimensions the encoders are unsure about count less;

  • cross-prediction error: mean squared error of RNA predicted from ATAC against the observed RNA, plus the same for ATAC predicted from RNA;

  • distance to reference: mean distance of the barcode’s joint (fused) embedding to its nearest training cells. This is not a disagreement score; it asks whether the barcode looks like any cell the model has seen, which is how doublets are usually found.

def posterior(adata, mod, batch=2048):
    mus, variances = [], []
    model.eval()
    with torch.no_grad():
        for s in range(0, adata.n_obs, batch):
            x = torch.as_tensor(dense(adata.X[s:s + batch]), dtype=torch.float32, device=device)
            mu, logvar = model.encode_modalities({mod: x})
            mus.append(mu[mod].cpu().numpy())
            variances.append(logvar[mod].exp().cpu().numpy())
    return np.concatenate(mus), np.concatenate(variances)


z_train = encode_fused_adata_pair(model, parts["train"], device=device, write_to_adatas=False)["mu"]
reference = NearestNeighbors(n_neighbors=N_KNN).fit(z_train)


def score_barcodes(pair):
    mu_r, var_r = posterior(pair["rna"], "rna")
    mu_a, var_a = posterior(pair["atac"], "atac")
    diff = mu_r - mu_a
    rna_hat = cross_modal_predict(model, pair["atac"], src_mod="atac", tgt_mod="rna", device=device)
    atac_hat = cross_modal_predict(model, pair["rna"], src_mod="rna", tgt_mod="atac", device=device)
    fused = encode_fused_adata_pair(model, pair, device=device, write_to_adatas=False)["mu"]
    return pd.DataFrame({
        "latent disagreement": np.linalg.norm(diff, axis=1),
        "standardized disagreement": (diff ** 2 / (var_r + var_a)).sum(1),
        "cross-prediction error": ((dense(pair["rna"].X) - rna_hat) ** 2).mean(1)
                                  + ((dense(pair["atac"].X) - atac_hat) ** 2).mean(1),
        "distance to reference": reference.kneighbors(fused)[0].mean(1),
    }, index=pair["rna"].obs_names)


cal_scores = score_barcodes(parts["cal"])
q_scores = score_barcodes(query)
assert QC_SCORE in q_scores.columns, f"QC_SCORE must be one of {list(q_scores.columns)}"
q_scores.groupby(kinds).median().round(3)
latent disagreement standardized disagreement cross-prediction error distance to reference
clean 1.058 1.431000 0.814 0.564
doublet 1.438 3.315000 1.085 0.955
swap, other type 4.208 32.077999 1.441 1.502
swap, same type 1.754 4.706000 0.986 0.582

4. Which score detects which problem?#

AUROC of each score for each kind of corrupted barcode against the clean test barcodes (0.5 = no better than chance).

auroc = pd.DataFrame({kind: {name: roc_auc_score(np.r_[np.zeros((kinds == "clean").sum()), np.ones((kinds == kind).sum())],
                                                 np.r_[q_scores[name][kinds == "clean"], q_scores[name][kinds == kind]])
                             for name in q_scores.columns}
                      for kind in ["swap, other type", "swap, same type", "doublet"]})
display(auroc.round(3))

fig, axes = plt.subplots(1, len(q_scores.columns), figsize=(3.4 * len(q_scores.columns), 2.8))
for ax, name in zip(axes, q_scores.columns):
    vals = np.log10(q_scores[name] + 1e-9)
    bins = np.histogram_bin_edges(vals, 40)
    for kind in ["clean", "swap, other type", "swap, same type", "doublet"]:
        ax.hist(vals[kinds == kind], bins=bins, density=True, histtype="step", lw=1.3, label=kind)
    ax.set(xlabel=f"log10 {name}", title=name)
axes[0].legend(frameon=False, fontsize=7)
plt.tight_layout()
plt.show()
swap, other type swap, same type doublet
latent disagreement 0.992 0.760 0.731
standardized disagreement 0.991 0.752 0.759
cross-prediction error 0.804 0.603 0.712
distance to reference 0.929 0.525 0.820
../../_images/5f2de4ba1d2dbad17cb4d1642d15d38ba718b5b8ce12fbe2e8746109366c82aa.png

Expect disagreement scores to find swaps between cell types and to struggle with swaps within a type (the two cells are genuinely similar), and expect a doublet, whose two modalities come from the same pair of cells, to be found by its distance to the reference rather than by disagreement. The table shows whether that holds here.

5. Flagging with a controlled false discovery rate#

A raw score has no natural cutoff. Conformal p-values give it one: for a barcode with score s, p = (1 + #{calibration barcodes with score ≥ s}) / (n_calibration + 1). If the barcode is clean and exchangeable with the calibration barcodes, this p-value is valid, and Benjamini–Hochberg applied to the p-values of many barcodes controls the false discovery rate (Bates et al., 2023, Annals of Statistics). Here the flagging score is QC_SCORE, fixed in the parameters cell; choosing the score after looking at the AUROC table would weaken the guarantee.

def conformal_p(cal, test):
    cal = np.sort(np.asarray(cal))
    n_ge = len(cal) - np.searchsorted(cal, np.asarray(test), side="left")
    return (1 + n_ge) / (len(cal) + 1)


def benjamini_hochberg(p, q):
    p = np.asarray(p)
    order = np.argsort(p)
    passed = p[order] <= q * np.arange(1, len(p) + 1) / len(p)
    k = passed.nonzero()[0].max() + 1 if passed.any() else 0
    flag = np.zeros(len(p), dtype=bool)
    flag[order[:k]] = True
    return flag


rows = {}
for name in q_scores.columns:
    flag = benjamini_hochberg(conformal_p(cal_scores[name], q_scores[name]), FDR_Q)
    row = {"flagged": int(flag.sum()),
           "realized FDR": (flag & (kinds == "clean")).sum() / max(flag.sum(), 1)}
    for kind in ["swap, other type", "swap, same type", "doublet"]:
        row[f"power: {kind}"] = flag[kinds == kind].mean()
    rows[name + (" (chosen)" if name == QC_SCORE else "")] = row
print(f"target FDR {FDR_Q:.0%}; {len(cal_scores)} calibration barcodes; "
      f"{(kinds != 'clean').sum()} corrupted among {len(kinds)} scored")
p_check = conformal_p(cal_scores[QC_SCORE], q_scores[QC_SCORE][kinds == "clean"])
print(f"validity check: {np.mean(p_check <= 0.1):.1%} of clean test barcodes have p <= 0.10 "
      f"(about 10% if they are exchangeable with the calibration barcodes)")
pd.DataFrame(rows).T.astype({"flagged": int}).round(3)

target FDR 10%; 965 calibration barcodes; 291 corrupted among 1268 scored
validity check: 9.2% of clean test barcodes have p <= 0.10 (about 10% if they are exchangeable with the calibration barcodes)
flagged realized FDR power: swap, other type power: swap, same type power: doublet
latent disagreement 0 0.000 0.000 0.00 0.000
standardized disagreement (chosen) 55 0.018 0.526 0.01 0.021
cross-prediction error 0 0.000 0.000 0.00 0.000
distance to reference 0 0.000 0.000 0.00 0.000

The realized FDR counts every flagged clean test barcode as a false discovery. Some of those may be real problems that were already in the data (the next section looks at them), so the realized FDR can sit slightly above the target without the procedure failing. The guarantee is for clean barcodes that are exchangeable with the calibration barcodes, so it does not extend to training cells (the model has fit them) or to a new sample processed differently.

Power is limited by the number of calibration barcodes. The smallest possible conformal p-value is 1 / (n_calibration + 1), and Benjamini–Hochberg needs many small p-values to flag anything, so a score that separates a problem well (high AUROC) can still flag nothing if few corrupted barcodes score above every calibration barcode. More calibration barcodes raise power.

6. What does it flag among real, uncorrupted cells?#

The same procedure applied to the clean test barcodes alone, with the chosen score. Flagged barcodes are compared with the others on sequencing depth and cell type.

clean = kinds == "clean"
p_clean = conformal_p(cal_scores[QC_SCORE], q_scores[QC_SCORE][clean])
flag_clean = benjamini_hochberg(p_clean, FDR_Q)
rna_counts = np.asarray(test_rna_raw.layers["counts"].sum(axis=1)).ravel()
atac_counts = np.asarray(test_atac_raw.layers["counts"].sum(axis=1)).ravel()
summary = pd.DataFrame({"flagged": flag_clean, "cell_type": test_labels,
                        "RNA counts": rna_counts, "ATAC counts": atac_counts})
print(f"{flag_clean.sum()} of {clean.sum()} clean test barcodes flagged at FDR {FDR_Q:.0%}")
display(summary.groupby("flagged")[["RNA counts", "ATAC counts"]].median().round(0))
by_type = summary.groupby("cell_type")["flagged"].agg(["sum", "size", "mean"])
by_type.columns = ["flagged", "cells", "fraction flagged"]
by_type.sort_values("fraction flagged", ascending=False).round(3)

0 of 977 clean test barcodes flagged at FDR 10%
RNA counts ATAC counts
flagged
False 3751.0 20465.0
flagged cells fraction flagged
cell_type
CD14 Mono 0 256 0.0
CD16 Mono 0 45 0.0
CD4 Naive 0 139 0.0
CD4 TCM 0 112 0.0
CD4 TEM 0 29 0.0
CD8 Naive 0 136 0.0
CD8 TEM_1 0 33 0.0
CD8 TEM_2 0 32 0.0
HSPC 0 3 0.0
Intermediate B 0 30 0.0
MAIT 0 13 0.0
Memory B 0 31 0.0
NK 0 41 0.0
Naive B 0 13 0.0
Plasma 0 2 0.0
Treg 0 17 0.0
cDC 0 19 0.0
gdT 0 15 0.0
pDC 0 11 0.0 
view = ad.AnnData(encode_fused_adata_pair(model, parts["test"], device=device, write_to_adatas=False)["mu"]
                  .astype(np.float32), obs=pd.DataFrame({"cell_type": pd.Categorical(test_labels),
                                                         "flagged": pd.Categorical(np.where(flag_clean, "flagged", "not flagged")),
                                                         f"log10 {QC_SCORE}": np.log10(q_scores[QC_SCORE][clean].to_numpy() + 1e-9)},
                                                        index=parts["test"]["rna"].obs_names))
sc.pp.neighbors(view, n_neighbors=30, use_rep="X")
sc.tl.umap(view, random_state=0)
sc.pl.umap(view, color=["cell_type", f"log10 {QC_SCORE}", "flagged"], wspace=0.45, legend_fontsize=6)
../../_images/2060496185e3f026141f2db0ed21a7b7a5227d0403227aae313e2e9376b26e57.png

Flagged barcodes concentrated in one cell type can mean either a population the model represents poorly (check the model’s alignment for that type) or a real artifact; unusually high counts point to multiplets, unusually low counts to damaged cells or empty droplets that passed filtering. None of this is confirmation, but it tells you which barcodes to inspect first.

Using this on your own data#

  • Score barcodes that the model did not train on, with calibration barcodes from the same sample and processing. To score every barcode of a dataset, use cross-fitting: train two models on halves of the data and score each half with the model trained on the other.

  • Pick the score before looking at the results, according to the problem you worry about: disagreement scores for mis-pairing, distance to the reference for doublets.

  • Report flagged barcodes as candidates with their conformal p-values rather than removing them silently.

Contents