Validation — the spike-in metric#

The question a UMI pipeline has to answer is not “can you remove errors”. Anything removes errors by discarding rare sequences. The question is:

can you remove PCR and sequencing error while keeping a real variant that is rarer than the error cloud around it?

Shugay et al. 2014 built the instrument that answers it. Three known IGH clones were spiked into a real B-cell library — EHEB and two deliberate variants at one and two substitutions from it:

clone

junction (48 nt)

subs

% plasmid

EHEB

TGTGCGAGAGATGATGGCGGGGG…GACTTTGG

0

~99

EHEB-V1

TGTGGGAGAGATGATGGCGGGGG…

1

~0.5

EHEB-V2

TGTGCGACACATGATGGCGGGGG…

2

~0.05

The metric is the ratio of a real spike-in to the worst error at the same substitution distance:

  • Err1 — the most abundant junction exactly 1 substitution from EHEB, excluding V1

  • Err2 — the most abundant junction exactly 2 substitutions from EHEB, excluding V2

quantity

raw reads

standard processing

MIGEC (UMI consensus)

EHEB / Err1

362

1041–1085

9007–24696

V1 / Err1

1.35

3.1–3.8

26.5–75.9

V2 / Err2

0.28

1.7–2.0

4.6–6.2

Read the second row of numbers carefully. V2 is less abundant than the worst 2-substitution PCR error — 0.28× — so no abundance threshold can keep V2 and drop that error, because the error is bigger. V1 is only 1.35× its worst competitor. This is the regime where abundance-based error correction is provably unsolvable, and it is precisely why molecular barcodes exist: consensus assembly moves V1/Err1 from ~1.4 to 26–76, which is a change in the evidence, not in a threshold.

Note

Measured independently on SRR1200517 by the arda benchmark project, which found the same ordering on a 1.95 M-read subsample and concluded that the missing capability was UMI consensus itself. Those numbers are a subsample and are not like-for-like against the paper’s; the ordering and the conclusion do not depend on that.

Running it#

python scripts/spikein_ratio.py reads.fq.gz                     # baseline, raw reads
python scripts/spikein_ratio.py consensus.fq.gz --label migec   # after assembly

The junction is located by its conserved 3′ anchor, so this runs on raw reads with no reference and no V/J calling, and it searches both orientations.

Warning

Anchor on the 3′ end only. Requiring the junction’s first twelve bases as well looks more robust and is catastrophically wrong here: V1 differs from EHEB at position 4 and V2 at positions 7–8, i.e. inside a 5′ anchor. Both variants then count as zero and the metric looks perfect. This is covered by a test.

Acceptance gate#

V1/Err1 in 26.5–75.9 and V2/Err2 in 4.6–6.2 after consensus assembly, on SRR1200517. That is the headline claim of the rewrite, and it is not yet met — consensus assembly lands in M1. Until then the script reports the raw baseline, which reproduces the published ~1.4 and ~0.3.

Data#

PRJNA239303, runs SRR1200517–SRR1200520. Only Experiment 2 is public; the 12-clonotype TRA/TRB truth of Supplementary Table 1a has no public raw reads and lives on the cluster. See SOURCES.md.

The shallow regime, on a real repertoire library#

Every claim migec makes about 1–3 reads per UMI — that the count ratio carries nothing there, that payload agreement is what remains, that singleton filtering is unaffordable — was measured on simulated libraries and on a deep amplicon. Neither is the shape the tool is actually used on. This is: one HiSeq lane of 5′-RACE bulk TCR β from an ageing cohort, ten donors multiplexed by a 4 nt sample tag, 149,588,907 read pairs, a 16 nt UMI, published in Britanova et al., J Immunol 2014;192(6):2689–98 (10.4049/jimmunol.1302064) — retrieved from PubMed.

Demultiplexing recovered 36.0% of the lane into the ten declared samples at 2.35–2.62 reads per UMI, which is the regime. Four donors then went through correction and consensus:

sample

reads

barcodes

merged

by payload

molecules

ε at depth

A2-i129

4,832,212

1,940,846

6.12%

10,819

1,822,070

8.1e-4

A2-i131

4,398,345

1,802,598

6.21%

10,144

1,690,723

1.1e-3

A2-i132

6,190,949

2,411,942

6.23%

13,716

2,261,640

1.6e-3

A2-i133

4,921,328

1,997,079

7.18%

13,083

1,853,650

1.4e-3

assets/shallow_repertoire.tsv is the full table.

What it says, in order of how much it matters:

The payload term is not decoration. Between 10,144 and 13,716 merges per donor — about 9% of all merges — were ones the count ratio alone would have refused. At 2.5 reads per UMI a parent with two reads and a child with one is not an asymmetry, and without the reads agreeing on the molecule those barcodes stay separate and inflate the count. See refine – correct the errors in the barcode.

