Cell Ranger: cell barcodes, cell calling, and what the cell count costs#

migec and Cell Ranger both read a 10x droplet library front to back: given R1 = 16 nt cell barcode + 10 nt UMI, which reads carry a real GEM barcode, which barcodes are cells, and how much of the library is in them. Those three questions are comparable and are compared here, on sc5p_v2_hs_PBMC_1k VDJ-T.

cat sc5p_v2_hs_PBMC_1k_t_S1_L00{1,2}_R1_001.fastq.gz > R1.fq.gz
cat sc5p_v2_hs_PBMC_1k_t_S1_L00{1,2}_R2_001.fastq.gz > R2.fq.gz
python scripts/compare_cellranger.py --r1 R1.fq.gz --r2 R2.fq.gz \
    --cellranger-dir cellranger_published/ --whitelist 737K-august-2016.txt \
    --out /tmp/cr --threads 8 --tsv assets/cellranger.tsv

Cell Ranger is run both ways. 10x publish their own Cell Ranger 5.0.0 output for this exact library, and cellranger vdj 10.1.0 was run here on the same reads (scripts/cellranger_vdj.sbatch, 16 cores) so the comparison has a cost axis and a second version to bound drift against.

Cell Ranger against itself, five major versions apart#

version

cells

contigs

shared

Jaccard

source

5.0.0

479

943

477

0.9938

published by 10x

10.1.0

478

939

477

0.9938

run here, 524.5 s / 936 MB

Version drift is negligible, so the migec-against-Cell-Ranger gap below is not a version artefact. That control is why 10.1.0 was run at all.

Note

5.0.0’s published ``Median TRB UMIs per Cell`` of 12.0 is not reproducible from its own deposited tables, under either denominator – both give 11.0. 10.1.0 reports 11.0 for TRB and 4.0 for TRA, and both are exactly the median over all cells in its own tables; 5.0.0’s TRA figure of 5.0 is the median over cells that have that chain. So 10.1.0 changed the denominator and is self-consistent, while 5.0.0’s TRB value matches neither. Derive these from the contig tables and say which denominator you used; do not quote metrics_summary for them.

The result#

tool

reads

valid barcode

cells

shared

molecules

reads in cells

Cell Ranger 10.1.0

6,301,573

90.60%

478

469

86.40%

migec, no whitelist

6,301,573

88.27%

890

469

496,373

84.26%

migec, 737K whitelist

6,301,573

88.90%

888

469

467,810

84.88%

What it costs#

stage

wall clock

peak RSS

what it produces

cellranger vdj 10.1.0, 16 cores

524.5 s

936 MB

cell calls and per-cell contigs, annotated

migec checkout + refine, 8 threads

35.0 s

679 MB

cell calls

+ assemble --contig + arda

~90 s total

679 MB

+ per-molecule consensus, annotated

Never: 35 s against 524 s is not like for like, and the table says so. Cell Ranger’s single invocation also assembles and annotates a contig per cell; the stages that answer the three axes above are 35 s of migec. The comparable end-to-end number is the third row – ~90 s against 524.5 s, still 5.8x – and it produces per-molecule resolution that Cell Ranger’s output does not have. Both numbers are on the same 6,301,573 read pairs, but on different hardware: Cell Ranger on a 16-core cluster node, migec on 8 threads of a laptop, so read the ratio and not the seconds.

migec calls 1.86x the cells and loses 1.5 points of reads-in-cells. That is the whole finding, and it says the extra barcodes are nearly empty: 419 barcodes beyond Cell Ranger’s set hold about one and a half percent of the library between them. A cell count is not an accuracy figure on its own.

The two cell sets are not nested and do not measure the same thing. migec’s gate is molecules of any sequence; Cell Ranger’s is “assembled a productive V(D)J contig”. A B cell or a monocyte with plenty of molecules is correctly a migec cell and correctly not a Cell Ranger VDJ cell. So the table carries five counts – 478, 888, 469 shared, 419 migec-only, 9 Cell-Ranger-only – and never a ratio. The 9 Cell Ranger cells migec misses sit below the OrdMag threshold on raw molecule counts.

