Variant calling: which caller, and what it can possibly see

migec is not a variant caller and will not become one. This page answers the question people actually arrive with – I want rare variants out of ctDNA or a tumour, what do I run after assemble? – and the question underneath it, which is the one that decides the answer: how many molecules does the caller get, and is that enough for the frequency you are chasing?

Two rules, then the numbers.

Compose or replace, never both

Downstream: what consumes the consensus draws the line for aligners and quantifiers: a tool that transports the barcode composes with migec, a tool that deduplicates on it replaces a stage of migec. Variant callers fall on both sides of that line, and running one from the wrong side is the most common way to lose molecules silently.

caller	after assemble?	why
Mutect2, LoFreq, FreeBayes, VarDict, bcftools	yes	they call from a BAM and never look at RX. One consensus record is one molecule, so their depth is a molecule count and their allele fraction is a molecule fraction
UMI-VarCal	no – it replaces assemble	it does its own UMI pileup and its own consensus. Feed it raw reads with the UMI, or feed migec’s consensus to a standard caller. Not both
UMIErrorCorrect	no – it replaces checkout + refine + assemble	it aligns first, then groups on (position, UMI) at edit distance <= 1 and consenses. It is an alternative pipeline, not a stage
DREAMS-vc, Shearwater	no as normally run	both fit a per-position error model across a panel of normals built from reads. Run on a consensus the model is fitted to a different noise process than the one it will see

Running a UMI-aware caller on a consensus counts each molecule once and then collapses the result again. Nothing errors; the molecule count just quietly drops.

Warning

Never: do not derive an allele fraction from a consensus BAM’s read depth and apply a UMI-aware caller’s family-size filter. The family size is 1 by construction after assemble – every record already is a family – so a --min-family-size 3 filter discards the entire library and reports zero variants without an error.

Warning

Never: turn a caller’s downsampling off before you run it on an amplicon consensus. Mutect2 defaults to --max-reads-per-alignment-start 50, and on a panel every molecule starts at the same coordinate – so on the certified ctDNA arms it saw 78-207 molecules where LoFreq saw 4,282-15,895 on the identical BAM. A 1% variant was then supported by one or two molecules and it reported about one call per sample, never the certified hotspot. That reads as “Mutect2 cannot do this” and is really “Mutect2 was shown 1.5% of the evidence”. Pass --max-reads-per-alignment-start 0. The general form of the trap: a default tuned for whole-genome reads, applied to records that are already molecules.

What to run

Warning

A standard caller on consensus reads needs a background model, not just a threshold. Run over the 0%-certified arm of a reference series, LoFreq on migec consensus returns 9-11 calls per sample at 0.4-1.4% VAF; 94% are -> G, eight positions recur in 3 of 3 replicates, and one of them is the PIK3CA H1047R hotspot itself. That is 2-colour chemistry’s dark-G bias, and no consensus removes it because it is systematic rather than random. See How low can you go: exome, ctDNA and MRD. What removes it is a per-position background built from normals – which is the real contribution of the UMI-aware and panel-of-normals callers below.

If you have UMIs and you are running migec, collapse first and use a standard caller:

migec checkout reads.fq.gz --bc-pattern '0:12' --sample S1 -o co/
migec refine   co/S1.fq.gz -o rf/
migec assemble rf/S1.fq.gz -o as/
minimap2 -ax sr -y ref.fa as/S1.consensus.fq.gz | samtools sort -o S1.bam
samtools index S1.bam

# then, in decreasing order of how much the choice matters:
lofreq call   -f ref.fa -o S1.vcf S1.bam                    # balanced, and it is the
                                                            # specificity end of the trade
gatk Mutect2  -R ref.fa -I S1.bam -O S1.vcf                 # more sensitive, needs
                                                            # FilterMutectCalls after it

The published comparison behind that ordering is Maruzani et al. 2024, which benchmarked six callers on ctDNA at 0.5-7.5% VAF (sources; full citation below). Their finding, in one line each:

caller	what the benchmark found
LoFreq	fewest putative false positives of any standard caller at every depth tested, and second only to Mutect2 on sensitivity. The balanced default
Mutect2	highest sensitivity; the most privately-called variants of any caller, which is the false-positive signature. Balanced once FilterMutectCalls runs – which their comparison deliberately did not run
UMI-VarCal	fewest putative false positives of all callers on UMI-encoded data. Also the fewest calls overall – 24 against 234-1,728 for the others on one sample
UMIErrorCorrect	most sensitive at the lowest VAF at every depth; false positives climb steeply with depth
FreeBayes	false-positive rate rises with allele frequency, unlike every other caller tested
bcftools	called nothing below ~8% VAF at any depth. Not a low-frequency caller and not written as one

Note

