cellqc: standardized quality control pipeline of single-cell RNA-Seq data

Cellqc standardizes the quality control of single-cell RNA-Seq (scRNA) data, turning Cell Ranger output into clean feature count matrices. It is implemented in Snakemake for reproducibility and scalability.

The pipeline starts from the Cell Ranger filtered matrix and, per sample:

1. Ambient RNA — SoupX (default) or DecontX estimates background contamination and subtracts it. Other methods can be run alongside for comparison without touching the counts.
2. Filtering — cells are removed on total UMI, detected genes and mitochondrial percentage, with every exclusion attributed to a specific criterion.
3. Doublets — DoubletFinder and/or scDblFinder. All callers score every cell; one configured caller decides removal.
4. Nuclear fraction — the intronic read fraction per cell, from the Cell Ranger BAM, computed when a BAM is present. Reported, not used for filtering.

Output is .h5ad matrices, a self-contained HTML report, and a presentation-ready PDF slide deck.

Cell calling is Cell Ranger EmptyDrops; cellqc does not re-call cells. Cell-type annotation is out of scope as of v0.2.0 — annotate downstream.

The diagram is generated from source: bash docs/make_figures.sh renders it from docs/workflow.dot.

Installation

From conda (recommended — this pulls the whole analysis stack):

mamba create -n cellqc -c conda-forge -c bioconda cellqc
conda activate cellqc

# DoubletFinder is not packaged for conda; see below
Rscript -e "remotes::install_github('chris-mcginnis-ucsf/DoubletFinder', upgrade=FALSE)"

From the environment file, if you want the exact development environment or are working from a clone:

mamba env create -n cellqc -f envs/cellqc.yaml
conda activate cellqc
Rscript -e "remotes::install_github('chris-mcginnis-ucsf/DoubletFinder', upgrade=FALSE)"
pip install -U cellqc          # or `pip install -e .` from a clone

pip install cellqc on its own installs the CLI and the workflow, but not the analysis stack: scanpy, pysam and the entire R side come from conda, because pip cannot install R packages. Use one of the two routes above.

If you would rather avoid the GitHub build entirely, set doublet.run: [scdblfinder] and doublet.decider: scdblfinder in the config; scDblFinder comes from bioconda.

v0.2.0 removed five of the six GitHub builds v0.1.0 needed (SeuratDisk, harmony, scPred, DropletQC and the lijinbio/DoubletFinder fork) and all four version pins (Seurat v4, r-matrix, pandas<2, anndata).

Dependent software:

To test the installation:

cellqc -h

Run the pipeline

cellqc requires a sample file and an optional configuration file.

• The sample file (e.g. samples.txt) is tab-delimited with headers sample, cellranger, and optionally nreaction.

  • sample is the sample ID.
  • cellranger is the Cell Ranger output directory. Relative paths are resolved against the sample file's directory.
  • nreaction is the number of reactions in the library prep, used to infer the expected doublet rate when one Cell Ranger run combines several reactions. Defaults to 1.

• The configuration file is YAML and optional. The defaults are:

seed: 42                  # every stochastic step is seeded; v0.1.0 seeded nothing
ambient:
  method: soupx           # soupx | decontx | none -- the ONE method applied to the counts
  compare: []             # e.g. [decontx] -- estimated and reported, never applied
nuclear_fraction:         # runs automatically when the sample has an indexed BAM
  numthreads: 12
  cbtag: CB
  retag: RE
  exontag: E
  introntag: N
filterbycount:
  mincount: 500
  minfeature: 300
  mito: 10
doublet:
  run: [doubletfinder, scdblfinder]   # callers to execute
  decider: doubletfinder              # the single caller whose call removes cells
  findpK: false
  numthreads: 5
  pK: 0.01
  rate: 0.1               # 10x multiplet rate at `capacity` cells recovered
  capacity: 13000

Inspection of configuration

1. ambient — ambient RNA correction

2. nuclear_fraction

Fraction of intronic reads per cell, intronic / (intronic + exonic), computed from the Cell Ranger BAM with pysam. There is no skip flag: the step runs for any sample with an indexed possorted_genome_bam.bam and is dropped for those without, so mixed cohorts work. The result is reported and plotted against log10(UMI) but is not used for filtering — DropletQC-style empty-drop and damaged-cell thresholds are sample- and tissue-dependent, so applying them automatically would be unreviewed auto-filtering.

3. filterbycount

4. doublet

There is no skip flag, for the same reason nuclear_fraction has none: what runs is the list of callers, and a caller you do not want is left out of doublet.run. Doublet detection itself always runs.

Why the expected doublet rate is linear in cell yield

Cells are loaded into GEMs at limiting dilution, so the number of cells per droplet is Poisson with mean λ = (cells loaded) / (number of GEMs). Among droplets that contain at least one cell, the fraction holding two or more is

P(≥2 | ≥1) = 1 − λ / (e^λ − 1)  ≈  λ/2      for small λ

