Resistance interpretation, internally

Three layers in crates/ligase-core/src/hivdb/, none of which knows about Python:

The ruleset is derived, not parsed as XML

Stanford’s algorithm is ASI2 XML, about 190 KB of it. Two files are committed per algorithm: that XML byte-for-byte, pinned by SHA-256 in rulesets/PROVENANCE.tsv, and a .ruleset derived from it by scripts/vendor_hivdb_algorithm.py.

The engine reads only the .ruleset, which means no XML parser in the core and no new dependency. The XML stays as the auditable thing you diff against upstream. tests/test_rulesets.py re-derives every .ruleset from its XML on every test run, so the two cannot drift apart unnoticed — a derived file that nothing re-derives is a file that can quietly stop matching its source.

The rulesets are then embedded with include_str!, so an installed wheel carries its own rules and a run cannot depend on where a file was unpacked. A test ties the embedded copy to the committed one.

The grammar is small enough to be strict about

Every HIVDB release from 9.8 to 10.2 uses one shape:

SCORE FROM ( <item>, <item>, ... )

where an item is <residue> => <score>, several residues joined by AND, or MAX( several such ). There is no ATLEAST, NOTMORETHAN, SELECT or EXCLUDE; all 27 drugs use USE_GLOBALRANGE; and the global range is the same in both versions:

(-INF TO 9 => 1, 10 TO 14 => 2, 15 TO 29 => 3, 30 TO 59 => 4, 60 TO INF => 5)

So both the derivation and the reader refuse anything they do not recognise, rather than skipping it. That is D21’s rule applied here: a reader that skipped what it did not understand would silently drop a drug, and a dropped drug reads as susceptible.

The derived format

#ligase-hivdb-ruleset 1
algorithm   HIVDB   9.8   2025-01-05
level       1       S     Susceptible
range       *       9     1
gene        RT      NRTI,NNRTI
class       NRTI    ABC,AZT,D4T,DDI,FTC,3TC,TDF
drug        ABC     abacavir        NRTI
rule        ABC     2       45      65R
rule        ABC     21      10      65RN+151M

A rule row is drug, group, score, residues. Rows sharing a group are MAX alternatives; residues joined by + are the rule language’s AND. A drug’s score is the sum over groups of the largest matching item in each, which is exactly what SCORE FROM with MAX means.

Two details that matter:

• A group with nothing matching contributes nothing, and the running maximum starts at None rather than at zero. Scores can be negative — 184V => -10 for zidovudine — and seeding at zero would silently turn a −10 into a 0.
• Nothing is clamped. Stanford’s own reports carry AZT -10.0 at level 1, and so does Ligase’s.

Amino acids are a bitset

Acids is a u32: bits 0–25 for A–Z, bit 26 for a deletion, bit 27 for an insertion, bit 28 for a stop codon. Matching a rule is then an intersection test, done once per rule position per drug per sample.

A stop codon can be observed but is never named by a rule, so it can never satisfy one — which is the right answer rather than a gap. A stop in the middle of protease is a defective virus, not a resistant one.

Coverage, and the three states of a position

Protein holds only the positions that were read. A position read and agreeing with the reference is in it, holding the reference amino acid; a position with no coverage is absent. That absence is the whole mechanism:

An item is undetermined only when every residue that was read matches and at least one was not. An item already refuted by a read position is settled, so its unread positions are not reported — otherwise every drug would list most of its rule positions and the warning would mean nothing.

A gene absent from Observed altogether was never supplied, and its drugs come back Not assessed with no score.

The frame the windows are stated in

hivdb::frame holds the one thing codon::GENES cannot say for itself: that its numbers are positions in HXB2 K03455 counted from the first base of the genome. Three constants and a verdict —

pub const FRAME: &str = "K03455";
pub const BASES: u64 = 9719;
pub const POL: u64 = 2085;
pub const LAST: u64 = /* the furthest base any window reaches: 5093 */;
 
pub fn mismatch(reference: &ReferenceId) -> Option<Mismatch>;

LAST is folded out of GENES in a const block rather than written down, so a window that moves takes the guard with it.

Mismatch has five shapes, and each owns its own sentence and its own instruction, for the reason Note does: the wording belongs to the kind of problem rather than to whichever surface prints it.

Checked in two places, for two different reasons

At check time, in graph.rs, keyed on the step’s output kind. A step annotated -> ResistanceReport has said in its signature what it reads, so the graph compares its Reference bindings against frame::expected() and reports a ReferenceFrame error with the rest of the pass. This is the one check that compares an input against the engine rather than against another input, which is why it does not live in check_step’s pairwise table: dimension checking needs two inputs that can both speak to a dimension, and only a reference can speak to a frame.