That benchmark ran every caller at default parameters and without base quality score recalibration, and skipped FilterMutectCalls on purpose, to keep the comparison even. It is a fair ranking of defaults, not a ceiling for any single tool. Mutect2 in particular is being scored without the filter its own authors require.

If you do not have UMIs, none of this applies and no consensus is possible: use LoFreq for a balanced call set or Mutect2 for sensitivity, and accept that the floor is the PCR error rate.

The number that decides it

Caller choice is a second-order effect. The first-order effect is whether the variant is present in enough molecules to be called by anything, and that is set by the wet lab:

molecules at the site  =  input DNA / 3.3 pg  x  strands recovered  x  efficiency
variant molecules      =  molecules at the site  x  VAF

A caller needs some minimum number of supporting molecules – three is typical – and the supporting count is a Poisson draw, not a guarantee. If the expectation is 3, a third of replicates see fewer than 3, and no threshold setting recovers a molecule that was never sampled.

This is why assemble reports molecules rather than reads, and why the RT floor (The RT/PCR floor — the quality no consensus can beat) is a per-molecule rate: an error made before amplification is in every read of that molecule and no consensus removes it. At a floor of 1e-4, 10,000 molecules at a site expect one false variant molecule from the chemistry alone – and that, not the sequencer’s Q30, is what a caller’s specificity is fighting.

Measured on cell-free DNA reference material

scripts/ctdna_titration.py runs the three migec stages over 100 runs of two SiMSen-Seq studies on commercial cfDNA reference material with certified mutant allele frequencies (sources; both are in SOURCES.md):

• PRJNA788522 – 0% (WT), 0.125%, 0.25% and 1% VAF, crossed with 5/20/80 ng input and 3.3/10/30x reads per UMI, three replicates each.

• PRJNA507366 – six polymerases on the same material, plus 0.031% and 0.0625% VAF.

The WT arms are true negatives: 33 runs and 6.1 million molecules in which the variant frequency is zero by construction. That is a false-positive floor measured on real chemistry rather than on a simulation.

python scripts/sra_fetch.py get SRR17220921 SRR17220924 ... -o simsen/
python scripts/ctdna_titration.py --reads simsen/ --out ctdna/

migec over both titrations, averaged over replicates

study	enzyme	input_ng	arm	depth	n	amplicons	reads	molecules_per_amplicon	reads_per_molecule	barcode_phred	variant_molecules	p_enough
PRJNA507366	Accuprime		WT		3	3	1.95274e+06	79751.7	8.09	30.8233	0	0
PRJNA507366	Accuprime_hifi		WT		3	3	2.25763e+06	93253.7	8.07	30.8367	0	0
PRJNA507366	Phusion		WT		3	3	1.84092e+06	95455.7	6.43667	31.3967	0	0
PRJNA507366	PlatRegular	80	WT		3	3	4.22022e+06	57728	24.3433	28.8667	0	0
PRJNA507366	PlatSuperfi	80	WT		3	3	3.91793e+06	134430	9.73667	29.37	0	0
PRJNA507366	Platinum		WT		3	3	2.15418e+06	86146.7	8.35667	31.0067	0	0
PRJNA507366	Platinum superfi	50	0.125%		1	3	3.8665e+06	73286	17.59	34.78	91.61	1
PRJNA507366	Platinum superfi	70	0.031%		3	3	4.07284e+06	89542	15.1367	35.1067	27.7567	1
PRJNA507366	Platinum superfi	70	0.0625%		3	3	3.60725e+06	86139	13.89	34.9567	53.8367	1
PRJNA507366	Platinum_Hifi		WT		3	3	1.81798e+06	80859.3	7.54	33.9867	0	0
PRJNA788522		20	0.125	10.0	3	5	1.06487e+06	11960.7	17.7367	28.6467	14.9533	0.999933
PRJNA788522		20	0.125	3.3	3	5	321953	7162.33	8.98667	27.65	8.95333	0.9935
PRJNA788522		20	0.25	10.0	3	5	1.15595e+06	13349.7	17.1167	29.02	33.3733	1
PRJNA788522		20	0.25	3.3	3	5	289023	8048.67	7.14333	27.7033	20.12	1
PRJNA788522		20	1	10.0	3	5	445961	8630	10.3233	27.27	86.3	1
PRJNA788522		20	1	3.3	3	5	131423	5970.67	4.40667	25.4833	59.7067	1
PRJNA788522		20	WT	10.0	3	5	1.07623e+06	13988.7	15.3567	29.6467	0	0
PRJNA788522		20	WT	3.3	3	5	334784	8479	7.89	27.3833	0	0
PRJNA788522		20	cell_line	10.0	3	5	813932	10299.3	15.7467	29.9133
PRJNA788522		20	cell_line	3.3	3	5	232375	6310	7.39333	26.27
PRJNA788522		20	cell_line	30.0	3	5	1.9048e+06	16808.7	22.6567	30.9867
PRJNA788522		5	0.125	10.0	3	5	434820	4530	20.02	29.03	5.66333	0.894
PRJNA788522		5	0.125	3.3	3	5	182150	2699	13.47	28.1	3.37	0.654933
PRJNA788522		5	0.25	10.0	3	5	909770	5899.33	28.9033	29.14	14.75	0.9994
PRJNA788522		5	0.25	3.3	3	5	182433	2945.67	12.32	28.17	7.36333	0.976367
PRJNA788522		5	1	10.0	3	5	501350	4129	24.2467	29.0267	41.29	1
PRJNA788522		5	1	3.3	3	5	167841	2888.33	11.5833	28.31	28.8833	1
PRJNA788522		5	WT	10.0	3	5	533736	4428	24.04	29.21	0	0
PRJNA788522		5	WT	3.3	3	5	146169	2868.67	10.1533	28.4567	0	0
PRJNA788522		5	cell_line	10.0	3	5	304544	4559.33	13.5433	27.53
PRJNA788522		5	cell_line	3.3	3	5	76465.3	2958.33	5.18	20.6033
PRJNA788522		5	cell_line	30.0	3	5	635122	6664.33	19.0733	29.2467
PRJNA788522		80	cell_line	10.0	3	5	1.49859e+06	26793	11.2	31.46
PRJNA788522		80	cell_line	3.3	3	5	582020	19310	6.02667	29.8167

