HIV-1 drug resistance

This page is the project’s own validation of its own pipeline. It is the acceptance test for v0.1’s last milestone, not an independent assessment, and should be read as that. 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.

What was compared

Two things were validated, separately, and the distinction matters.

The scoring engine was checked against Stanford’s own answers: thirty clinical samples at four ambiguity thresholds, 3,348 expectations, produced by Sierra 3.5.3 against HIVDB 9.8 and published as supplementary data by the HIVGenoPipe study. All 3,348 are reproduced exactly — every drug, every score, every level.

That oracle was chosen because the obvious alternative is circular. Checking the scorer against expectations derived from the same XML rule file 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 was checked by running it on reads whose published answers are known — four runs from PRJNA1218827, the HIVGenoPipe study’s MiSeq data, with SRA metadata mapping each run to the sample identifier the paper uses, so that 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

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 gaining L100LI and high-level resistance to four NNRTIs when the threshold drops from 15% to 5%, on the same reads and the same alignment; Sample 4’s T97A + E138T + G140S + Q148H giving high-level resistance to all five INSTIs; Sample 30 coming back susceptible on all twenty-seven with no false positives; and the no-template control returning Not assessed on all twenty-seven with no score.

What disagrees

Three of the 243 drug calls differ, all at the 2% threshold, and all three are the same missed mutation counted once per drug. Each is listed with its cause in KNOWN_AT_TWO_PERCENT in the test, 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. Six were cosmetic: Ligase read K219KEGR where Stanford reads K219KE, and A105AT where Stanford reads A105T, with every score and every level identical. The cause was reading codons from a VCF, which says what is at each position and not what travels together, so enumerating two ambiguous bases yielded amino acids no single read carried. Codons are now read from the reads spanning them, and those six went with the change (D31).

The three that remain are one missed mutation, and it is a real one. Stanford reports V75VA for Sample 30 at 2%. The caller credits that base in 4 of the 239 reads it counted — 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.

That 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.

What validating it changed

This is the part worth reading, and it is not a success story.

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, because the majority allele is the reference and there was therefore 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.

Reading further

The full comparison, with every discrepancy and its cause, is in validation. How the HIVdb rules are read and scored is in resistance interpretation. How to write a resistance step and choose a threshold is in interpreting resistance.