citeqc

Measure how well each antibody in a CITE-seq experiment is predicted by its own transcript, and check your analysis for artefacts that produce confident, biologically plausible false positives.

Documentation: https://elmasnuryilmaz.github.io/citeqc/

The short version of why this exists: a third to a half of surface markers cannot be read off the transcriptome, and which ones is predictable from molecular mechanism rather than random. CD45RA is an alternative splice product of PTPRC; an antibody distinguishes CD45RA from CD45RO, total PTPRC mRNA cannot. CTLA-4 and CD40L sit in intracellular vesicles until the cell is triggered, so their surface level is set by trafficking, not transcription.

Across four public CITE-seq datasets (~410,000 cells, 118 markers, 2017–2024), CD45RA was the lowest-coupled marker in every dataset that measured it (rank 1 of 13, 1 of 113, 1 of 29).

Install

pip install citeqc

Optional .h5mu support:

pip install 'citeqc[anndata]'

Python 3.10 or newer.

Quick start

citeqc ships a small synthetic dataset (simulated, not real data) in which every antibody's behaviour is known by construction, so you can try the whole workflow in seconds:

citeqc example citeqc_example
citeqc couple --adt citeqc_example/adt.csv --rna citeqc_example/rna.csv \
  --meta citeqc_example/meta.csv --normalise
citeqc diagnose --adt citeqc_example/adt.csv --meta citeqc_example/meta.csv \
  --groups citeqc_example/groups.csv

The quick start walks through the output, and examples/quickstart.py does the same from Python.

Measure coupling

import pandas as pd
from citeqc import coupling, mapping

adt  = ...   # cells x antibodies, CLR-normalised
rna  = ...   # cells x genes, log-normalised
meta = ...   # cells x [donor, celltype, nCount_ADT, nCount_RNA]

pairs, unmapped = mapping.map_antibodies(adt.columns, rna.columns)
result = coupling.measure(adt, rna, meta, pairs)
result["pct_rank"] = coupling.rank_percentile(result)

result has one row per antibody: median_rho is the headline number, with q25/q75 across strata and draws.

Same thing from the shell:

citeqc couple --adt adt.csv --rna rna.csv --meta meta.csv --out coupling.csv

Add --normalise if your matrices are raw counts, and --transpose if they are features × cells.

measure(..., return_raw=True) also returns every stratum × draw observation, which is what a donor-level bootstrap needs (resample donors, recompute the per-marker medians, and ask how often a marker keeps its rank).

Is a low value specific to that marker?

A correlation of 0.02 is hard to read on its own: it may be a decoupled marker, or a marker the assay cannot measure. specificity correlates every antibody with every gene in a background set, using the same stratification, grouping and depth adjustment, and reports where the cognate gene falls in that antibody's row.

matrix, summary = coupling.specificity(adt, rna, meta, pairs)
summary[["marker", "cognate_rho", "spec_percentile", "rank_of_cognate"]]

spec_percentile of 1.0 means the protein tracks its own transcript better than any other gene; about 0.5 means the cognate gene is indistinguishable from the background. The background is descriptive rather than a null (genes that truly co-vary with the antibody sit in it), and a marker can be low on measure yet specific here: weak coupling is not the same as no coupling. In our own data CD45RA's transcript, PTPRC, ranked 110th of 111 genes in the CD45RA row, while the same gene ranked first for CD45RO.

How to read the result

Absolute correlations depend on how finely you split cell types, so they are not comparable between datasets; compare pct_rank instead. Within a dataset, calibrate against a marker that should couple well in the cells you are looking at (CD8/CD8A in CD8⁺ T cells, CD19/CD19 in B cells). That value is your ceiling; markers far below it are not readable from RNA.

Three constraints the function enforces for you

1. Stratification by donor × cell type. Cell-type composition moves protein and mRNA together, so an unstratified correlation largely measures composition. This is the single biggest source of inflated coupling.
2. Equal-sized pseudo-groups. A mean over more cells is less noisy, so unequal group sizes create a size-driven gradient across every feature.
3. Library depth regressed out in both modalities.

Check your analysis

Five checks, each built on a negative control the data already contains. They answer one question: what would this analysis report if there were no signal? A check flags a problem; it does not say where the problem comes from.

from citeqc import diagnostics as dx

dx.panel_consistency(adt, meta, by="donor")   # is one panel measured everywhere?
dx.control_check(my_zscores, control_markers=["CD19", "CD20"])
dx.size_confounding(pseudobulk, sizes=n_cells_per_group)
dx.pvalue_sanity(my_pvalues)
dx.saturation_check(psi_matrix)               # bounded values: PSI, fractions

Each returns a Finding with .passed, a .summary you can print, and .detail with the numbers.

citeqc diagnose --adt adt.csv --meta meta.csv --pvalues p.csv

Exit code is non-zero if any check flags, so it drops into CI.

size_confounding deserves emphasis. In our own analysis, before correction it produced 8,881 of 11,364 genes at FDR<0.05: 8,797 "up" against 84 "down". After resampling to a fixed number of cells per group, 2 genes survived (both "down"). diagnostics.equal_n_resample does the resampling.

control_check and size_confounding are meant to be run in that order, and the second one first if the first one flags. In that analysis the four control markers scored a mean |Z| of 2.56 against 3.66 for genuine T-cell markers; equal-N resampling took the controls to 0.39, whereas regressing out an ambient score built from the other controls left them at 2.40. We first attributed the gradient to ambient signal, and that attribution was wrong: an ambient correction should be judged on control markers held out of the correction. ambient_check remains as an alias of control_check.

What each check guards against and how it decides is explained in the diagnostics guide.

The antibody → gene table

145 entries covering common TotalSeq naming, with mechanism annotations where a documented post-transcriptional mechanism exists.

citeqc map --query cd45
citeqc map --mechanism vesicular

map_antibodies tries the raw name, a normalised form (CD3--UCHT1-TSA → CD3), then a form with a trailing replicate index removed (CD38-1 → CD38). Antibodies it cannot place are returned as unmapped, never guessed. Extend it with map_antibodies(..., extra={"MyAb": "GENE"}).

Three caveats: the mapping is one-to-one, so it cannot express multi-gene or post-translationally modified epitopes; HLA-DR is mapped to HLA-DRA alone; and CD45RA and CD45RO both map to PTPRC, which is the point rather than a limitation.

What it does not do

It does not normalise for you beyond the helpers in coupling, denoise ADT (see totalVI or dsb), or correct ambient RNA (see SoupX, CellBender; note that neither addresses the antibody signal). It measures and it checks. Run the measurement on the normalisation you trust and compare ranks between them; the accompanying paper does this for CLR, log-CP10K, isotype subtraction, dsb and totalVI (under totalVI the headline marker's rank weakens, which is the kind of result this comparison is for).

Development

pip install -e '.[test,anndata]'
pytest

The tests include synthetic data where the coupled and decoupled markers are known by construction. tests/validate_on_real_data.py (run separately; it needs the analysis repository's processed data) checks that the package reproduces the prototype's published results on real data, to Spearman ρ = 1.000.

Release steps for maintainers are in the releasing guide; changes are listed in CHANGELOG.md.

Citing

If you use this, please cite it using CITATION.cff, the accompanying paper (in preparation), and the datasets you analyse.

License

MIT