molecules_per_amplicon is the evidence a caller gets at one site; variant_molecules is that times the certified VAF; p_enough is the probability that at least three of them are really there. assets/ctdna_titration_runs.tsv beside it has all 100 runs individually.

Where detection stops being about the caller

Calling variants on the consensus and scoring against the certified frequencies gives the honest performance of this panel and pipeline. The locus is PIK3CA H1047R (3:179234297 A>G), twelve runs per arm:

certified VAF	detected	sensitivity	mean measured	verdict
undiluted (~3.6%)	24/24	100%	3.61%	reliable
1%	11/12	92%	0.93%	reliable, and accurate
0.25%	4/12	33%	0.23%	accurate when seen, mostly missed
0.125%	1/12	8%	0.04%	effectively undetectable
0% (true negative)	3/12	– (25% false positive)	0.16%	calls a variant that is not there

Read the last two rows together. Below 1% VAF sensitivity collapses faster than specificity does: at 0.125% the assay finds 8% of real variants while still calling 25% of the true negatives. That is the worst possible shape, and it is not fixed by choosing a different caller – How low can you go: exome, ctDNA and MRD shows the artifact is systematic, reproducible and -> G biased.

Never: that table is ``–min-reads 1``, which is the wrong setting for variant calling, and it is why the numbers looked so poor. Re-run at --min-reads 3 (How low can you go: exome, ctDNA and MRD), three replicates per arm at 20 ng / 10x:

truth	mr=1 called	mr=1 VAF	mr=3 called	mr=3 VAF
0%	3/3 wrong	0.66%	0/3 correct	–
0.125%	1/3	0.17%	1/3	0.05%
0.25%	3/3	0.79% (3.2x too high)	3/3	0.22%
1%	3/3	1.01%	3/3	1.02%

Note: the artifact was inflating true positives, not only inventing false ones. At 0.25% it read 0.79%, and 0.25% + the 0.57% artifact floor is 0.82% – the contamination is additive, so a quantitative result was wrong by threefold at a frequency where the call itself looked fine.

The reliable limit is 0.25%, quantified accurately, once singleton molecules are excluded. That is fourfold better than the same pipeline at --min-reads 1 and in the range Illumina specify for TruSight Oncology 500 ctDNA v2 (0.2% for SNVs at 20 ng). 0.125% remains out of reach here at 1 of 3 replicates – that arm is molecule-limited, and the published <0.1% claims for this chemistry assume both more input and a per-position background model.

Note: depth does buy molecules, until it does not. Deeper sequencing recovers more of the molecules that are in the tube – 20 ng of undiluted material gives 6,310 / 10,299 / 16,809 molecules per amplicon at 3.3 / 10 / 30x – but the ceiling is the number of input molecules, and past it further reads only raise reads-per-molecule. The honest summary is that depth and input are both worth spending on until the molecule count stops rising.

Never: the molecule total of a multiplex panel is not the count at a site, and dividing by the panel size is not enough either. Aligning to GRCh38 and counting per target shows the weakest one holds 0.09-0.64x the panel mean (median 0.36) – so an average overstates the weakest target by up to elevenfold. The panel is TP53 x2, PIK3CA x2 and KIT, inferred from coverage; How low can you go: exome, ctDNA and MRD has the method and the off-target share that comes with it.

The barcode error tracks polymerase fidelity

