Validation

The milestones sign 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.

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.

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

Three of the 243 drug calls differ, all of them at the 2% threshold, and all three are the same missed mutation counted once per drug. They are listed with their 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”.

It was nine until codons were read from the reads rather than from the combinations of their bases. Six were a phasing artefact, and what fixed them is below.

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.

The caller credits the V75A base in 4 of the 239 reads it counted at HXB2 2773 — 1.7%, below the 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 disagreement left that would change what anyone did.

Phasing does not reach it, and must not. The raw read bases at 2773 carry the alanine base in 3.2% of themselves, where bcftools mpileup credits 1.7% once base quality is taken into account. Counting raw bases would “find” V75A at a 2% threshold and would agree with the published report by accident, on evidence the caller had already judged too poor to report. Which alleles are present stays the caller’s judgement; phasing only decides which combinations of called alleles are real (D31).

What this exercise changed about the pipeline

Validating is what found the biggest problem, and it was not in the interpreter. It also found the second, which was.

Codons are read from the reads, not from the combinations

Sample 3’s five NRTI rows used to read K219KEGR where Stanford reads K219KE, and Sample 4’s lenacapavir row A105AT where Stanford reads A105T: six discrepancies, every score and level identical, only the mutation string wrong.

A codon was assembled from the VCF — the alleles credited at each of the three bases, then every combination of them. The caller credits G at HXB2 3204 in 43% of the reads it counted and at 3205 in 2.4%, so at a 2% threshold both pass and the combinations give AAA, AGA, GAA, GGA — K, R, E and G. Of the 4,347 reads spanning the whole codon, 54% carry AAA and 45% carry GAA; AGA and GGA are under two tenths of a per cent each. Fewer than two reads in a thousand carried both substitutions. A VCF says what is at each position, not what travels together.

So codons are now read from the alignment, which was already an input for coverage: reads spanning all three bases are counted as the triplets they carry, and the codon holds the amino acids whose share of those reads reaches the threshold. Six of the nine discrepancies went with it, and the agreement at 15% and 5% is unchanged, because a majority allele is carried by the reads that support it.

A consensus caller could not see the minority virus (M6)

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

The milestones also name 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.