Validation

docs/milestones.md signs M6 off when "resistance calls match published results ... with discrepancies explained". This is what was compared, what agreed, and what did not.

Everything here is reproduced by tests/test_validation.py, which runs the real pipeline over real reads on every CI run of the demo job.

Two things are validated, separately

The scoring engine is checked against Stanford's own answers. crates/ligase-core/tests/fixtures/hivdb/ holds the Stanford HIVdb reports published as supplementary data by the HIVGenoPipe study — thirty clinical samples at four ambiguity thresholds, 3,348 expectations — produced by Sierra 3.5.3 against HIVDB 9.8.

All 3,348 are reproduced exactly: every drug, every score, every level.

That oracle matters because the obvious alternative is circular. Checking the scorer against expectations derived from the same XML it reads would test the parser against itself. Sierra is a different implementation, in a different language, by the algorithm's own authors.

The pipeline end to end is checked by running it on reads whose published answers are known, which is the rest of this page.

The samples

PRJNA1218827 is the HIVGenoPipe study's MiSeq data. Four runs are fetched, and the SRA metadata maps each to the sample identifier the paper uses — so nothing rests on reading a figure.

Run	Sample	What it is
SRR32217726	Sample 3	a mixture: K219KE, and L100LI below it
SRR32217715	Sample 4	multi-class positive: NRTI, INSTI and capsid
SRR32217716	Sample 30	fully susceptible — a clean negative
SRR32217714	—	the study's no-template control

A mixture, a strong positive, a true negative and a blank. Between them they exercise every judgement the interpreter makes.

What agrees

The pipeline's mixture parameter answers the same question as the paper's Amb2/Amb5/Amb15 reports, so one set of alignments validates all three thresholds.

Threshold	Drug calls compared	Agree
15%	81	81
5%	81	81
2%	81	72

At 15% and 5% the agreement is exact — every score, every level, every mutation string, for all twenty-seven drugs on all three samples.

That includes the cases worth naming:

• Sample 3 at 15%: K219KE, potential low-level resistance to zidovudine and stavudine. A mixture, scored as both amino acids.
• Sample 3 at 5%: adds L100LI — high-level resistance to four NNRTIs. Same reads, same alignment; a 5% minority virus changes the regimen.
• Sample 4: M184V, T97A + E138T + G140S + Q148H giving high-level resistance to all five INSTIs, and M66I + Q67H + N74D + A105T in capsid giving high-level lenacapavir resistance at a score of 160.
• Sample 30: susceptible on all twenty-seven, with no false positives.
• The no-template control: all twenty-seven drugs Not assessed, with no score. Sixty of its thirty-two thousand reads map to HXB2.

What disagrees, and why

Nine of the 243 drug calls differ, all of them at the 2% threshold. Each is listed with its cause in KNOWN_AT_TWO_PERCENT in tests/test_validation.py, and the test fails if one of them changes shape or a new one appears — which is what keeps "discrepancies explained" from decaying into "discrepancies tolerated".

Six are cosmetic: phasing a VCF cannot carry

Sample 3's five NRTI rows read K219KEGR where Stanford reads K219KE, and Sample 4's lenacapavir row reads A105AT where Stanford reads A105T. Every score and every level is identical.

The cause is the real limitation of reading codons from a VCF. At HXB2 3204 the reads are 44.65% G; at 3205 they are 2.75% G. At a 2% threshold both pass, and enumerating the combinations of the two ambiguous bases gives AAA, AGA, GAA, GGA — K, R, E and G. No single read need carry both substitutions. A VCF says what is at each position, not what travels together.

Neither extra amino acid matches a rule, so nothing clinical follows. The fix, if it is wanted, is to count codon haplotypes from reads spanning the codon, which the alignment already makes available.

Three are one missed mutation, and it is a real one

Stanford reports V75VA for Sample 30 at 2%, reading potential low-level resistance to zidovudine, intermediate to stavudine and low-level to didanosine. Ligase reports susceptible to all three.

This alignment supports the V75A base in 3 of 392 reads at HXB2 2773 — 0.77%, below any threshold that would report it. HIVGenoPipe realigns against a sample-specific reference, which is that paper's whole method, and places reads that an alignment to HXB2 does not.

This is a sensitivity limit, not a disagreement about what the reads mean, and it is the honest cost of the demo's approach. It is also the only one of the nine that would change what anyone did.

What this exercise changed about the pipeline

Validating is what found the biggest problem, and it was not in the interpreter.

The demo's call step used to run bcftools call -mv --ploidy 1. On Sample 3, the reads at HXB2 3204 are A 238, G 192 — a 45% K219E — and bcftools call emitted no record for it at all. The majority allele is the reference, so there was no variant to call. The report came back susceptible to all twenty-seven drugs.

It was right to, given that input: it matches Stanford's own majority-consensus report for the same sample exactly, which is itself useful evidence that the interpreter is sound. But Stanford's 5% report for that sample finds L100I as well and reads high-level resistance to four NNRTIs. A pipeline whose answer flips from "fully susceptible" to "resistant to most of the NNRTI class" on a threshold it never exposed is not doing what it says it does.

So the caller is gone, the pileup's own allele depths are the evidence, and min_freq now means what its name says.

PRJNA517147 is not used for this

docs/milestones.md also names PRJNA517147, the Washington DC cohort study (Jair et al., PLoS ONE 2019;14(4):e0214820). Its per-sample resistance calls are published only as colour in a heatmap — Fig 4, reverse transcriptase only, with protease and integrase "not shown" and no supporting tables.

Of the four runs fetched from it, three are blank rows in that figure and the fourth cannot be separated from its neighbour by eye. The mutations the paper quantifies in text are attached to PL- identifiers that are never mapped to the specimen numbers the figure uses.

That project remains the demo's default dataset, and M3, M4 and M5 are still signed off on it. It is not used as a source of expected resistance answers, because a number transcribed from a figure is not a citable value. PRJNA1218827 publishes machine-readable per-sample reports with the algorithm version stamped in every row, and that is what M6 is checked against.