The distance-1 estimator correctly refused to answer. It read 1.3e-2, sixteen times the children-based estimate, and flagged itself: 10% of every barcode’s 1-substitution neighbourhood is occupied at these densities, so the excess it measures is a small difference of two large numbers. The estimate from the children of well-covered parents — 8.1e-4 to 1.6e-3, Q28–Q31 — is the one to read, and err_unreliable is what says so. See Barcode error against depth.

A 16 nt UMI was worth 12.1–13.0. The effective length is three to four bases below the nominal one, which is a 100–1000× smaller barcode space than the sheet implies, and it is why the collision correction declined: at 8.7–11.9% occupancy the space is saturated and the inversion would collapse onto the observed count. See Barcode space, collisions and the error budget.

The residual FDR is 5.0–5.7% of one-read molecules, and at ≥2 reads it is inside the 5% target. Reported, never applied.

Note: the numbers above are from the published wheel of the previous release, run on the cluster that holds the reads — so they are also a check that the artefact on PyPI works on a machine nobody developed on. Peak RSS was 932 MB for the demultiplex, of which 568 MB was the UMI counters, on a lane that is 21 GB compressed: the allocation that Speed and memory bounds, measured on the library that produces it.

MIGEC 1.2.9, the implementation this replaced#

The rewrite’s claim is that the algorithms are the specification and the code is not, so the comparison worth running is against the Groovy MIGEC this repo replaced — the same barcode dialect, the same sheet, the same simulated library, both pipelines end to end.

gh release download 1.2.9 --repo antigenomics/migec -p 'migec-1.2.9.zip'
python scripts/compare_migec_v1.py --out /tmp/v1 --jar migec-1.2.9.jar \
    --molecules 20000 --clones 200 --coverage 8 --min-count 1

20,000 molecules over 200 clones, 149,103 reads, a 12 nt barcode at a 2·10⁻³ per-base error rate. v1 is Checkout → Assemble --filter-collisions; migec 2 is checkout → refine → assemble. assets/migec_v1.tsv:

min reads

tool

consensuses

exactly a template

precision

seconds

peak RSS

1

MIGEC 1.2.9

22,717

21,338

0.9393

2.95

974 MB

1

migec 2

20,017

19,977

0.9980

0.31

267 MB

5

MIGEC 1.2.9

16,292

15,490

0.9508

3.54

941 MB

5

migec 2

16,636

16,635

0.9999

0.30

266 MB

  • 9.4–11.9× the wall clock, against an M5 gate of 3×, and 3.5–3.7× less memory. Single core in both cases; migec 2 threads and v1 does not, so this understates it.

  • The molecule count is right. 20,017 consensuses for 20,000 molecules — 0.09% over. v1 emits 22,717, 13.6% over, which is the barcode errors --filter-collisions did not catch: its rule is a count ratio, and refine – correct the errors in the barcode is where the measurement lives that says a count ratio carries nothing below ~3 reads per UMI.

  • 99.8–99.99% of consensuses are exactly a template, against 93.9–95.1%.

Warning

--min-count is applied to both. v1 defaults to 5 and migec 2 to 1 — v1 names its output .t5. for exactly this reason — so leaving each at its own default compares defaults rather than implementations, and would have credited us with recovering molecules v1 was told to throw away.

MAGERI 1.1.1, the other descendant#

MAGERI is the closer comparator of the two: same author, same UMI model, and the assembler this repo’s consensus is a rewrite of. It is also a superset — it assembles, maps and calls variants in one run — so what is compared is the part they share, and the extra work is named rather than divided out.

gh release download 1.1.1 --repo mikessh/mageri -p mageri.zip
python scripts/compare_mageri.py --out /tmp/mageri --jar mageri.jar \
    --molecules 20000 --clones 200 --min-count 1

The same library as above: 20,000 molecules over 200 clones, 149,103 reads, a 12 nt barcode at a 2·10⁻³ per-base error rate. assets/mageri.tsv:

min reads

tool

consensuses

exactly a template

precision

seconds

peak RSS

1

MAGERI 1.1.1

22,707

22,377

0.9855

2.02

2,246 MB

1

migec 2 + minimap2

20,017

19,977

0.9980

0.41

280 MB

2

MAGERI 1.1.1

19,966

19,930

0.9982

1.98

1,171 MB

2

migec 2 + minimap2

19,974

19,940

0.9983

0.39

279 MB

5

MAGERI 1.1.1

16,282

16,282

1.0000

1.98

1,146 MB

5

migec 2 + minimap2

16,636

16,635

0.9999

0.39