Note: the whitelist is worth 0.63 points of validity and two cells. Snapping off-list barcodes onto the list lifts read validity from 88.27% to 88.90% against Cell Ranger’s 90.60%, and moves the cell count by two. The remaining 1.7-point gap is reads whose R1 is not a barcode at all: the deepest off-list 26-mers are phase shifts of one fixed sequence (GGTCCGTCTTGCGCCG and its rotations), which refine’s pass 0 remaps rather than drops.

Per-cell receptor chains, without a per-cell assembler#

Cell Ranger assembles a contig per cell from every read of a barcode, then annotates it. migec assembles a consensus per molecule, annotates each one with arda, and lets the cell’s chain be a vote over its molecules. Same question – which chain is in which cell – reached with and without an assembler.

migec assemble ref/PBMC.fq.gz -o asm/ --contig --min-reads 30
python scripts/compare_cellranger_chains.py --consensus asm/PBMC.consensus.fq.gz \
    --cellranger-dir cellranger_published/ --min-reads 30 --out /tmp/chains \
    --tsv assets/cellranger_chains.tsv

locus

Cell Ranger chains

migec chains

shared

chain recall

junction agreement

TRA, against 5.0.0

426

451

426

1.0000

0.9507

TRB, against 5.0.0

469

474

468

0.9979

0.9915

TRA, against 10.1.0

424

450

424

1.0000

0.9505

TRB, against 10.1.0

469

473

468

0.9979

0.9915

Every TRA chain and all but one TRB chain Cell Ranger found is also found by migec plus arda, in 22 seconds over 47,584 consensuses, with no per-cell assembly step anywhere – and the result is the same against both Cell Ranger versions, which is the point of scoring it twice. Recall is the metric that leads because a missed chain is the unrecoverable error; junction agreement is secondary and is scored only over the chains both tools called, which is a self-selecting denominator.

migec calls slightly more chains than Cell Ranger on both loci (451 against 426, 474 against 469). Those extra calls are not scored here – Cell Ranger’s set is the reference, not the truth, and a second productive TRA is allelic inclusion rather than an error.

Warning

Depth does not buy junction coverage on this chemistry, and the obvious arithmetic says it does. The tempting model places each read uniformly over the ~508 nt amplicon: a 90 nt read then spans the median 42 nt junction with probability 0.114, and 30 reads give \(1 - 0.886^{30} = 0.975\). Measured, that is wrong. Reads of one (CB, UMI) are co-terminal in 92% of 10x groups (Fragmented libraries — why 10x needs a different consensus), so a molecule is a pile at one position, not a tiling, and its consensus covers one window however deep it is. At --min-reads 30 the mean consensus is 204 nt, not 508, and 7,855 of 47,584 molecules (16.5%) carry a cell, a locus and a junction – close to the single-window 0.32 the geometry predicts, nowhere near 0.975. The depth cut is still right, because a deep pile gives a clean consensus; it just does not extend one.

Per-cell contigs, reference-free#

The window a molecule covers is one window, but different molecules start at different positions, so a cell’s molecules tile its transcript. That is measurable before any assembler exists, and it is the ceiling everything below is scored against: 99.90% of the 25-mers of Cell Ranger’s 943 filtered contigs are already present in their own cell’s raw migec molecules, and 942 of 943 CDR3 nucleotide sequences appear verbatim in one of them. The contig is in the data.

arda cells assembles it – adapter trim, 25-mer seeds, verified overlaps, union-find layout, haplotype phasing, weighted column consensus – with no germline reference until the finished contigs are annotated.

migec assemble ref/PBMC.fq.gz -o asm_all/ --contig --min-reads 1
python scripts/compare_cellranger_contigs.py --consensus asm_all/PBMC.consensus.fq.gz \
    --cellranger-dir cellranger_published/ --out /tmp/contigs \
    --tsv assets/cellranger_contigs.tsv

