Worked example: the calibration computed by hand

The main text derives the method; this supplement follows one data set through every step of it, by hand. The example is deliberately small, but it is constructed so that each step has an effect: the bound τ must widen once for lack of pairs and is once stopped by biological drift, and the three laboratories do not agree perfectly, so that the least-squares solution is a genuine compromise between them. Equations of the main text are cited as Eq. (1)–(9); those of this supplement are numbered (S1), (S2), ….

S1The data

Three laboratories (labA, labB, labC) measured a phenotype (an MIC, in mg/L) on 24 isolates of 12 strains. Each strain was measured by two of the laboratories.

The file examples/worked_phenotypes.tsv has one row per isolate, S1_A  1  labA and so on; the isolate name is the strain followed by the laboratory.

The genomes are given as cgMLST allelic profiles, 12 loci per isolate, as chewBBACA [1] writes them (examples/worked_profiles.tsv). Each number is the id of the allele found at that locus. Highlighted: INF- alleles (new alleles inferred by chewBBACA), missing calls and alleles that differ from the rest of the strain.

S2From allelic profiles to distances

The distance between two isolates is the number of loci where both have a call and the two alleles differ. Each call is read as follows:

• INF-1 is allele 1: chewBBACA adds the prefix when the allele is new, but it is an allele like any other.
• LNF, ASM, NIPH, PLOT3, … and 0 or an empty cell are missing calls. A locus missing in either isolate is skipped, so a missing call never counts as a difference.

Isolates S7_A and S7_C, locus by locus:

Eleven loci are compared and two differ, so \(d(\mathrm{S7\_A}, \mathrm{S7\_C}) = 2\). For S1_A and S1_B the only unusual call is INF-1, which is allele 1, so \(d = 0\): the two isolates have identical profiles. Isolates of different strains differ at every locus they share, so their distance is 11 or 12.

This is the distance of Eq. (3), computed in the same way by chewBBACA, cgmlst-dists [2] and ReporTree [3]. phenocal computes it only for the isolates that have a phenotype, and only between isolates of different laboratories, since those are the only pairs it uses.

S3Cross-cohort pairs and log-ratios

Values are compared on the \(\log_{10}\) scale, where a constant factor between laboratories becomes a constant difference:

\(\mu_i\) is the true value of strain \(i\), \(\delta_c\) the offset of laboratory \(c\) (what we want), \(\varepsilon\) the measurement noise. For a pair of isolates from laboratories \(a\) and \(b\) the log-ratio is

If the two isolates are the same strain, \(\mu\) cancels and \(s\) measures \(\delta_a - \delta_b\) plus noise. For example S1_A = 1 and S1_B = 2: \(s = \log_{10}(1/2) = -0.301\), so on this strain labB reads twice as high as labA.

The 9 isolates of labA, 6 of labB and 9 of labC give \(9\times 6 + 9\times 9 + 6\times 9 = 189\) cross-cohort pairs, grouped by the two laboratories they connect. Pairs within one laboratory say nothing about the offsets and are not used.

S4One edge per pair of laboratories: τ, Δ, SE

For each pair of laboratories phenocal decides which pairs to trust, summarises them in a single value \(\Delta\), and attaches a standard error. Three quantities are fixed by the options of the run (here the defaults, with --min-support 3):

• \(\sigma = \log_{10}(2)/\sqrt 2 = 0.2129\), the measurement resolution of a laboratory working with two-fold dilutions (--sigma fixed);
• \(\lambda = \log_{10}(2)^2 = 0.0906\), how much variance one allele of distance adds (--lambda);
• \(\kappa = 1\), how far the median may drift, in units of \(\bar\sigma = \sqrt{(\sigma_a^2+\sigma_b^2)/2} = 0.2129\) (--max-drift).

The rule for τ. Sort the pairs by distance. Let \(m_0\) be the median of \(s\) over the closest pairs (distance \(d_{\min}\)). Try the distinct distances in increasing order; τ is the first one for which

The first condition asks for enough pairs; the second stops τ when more distant pairs start to pull the median away, the signature of biology rather than protocol. The edge value and its uncertainty are then

A median of an even number of values is the mean of the two middle ones.

labA – labB

One strain (S1) is identical in the two laboratories; two more (S2, S3) differ by one allele. At τ = 0 there is a single pair, fewer than 3, so τ widens to 1. The median does not move, so τ = 1 is accepted. Note that S3 reads the same in both laboratories: one noisy pair does not move a median of three.

The other 51 pairs of these two laboratories are different strains, at distance 11 or more.

a factor of \(10^{-0.301} = 0.50\), and

labA – labC

The interesting edge. S4 and S5 are identical in the two laboratories and both say that labC reads four times higher (\(s = -0.602\)). That is two pairs, not enough. At distance 2 come S6 and S7, which read the same in both laboratories: they are not quite the same strain, and their biology differs. With them the median moves by a whole dilution (0.301), 1.41 times σ̄, more than κ = 1 allows, so τ = 2 is rejected. At τ = 3 two more pairs bring the median back, and τ = 3 is accepted with six pairs.

The other 75 pairs of these two laboratories are different strains, at distance 11 or more.

a factor of \(10^{-0.602} = 0.25\), and

labB – labC