At read time, in codon::read, but only for TooShort. A declaration is taken at its word because there is nothing in a FASTA to check it against, so a pol-relative reference declared as a whole genome passes the check; the reading can still see that the windows reach 5093 and the file does not. It returns LoadError::Frame, which write_report keeps apart from every other failure as ReportError::Frame, which the bindings raise as _core.FrameError, which lg.hivdb.interpret turns into the same ReferenceFrameError a check would have raised. Four layers, one class, because what a caller branches on should not depend on which layer caught it.

This replaced a Note::GenePastReference, which said per gene that a gene ran off the end of the reference. It was not enough: on a pol-relative reference the genes that do fit are read out of the wrong genes, and the note said nothing about those. A reference the windows cannot be read against is now refused rather than annotated, which makes that note unreachable, and it is gone.

Codons

codon::GENES holds the HXB2 landmarks:

Two of these were wrong in the demo before M6, and both errors were silent. Reverse transcriptase ran to 440 codons — the p51 subunit — where HIVdb numbers RT 1–560 and its published mutation lists reach 558. Integrase ran to 5096, which is where pol ends and includes its stop codon; the protein is 288 amino acids. hivdb_codons.rs now checks the windows against the rules, because a rule about a position outside the window would never fire and would do so quietly.

A gene has to be read before it is assessed

gene_coverage is the share of a gene’s scored codons that must be read for the gene to go into Observed at all, defaulting to 0.8. codon::needed turns it into a count, rounded up, so the comparison is in integers and the note can say 5 of 26, where 21 were needed rather than a percentage the reader has to do arithmetic on.

The denominator is Ruleset::scored(gene) — every position some drug scored against that gene names, restricted to the gene’s own window — and the numerator counts those of them Protein::at returns something for. codon::read therefore takes the &Ruleset, which is why write_report loads the algorithm before it reads the calls. The reading still does no scoring; it needs the positions to count over.

Counting over the gene instead was close to unsatisfiable. RT is 560 codons and HIVDB 9.8 scores 40 of them, the highest at 348, so four fifths of the gene asked for 448 codons — a floor a partial-pol genotype cannot clear however much of what the rules consult it covered (D34). It cuts the other way too: a dropout over RT 40–120 leaves 86% of the gene read and takes 21 of the 40 scored positions, leaving 19 against a floor of 32 — which the old share waved through and this one refuses.

A gene below it is left out, which makes its drugs not assessed — the same state as a gene the pipeline never supplied. The reading emits Note::GeneBarelyRead, whose counts are of scored codons.

That it is one verdict for a whole gene is its real limit, and a consumer running a per-drug policy measured what the limit costs. An RT 40–120 dropout with M184V read — the mutation in covered sequence — keeps four RT rows at gene_coverage=0.0, a high-level 3TC call among them, each marked with its unread positions, while nine susceptible calls are declined; at the 0.8 default the gene is skipped and the 3TC call goes with it. Their policy is asymmetric on purpose: a resistant call rests on a mutation that was read and unread positions could only add to it, while a susceptible call leaning on unread codons is absence of evidence. A share over a gene cannot say that, so gene_coverage=0.0 plus the per-drug unread column is how such a policy is expressed here — a supported configuration rather than a workaround for the old denominator. Whether a per-drug floor should be the shape of the rule is D34’s open question, now with a measurement behind it.

A position a ruleset names past the end of a gene is left out of the count: no reading could cover it, and a floor nothing can clear is not a floor. hivdb_codons.rs checks the windows against the rules for the same reason — a rule about a position outside the window never fires, and did so quietly.

It is enforced here, in the reading, rather than in score, and that placement is the whole reason the rule is a share over a gene rather than the sharper per-drug one. hivdb_scoring.rs builds Observed directly from 3,348 published Stanford expectations whose fixtures name only the positions that differ from consensus B; a coverage rule in the scoring layer would see almost every position as unread and turn the oracle’s own answers indeterminate. codon::read is the layer the oracle does not go through. D34 has the argument and what it would take to do the sharper thing.

Reading a codon:

1. Every base of the codon must be covered to depth, or the codon is left out. Two bases out of three is not two thirds of an amino acid.
2. Each base’s alleles are taken from the VCF’s AD, keeping those at or above mixture. Every combination of them is a candidate.
3. The reads decide which candidates are real. Reads spanning all three bases of the codon are counted as the triplets they carry, and the codon holds the amino acids whose share of those reads reaches mixture.

