Single-Cell Quality Control & Preprocessing

Canonical QC Order

The dependency graph

The pipeline is not strictly serial. Background-info-based ambient correction (SoupX-style, which needs the empty-droplet soup profile) and the cell-QC chain (cell calling → demultiplex → Gate 1 → Gate 2) are two parallel inputs that only converge at the ambient-correction step — treating the whole thing as one straight line is what makes "just run steps in card order" feel plausible when it isn't.

                  raw droplets
                       │
              ┌────────┴────────┐
              │                 │
       background info       cell calling
       (empty-droplet /         │
        soup profile)      demultiplex
              │                 │
              │              Gate 1
              │                 │
              │              Gate 2
              │                 │
              └──────┐   ┌──────┘
                     ↓   ↓
               ambient correction
                     │
               final annotation
                     │
               lineage extraction
                     │
        final analysis / benchmark freeze
                     │
          ┌──────────┼──────────┐
          ↓          ↓          ↓
         PCA         FM         DE

(There is an earlier, narrower freeze too — right after Gate 2, before ambient correction — that only decides which cells are eligible for artifact modeling. See "Two freezes, not one" below; the graph above collapses both into one node for brevity.)

(CellBender is deliberately not on this graph — see "The CellBender branch" below; it replaces conventional cell calling plus the "background info"-based ambient-correction branch with one joint step. Demultiplexing, Gate 1, and Gate 2 are unaffected — they still run downstream on CellBender's output.)

Same graph, as a linear walk-through

Useful for reading top-to-bottom in plain language, but remember: this reads as one line only because it picks a single path through the branching graph above — it is not evidence that every step is a hard serial dependency on the one before it.

RAW DATA
│
├─ FASTQ / raw droplet counts
│
▼
PROVENANCE
capture · channel · donor · sample
│
▼
CELL CALLING
Is the barcode a cell?
│
▼
DEMULTIPLEXING
SNG · DBL · AMB
│
▼
GATE 1
mild low-quality-cell QC
counts ↓ · genes ↓ · mito ↑
│
▼
GATE 2
per-capture doublet scoring
raw counts · full lineage context
│
▼
ARTIFACT-MODELING POPULATION FREEZE
confirmed doublet → remove
unresolved → retain + flag
(eligibility for ambient correction only —
 NOT the final analysis population)
│
▼
AMBIENT CORRECTION
per capture · full lineage context
raw counts preserved
│
▼
FINAL ANNOTATION
harmonize lineage labels
│
▼
LINEAGE EXTRACTION
Platelet / Monocyte / target population
│
▼
FINAL ANALYSIS / BENCHMARK POPULATION FREEZE
same cells across models
│
├──────── PCA preprocessing
├──────── FM A preprocessing
├──────── FM B preprocessing
└──────── statistical-analysis preprocessing

Gate 1 stays deliberately mild for two reasons that pull in opposite directions: too little Gate 1 lets extreme low-information cells distort the PCA/kNN geometry the doublet detector relies on; too much Gate 1 (especially on the upper tail) pre-deletes the very doublets Gate 2 is meant to catch. High counts/high genes are Gate-2 diagnostic evidence, not a Gate-1 deletion rule.

The CellBender branch

CellBender does not slot into the same pipeline position as DecontX/SoupX — it is a clearly separate alternate branch off "raw droplets" in the dependency graph above, not a stage merged into the main chain. It replaces conventional cell calling and ambient correction only — not demultiplexing or Gate 1, which still run downstream of it:

versus the conventional path: Cell Ranger filtered matrix → demultiplexing → Gate 1 → Gate 2 → DecontX/SoupX. Treat these as two alternative branches, not interchangeable steps in one universal order — running CellBender and then DecontX on its output double-corrects for ambient RNA.

Why ambient correction follows doublet handling, not the reverse

A doublet and an ambient-contaminated singlet are different generative models, and conflating them motivates the wrong step order:

• Doublet: Xi ≈ XA + XB — both components are real cells captured together.
• Ambient contamination: Xi ≈ (1−ε)XA + εXsoup — the second component is extracellular background, not a second cell.

If DecontX (or SoupX) runs before doublet detection, a real doublet's second-cell signal can be misattributed to background contamination, and the corrected counts handed to the doublet detector no longer reflect the original capture geometry. For the DecontX-oriented workflow used here, doublet decisions are frozen before decontamination; SoupX and CellBender have method-specific dependencies described separately — SoupX's soup estimation draws on the raw, unfiltered droplet pool rather than obeying this same doublets-then-ambient ordering.

What matrix, cells, and grouping go where

SoupX uses raw droplets to estimate the soup profile, then applies that estimate to retained cells — the table lists both inputs on one row because both feed the same correction step, not two separate pipeline stages.

Two rules follow directly from this table: do not subset to the target lineage before Scrublet or DecontX unless the method specifically requires it (a monocyte is not a QC failure just because the final benchmark is platelet-only — it is valid background context for artifact modeling), and the unit of QC (capture/channel) and the unit of inference (donor/patient) can differ — cell QC determines which measurements are trustworthy, it does not turn cells into independent biological replicates.

Data & Provenance

Stage 0 — identity mapping. Before any QC decision, every barcode needs a resolved chain of custody. At minimum, store:

barcode → capture_id → channel_id → sample_id → donor_id → multiplexing_status → source_matrix → count_layer

This is especially important for public datasets, where several of these stages may already be baked into the deposited matrix and cannot be reconstructed after the fact.

Public-cohort provenance audit

For every external cohort, audit what was already done upstream before deciding what you need to (re)do. The governing principle: start from the earliest recoverable stage; never pretend an unavailable upstream stage was reconstructed.

If cell calling was performed upstream and only a filtered matrix was deposited, say so explicitly rather than treating the filtered matrix as though it were raw droplets. Re-running doublet detection on a public deposit that has already had doublets removed is a distinct methodological decision — not a routine rerun — because the residual doublet population and its expected rate are no longer the same as an unfiltered capture.

Cell Calling

Cell calling answers one question: barcode → cell or not-cell. It precedes cell-level QC — a barcode with almost no transcripts is not a low-quality cell, it is an empty droplet that never held a cell (see Guest 10 · Empty Droplets & Cell Calling in the parable below). Conventional cell calling (CellRanger's default, EmptyDrops-style methods) operates on the raw droplet matrix and outputs the filtered/called matrix that Gate 1 then consumes. CellBender replaces this step with a joint cell-calling + ambient-correction model — see the CellBender branch in the Canonical QC Workflow above; the two are alternative branches, not sequential steps.

Demultiplexing

For genetically or hashtag-multiplexed captures, demultiplexing resolves each called cell to one of three states before any transcriptome-based doublet detection runs:

called cells → donor assignment → SNG (singlet) · DBL (cross-donor doublet) · AMB (ambiguous)

This step matters for Scrublet/scDblFinder calibration, not just sample bookkeeping: if genotype/HTO demultiplexing has already removed cross-donor doublets, the residual doublet prevalence a transcriptome-based detector should expect is no longer the raw loading-chart rate —

P(doublet | survived demultiplexing) ≠ P(doublet)

Scrublet should not be taught, or run, as though it always sees the original, unfiltered doublet population.

And Then There Were Some · Ten Guests

Agatha Christie's And Then There Were None (1939): ten guests on Soldier Island, judged one by one against a nursery rhyme. Here each guest stands for one quality-control theme — a failure mode, not a step in a pipeline. The ten guests are failure modes, not workflow order: see the Canonical QC Order section above for the actual execution sequence. (Classic-novel spoilers inside.)

"And then there were none" is the warning, not the goal: ten guests, zero executions — the matrix that enters the model is the witnesses who survived cross-examination.

Gate 1 — Cell Quality & QC Pipelines

Gate 1 is the mild, capture-aware cell-quality pass described in the Canonical QC Workflow: lower-tail counts/genes, upper-tail mito — diagnostics first, a named gate1_pass decision second. The frameworks below implement Gate-1-style filtering (and more) end-to-end.

Gate 2 — Doublet Detection

Ambient RNA & Contamination Removal

DecontX applies directly to the retained cell matrix after doublet handling, in this workflow's doublet-then-ambient order — see why ambient correction follows doublet handling. SoupX instead estimates its ambient/"soup" profile from the raw, unfiltered droplet pool, then uses that estimate to correct the doublet-cleared, full-lineage cell population — its soup estimation step does not itself obey the doublets-then-ambient cell-population ordering, even though the correction it applies does. CellBender is a separate branch that jointly performs cell calling and ambient correction on the raw droplet matrix; it is not interchangeable with running DecontX (or SoupX) on an already cell-called matrix.

Feature Selection & Dimensionality Reduction

Layer III — downstream preprocessing, not cell QC. These methods select or transform features for a QC-retained population; they are not a substitute for the identity/QC steps in the Canonical QC Order. See also Metacell & Statistical Methods below, which sits in the same downstream layer.

Foundation Models for Single-Cell Analysis

Layer III — downstream preprocessing, not a QC stage. Foundation models consume the same QC-retained biological cells that any other analysis does; a model's own input contract (tokenization, normalization) runs after the canonical identity/QC backbone, never before it decides which cells pass QC.

Metacell & Statistical Methods

Data Repositories & Standardization

Part of the Method Library — for a full index of every method on this page including the ones above (doublets, ambient RNA, features, foundation models), see the Quick Reference: Choosing the Right Tool section further down the page.

Perturb-seq QC — CRISPR Screen-Specific

STANDARD scRNA QC
       ↓
QC-retained cells
       ↓
PERTURBATION-SPECIFIC QC
guide identity · coverage · efficacy · escape status · control quality
       ↓
perturbation inference

Why standard scRNA-seq QC is necessary but not sufficient

Perturbation-specific QC extends the standard pipeline — it does not replace or reorder it. The same QC principles apply, but thresholds must remain experiment-aware and must not erase genuine perturbation-induced phenotypes: a perturbation can legitimately alter RNA content, mitochondrial programs, and cell state, so a threshold tuned on unperturbed data can wrongly delete a real perturbation effect. Every barcode still needs to pass Gate 1/Gate 2/ambient checks first. What changes is that each surviving cell additionally carries a perturbation identity: which guide RNA(s) it received. Poor perturbation QC inflates both false negatives (real effects missed because escape cells dilute the signal) and false positives (guide-unrelated technical variation misread as a perturbation effect). The five challenges below are unique to CRISPR-based designs and require targeted QC steps that standard pipelines do not implement.

Perturb-seq QC Metrics at a Glance

Minimal code walkthrough — Perturb-seq guide QC (Python / Scanpy)

import scanpy as sc
import pandas as pd
import numpy as np

# ── Assumes adata has already passed standard cell QC + doublet removal ──────────
# adata.obs will receive guide assignment columns.
# Cell Ranger multi output directory structure:
#   /outs/per_sample_outs/<sample>/count/crispr_analysis/protospacer_calls_per_cell.csv

# ── 1. Load guide assignments from Cell Ranger protospacer calls ─────────────────
guide_calls = pd.read_csv(
    "crispr_analysis/protospacer_calls_per_cell.csv",
    index_col="cell_barcode"
)
# Columns: cell_barcode | num_features | feature_call | num_umis
adata.obs = adata.obs.join(guide_calls[["feature_call", "num_umis"]], how="left")
adata.obs["guide_id"] = adata.obs["feature_call"].fillna("unassigned")

# ── 2. Guide assignment QC ────────────────────────────────────────────────────────
n_cells = adata.n_obs
assignment_rate = (adata.obs["guide_id"] != "unassigned").mean()
multi_guide_rate = adata.obs["feature_call"].str.contains("|", na=False, regex=False).mean()
# Thresholds below are illustrative starting points, not universal cutoffs —
# use assay/caller-specific values for your library and chemistry.
print(f"Guide assignment rate: {assignment_rate:.1%}  (illustrative reference: > 85 %)")
print(f"Multi-guide rate:      {multi_guide_rate:.1%}  (illustrative reference: < 5 % in low-MOI)")

# ── 3. Perturbation coverage ──────────────────────────────────────────────────────
coverage = adata.obs["guide_id"].value_counts()
low_coverage = coverage[coverage < 10]
print(f"
{len(low_coverage)} guides with < 10 cells (will be flagged):")
print(low_coverage.to_string())
adata.obs["low_coverage_guide"] = adata.obs["guide_id"].isin(low_coverage.index)

# ── 4. Target gene knockdown score ────────────────────────────────────────────────
# guide_target_map: dict {guide_id: "TARGET_GENE"} — from your library design file
# Assumes adata.X is log1p-normalised

nt_mask = adata.obs["guide_id"].str.startswith("NT")
kd_scores = {}
for guide, target_gene in guide_target_map.items():
    if target_gene not in adata.var_names:
        continue
    guide_mask = adata.obs["guide_id"] == guide
    if guide_mask.sum() < 5:
        continue
    gene_idx = adata.var_names.get_loc(target_gene)
    nt_vals   = np.asarray(adata.X[nt_mask,    gene_idx]).flatten()
    guide_vals = np.asarray(adata.X[guide_mask, gene_idx]).flatten()
    nt_med, nt_sd = np.median(nt_vals), np.std(nt_vals)
    kd_scores[guide] = float((nt_med - np.median(guide_vals)) / (nt_sd + 1e-6))

# kd_scores are standardized effect sizes in SD units (NOT percent reduction)
kd_df = pd.Series(kd_scores, name="KD_score_sd").sort_values()
low_kd = kd_df[kd_df < 0.5].index.tolist()
print(f"
Guides with KD score < 0.5 SD (consider Mixscape): {len(low_kd)}")

adata.obs["kd_score"]   = adata.obs["guide_id"].map(kd_scores)

# ── 5. NT control UMAP sanity check ─────────────────────────────────────────────
# If UMAP is already computed upstream:
nt_adata = adata[nt_mask].copy()
# Inspect for a batch association; mixing here is not proof of correction,
# and separation is not proof it failed — pair with a quantitative diagnostic.
sc.pl.umap(nt_adata, color="batch",
           title="NT controls by batch — inspect for batch association")

# ── 6. Save QC-annotated object ──────────────────────────────────────────────────
adata.write_h5ad("adata_perturbseq_qc.h5ad")
print("Saved: adata_perturbseq_qc.h5ad")

⚖️ Side-by-Side: Ambient RNA Removal — Subtraction vs Bayesian Mixture

SoupX and DecontX both correct ambient RNA contamination in droplet-based scRNA-seq, but they diverge on both their input requirements and the fundamental statistical question: should contamination be estimated via cluster-driven marker subtraction using the empty-droplet soup profile (SoupX), or inferred as a per-cell Bayesian mixture directly from the cell count matrix (DecontX)? Neither is interchangeable with CellBender, which instead operates on the raw, unfiltered droplet matrix and jointly performs cell calling — see the CellBender branch above.

Annotation & Lineage Extraction

Two different annotation layers can coexist without contradiction — a preliminary layer used to calibrate QC, and a final layer used to define the analysis population:

Default: capture-aware QC. Gate 1 stratifies by capture alone unless there is strong, independent prior justification for a lineage-specific pass. Optional refinement: broad-lineage-aware thresholds — worth adding only when a lineage is known to have a systematically different RNA-content distribution (platelets, for instance, where a sample-wide lower-tail count threshold can misclassify an entire low-RNA lineage as poor quality). Any lineage-aware refinement should use coarse labels that are independently stable (not derived from the same metric Gate 1 is about to threshold), and its sensitivity should be audited against the capture-only default — otherwise expression structure informs the preliminary annotation, which then shapes which expression structure survives QC, a circularity worth naming rather than assuming away. Regardless of which is used, the population used to estimate doublets/ambient contamination should stay full-lineage; only the final, harmonized annotation is used to extract the target population (see Off-target / Non-target Lineage Cells), and that extraction happens after Gate 2 and ambient correction, not before.

Benchmark / Analysis Assembly

Two freezes, not one

"Freeze the population" collapses two genuinely different decisions into one word. Keeping them apart resolves an inconsistency early drafts of this page had between the prose and the hands-on code:

• Artifact-modeling population freeze — right after Gate 2 (confirmed doublets removed, unresolved cells retained + flagged), before ambient correction. This freeze only decides which cells are eligible to enter DecontX/SoupX — a monocyte belongs here even in a platelet-only study, because it is valid background context for the ambient model. In the hands-on walkthrough below, this is the population returned by module_08.
• Final analysis / benchmark population freeze — after ambient correction, final harmonized annotation, and lineage extraction. This is the population that actually goes into PCA/FM/DE preprocessing, and the one that must stay identical across every model being compared. In the walkthrough below, this is the population returned by module_11 (adata_platelets).

Record a ledger spanning both freezes — gate1_pass, doublet_call, doublet_method_status, decontx_contamination, artifact_modeling_pass, final_lineage, final_analysis_pass — so "why is cell X absent?" can be answered without rerunning the pipeline, and so it's always clear which of the two freezes a given cell dropped out at.

Quick Reference: Choosing the Right Tool

These are independent method-selection questions, not a pipeline — the arrows below name the recommended tool for each problem, not an execution order (see Canonical QC Order for that).

All methods at a glance

🛠️ Hands-On Practice

The steps below walk through a minimal but realistic single-cell QC run using Python/Scanpy — from a cell-called 10x Genomics matrix to a QC-annotated, policy-filtered AnnData object ready for clustering. All code runs in a standard conda environment and is designed to be self-contained.

Environment & packages

Install the core Python QC stack into a fresh environment. Scrublet handles doublet detection; miQC (via rpy2) or adaptive MAD thresholds replace hard cutoffs for cell filtering.

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

pip install scanpy anndata scrublet scvi-tools doubletdetection
# optional ambient-RNA correction (R bridge)
# pip install rpy2
# R: install.packages(c("SoupX", "celda"))   # celda for decontX

Hardware. A standard laptop (8 GB RAM) handles datasets up to ~50k cells; for 100k+ cells use a compute node with 32–64 GB RAM. GPU is only required if running scVI-SOLO doublet detection.

Data structures & formats

• AnnData — central object: adata.X (current working expression matrix), adata.obs (per-cell metadata), adata.var (per-gene metadata)
• Per-cell QC metrics computed by sc.pp.calculate_qc_metrics(): n_genes_by_counts, total_counts, pct_counts_mt
• Raw vs filtered matrix — 10x output contains raw_feature_bc_matrix/ (all barcodes) and filtered_feature_bc_matrix/ (cell-called barcodes). The raw matrix also matters beyond ambient correction: cell calling, EmptyDrops-style approaches, CellBender, and auditing an upstream call all consume it. If only an author-filtered matrix is available, say explicitly that cell calling was performed upstream and cannot be reconstructed from the deposit.
• Representation contract — adata.X is a mutable current working matrix: it gets reassigned as you normalize/transform. adata.layers["raw_counts"] is the immutable raw-count representation, set once and never overwritten. adata.raw is not a reliable stand-in for this — it is just Scanpy's snapshot of whatever was in .X at the time you called it, not a guarantee of untouched sequencer counts.
• Doublet scores — stored in adata.obs["doublet_score"] and adata.obs["predicted_doublet"] after Scrublet; a diagnostic, not yet a filtering decision (see the anti-patterns box below)
• h5ad — HDF5-backed AnnData format; save with adata.write_h5ad() for downstream sharing

Minimal code walkthrough — as a module graph, not a recipe

This walkthrough follows the Canonical QC Workflow, split into named modules. Each module's comment states its INPUT, OUTPUT, and — the part a linear recipe hides — whether the cell population can change. Only 02/03/08/11 are allowed to change which cells are present; 05 (Gate 1) only ever adds a gate1_pass column — filtering is deferred to the caller so the full ledger survives (see master_adata below). Filtering always happens on a decision column, never inline on a raw metric.

import scanpy as sc
import numpy as np
import pandas as pd
from scipy.stats import median_abs_deviation

def _mad_bounds(metric, n_mads=5):
    med, mad = np.median(metric), median_abs_deviation(metric)
    return med - n_mads * mad, med + n_mads * mad

USE_LINEAGE_AWARE_GATE1 = False        # opt-in only — see "Default: capture-aware QC" above; flip to True
                                       # only with strong, independent prior justification for a lineage-
                                       # specific RNA-content distribution (e.g. platelets), and audit
                                       # sensitivity against the capture-only default before trusting it.
MIN_STRATUM_SIZE = 30                 # capture x lineage strata smaller than this fall back to capture-only MAD
                                       # (only relevant when USE_LINEAGE_AWARE_GATE1 is True)
DIAGNOSTIC_MIN_GENES_REFERENCE = 200  # illustrative reference only — never a deletion rule
EXPECTED_DOUBLET_RATE_BY_CAPTURE = {
    # Illustrative starting points, one per capture's own loading/recovery target —
    # NOT a single global constant. If demultiplexing (03) already removed
    # cross-donor doublets for a capture, lower its value: residual prevalence
    # after demux is not the same as the raw loading-chart rate.
    "capture_1": 0.06,
}
DEFAULT_EXPECTED_DOUBLET_RATE = 0.06  # fallback only for captures missing from the dict above


def module_00_validate_input(path):
    # 00_validate_input — INPUT: raw file(s) / OUTPUT: validated raw matrix / NO CELLS REMOVED
    adata = sc.read_10x_mtx(path, var_names="gene_symbols", cache=True)
    adata.var_names_make_unique()
    # Representation contract: raw_counts is set ONCE, here, and never overwritten.
    # adata.X is a mutable "current working matrix"; adata.raw is NOT a reliable
    # stand-in for this — it only snapshots whatever was in .X when called.
    adata.layers["raw_counts"] = adata.X.copy()
    raw = adata.layers["raw_counts"].data if hasattr(adata.layers["raw_counts"], "data") else adata.layers["raw_counts"]
    assert (raw >= 0).all(), "raw_counts must be non-negative"
    assert np.allclose(raw, np.round(raw)), "raw_counts must be integer-like"
    # Gene filtering is a feature-availability filter — conceptually separate from
    # cell QC; it changes which genes are represented, not which cells survive.
    sc.pp.filter_genes(adata, min_cells=3)
    return adata


def module_01_build_provenance(adata):
    # 01_build_provenance — INPUT: validated matrix / OUTPUT: + capture/donor/sample cols / NO CELLS REMOVED
    # Real run: pull these from sample metadata / CellRanger multi manifests.
    adata.obs["capture_id"] = adata.obs.get("capture_id", "capture_1")
    adata.obs["donor_id"]   = adata.obs.get("donor_id", "unknown")
    adata.obs["sample_id"]  = adata.obs.get("sample_id", "sample_1")
    # Preliminary broad-lineage label for Gate-1 calibration ONLY (e.g. from a
    # coarse marker-based pass upstream) — distinct from the final, post-decontam
    # lineage assigned in module_10. Placeholder here if none is supplied.
    adata.obs["broad_lineage"] = adata.obs.get("broad_lineage", "unassigned")
    return adata


def module_02_cell_call_or_record_upstream(adata, called_upstream=True):
    # 02_cell_call_or_record_upstream — INPUT: provenance-tagged matrix /
    # OUTPUT: called_cell flag (or upstream-provenance note) / CELL POPULATION MAY CHANGE
    if called_upstream:
        # filtered_feature_bc_matrix was already cell-called by CellRanger — say so,
        # don't pretend this pipeline reconstructed cell calling from raw droplets.
        adata.obs["called_cell"] = True
        adata.uns["cell_calling_provenance"] = "upstream (CellRanger filtered matrix)"
    else:
        adata.obs["called_cell"] = True  # a real EmptyDrops-style caller would run here
        adata.uns["cell_calling_provenance"] = "this pipeline"
    return adata[adata.obs["called_cell"]].copy()


def module_03_demultiplex_or_record_upstream(adata, demuxed_upstream=False):
    # 03_demultiplex_or_record_upstream — INPUT: called cells /
    # OUTPUT: SNG/DBL/AMB status (or upstream note) / CELL POPULATION MAY CHANGE (AMB/DBL removed if applicable)
    if not demuxed_upstream:
        adata.obs["demux_status"] = "SNG"  # single-sample capture — nothing to resolve
        adata.uns["demultiplex_provenance"] = "not applicable"
        return adata
    # Genotype/HTO demultiplexing would populate SNG/DBL/AMB here. Once cross-donor
    # doublets are removed, P(doublet | survived demux) != P(doublet) — module_06's
    # per-capture expected rate should reflect that, not the raw loading-chart number.
    return adata[adata.obs["demux_status"] == "SNG"].copy()


def module_04_compute_gate1_metrics(adata):
    # 04_compute_gate1_metrics — INPUT: called/demuxed cells + raw counts /
    # OUTPUT: same cells + QC metrics (diagnostics) / NO CELLS REMOVED
    adata.var["mt"] = adata.var_names.str.startswith("MT-")
    sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], percent_top=None, log1p=True, inplace=True)
    sc.pl.violin(adata, ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
                 jitter=0.4, multi_panel=True)
    return adata


def module_05_apply_gate1_policy(adata):
    # 05_apply_gate1_policy — INPUT: metrics + frozen config / OUTPUT: gate1_pass column on ALL called cells /
    # NO CELLS REMOVED (filtering is deferred to the caller — see master_adata pattern below)
    # DEFAULT is capture-only MAD (USE_LINEAGE_AWARE_GATE1 = False). Stratifying further
    # by broad_lineage is an explicit OPT-IN refinement, justified only when a lineage is
    # known to have a systematically different RNA-content distribution (e.g. platelets)
    # AND the preliminary lineage label is independently stable (not itself derived from
    # the metric being thresholded) — otherwise expression -> preliminary annotation ->
    # lineage-specific QC -> retained expression becomes circular. Even with the switch
    # on, a stratum below MIN_STRATUM_SIZE falls back to capture-only bounds.
    low_counts = pd.Series(False, index=adata.obs.index)
    low_genes  = pd.Series(False, index=adata.obs.index)
    high_mito  = pd.Series(False, index=adata.obs.index)

    for cap, cap_group in adata.obs.groupby("capture_id"):
        cap_counts_bounds = _mad_bounds(cap_group["log1p_total_counts"])
        cap_genes_bounds  = _mad_bounds(cap_group["log1p_n_genes_by_counts"])
        cap_mito_bounds   = _mad_bounds(cap_group["pct_counts_mt"])
        # Only split into per-lineage strata when the opt-in switch is on; otherwise
        # every cell in this capture is one stratum, i.e. genuinely capture-only.
        strata = cap_group.groupby("broad_lineage") if USE_LINEAGE_AWARE_GATE1 else [(None, cap_group)]
        for lineage, stratum in strata:
            if USE_LINEAGE_AWARE_GATE1 and len(stratum) >= MIN_STRATUM_SIZE:
                counts_lo, counts_hi = _mad_bounds(stratum["log1p_total_counts"])
                genes_lo,  genes_hi  = _mad_bounds(stratum["log1p_n_genes_by_counts"])
                _,         mito_hi   = _mad_bounds(stratum["pct_counts_mt"])
            else:
                # Lineage-aware refinement is off (capture-only default), or this
                # capture x lineage stratum is too small for a stable per-lineage MAD —
                # either way, fall back to capture-only bounds.
                counts_lo, counts_hi = cap_counts_bounds
                genes_lo,  genes_hi  = cap_genes_bounds
                _,         mito_hi   = cap_mito_bounds
            idx = stratum.index
            low_counts.loc[idx] = stratum["log1p_total_counts"] < counts_lo
            low_genes.loc[idx]  = stratum["log1p_n_genes_by_counts"] < genes_lo
            high_mito.loc[idx]  = stratum["pct_counts_mt"] > mito_hi
            # High counts/genes are Gate-2 diagnostic evidence — NOT filtered here.

    # Optional extreme safety net — NOT part of the canonical policy. Opt in only
    # for a specific, stated reason (e.g. known multiplet-prone capture):
    # MITO_HARD_SAFETY_CAP = 50
    # high_mito |= (adata.obs["pct_counts_mt"] > MITO_HARD_SAFETY_CAP)

    adata.obs["gate1_low_counts"] = low_counts
    adata.obs["gate1_low_genes"]  = low_genes
    adata.obs["gate1_high_mito"]  = high_mito
    # DECISION: an explicit policy built by OR-ing named diagnostics — never an
    # opaque single outlier() call applied straight to a raw metric.
    adata.obs["gate1_pass"] = ~(
        adata.obs["gate1_low_counts"] | adata.obs["gate1_low_genes"] | adata.obs["gate1_high_mito"]
    )
    print("Gate 1 retention by capture (audit habit — check every grouping variable):")
    print(adata.obs.groupby("capture_id")["gate1_pass"].mean())
    return adata  # FULL population retained — gate1_pass is a column, not a subset


def module_06_run_doublet_scoring(adata, expected_rate_by_capture=None):
    # 06_run_doublet_scoring — INPUT: Gate-1-retained cells, raw counts, per capture /
    # OUTPUT: doublet_score, doublet_method_status / NO CELLS REMOVED
    expected_rate_by_capture = expected_rate_by_capture or EXPECTED_DOUBLET_RATE_BY_CAPTURE
    raw = adata.layers["raw_counts"]
    raw_values = raw.data if hasattr(raw, "data") else raw
    assert (raw_values >= 0).all() and np.allclose(raw_values, np.round(raw_values)), \
        "Scrublet requires raw, non-negative, integer-like counts"

    adata.obs["doublet_score"] = np.nan
    # Nullable boolean: an unresolved (too_small / runtime_failed) capture must
    # read as "no call made", not as a genuine negative — plain False would be
    # indistinguishable from an actual singlet call.
    adata.obs["predicted_doublet"] = pd.Series(pd.NA, index=adata.obs_names, dtype="boolean")
    adata.obs["doublet_method_status"] = "not_run"  # reliable | too_small | runtime_failed

    for cap in adata.obs["capture_id"].unique():
        sub = adata[adata.obs["capture_id"] == cap].copy()
        sub.X = sub.layers["raw_counts"]
        if sub.n_obs < 50:
            # Too few cells for a reliable simulated-doublet model. Do NOT
            # silently coerce this into "singlet" — leave it explicit/unresolved.
            adata.obs.loc[sub.obs_names, "doublet_method_status"] = "too_small"
            continue
        try:
            rate = expected_rate_by_capture.get(cap, DEFAULT_EXPECTED_DOUBLET_RATE)
            sc.pp.scrublet(sub, expected_doublet_rate=rate)
            adata.obs.loc[sub.obs_names, "doublet_score"] = sub.obs["doublet_score"].values
            adata.obs.loc[sub.obs_names, "predicted_doublet"] = sub.obs["predicted_doublet"].values.astype("boolean")
            adata.obs.loc[sub.obs_names, "doublet_method_status"] = "reliable"
        except Exception:
            # No silent numeric-threshold fallback if Scrublet's automatic call fails.
            adata.obs.loc[sub.obs_names, "doublet_method_status"] = "runtime_failed"
    return adata


def module_07_review_doublet_status(adata):
    # 07_review_doublet_status — INPUT: doublet scores/status /
    # OUTPUT: confirmed_doublet, unresolved flags / NO CELLS REMOVED (decision only)
    adata.obs["confirmed_doublet"] = (
        (adata.obs["doublet_method_status"] == "reliable") & adata.obs["predicted_doublet"]
    ).astype("boolean")
    adata.obs["doublet_unresolved"] = adata.obs["doublet_method_status"].isin(["too_small", "runtime_failed"])
    return adata


def module_08_freeze_artifact_modeling_population(adata):
    # 08_freeze_artifact_modeling_population — INPUT: master ledger (gate1_pass + confirmed_doublet, ALL called cells) /
    # OUTPUT: artifact_modeling_pass column on the ledger, plus a derived frozen population /
    # CELL POPULATION CHANGES (freeze point) — only in the returned derived view; the input ledger is untouched
    # This is the ARTIFACT-MODELING population freeze (eligibility for DecontX/SoupX),
    # NOT the final analysis/benchmark population — that freeze happens later, in
    # module_11, after ambient correction + harmonized annotation + lineage extraction.
    # Outcome-blind: nothing above this line may be tuned using downstream labels
    # (e.g. fatal/survived) or model performance — that is leakage, not QC.
    adata.obs["artifact_modeling_pass"] = adata.obs["gate1_pass"] & ~adata.obs["confirmed_doublet"].fillna(False)
    print("Artifact-modeling population retention by capture (NOT the final analysis population):")
    print(adata.obs.groupby("capture_id")["artifact_modeling_pass"].mean())
    frozen = adata[adata.obs["artifact_modeling_pass"]].copy()
    frozen.uns["qc_freeze_config"] = {                       # frozen QC configuration
        "gate1_n_mads": 5, "gate1_lineage_aware": USE_LINEAGE_AWARE_GATE1,
        "gate1_min_stratum_size": MIN_STRATUM_SIZE,
        "scrublet_expected_rate_by_capture": EXPECTED_DOUBLET_RATE_BY_CAPTURE,
        "min_capture_size_for_doublet_calling": 50,
    }
    return frozen


def module_09_run_decontamination(adata):
    # 09_run_decontamination — INPUT: frozen population, full lineage, raw counts /
    # OUTPUT: decontx_counts (new layer), decontx_contamination / NO CELLS REMOVED, raw_counts untouched
    # adata.layers["decontx_counts"], contamination = run_decontx(adata.layers["raw_counts"])
    adata.obs["decontx_contamination"] = np.nan  # populate from DecontX/SoupX output
    return adata


def module_10_harmonize_lineage(adata):
    # 10_harmonize_lineage — INPUT: frozen population / OUTPUT: broad_lineage / final_lineage cols /
    # NO CELLS REMOVED
    adata.obs["final_lineage"] = adata.obs.get("final_lineage", adata.obs["broad_lineage"])  # post-decontam label
    return adata


def module_11_build_analysis_population(adata, target_lineage=None):
    # 11_build_analysis_population — INPUT: harmonized population /
    # OUTPUT: target-lineage-extracted analysis object / CELL POPULATION CHANGES (lineage extraction)
    # This IS the FINAL ANALYSIS / BENCHMARK POPULATION FREEZE (see "Two freezes, not
    # one" in Benchmark / Analysis Assembly above) — distinct from module_08's
    # artifact_modeling_pass, which only gated entry into ambient correction.
    if target_lineage is None:
        adata.obs["final_analysis_pass"] = True
        return adata
    adata.obs["final_analysis_pass"] = adata.obs["final_lineage"] == target_lineage
    return adata[adata.obs["final_analysis_pass"]].copy()


# ── run the module graph ──────────────────────────────────────────────────────
master_adata = module_00_validate_input("filtered_feature_bc_matrix/")
master_adata = module_01_build_provenance(master_adata)
master_adata = module_02_cell_call_or_record_upstream(master_adata, called_upstream=True)
master_adata = module_03_demultiplex_or_record_upstream(master_adata, demuxed_upstream=False)
master_adata = module_04_compute_gate1_metrics(master_adata)
master_adata = module_05_apply_gate1_policy(master_adata)
# master_adata now holds ALL called/demuxed cells with a gate1_pass column —
# no rows dropped yet. "Why was cell X absent?" is answerable from this ledger
# all the way to the end of the pipeline.

# Derived VIEW: only Gate-1-retained cells continue into Gate 2 — master_adata
# itself is not replaced.
gate1_adata = master_adata[master_adata.obs["gate1_pass"]].copy()
gate1_adata = module_06_run_doublet_scoring(gate1_adata, expected_rate_by_capture=EXPECTED_DOUBLET_RATE_BY_CAPTURE)
gate1_adata = module_07_review_doublet_status(gate1_adata)

# Write Gate-2 columns back onto the master ledger. Gate-1 failures never ran
# through Gate 2 — they get an explicit "gate1_excluded" status, not a blank.
doublet_cols = ["doublet_score", "predicted_doublet", "doublet_method_status",
                "confirmed_doublet", "doublet_unresolved"]
master_adata.obs = master_adata.obs.join(gate1_adata.obs[doublet_cols], how="left")
master_adata.obs["doublet_method_status"] = master_adata.obs["doublet_method_status"].fillna("gate1_excluded")

adata_frozen = module_08_freeze_artifact_modeling_population(master_adata)  # artifact-modeling freeze
# module_08 sets artifact_modeling_pass on master_adata.obs in place (adata_frozen is a
# .copy() taken AFTER that assignment) -- so it is already on the master ledger below.
adata_frozen = module_09_run_decontamination(adata_frozen)
adata_frozen = module_10_harmonize_lineage(adata_frozen)
adata_platelets = module_11_build_analysis_population(adata_frozen, target_lineage="Platelet")  # final analysis/benchmark freeze
# module_11 likewise sets final_analysis_pass on adata_frozen.obs in place before
# subsetting into adata_platelets, so adata_frozen.obs now also carries
# decontx_contamination (09), final_lineage (10), and final_analysis_pass (11).

# Join those post-artifact-modeling-freeze columns back onto the master ledger so
# qc_ledger_full.csv genuinely spans BOTH freezes, not just Gate 1/Gate 2. Cells
# excluded before the artifact-modeling freeze (gate1 fail / confirmed doublet) never
# entered adata_frozen and correctly land as NaN here -- "did not reach this stage",
# not a false/failed decision.
post_freeze_cols = ["decontx_contamination", "final_lineage", "final_analysis_pass"]
master_adata.obs = master_adata.obs.join(adata_frozen.obs[post_freeze_cols], how="left")

adata_frozen.write_h5ad("adata_qc.h5ad")
master_adata.obs.to_csv("qc_ledger_full.csv")  # complete audit trail, every called cell, spanning both freezes
print(f"Saved adata_qc.h5ad — {adata_frozen.n_obs} cells (master ledger: {master_adata.n_obs} called cells)")

Diagnostics vs. decisions in this walkthrough

Common pitfalls & tips

• One-size-fits-all thresholds hurt. A hard cutoff of >200 genes / <5% mito fails for hepatocytes (high mito) or RBCs (few genes). Use asymmetric, MAD-based adaptive thresholds per capture; add broad-lineage stratification only as a justified, sensitivity-audited refinement (see USE_LINEAGE_AWARE_GATE1 above), not a default.
• Ambient correction follows doublet handling, not the reverse. DecontX applies directly to the doublet-cleared, full-lineage population — see why ambient correction follows doublet handling. SoupX estimates its soup profile from the raw, unfiltered droplet pool, but still applies its correction to that same doublet-cleared population — don't run either tool before doublet decisions are frozen. (CellBender is the exception: it is a joint cell-calling + background-removal step on the raw, unfiltered matrix — a separate branch, not a substitute for this order.)
• Doublet rate depends on your loading design, not a global default. Derive expected_doublet_rate from your own 10x loading chart; do not force the predicted-doublet fraction to match it exactly.
• Run doublet detection per capture/channel, never pooled. Doublets arise during physical loading — pooling unrelated captures into one Scrublet model teaches the wrong default.
• Do not remove biologically high-mito cells blindly. Cardiomyocytes, platelets, and cells under oxidative stress have genuinely high mitochondrial fractions; inspect marker genes before discarding, and keep any hard mito cap as an extreme safety net only.
• Scrublet needs unnormalised, raw integer counts. Always pass adata.layers["raw_counts"], never normalised/log/scaled/decontaminated data — and assert the non-negative, integer-like contract before calling it.