Three strains (S10, S11, S12) are identical in the two laboratories, so τ = 0. Two say that labC reads four times higher than labB, one says twice as high; the median is \(\log_{10}(1/4)\). With τ = 0 and three pairs this is the most precise edge, with twice the weight of the others.

The other 51 pairs of these two laboratories are different strains, at distance 10 or more.

a factor of \(10^{-0.602} = 0.25\), and

S5Putting the edges together: weighted least squares

Each edge is an equation in the offsets, \(\delta_a - \delta_b \approx \Delta_{ab}\). Fixing the anchor \(\delta_A = 0\) (any laboratory can be the anchor; it only sets the reference scale):

Three equations, two unknowns, and they disagree: going round the triangle, the edge values should add up to zero, but

The labA–labB edge says labB reads twice as high as labA; the other two together say it reads the same. Weighted least squares finds the offsets that minimise the weighted sum of squared disagreements:

Setting both derivatives to zero gives the normal equations \((A^{\top} W A)\,\delta = A^{\top} W \Delta\):

The matrix on the left is the weighted graph Laplacian with the anchor removed. Its determinant is \(49.66^2 - 33.11^2 = 1370.0\), and inverting it,

As factors: labB reads \(10^{0.1806} = 1.52\times\) and labC \(10^{0.7225} = 5.28\times\) higher than labA. (Exactly, \(\hat\delta_B = 0.6\log_{10}2\) and \(\hat\delta_C = 2.4\log_{10}2\).) labB ends up between the 2× of its direct edge and the 1× implied by the other two, closer to the latter because the labB–labC edge weighs twice as much.

What is left over on each edge:

A residual of a third of a dilution is large for real data; here it was put in on purpose. phenocal writes it at the end of every run, as a check on how well the laboratories agree.

S6Uncertainty

The edge values are uncertain, with standard errors SE. phenocal propagates this by a parametric bootstrap [4]: it draws every edge from \(\mathcal N(\Delta, \mathrm{SE}^2)\), solves the least squares again, repeats this 10,000 times, and reports the median of the offsets and their 2.5th–97.5th percentiles.

Because the solution is a linear function of the edge values with weights \(1/\mathrm{SE}^2\), the same answer can be written down exactly (Proposition 2 of the main text): the offsets are normal, with covariance \((A^{\top} W A)^{-1}\), the matrix inverted above. So

phenocal's bootstrap gives:

The differences from the exact values are Monte Carlo error: 10,000 random draws give each number with some noise, which shrinks with more draws (--bootstrap). With the default seed the upper bound for labC is 12.07 instead of 12.46, a rather unlucky draw: over the seeds 1 to 30 the same bound averages 12.45, with a standard deviation of 0.0042 in \(\log_{10}\) (1%). The intervals are wide because each edge rests on three to six pairs; with real collections, edges have tens or hundreds of pairs.

Why a bootstrap at all, if the exact answer is available? Because phenocal also offers a robust bootstrap with Student-\(t\) draws (--robust) and one that resamples clonal clusters (--clusters), which have no closed form.

S7Harmonised values and probabilistic labels

Each value is brought to the scale of the anchor by removing its laboratory's offset, using the reported offsets:

Whether an isolate is above a threshold \(t\) (say 1.25 mg/L) depends on the offset, which is uncertain. phenocal treats the offset of laboratory \(c\) as normal, centred on \(\hat\delta_c\) with standard deviation \(s_c\) read off its 95% interval, \(s_c = (\log_{10}\mathrm{hi} - \log_{10}\mathrm{lo})/(2 \times 1.96)\), and gives the probability

For S11_B and \(t = 1.25\): \(s_{\mathrm{B}} = (0.5514 - (-0.1949))/3.92 = 0.1904\),

phenocal estimates the same probability by drawing the offset 10,000 times and gets 0.562. The isolate's harmonised value, 1.33, is just above 1.25, but the offset of labB is so uncertain that the label is close to a coin toss. That is the point of probabilistic labels: they separate isolates whose class is clear from those whose class depends on the calibration.

All isolates (labels.tsv; isolates of the anchor have a certain offset, so their probabilities are 0 or 1):

S8Run it yourself

The two files are in the repository and in the Docker image:

phenocal --phenotypes examples/worked_phenotypes.tsv \
         --profiles examples/worked_profiles.tsv \
         --anchor labA --min-support 3 --out worked \
         --trace worked.md --dashboard worked.html

worked.md repeats these calculations for your own data, with the real numbers of the run; worked.html shows them interactively.

References

1. Silva M, Machado MP, Silva DN, Rossi M, Moran-Gilad J, Santos S, Ramirez M, Carriço JA. chewBBACA: a complete suite for gene-by-gene schema creation and strain identification. Microbial Genomics. 2018;4(3):e000166. doi:10.1099/mgen.0.000166.
2. Seemann T. cgmlst-dists: calculate distance matrix from cgMLST allele call tables. Software. github.com/tseemann/cgmlst-dists.
3. Mixão V, Pinto M, Sobral D, Di Pasquale A, Gomes J, Borges V. ReporTree: a surveillance-oriented tool to strengthen the linkage between pathogen genetic clusters and epidemiological data. Genome Medicine. 2023;15(1):43. doi:10.1186/s13073-023-01196-1.
4. Efron B. Bootstrap methods: another look at the jackknife. The Annals of Statistics. 1979;7(1):1–26. doi:10.1214/aos/1176344552.