281 MB

  • 4.9–5.1× the wall clock, with the alignment MAGERI does folded into migec’s row: minimap2 -ax sr -y onto the same reference plus a samtools sort costs 0.08 s of the 0.39. MAGERI’s clock still also covers variant calling, which migec does not do at all.

  • 4.1–8.0× less memory. The range is the JVM’s, not ours: MAGERI’s peak RSS moves between 1,146 and 2,246 MB on the same input across runs, where migec 2 sits at 279–281 MB across all six.

  • At one read per MIG the molecule count separates them. MAGERI emits 22,707 consensuses for 20,000 molecules — 13.5% over, within a tenth of a point of MIGEC 1.2.9’s 13.6%, which is what a shared count-ratio correction lineage looks like. migec 2 is 0.09% over. At --min-count 2 the two agree to 0.04%, because a threshold of 2 discards most of what a count ratio cannot correct.

  • At five reads per MIG migec keeps 2.2% more molecules at identical recall (16,636 against 16,282, both at 0.9873 of the recoverable templates), and one of its 16,636 is not exactly a template against none of MAGERI’s 16,282. That is the trade the whole of refine – correct the errors in the barcode is about: err on precision when merging, and do not throw a molecule away to buy a fourth decimal place.

Warning

The MIG size threshold is applied to both, and it is checked afterwards rather than assumed. MAGERI’s preset carries forceOverseq=true with defaultOverseq=5, so out of the box it assembles only MIGs of five reads or more; compare_mageri.py rewrites that value from the exported preset and then reads back the threshold MAGERI reports it used, refusing to score the run if the two differ.

Note

MAGERI’s sub-clustering test runs at pcrMinorTestPValue = 0.01 — the nominal threshold Three permutation nulls measured as over-calling by 19× once both margins of the reads × positions matrix are preserved. It is left at its default here on purpose: a head to head measures what the other tool does, not what it would do if it were this one.

The same comparison at the variant level#

MAGERI calls variants; migec stops at the consensus and hands it to a caller. Stopping the comparison at the consensus therefore stops it one step short of what either tool is for. The simulator’s variant_af turns the clone set into one reference plus five point variants of it at a known allele fraction — the ctDNA shape, and the only shape in which a call set has an answer — and the two pipelines are then scored on the calls.

python scripts/compare_mageri.py --out /tmp/mageri --jar mageri.jar \
    --molecules 20000 --clones 6 --coverage 1.5 --min-count 1 \
    --variant-af 0.01 --caller lofreq

assets/mageri_variants.tsv, 20,000 molecules, 5 variants, one replicate per cell:

reads/UMI

pipeline

variant AF

called

true

false

sensitivity

precision

8

MAGERI

1%

5

5

0

1.000

1.000

8

migec 2 + LoFreq

1%

5

5

0

1.000

1.000

8

MAGERI

0.1%

5

5

0

1.000

1.000

8

migec 2 + LoFreq

0.1%

5

5

0

1.000

1.000

1.5

MAGERI

1%

142

5

137

1.000

0.035

1.5

migec 2 + LoFreq

1%

5

5

0

1.000

1.000

1.5

MAGERI

0.1%

142

5

137

1.000

0.035

1.5

migec 2 + LoFreq

0.1%

1

1

0

0.200

1.000

  • At 8 reads per UMI the two are indistinguishable: 5 of 5 with no false positive, at 1% and at 0.1%, and the called allele fraction is within 1–4·10⁻⁴ of the injected one for both. A consensus over 8 reads removes essentially all sequencing error, and what is left is set by the molecule count, which is the same number for both tools.

  • At 1.5 reads per UMI they separate by 28×, and it is the calling that separates them, not the collapsing. Both tools emit the same 20,170 consensuses at the same accuracy — 0.8990 of MAGERI’s are exactly a template against 0.8985 of migec’s — so the input to the two callers is matched to within 5·10⁻⁴. On that input MAGERI reports 142 variants of which 137 are at positions nothing was injected at; LoFreq on migec’s consensus reports 5 and is right about all of them.

  • At 0.1% and 1.5 reads per UMI the trade is explicit: MAGERI buys 4 extra true positives with 137 false ones, a positive predictive value of 3.5%. migec 2 finds 1 of 5 and is right about it. Which of those a study wants is the study’s choice, and How low can you go: exome, ctDNA and MRD is where the arithmetic for making it lives — but a caller that returns 137 false positives per sample cannot be run without a per-position background model, which is exactly what that page concludes.

Note

The variant arm ran on the cluster rather than on this laptop, because LoFreq is what migec would hand its consensus to and that is where it is installed (scripts/compare_mageri.sbatch, SLURM, 8 cores). Its wall clocks are therefore not comparable with the consensus table above, which was measured on one machine end to end; on the cluster the ratio is 2.3–3.7× including LoFreq, which migec’s row pays and MAGERI’s does not.