λ is proportional to how many cells were loaded, and the cells recovered are proportional to λ as well, so over the loading range the instrument supports, the multiplet fraction is proportional to the number of cells recovered. That is why the multiplet rate is quoted as a rate per thousand cells rather than as a single number: 10x Genomics user guides give ≈0.8% multiplets per 1,000 cells recovered (≈8% at 10,000 cells), and scDblFinder's default dbr uses the same rule of thumb at ≈1% per 1,000 cells captured. Bloom (2018) derives the Poisson treatment exactly, including the correction needed when the mixed cell types are not in equal proportion.

doublet.rate and doublet.capacity are the two ends of that line: rate multiplets at capacity cells recovered. The defaults (0.1 at 13,000) give 0.77% per 1,000 cells, i.e. the 10x specification, and reproduce v0.1.0's hard-coded constants exactly. To use scDblFinder's 1% per 1,000 instead, set rate: 0.1, capacity: 10000.

Two limits are worth knowing. The linear form is the small-λ limit: the exact Poisson expression bends below the line as loading increases (at λ = 0.2 it is 9.7% rather than 10%), so the linear rule slightly over-estimates at high yields. And nreaction divides the fraction because pooled reactions are separate emulsions — a cell from one reaction cannot share a droplet with a cell from another.

References:

• Bloom JD (2018) Estimating the frequency of multiplets in single-cell RNA sequencing from cell-mixing experiments. PeerJ 6:e5578. https://peerj.com/articles/5578/
• 10x Genomics Chromium Single Cell reagent user guides / technical notes, multiplet rate vs targeted cell recovery (e.g. CG000422).
• McGinnis CS, Murrow LM, Gartner ZJ (2019) DoubletFinder. Cell Systems 8:329–337 — takes nExp from the 10x multiplet-rate table. https://doi.org/10.1016/j.cels.2019.03.003
• Germain P-L et al. (2021) Doublet identification in single-cell sequencing data using scDblFinder. F1000Research 10:979 — "roughly 1% per 1000 cells captured". https://f1000research.com/articles/10-979/v2

Both callers are given the same expected doublet rate, so a difference between them reflects the methods rather than differing priors. Their concordance (2×2 table and Cohen's κ) is reported. Concordance is a consistency measure, not an accuracy measure — with no ground-truth doublets, neither caller can be shown superior on real data.

Note that homotypic doublets are not modelled (modelHomotypic is deliberately not called), so the expected count over-estimates the detectable doublet count and the step removes slightly more cells than the true heterotypic count. The bias direction is known, constant, and stated in every report.

Result files

Per-stage outputs (ambient/, barcoderank/, nuclear_fraction/, filterbycount/, doubletfinder/, scdblfinder/) keep the statistics tables and figures. Every figure is written as a vector PDF with editable text alongside a 300 dpi PNG for the HTML report.

The intermediate matrices (filterbycount/{sample}.h5ad, filterdoublet/{sample}.h5ad) are working files. filterdoublet/'s is marked temp and deleted once result/{sample}.h5ad is written: it held the same cells and the same counts, differing only in the barcode prefix and the nuclear-fraction columns, so keeping it wrote every count matrix to disk twice. To keep it, run the workflow through Snakemake directly with --notemp — the cellqc CLI does not pass Snakemake flags through.

An example

One sample

No sample file needed — -D writes one for you:

cellqc -d out -t 8 \
  -D sample:=:S1 \
  -D cellranger:=:/path/to/cellranger/S1/outs

The cellranger path must be absolute here: -D writes out/samples_<timestamp>.txt, and relative paths in a sample file are resolved against that file's directory, which is the outdir. Add -D nreaction:=:2 if the run pooled more than one 10x reaction, and -c config.yaml to change any threshold. That gives:

out/result/S1.h5ad            the final QC'd matrix
out/result/S1_obs.txt.gz      per-cell QC metrics and doublet scores, indexed by barcode
out/result/S1_var.txt.gz      the feature table, indexed by gene
out/result/report.html        self-contained QC report
out/result/report_slides.pdf  slide deck

Equivalently, with a one-line sample file — this is the form to prefer, because the file is a record of what was run and relative paths work in it:

sample	cellranger	nreaction
S1	/path/to/cellranger/S1/outs	1

cellqc -d out -t 8 -- samples.txt

A cohort

A sample file (e.g. samples.txt) for two samples:

sample	cellranger	nreaction
AMD1	/path/to/cellranger/AMD1/outs	1
AMD2	/path/to/cellranger/AMD2/outs	1

Run it with the installed entry point:

cellqc -d out -t 8 -- samples.txt                 # default parameters
cellqc -d out -t 8 -c config.yaml -- samples.txt  # customized parameters
cellqc -d out -t 8 -n -- samples.txt              # dry run; writes out/config_<timestamp>.yaml

The dry run writes the fully resolved configuration, defaults included, to outdir/config_<timestamp>.yaml — copy that file, edit it, and pass it back with -c.

To see the jobs Snakemake will run before running them, use the dry run above; snakemake --dag renders the graph itself if you want a picture of a particular cohort.

Example outputs from the reference run (GSE188280, 13,559 cells) are in docs/tests/: report.html, report_slides.pdf and metrics.csv.