variant

k-mer coverage

contigs at >=90%

CDR3 exact

N50

chain recall

doublets

molecules, no assembly

0.9990

943

942 (0.9989)

–

–

–

arda cells, default

0.9759

892

933 (0.9894)

536

0.9777

17

no phasing

0.9663

868

926 (0.9820)

460

0.9714

12

no adapter trim

0.9525

879

907 (0.9618)

580

0.9491

16

Per chain at the default: TRA 454/464 (0.9784), TRB 479/479 (1.0000). 479 cells and 249,635 molecules in 23 s.

Warning

``–min-reads`` throws away the tiling. A one-read molecule is a poor consensus and one more window of the transcript, and the second thing is what the assembly needs. This axis runs at --min-reads 1 on purpose; the depth cut belongs to the per-molecule route above, not to this one.

Warning

Contig N50 is a description, never a score. The no-adapter-trim row has the highest N50 in the table, 580 against 536, and the lowest value in every other column. A contig built across an adapter is a longer contig and a wronger one.

Doublets: two chains of the same locus#

One TRA and one TRB in a droplet is a paired T cell. Two TRB is a doublet – and until the phasing landed it was invisible by construction: the two chains share their constant region, so the overlap layout puts them in one component and the column consensus averages their junctions into a sequence that is neither. On a synthetic cell built from two TRB receptors that produced one 918 nt contig with no callable junction; with the phasing on, both true junctions come back. On this library it moves chain recall 0.9714 to 0.9777 and doublet candidates 12 to 17.

What separates a real second chain from ambient RNA is productivity first and count second: of the extra chains Cell Ranger agrees with, 60/60 are productive; of those it does not, 20/128 are. An extra chain carried by exactly one molecule is contamination 96-97% of the time. The sweep that fixes the thresholds, the per-cell table and the QC panels are documented in arda’s single-cell page.

Measured, and not comparable#

Each of these was computed and then deliberately kept out of the table, because the number would have been read as a comparison and is not one.

  • Molecules per cell. Cell Ranger’s umis counts UMIs incorporated into a filtered contig – 7,623 over 479 cells, median 15 – while migec counts every molecule of any sequence, a median of 178 on the same barcodes. That is a different population, not a 12x over-count, so no ratio of the two appears.

  • migec’s “100% assigned” against “Valid Barcodes 90.6%”. Not the same measurement. On ^XXXXXXXXXXXXXXXXNNNNNNNNNN there is nothing to score, so the assigned rate is 100% by construction. The comparable quantity is the read share on the whitelist, which is what the table reports.

  • Cell Ranger’s per-contig ``reads`` as a reads-in-cells figure. Summing it gives 63.61%, which is not any published metric – it excludes reads in cells that went into no contig. The 86.80% comes from metrics_summary.csv and is not recoverable from the per-contig tables.

  • ``Median TRB UMIs per Cell``. Not cut but corrected – see the note at the top of the page. 5.0.0’s published 12.0 matches neither denominator in its own tables; 10.1.0 reports 11.0 and is self-consistent.

Two traps worth naming#

Warning

Both lanes, or neither. L001 is 3,155,166 read pairs and L002 is 3,146,407; only their sum, 6,301,573, is what Cell Ranger reports. checkout takes one R1 and one R2, so the lanes are concatenated first – and in the same order for both mates. assemble matches mates by position and refuses only on a length mismatch, so L001+L002 against L002+L001 has the right length and mis-mates every pair with nothing flagging it. The script asserts the read total before scoring anything.

Warning

Cell Ranger’s AIRR export means the opposite of arda’s by the same column name. In _airr_rearrangement.tsv, consensus_count is READS and duplicate_count is UMIs – verified on all 943 rows against the contig CSV, and the inverse of what arda writes. A join on the column names silently swaps reads for molecules, which on this library is a factor of 526. The script reads counts only from _filtered_contig_annotations.csv and asserts the identity at runtime rather than trusting it.