Step 3 is what D31 added, and step 2 is why it is needed: a VCF says what is at each position, not what travels together, so the candidates for a codon with two ambiguous bases include codons no read carries. All four used to be reported, which was six of M6’s nine published discrepancies.

Three properties of step 3 are deliberate, and D31 argues each of them:

• It only ever removes a candidate. Which alleles are present at a position stays the caller’s quality-aware judgement.
• The threshold is applied to the amino acid, summed over the triplets that spell it, not to each triplet on its own.
• A codon the reads cannot settle — no read spanning it depth times, or no candidate among the reads at all — keeps its per-base reading, and Reading::notes says how many codons of which gene did that.

Reading::notes is a Vec<Note>, and the step hands it back to the engine rather than keeping it: write_report returns the notes, interpret_resistance returns them as the code-and-value pairs lg.write expects, and lg.hivdb.interpret’s body returns those. From there the engine says them wherever it says what a run did (Notes). A reading made without an alignment, a gene nobody sequenced deeply enough and a codon no read spans are each a hedge with a clinical direction that the report file itself cannot state (D32).

What it costs

codon::scan makes one pass over the BAM and answers both questions from the same bases: a dense u32 counter per gene position for coverage, and a dense [u32; 64] of triplet counts per codon, indexed by the triplet’s place in the codon table. Dense because it is written once per codon per read — a few million times on a deep amplicon — and the whole table is a megabyte, which is cheaper to allocate once than to look a key up that often.

Reading sequences rather than spans roughly doubles the step: on the deepest run of the validation, 95,000 reads, resistance goes from 0.6s to 1.1s. It stays the cheapest step in the pipeline — half a second against an alignment that takes tens of them, and less than the demo’s own run-to-run variation on the same code.

How the rules reach the cache key

A callback step’s recipe describes as the constant "callback", and code is the step function’s source, not what a library it calls reads. So the rules had no home in a key at all.

They go in through tool=. That field already means “something whose version is part of the result”, and its version string is already free-form — the whole of <tool> --version. exec::BUILT_IN, one entry long, makes the exception visible rather than hiding it inside the probe, and hivdb resolves to:

HIVdb rules as shipped: 9.8 2025-01-05 blake3:a1139c…; 10.2 2026-04-26 blake3:3268c9…; read by ligase 0.4.0

Four consequences, all deliberate:

• By content, not by name. Upstream publishes corrections under near-identical names, so a key that moved only with the version number would serve a result computed under rules that have since changed.
• Every shipped algorithm, not only the one in use, because the probe runs before any step body has said which it wants. Editing 10.2 reruns a pipeline pinned to 9.8. That is over-conservative in the direction D20 already commits to: a cache that is too cold costs time, one that is too warm is wrong.
• “as shipped”, because a methods paragraph renders a tool version after “ran with”, and a step that named 9.8 must not appear to have run with 10.2. Which one a step used is its algorithm parameter, reported beside it.
• “read by ligase <version>”, because the rules are only half of what a hivdb step is. The other half is this engine’s judgements about reading them — which codons count as read, how much of a gene has to be, what a mixture satisfies — and those move between releases. Moving gene_coverage onto the scored codons changed what the engine answers for the same calls, the same reference, the same parameters and the same rules; without this a cache would have gone on serving the old refusals. Over-conservative in the same direction: a release that touched none of it still reruns every hivdb step.
• The build, not the release, which is a distinction the first version of this line did not draw. CARGO_PKG_VERSION is read from the manifest, and the manifest is bumped by a release commit — so every commit between one release and the next reported the number of the release it followed. A consumer measured the cost on a source build of the commit carrying the scored-codon floor: it reported read by ligase 0.3.0 while answering differently from the published 0.3.0, so a shared store would serve one engine’s answer to the other. A build that is not a release now says so — read by ligase 0.4.0+13.g0c5980e — and the key follows the source for free (D62).

Doing it this way also picks up provenance, the methods sentence and D24’s citation table for free — and HIVdb genuinely needs a citation, which a plain parameter gave no home.

What is not here

Comments. Stanford’s reports carry explanatory prose per mutation and per result. They are in the vendored XML, the derivation drops them, and Ligase reports levels, scores and mutations. They are the obvious addition if a report ever wants prose.

Any network call. Stanford’s Sierra service appears in this repository exactly once: as the source of tests/fixtures/hivdb/expected.tsv, which was built once from published supplementary data and committed. Nothing queries it at run time or in CI.