Abstract
Phenotypes such as minimum inhibitory concentrations are routinely pooled from studies that measured them with different protocols, and a protocol that shifts the whole measurement scale confounds the phenotype with its source before any analysis begins. phenocal estimates these shifts without shared reference strains. It treats pairs of genetically identical or near-identical isolates measured in different laboratories as natural controls, in which the biology cancels and only the difference between the two laboratory offsets remains. For every pair of laboratories an adaptive bound on the allelic distance selects the closest pairs whose median log-ratio is stable, the resulting edge values are combined by weighted least squares on the graph of laboratories, and the uncertainty of the offsets is propagated to harmonised values and to isolate-level probabilities of exceeding each threshold. The estimator uses only the phenotypes, the cohort labels and the allelic profiles, never the features of a downstream model, and is described in the AIS 2026 paper [1]; this page gives the method in full and the supplementary material works every calculation by hand.
Keywords: label calibration · batch effect · natural controls · cgMLST · wgMLST · weighted least squares · probabilistic labels
1Introduction
A phenotype such as a minimum inhibitory concentration (MIC) is often assembled from several studies or laboratories, each with its own protocol. Differences in medium, inoculum or reading do not add noise so much as move the whole scale, so that one laboratory may read every strain two dilutions higher than another. When the raw values are pooled, biology and protocol become inseparable: a classifier trained on the pooled data can learn which laboratory produced a value rather than the trait itself, and a threshold such as “tolerant if MIC ≥ 1.25 mg/L” takes a different meaning in each cohort (Figure 1).
examples/toy_phenotypes.tsv and of the
phenocal output.The scale is normally fixed with shared reference material, the same control strains measured by every laboratory. Data collected from published studies rarely include such controls, but they often contain something equivalent. Surveillance collections sample the same clones repeatedly, so the same genome is frequently measured by more than one laboratory, and such a pair behaves as a control whose true phenotype is unknown but identical on both sides (Figure 2). phenocal reads these natural controls from the genomes, through cgMLST or wgMLST allelic distances, and estimates one offset per laboratory together with its uncertainty.
2Model and estimator
The procedure has five steps, summarised in Figure 3: collecting the cross-cohort pairs with their allelic distance, choosing for each pair of cohorts the pairs close enough to act as controls, summarising them in an edge value with a standard error, solving for the offsets on the graph of cohorts, and propagating the uncertainty to the labels.
2.1Measurement model
Let \(y_i = \log_{10} m_i\) be the logarithm of the value \(m_i\) measured on isolate \(i\), and \(c(i)\) the cohort that measured it. We assume
$$ y_i = \mu_i + \delta_{c(i)} + \varepsilon_i, \tag{1} $$where \(\mu_i\) is the unknown biological value of the isolate, \(\delta_c\) the offset of cohort \(c\), and \(\varepsilon_i\) measurement noise with median zero. An offset on the logarithmic scale is a constant factor \(10^{\delta_c}\) on the original scale, which is how a shift of the dilution series or of the reading behaves.
Definition 1 (log-ratio). For isolates \(i\) in cohort \(a\) and \(j\) in cohort \(b \neq a\), the log-ratio is $$ s_{ij} = y_i - y_j = (\delta_a - \delta_b) + (\mu_i - \mu_j) + (\varepsilon_i - \varepsilon_j). \tag{2} $$
If the two isolates have the same genome, \(\mu_i = \mu_j\) and the log-ratio measures the contrast \(\delta_a - \delta_b\) up to noise. Pairs within one cohort carry no information on the offsets and are not used.
2.2Genomic distance
Relatedness is measured on cgMLST or wgMLST allelic profiles, which phenocal reads as
written by chewBBACA [2] (--profiles).
Definition 2 (allelic distance). Let \(a_{il}\) be the allele of isolate \(i\) at locus \(l\), or \(\varnothing\) if the call is missing. The distance between isolates \(i\) and \(j\) is the number of loci called in both whose alleles differ, $$ d_{ij} = \sum_{l} \mathbf 1\big[a_{il} \neq \varnothing,\; a_{jl} \neq \varnothing,\; a_{il} \neq a_{jl}\big]. \tag{3} $$
This is the distance computed by chewBBACA, cgmlst-dists [3] and
ReporTree [4], so a missing call never counts as a difference. The chewBBACA
prefix INF- marks a newly inferred allele and is removed; the chewBBACA
classes LNF, ASM, ALM, NIPH,
PLOT3, … and the value 0 are missing calls; integer alleles are compared by
value, so that 32-bit CRC hashes are handled at any size, and any other text is compared
as text. Distances computed elsewhere, a distance matrix or a table of pairs on any other
ruler such as SNP distances, can be supplied instead (--pairs).
2.3Adaptive relatedness bound
Pairs at distance zero are the cleanest controls, but there may be too few of them, while pairs a few alleles apart are almost as informative and distant pairs differ in their biology. The bound is therefore chosen from the data, separately for every pair of cohorts.
Definition 3 (adaptive bound). Let \(P_{ab}\) be the cross-cohort pairs of cohorts \(a\) and \(b\), \(d_{\min}\) their smallest distance, \(m_0\) the median log-ratio of the pairs at \(d_{\min}\), and \(\bar\sigma_{ab} = \sqrt{(\sigma_a^2 + \sigma_b^2)/2}\) the resolution of the two cohorts. The bound \(\tau_{ab}\) is the smallest distance \(\tau\) observed in \(P_{ab}\) such that $$ n(\tau) = \#\{(i,j) \in P_{ab} : d_{ij} \le \tau\} \ge n_{\min} \quad\text{and}\quad \big|\operatorname{median}\{s_{ij} : d_{ij} \le \tau\} - m_0\big| \le \kappa\,\bar\sigma_{ab}. \tag{4} $$ If no distance satisfies both conditions, all pairs of \(P_{ab}\) are used.
The first condition requires enough pairs (--min-support, \(n_{\min}=5\)
by default); the second stops the bound as soon as more distant pairs move the median by
more than \(\kappa\) units of measurement resolution (--max-drift,
\(\kappa = 1\)), the signature of biology rather than protocol (Figure 4). Taking the
smallest admissible bound keeps the evidence as close to identical genomes as the data
allow, and because the bound is chosen per edge there is no global distance threshold to
tune, which matters since cgMLST and wgMLST distances live on different scales. When no
bound is admissible, the large \(\tau\) and small \(n\) give the edge a low weight through
Eq. (6).
2.4Edge estimates
Each pair of cohorts with supporting pairs becomes an edge of the comparison graph, with value
$$ \Delta_{ab} = \operatorname{median}\{\, s_{ij} : (i,j) \in P_{ab},\ d_{ij} \le \tau_{ab} \,\} \;\approx\; \delta_a - \delta_b, \tag{5} $$a median rather than a mean, so that a few discordant pairs do not move it and only the symmetry of the residuals is required, not a parametric error law. Its standard error combines the resolution of the two cohorts, a term for the biology that may still differ between isolates up to distance \(\tau\), and the number of pairs \(n_{ab} = n(\tau_{ab})\):
$$ \mathrm{SE}_{ab} = \sqrt{\frac{\sigma_a^2 + \sigma_b^2 + \lambda\,\tau_{ab}}{n_{ab}}}, \qquad w_{ab} = \frac{1}{\mathrm{SE}_{ab}^2}. \tag{6} $$By default \(\sigma_c = \log_{10} 2 / \sqrt 2\), the resolution of a two-fold dilution
series, and \(\lambda = (\log_{10} 2)^2\) per allele; --sigma empirical estimates
each \(\sigma_c\) from the spacing of the cohort's own measurement grid. With the defaults
Eq. (6) reduces to \(\log_{10} 2 \cdot \sqrt{(1 + \tau)/n}\).
2.5Calibration on the comparison graph
The cohorts are the nodes of a graph whose edges carry the values \(\Delta_{ab}\) and
weights \(w_{ab}\) (Figure 5). Since only contrasts are identified, one cohort, the anchor,
is fixed at zero (--anchor, the largest cohort by default), and the offsets
solve the weighted least-squares problem
whose normal equations
$$ A^{\top} W A\, \hat\delta = A^{\top} W \Delta \tag{8} $$involve the signed edge–node incidence matrix \(A\) without the anchor column and \(W = \operatorname{diag}(w_{ab})\); \(A^{\top} W A\) is the weighted graph Laplacian with the anchor row and column removed. The system has a unique solution if and only if every cohort is connected to the anchor through the graph, possibly through other cohorts; phenocal stops with an error otherwise.
Proposition 1 (anchor invariance; Lemma 1 of [1]). Let \(\hat\delta\) solve Eq. (7) with anchor \(a\) and \(\hat\delta'\) with anchor \(b\). Then \(\hat\delta' = \hat\delta - \hat\delta_b \mathbf 1\). Consequently every contrast \(\hat\delta_a - \hat\delta_b\), and every label obtained by comparing a harmonised value with a threshold expressed on the same anchored scale, does not depend on the choice of the anchor.
The objective depends on \(\delta\) only through differences, hence is invariant under \(\delta \mapsto \delta + k\mathbf 1\). For a connected graph the weighted Laplacian has a one-dimensional null space spanned by \(\mathbf 1\), so each anchor selects the unique minimiser on that line, and two anchors give solutions that differ by a constant. Re-anchoring multiplies harmonised values and thresholds by the same factor, so no label changes.
2.6Uncertainty
phenocal propagates the uncertainty of the edges by a parametric bootstrap
[5]: each edge value is drawn from \(\mathcal N(\Delta_{ab}, \mathrm{SE}_{ab}^2)\),
Eq. (7) is solved again, and this is repeated \(B\) times (--bootstrap,
\(B = 10\,000\)). The reported offset is the median of the draws, and the 95% interval
spans their 2.5th and 97.5th percentiles. For this Gaussian scheme the result is also
available in closed form.
Proposition 2 (exact distribution). If the edge values are perturbed independently as \(\Delta^{*}_{ab} \sim \mathcal N(\Delta_{ab}, \mathrm{SE}_{ab}^2)\), the bootstrap offsets are exactly normal, \(\hat\delta^{*} \sim \mathcal N\big(\hat\delta,\ (A^{\top} W A)^{-1}\big)\). The bootstrap median therefore estimates the least-squares solution \(\hat\delta\), and the interval estimates \(\hat\delta_c \pm 1.96\,[(A^{\top} W A)^{-1}]_{cc}^{1/2}\), up to Monte Carlo error.
By Eq. (8), \(\hat\delta^{*} = (A^{\top} W A)^{-1} A^{\top} W \Delta^{*}\) is linear in \(\Delta^{*}\), whose covariance is \(W^{-1}\). Hence \(\hat\delta^{*}\) is normal with mean \(\hat\delta\) and covariance \((A^{\top} W A)^{-1} A^{\top} W\, W^{-1}\, W A\,(A^{\top} W A)^{-1} = (A^{\top} W A)^{-1}\).
The bootstrap is kept because two variants have no closed form: --robust
draws the perturbations from a Student-\(t\) distribution, and --clusters
resamples whole clonal clusters. The latter matters because near-identical pairs are
often nested in clusters, and pairs from the same cluster are not independent; given a
table assigning isolates to clusters, phenocal reports a cluster-weighted estimate and an
interval that resamples clusters across the graph, the preferred summary when cluster
identifiers are available. Proposition 2 is verified numerically in Section 6 of the
supplementary material.
2.7Harmonised values and probabilistic labels
The harmonised value of isolate \(i\) is its value on the anchor's scale,
\(m_i / 10^{\hat\delta_{c(i)}}\). Because \(\hat\delta_c\) is uncertain, so is the class of an
isolate close to a threshold. For each threshold \(t\) (--thresholds) phenocal
reports
where \(s_c\) is the half-width of the 95% interval of \(\delta_c\) divided by 1.96 and \(\Phi\) the standard normal distribution function; phenocal evaluates Eq. (9) by Monte Carlo. Isolates of the anchor cohort have no offset uncertainty and receive probabilities 0 or 1 (Figure 6).
labels.tsv is 0.256).A downstream model can use these probabilities in place of a hard label that might flip under a small change of the offset, or leave out the isolates whose label is uncertain. Finally, phenocal estimates a data-driven cutoff on the harmonised scale by two methods that bracket it, the antimode of a kernel density estimate [6] and the crossover of a two-component Gaussian mixture, with a bootstrap interval.
3Calibration before supervised learning
phenocal is a pre-processing step for supervised learning on pooled cohorts: it corrects the target, the phenotype, before any model is trained. Its inputs are the phenotypes, the cohort labels and the genomes, never the features or the model, so it can precede any classifier or regressor without leaking information from the evaluation (Figure 7).
3.1Classification
A classifier needs a label, usually obtained by cutting the phenotype at a threshold. With raw values, a cohort whose protocol reads high falls entirely above the threshold and one that reads low falls mostly below it, and two problems follow. First, since cohorts usually sample different populations, the genome reveals the cohort, and a classifier trained on pooled raw labels can appear accurate by recognising the cohort rather than the trait, which a random cross-validation does not reveal. Second, the question of interest, whether the model transfers to cohorts it was not trained on, may become impossible to answer: if the held-out cohorts contain a single class, balanced accuracy and ROC-AUC are undefined.
On the restored endpoint the paper reports, for a gradient-boosted tree classifier on allelic profiles, a forward ROC-AUC of 0.999 with wgMLST features and 0.748 with cgMLST features, and 0.880 in the reverse direction. The endpoint remains evaluable across the pre-declared thresholds 0.6–2.0 mg/L, with the AUC changing only mildly with the cutoff. The result does not arise from leakage: it falls to chance with permuted training labels and with a cohort-only baseline (AUC 0.50), and it is unchanged when the isolates used as calibration controls are removed from training. These figures show that the endpoint has become evaluable; they do not show that the calibrated MICs are the biological truth.
The probabilities of Eq. (9) offer two further options beyond the hard label: they can serve as soft targets or sample weights, or training can be restricted to isolates whose class is stable (for instance \(P_i(t) < 0.1\) or \(> 0.9\)).
3.2Regression
When the phenotype itself is predicted, the target is \(y_i = \mu_i + \delta_{c(i)} + \varepsilon_i\). Trained on raw values, a regressor must account for the offsets as well, and because cohort and genome are correlated it attributes them to the genome; its predictions then mix biology with protocol, refer to no particular laboratory's scale, and carry the offset of whichever cohort a new isolate resembles. Trained on the harmonised values \(y_i - \hat\delta_{c(i)}\), it sees a single scale, the anchor's, and only biology and noise remain to be explained. The calibration uncertainty can be carried along, since \(s_c^2\) adds to the noise of every target in cohort \(c\) and weights such as \(1/(\sigma_c^2 + s_c^2)\) let well-calibrated cohorts count more. The paper evaluates the classification endpoint; the argument for regression is the same, but no regression results are reported.
3.3Relation to batch-effect correction and label noise
Batch-effect methods such as ComBat [12] and the linear-model correction of limma [13] address shifts in the features, which are well documented in high-throughput data [14]; covariate-shift corrections likewise reweight the feature distribution [15]. Here the shift is in the label. Applied to the phenotype as a per-cohort centring, a batch correction removes every difference between cohort means, including the biological ones when cohorts sample different populations, whereas phenocal estimates the offset only from isolates whose biology is held fixed by the genome. Nor is the problem one of random label noise [16]: the error is systematic and shared by a whole cohort, so model capacity does not remove it, and in the paper's simulations larger neural models trained on uncalibrated labels did no better than logistic regression, while calibration improved every learner tested.
4Illustrative example
The toy example of Figure 1 is included in the repository and runs in a fraction of a second:
phenocal --phenotypes examples/toy_phenotypes.tsv \
--profiles examples/toy_profiles.tsv \
--anchor labA --min-support 2 --out toy
Three strains are each measured by three laboratories, so every pair of cohorts has three pairs at distance zero. For labA and labB the three log-ratios all equal \(\log_{10} 0.5\), so \(\Delta_{AB} = -0.301\) with \(\mathrm{SE} = \sqrt{(0.213^2 + 0.213^2)/3} = 0.174\). The three edges are mutually consistent (\(2 \times 2.5 = 5\)), and Eq. (7) gives \(\delta_B = \log_{10} 2\) and \(\delta_C = \log_{10} 5\). phenocal reports labB 1.99× [1.05, 3.78] and labC 5.01× [2.64, 9.26] relative to labA, centred on the true factors and wide because each edge rests on only three pairs. After calibration each strain has the same value in every laboratory to within 1% (Table 1).
Table 1. Raw and harmonised values of the toy example
(harmonised.tsv).
| sample | cohort | raw value | harmonised |
|---|---|---|---|
| s1_A | labA | 0.5 | 0.500 |
| s1_B | labB | 1 | 0.504 |
| s1_C | labC | 2.5 | 0.499 |
| s2_A | labA | 1 | 1.000 |
| s2_B | labB | 2 | 1.007 |
| s2_C | labC | 5 | 0.999 |
| s3_A | labA | 2 | 2.000 |
| s3_B | labB | 4 | 2.015 |
| s3_C | labC | 10 | 1.998 |
The same results are obtained from the distance matrix
examples/toy_distance_matrix.tsv and from the pair table
examples/toy_pairs.tsv given with --pairs; the test suite checks
that the three inputs give identical outputs. Because the toy is noiseless, it exercises
neither the adaptive bound nor the least-squares compromise. The
supplementary material therefore follows a small example with
noise through every step by hand, from the allelic profiles to the probabilistic labels,
and checks each number against the output of phenocal.
5Assumptions and limitations
The method rests on three assumptions. Isolates at small allelic distance are assumed to share their true phenotype; the adaptive bound and the \(\lambda\tau\) term limit the influence of the biology that remains between near-identical isolates, and pairs at distance zero are the cleanest evidence. Each laboratory is assumed to scale all its values by one factor, a constant shift on the logarithmic scale; effects that depend on the value itself are not modelled. Finally, the comparison graph must be connected through genetically linked isolates, and cohorts that share no clone with any other cannot be calibrated.
Remark (identifiability). A biological difference that affects every isolate of one cohort in the same way, including its clones shared with other cohorts, is indistinguishable from a protocol offset and would be calibrated away. Biology specific to some strains does not have this problem, since it makes the log-ratios grow with genomic distance, which is what the adaptive bound detects.
6Software availability
phenocal is written in Rust without external dependencies and released under the MIT
licence at github.com/genpat-it/phenocal.
It is distributed as static binaries for Linux, macOS and Windows (installable with a
one-line script), as the Docker image ghcr.io/genpat-it/phenocal, and on
crates.io (cargo install phenocal). Each run writes the offsets, the per-cohort
resolutions, the edges, the harmonised values, the probabilistic labels and the cutoff as
tab-separated tables, and optionally an interactive report (--dashboard) and a
step-by-step trace of the calculation with the numbers of the run (--trace).
The code as used in the paper is frozen at the tag
ais2026; later
releases add input formats and packaging and give identical results on the same input.
How to cite
If you use phenocal, please cite [1]:
@inproceedings{deruvo2026phenocal,
title = {{PhenoCal}: Natural-Control Label Calibration for
Cross-Cohort Machine Learning},
author = {de Ruvo, Andrea and Castelli, Pierluigi and
Di Pasquale, Adriano and Radomski, Nicolas},
booktitle = {12th International Artificial Intelligence Symposium (AIS 2026)},
series = {Lecture Notes in Computer Science},
publisher = {Springer Nature},
year = {2026}
}
References
- de Ruvo A, Castelli P, Di Pasquale A, Radomski N. PhenoCal: natural-control label calibration for cross-cohort machine learning. In: 12th International Artificial Intelligence Symposium (AIS 2026). Lecture Notes in Computer Science. Springer Nature; 2026. Reproducibility archive: doi:10.5281/zenodo.21461769.
- 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.
- Seemann T. cgmlst-dists: calculate distance matrix from cgMLST allele call tables. Software. github.com/tseemann/cgmlst-dists.
- 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.
- Efron B. Bootstrap methods: another look at the jackknife. The Annals of Statistics. 1979;7(1):1–26. doi:10.1214/aos/1176344552.
- Silverman BW. Density Estimation for Statistics and Data Analysis. London: Chapman & Hall; 1986.
- Palma F, Radomski N, Guérin A, et al. Genomic elements located in the accessory repertoire drive the adaptation to biocides in Listeria monocytogenes strains from different ecological niches. Food Microbiology. 2022;106:103757. doi:10.1016/j.fm.2021.103757.
- Kragh ML, Ivanova M, Dalgaard P, Leekitcharoenphon P, Truelstrup Hansen L. Sensitivity of 240 Listeria monocytogenes isolates to common industrial biocides is more dependent on the presence of residual organic matter or biofilm than on genetic determinants. Food Control. 2024;158:110244. doi:10.1016/j.foodcont.2023.110244.
- Ivanova M, Kragh ML, Szarvas J, et al. Large-scale phenotypic and genomic analysis of Listeria monocytogenes reveals diversity in the sensitivity to quaternary ammonium compounds but not to peracetic acid. Applied and Environmental Microbiology. 2025;91:e01829-24. doi:10.1128/aem.01829-24.
- Cooper AL, Carrillo CD, Deschênes M, Blais BW. Genomic markers for quaternary ammonium compound resistance as a persistence indicator for Listeria monocytogenes contamination in food manufacturing environments. Journal of Food Protection. 2021;84(3):389–398. doi:10.4315/JFP-20-328.
- He Y, Xu T, Li S, Mann DA, Britton B, Oliver HF, den Bakker HC, Deng X. Integrative assessment of reduced Listeria monocytogenes susceptibility to benzalkonium chloride in produce processing environments. Applied and Environmental Microbiology. 2022;88(21):e01269-22. doi:10.1128/aem.01269-22.
- Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007;8(1):118–127. doi:10.1093/biostatistics/kxj037.
- Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research. 2015;43(7):e47. doi:10.1093/nar/gkv007.
- Leek JT, Scharpf RB, Bravo HC, et al. Tackling the widespread and critical impact of batch effects in high-throughput data. Nature Reviews Genetics. 2010;11(10):733–739. doi:10.1038/nrg2825.
- Shimodaira H. Improving predictive inference under covariate shift by weighting the log-likelihood function. Journal of Statistical Planning and Inference. 2000;90(2):227–244. doi:10.1016/S0378-3758(00)00115-4.
- Frénay B, Verleysen M. Classification in the presence of label noise: a survey. IEEE Transactions on Neural Networks and Learning Systems. 2014;25(5):845–869. doi:10.1109/TNNLS.2013.2292894.