Barcode error against depth#

How often does a barcode base get called wrong? refine answers it twice, from two different things in the same data, and reports both — because the two fail in opposite directions and a single number would hide which failure you are looking at.

The first is the distance-1 excess: count pairs of observed barcodes one substitution apart, subtract what independent draws would produce, and invert the remainder. That is estimated_error, and it is described in UMI statistics.

The second is on this page. It reads the errors off the children correction actually found.

Two estimators, one number#

A parent barcode carrying \(c\) reads offered the instrument and the polymerase \(c L\) barcode bases to get wrong. What they got wrong is sitting in the table as that parent’s error children, so the same \(\varepsilon\) falls out of the row two ways:

\[\begin{split}u(c) &= 3L \left( 1 - e^{-c\varepsilon/3} \right) &&\Rightarrow\quad \varepsilon = -\frac{3}{c}\ln\!\left(1 - \frac{u}{3L}\right) \\ r(c) &= c\,L\,\varepsilon &&\Rightarrow\quad \varepsilon = \frac{r}{c\,L}\end{split}\]

where \(u\) is the distinct child barcodes per parent and \(r\) is the reads in those children per parent. Both are columns of <sample>.umi_errors.tsv, one row per exact parent depth.

The two differ in one respect that matters. A barcode has exactly \(3L\) neighbours one substitution away, so \(u\) saturates: past the point where a parent has spawned all of them, more errors add no new children and the estimate bends down. \(r\) counts reads, of which there is no ceiling, so it does not bend. Where the two curves separate on the figure is where this library’s barcode neighbourhood filled up — measured, not predicted.

Note

error_from_variants is left blank past saturation rather than reported as a small number. Inverting a full neighbourhood returns “no errors” for the most error-ridden library there can be, and a blank cell is harder to misread than a zero.

What it is worth, and where it stops#

Both estimators are bounded by the merges correct_umis actually made. Neither is saturation-free, and this table must not be read as though it were. Measured against a known injected rate on simulated libraries, as a fraction of the truth:

occupancy	0.2%	2.3%	9.8%	33%	100%
distance-1 excess	0.97	0.96	0.76	0.45	0.001
from the children	0.99	0.95	0.88	0.62	0.00

So the children estimate is the better of the two wherever either works, and at 0.2% occupancy on a library sequenced 25 deep it lands within 1% of the injected rate. tests/synthetic/ test_umi_errors.py is what holds those numbers; the bounds there are loose on purpose, because the claim is the trend and not a third decimal place.

Warning

At 100% occupancy both go to zero, and for the same reason: on a full barcode space correct_umis refuses to merge — correctly, because a distance-1 neighbour there is more likely a real molecule than an error child. The saturated flag in the report is what says the answer is a floor. Read the flag. Do not read this table instead of it.

Read it at depth#

A child whose parent was never sequenced cannot be merged into anything, so it is never counted. At 1–3 reads per UMI that is about 80% of all barcode errors (Grouping accuracy: Calib, UMI-tools, fgbio), which makes the shallow end of this table a lower bound rather than a measurement.

That is why the report gives two figures:

error_at_depth — restricted to parents seen at least error_depth (10) times, where correction is close to complete. This is the number to quote, and error_phred is the same thing as a Phred so it can be put beside the barcode’s own reported quality.
error_from_children — every depth, stated as the lower bound it is.

barcode error   9.73e-04 per base, estimated from the distance-1 excess
                9.89e-04 from the reads in the children of molecules seen >= 10 times -- Q30
                9.98e-04 over all depths, which is a LOWER bound: a child with no sequenced
                parent cannot be found

The library there had 1e-3 injected. Three routes to the same number, agreeing to within 3%, is what makes any of them believable — no one of them is the reference for the others.

The figures#

migec plot draws two panels from this table.

umi_error_children puts the distinct children and the reads in them against the parent’s depth, with \(3L\) as a dashed ceiling. Both series climb with \(c\); only one of them can climb forever.

umi_error_rate puts the two implied error rates against the same axis, with what refine reports drawn across them. The y axis is a log error rate, so one decade is exactly ten Phred and 1e-3 is Q30 — which is the comparison the panel exists for.

Note

Both panels are drawn with points, never lines. One row is one exact depth, and past the head of the distribution most depths hold a handful of parents, so a line would render integer quantisation as structure and bridge gaps in the support where nothing was observed at all. The same correction was made to mig_size_spectrum, where a line through one-molecule-per-size rows was drawing the \(y = x\) diagonal as the most prominent feature of the figure.

The table#

<sample>.umi_errors.tsv, one row per distinct parent depth:

column	meaning
`parent_reads`	reads carried by the parent, exact rather than binned
`parents`	surviving barcodes at that depth. Weight by this, or filter on it
`child_barcodes`	distinct barcodes merged into those parents
`child_reads`	reads carried by those children
`children_per_parent`	\(u\)
`reads_per_parent`	\(r\)
`neighbours`	\(3L\), the saturation ceiling. Constant down the column
`error_from_variants`	\(\varepsilon\) from \(u\); `.` past saturation
`error_from_reads`	\(\varepsilon\) from \(r\)
`phred_from_reads`	the same, as a Phred
`estimate`	`error_at_depth`, repeated so the panel can draw it. Constant down the column

neighbours and estimate are constant on every row deliberately. The panels draw them as reference lines, and a figure that needs a value its own table does not carry is a figure that will one day disagree with the report.