PRJNA507366 varies the polymerase while holding template and protocol fixed, which makes it an independent check on refine’s barcode-error estimate (Barcode error against depth): a higher-fidelity enzyme should miscall fewer barcode bases, and the estimator should say so without being told which enzyme it is looking at.

polymerase	barcode Phred	error per base	runs
Accuprime	30.82	8.3e-4	3
Accuprime HiFi	30.84	8.2e-4	3
Platinum	31.01	7.9e-4	3
Phusion	31.40	7.2e-4	3
Platinum HiFi	33.99	4.0e-4	3
Platinum SuperFi	34.78-35.11	3.1-3.3e-4	7

The estimator separates the high-fidelity enzymes from the standard ones by 3 Phred, a factor of two in error rate, and puts Platinum HiFi and Platinum SuperFi together at the top – from the barcodes alone, with no reference, no alignment and no knowledge of the experiment. Note: these are within-study comparisons; the 80 ng paired-enzyme arm used a different library prep and its absolute numbers are not comparable with the rows above.

Why this dataset and not the published ctDNA benchmark

Maruzani et al. benchmarked on SRR10296599 and eight metastatic breast cancer runs, and had to generate the UMIs in silico – 9 nt, Phred fixed at 37, assigned by Poisson to reads sharing start and end positions – because neither deposit kept its real ones. We confirmed that independently: both runs report nreads=2 with no index read, migec suggest finds no barcode in either mate, and vdb-dump shows an empty linkage group. scripts/sra_fetch.py probe is the one-line version of that check.

Two consequences for how far their UMI results carry:

• Their synthetic UMIs cannot be wrong. A random 12-mer with a fixed Phred has no error rate, no collisions beyond chance, and no ambiguous parent. On the real reference material above, migec measures the barcode error at Q30 – so roughly 1.3% of 12 nt barcodes carry at least one miscalled base, and correcting them is most of what separates one UMI-aware caller from another. The benchmark removed the variable it was measuring.

• Their assignment rule assumes reads of a molecule are co-terminal, which X1 measured as false (7.8% of 10x groups overall). cfDNA has preferred cut sites so it is less wrong for a capture panel than for 3’ GEX, but it is not free.

This is not a criticism of the ranking – it is the best available comparison of caller defaults, and we quote it above. It is a statement about which half of it transfers: the standard-caller ordering rests on real reads and carries over; the UMI-aware arm rests on UMIs that were invented after the fact.

Note

Public ctDNA data with real, recoverable UMIs does exist – it just was not the data that benchmark used. PRJNA507366 and PRJNA788522 both carry a 12 nt inline UMI that survived deposition, and migec suggest recovers it from base composition alone with no prior knowledge of the protocol. SOURCES.md records both.

Citations

• Maruzani R, Brierley L, Jorgensen A, Fowler A. Benchmarking UMI-aware and standard variant callers for low frequency ctDNA variant detection. BMC Genomics 2024;25(1):827. doi:10.1186/s12864-024-10737-w, PMID 39227777.

• Sater V, Viailly PJ, Lecroq T, Prieur-Gaston E, Bohers E, Viennot M, Ruminy P, Dauchel H, Vera P, Jardin F. UMI-VarCal: a new UMI-based variant caller that efficiently improves low-frequency variant detection in paired-end sequencing NGS libraries. Bioinformatics 2020;36(9):2718-2724. doi:10.1093/bioinformatics/btaa053, PMID 31985795.

• Sater V, Viailly PJ, Lecroq T, Ruminy P, Berard C, Prieur-Gaston E, Jardin F. UMI-Gen: A UMI-based read simulator for variant calling evaluation in paired-end sequencing NGS libraries. Comput Struct Biotechnol J 2020;18:2270-2280. doi:10.1016/j.csbj.2020.08.011, PMID 32952940.

• Osterlund T, Filges S, Johansson G, Stahlberg A. UMIErrorCorrect and UMIAnalyzer: Software for Consensus Read Generation, Error Correction, and Visualization Using Unique Molecular Identifiers. Clin Chem 2022;68(11):1425-1435. doi:10.1093/clinchem/hvac136, PMID 36031761. The source of PRJNA788522, the titration above.

• Filges S, Yamada E, Stahlberg A, Godfrey TE. Impact of Polymerase Fidelity on Background Error Rates in Next-Generation Sequencing with Unique Molecular Identifiers/Barcodes. Sci Rep 2019;9(1):3503. doi:10.1038/s41598-019-39762-6, PMID 30837525. The source of PRJNA507366.

• BesenbacherLab/UMIseq_variant_calling – a gwf workflow comparing Shearwater, Mutect2, VarScan2 and DREAMS-vc on cfDNA, with a panel of normals built first. Read as the design reference for the panel-of-normals approach; not run here.