Show target codons, vaccine mismatches, and recovery by RNA support

README.md

@@ -79,6 +79,29 @@ Benchmark-matched Sid counts and Exacto protein/RNA-only/no-call outcomes stay v
 below the sequences; other source libraries are available in a separate disclosure.
 Missing counts, zero coverage, and covered loci with zero mutant reads remain
 different states. RNA allele support does not establish full peptide recovery.
+The report separates 2+ mutant-read evidence from single-read evidence, observed
+zero, and unknown counts. This uses the maximum of Sid’s count and distinct
+Exacto raw-read allele-call inputs per benchmark library at any timepoint, never
+a sum across libraries or methods. Assembly and corrected counts do not raise
+this tier; per-sample evidence retains both underlying counts. The main
+page retains the overall score and shows recovery within each evidence group;
+amino-acid failures name the expected and observed residues inline.
+
+Blue underlines locate each translated target codon using the archived RNA-call
+position. Deletion markers mean the codon at/after the junction. A dotted tail
+means downstream of a frameshift, not confirmed novel sequence throughout; fusion
+breakpoints and novel-ORF boundaries are not annotated in this panel. Reference
+markers require an exact local anchor with an unambiguous projected position.
+Rose marks vaccine-versus-reconstruction differences on both rows, independently
+of the amber comparison against W1. Missing ends are coverage, not substitutions.
+Stars mark every RNA-support leader (including ties) per sample and method.
+Exacto itself emits unranked candidates; this report ranks displayed windows by
+distinct input RNAs without using vaccine agreement. The top-window filter finds
+covered-position disagreements, including leaders that contain another peptide;
+details distinguish incomplete, unaligned, and allele-unconfirmed sequences.
+The recorded CD109 peptide ends in KKK, which may be a synthetic modification but
+is not documented as such by its source. Full recorded sequences remain the
+scoring reference; no terminal residues are silently removed.
 Its metric is **vaccine sequence contained in reconstruction**: the recorded
 vaccine peptide must appear as a contiguous substring of a reconstructed protein,
 which can extend on either side. Every counted occurrence is tied to an RNA call

pipeline/build_site.py

@@ -53,14 +53,52 @@ def current_peptide_report(catalogue: dict, results: dict | None) -> dict | None
     return report
 
 
+def source_rna_evidence(variant: dict) -> list[dict]:
+    samples = {sample.bam_url.rsplit("/", 1)[-1]: sample.name for sample in SAMPLES}
+    return [
+        {**entry, "benchmark_sample": samples.get((entry.get("bam_file") or "").rsplit("/", 1)[-1])}
+        for entry in variant.get("assay_support", [])
+        if entry.get("tissue") == "tumor"
+        and assays.ASSAY_META.get(entry.get("assay_type"), {}).get("kind") == "rna"
+    ]
+
+
+def rna_support_tier(evidence: list[dict], recovery: dict | None = None) -> dict:
+    """Maximum observed mutant count, never pooled across libraries or methods."""
+    def classify(entries):
+        counts = [entry["alt_reads"] for entry in entries
+                  if isinstance(entry.get("alt_reads"), int) and entry["alt_reads"] >= 0]
+        maximum = max(counts, default=None)
+        return {"max_alt_reads": maximum, "category": "unknown" if maximum is None else
+                "multiple" if maximum >= 2 else "single" if maximum == 1 else "zero"}
+
+    by_sample = {}
+    for sample in SAMPLES:
+        sid = classify([entry for entry in evidence if entry.get("benchmark_sample") == sample.name])
+        raw = (recovery or {}).get("samples", {}).get(sample.name, {}).get("arms", {}).get("reads", {})
+        calls = raw.get("rna_variant_calls")
+        # Only the raw-read arm is one input RNA per transcript model. Corrected
+        # and assembled inputs reuse those reads and must not raise this tier.
+        exacto_count = len({call["transcript_model_id"] for call in calls}) if calls is not None else None
+        by_sample[sample.name] = {
+            **classify([{"alt_reads": sid["max_alt_reads"]}, {"alt_reads": exacto_count}]),
+            "sid_alt_reads": sid["max_alt_reads"], "exacto_raw_rnas": exacto_count,
+        }
+    return {
+        **classify([{"alt_reads": counts["max_alt_reads"]} for counts in by_sample.values()]),
+        "by_sample": by_sample,
+        "missing_samples": [name for name, counts in by_sample.items() if counts["category"] == "unknown"],
+        "basis": "Maximum of Sid mutant reads and distinct Exacto raw-read allele-call inputs per benchmark library; never summed",
+    }
+
+
 def peptide_report_with_rna_evidence(report: dict, catalogue: dict, results: dict | None = None) -> dict:
     """Attach source counts without pooling libraries or inferring missing counts.
 
     These are allele counts at the locus, not vaccine peptide observations.
     Match benchmark samples by the actual BAM basename, not timepoint alone.
     """
     by_id = {variant["variant_id"]: variant for variant in catalogue["variants"]}
-    samples = {sample.bam_url.rsplit("/", 1)[-1]: sample.name for sample in SAMPLES}
     report = json.loads(json.dumps(report))
     scored = {variant["variant_id"]: variant for variant in (results or {}).get("variants", [])}
     for variant in report["variants"]:
@@ -73,12 +111,8 @@ def peptide_report_with_rna_evidence(report: dict, catalogue: dict, results: dic
                 for sample, values in scored.get(variant["variant_id"], {}).get("samples", {}).items()
                 for arm, entry in values.get("arms", {}).items()
             ]
-        variant["source_rna_support"] = [
-            {**entry, "benchmark_sample": samples.get((entry.get("bam_file") or "").rsplit("/", 1)[-1])}
-            for entry in by_id[variant["variant_id"]].get("assay_support", [])
-            if entry.get("tissue") == "tumor"
-            and assays.ASSAY_META.get(entry.get("assay_type"), {}).get("kind") == "rna"
-        ]
+        variant["source_rna_support"] = source_rna_evidence(by_id[variant["variant_id"]])
+        variant["rna_support"] = rna_support_tier(variant["source_rna_support"], scored.get(variant["variant_id"]))
     return report
 
 
@@ -286,6 +320,7 @@ def build_payload() -> dict:
                 "assay_matrix": assays.matrix(variant),
                 "germline_matrix": assays.germline(variant),
                 "recovery": recovery,
+                "rna_support": rna_support_tier(source_rna_evidence(variant), recovery),
             }
         )
 

tests/test_vaccine_peptide_analysis.py

@@ -146,3 +146,39 @@ def test_source_rna_counts_match_the_actual_bam_and_do_not_pool_libraries():
     assert evidence[0]["alt_reads"] == 4
     assert evidence[1]["benchmark_sample"] is None
     assert "source_rna_support" not in report["variants"][0]
+
+
+def test_rna_support_tiers_use_maximum_without_pooling_or_imputing_missing_samples():
+    rows = [
+        {"benchmark_sample": "T1-ONT", "alt_reads": 1},
+        {"benchmark_sample": "T2-ONT", "alt_reads": 1},
+        {"benchmark_sample": "T2-ONT", "alt_reads": 1},
+        {"benchmark_sample": None, "alt_reads": 100},
+    ]
+    tier = build_site.rna_support_tier(rows)
+    assert tier["category"] == "single"
+    assert tier["max_alt_reads"] == 1
+    assert tier["by_sample"]["T2-ONT"]["category"] == "single"
+    assert "T1-PacBio" in tier["missing_samples"]
+    rows.append({"benchmark_sample": "T3-ONT", "alt_reads": 2})
+    assert build_site.rna_support_tier(rows)["category"] == "multiple"
+    assert build_site.rna_support_tier([])["category"] == "unknown"
+    assert build_site.rna_support_tier([{"benchmark_sample": "T1-ONT", "alt_reads": None}])["category"] == "unknown"
+    assert build_site.rna_support_tier([{"benchmark_sample": "T1-ONT", "alt_reads": 0}])["category"] == "zero"
+
+
+def test_raw_rna_evidence_promotes_a_tier_without_double_counting_calls_or_assemblies():
+    evidence = [{"benchmark_sample": "T2-ONT", "alt_reads": 1}]
+    calls = [{"transcript_model_id": "read-a"}, {"transcript_model_id": "read-a"}]
+    recovery = {"samples": {"T1-PacBio": {"arms": {
+        "reads": {"rna_variant_calls": calls},
+        "corrected": {"rna_variant_calls": [{"transcript_model_id": str(i)} for i in range(20)]},
+    }}}}
+    assert build_site.rna_support_tier(evidence, recovery)["category"] == "single"
+    calls.append({"transcript_model_id": "read-b"})
+    tier = build_site.rna_support_tier(evidence, recovery)
+    assert tier["category"] == "multiple"
+    assert tier["max_alt_reads"] == 2
+    assert tier["by_sample"]["T1-PacBio"]["sid_alt_reads"] is None
+    assert tier["by_sample"]["T1-PacBio"]["exacto_raw_rnas"] == 2
+    assert tier["by_sample"]["T2-ONT"]["category"] == "single"

tests/test_web_peptides.js

@@ -3,7 +3,8 @@ const fs = require("node:fs");
 const path = require("node:path");
 const test = require("node:test");
 const { sequenceRows, highlightedPieces, sharedOffset, referenceGroups,
-  sequenceComparison, alignmentSlice, reconstructionWindows, windowDifferences, windowSlice, rnaSupport, sequenceStatus, sourceRnaState } = require("../web/peptides.js");
+  sequenceComparison, alignmentSlice, reconstructionWindows, windowDifferences, windowSlice, rnaSupport, sequenceStatus, sourceRnaState,
+  vaccineComparison, vaccineDifferences, windowEvents, referenceEvents, prepareComparisons, topWindows } = require("../web/peptides.js");
 
 const reportPath = path.join(__dirname, "../results/vaccine_peptide_analysis.json");
 const report = fs.existsSync(reportPath) ? JSON.parse(fs.readFileSync(reportPath, "utf8")) : null;
@@ -338,3 +339,127 @@ test("RNA support counts input RNAs once per sample and method, without counting
   assert.deepEqual(rnaSupport(data, hits, { sample: "T1-ONT", method: "reads" })[0].methods,
     [{ method: "reads", count: 2 }]);
 });
+
+function annotatedFixture() {
+  const sequence = 'SFSGPGMSGMALMEVNLLSGKKK';
+  const variant = { variant_id: 'CD109', gene: 'CD109', ref: 'G', alt: 'T', protein_change: 'p.Arg1310Met',
+    reconstructed_candidates: [], published_vaccine_peptides: [{ peptide_id: 'p1', sequence, in_vaccines: ['V3'], matches: [] }] };
+  const data = { analysis: { artifacts: [{ sample: 'T1', arm: 'reads' }, { sample: 'T1', arm: 'corrected' }] },
+    variants: [variant], reconstructions: {} };
+  function add(id, sequence, position, matched = false, arm = 'reads') {
+    const candidate = { sample: 'T1', arm, reconstruction_id: `${arm}/${id}`, protein_id: id,
+      rna_call_id: id, variant_amino_acid_start: position };
+    variant.reconstructed_candidates.push(candidate);
+    data.reconstructions[candidate.reconstruction_id] = sequence;
+    if (matched) variant.published_vaccine_peptides[0].matches.push({ ...candidate, amino_acid_starts: [1] });
+  }
+  return { data, variant, add, sequence };
+}
+
+test('CD109 vaccine mismatches highlight KKK versus FMV on both rows, even in W1', () => {
+  const { data, add } = annotatedFixture();
+  add('short|orf_0-59', 'MSGMALMEVNLLSGFMVPSEA', 1);
+  const rows = sequenceRows(data), [group] = prepareComparisons(data, rows);
+  const entry = group.references[0], window = group.windows[0];
+  const diff = vaccineComparison(window, entry);
+  assert.equal(diff.missing, 6); // absent prefix is coverage, not a substitution
+  assert.deepEqual(diff.substitutions.map(x => [x.position, x.reference, x.residue]),
+    [[21, 'K', 'F'], [22, 'K', 'M'], [23, 'K', 'V']]);
+  const highlights = vaccineDifferences(group);
+  assert.equal(highlights.windows.get(window).size, 3);
+  assert.deepEqual([...highlights.references.get(entry).keys()], [...highlights.windows.get(window).keys()]);
+  assert.equal(group.comparisons.get(window).substitutions.length, 0);
+  // A protein starting at the mutation has only four local anchor residues;
+  // this intentionally does not infer a vaccine-row target position.
+  assert.equal(referenceEvents(entry, group.windows, group.events).size, 0);
+});
+
+test('target codon uses RNA-call coordinates, independently of the vaccine mismatch', () => {
+  const { data, add } = annotatedFixture();
+  add('full|orf_0-68', 'SFSGPGMSGMALMEVNLLSGFMV', 7);
+  const rows = sequenceRows(data), [group] = prepareComparisons(data, rows);
+  const event = group.events.get(group.windows[0]);
+  assert.equal(event.size, 1);
+  assert.equal(event.get(6).kind, 'mutation');
+  assert.equal(group.vaccineDifferences.windows.get(group.windows[0]).has(6), false);
+  assert.deepEqual([...referenceEvents(group.references[0], group.windows, group.events).keys()], [6]);
+});
+
+test('grouped target locations preserve ambiguity and do not invent a reference annotation', () => {
+  const { data, add, sequence } = annotatedFixture();
+  add('a|orf_0-68', sequence, 7, true);
+  add('b|orf_0-68', sequence, 9, true);
+  const [group] = prepareComparisons(data, sequenceRows(data));
+  const events = group.events.get(group.windows[0]);
+  assert.deepEqual([...events.keys()].sort((a,b) => a-b), [6, 8]);
+  assert.ok([...events.values()].every(x => x.kind === 'ambiguous'));
+  assert.equal(referenceEvents(group.references[0], group.windows, group.events).size, 0);
+});
+
+test('frameshift tails survive cropping without claiming novel ORF boundaries', () => {
+  const { data, variant, add } = annotatedFixture();
+  variant.ref = 'AG'; variant.alt = 'A';
+  add('a|orf_0-68', 'SFSGPGMSGMALMEVNLLSGFMV', 1);
+  const [group] = prepareComparisons(data, sequenceRows(data));
+  const window = group.windows[0];
+  const cropped = { ...window, sequence: window.sequence.slice(4), offset: window.offset + 4 };
+  const events = windowEvents(cropped, variant);
+  assert.equal(events.get(0).kind, 'mutation');
+  assert.equal(events.get(4).kind, 'tail');
+  assert.match(events.get(4).title, /not annotated/);
+});
+
+test('top windows rank RNA inputs, not vaccine containment or row order, and retain ties', () => {
+  const { data, add, sequence } = annotatedFixture();
+  add('match|orf_0-68', sequence, 7, true);
+  add('mismatch1|orf_0-68', 'SFSGPGMSGMALMEVNLLSGFMV', 7);
+  add('mismatch2|orf_0-68', 'SFSGPGMSGMALMEVNLLSGFMV', 7);
+  add('match|orf_3-71', sequence, 7, true); // same input, another ORF: still one RNA
+  add('corrected|orf_0-68', sequence, 7, true, 'corrected');
+  let rows = sequenceRows(data), groups = prepareComparisons(data, rows);
+  let top = topWindows(data, groups, rows);
+  const raw = top.find(x => x.method === 'reads'), corrected = top.find(x => x.method === 'corrected');
+  assert.equal(raw.count, 2);
+  assert.equal(raw.leaders.length, 1);
+  assert.notEqual(raw.leaders[0].window.number, 1);
+  assert.equal(raw.disagreement, true);
+  assert.equal(corrected.count, 1);
+  assert.equal(corrected.disagreement, false);
+  add('match2|orf_0-68', sequence, 7, true);
+  rows = sequenceRows(data); groups = prepareComparisons(data, rows);
+  top = topWindows(data, groups, rows, { method: 'reads' });
+  assert.equal(top.length, 1);
+  assert.equal(top[0].leaders.length, 2);
+  assert.deepEqual(top[0].leaders.map(x => x.peptides[0].status).sort(), ['contained', 'disagreement']);
+});
+
+test('top-window missing coverage, unaligned output and no output stay separate from disagreement', () => {
+  const { data, add } = annotatedFixture();
+  add('short|orf_0-50', 'MSGMALMEVNLLSGKKK', 1);
+  add('unrelated|orf_0-29', 'AAAAAAAAAA', 3, false, 'corrected');
+  data.analysis.artifacts.push({ sample: 'T2', arm: 'reads' });
+  const rows = sequenceRows(data), groups = prepareComparisons(data, rows);
+  const tops = topWindows(data, groups, rows);
+  assert.equal(tops[0].leaders[0].peptides[0].status, 'incomplete');
+  assert.equal(tops[1].leaders[0].peptides[0].status, 'unaligned');
+  assert.deepEqual(tops[2].leaders, []);
+  assert.ok(tops.every(x => !x.disagreement));
+});
+
+test('RNA-linked containment cannot leak into a leader from a different method', () => {
+  const { data, add, sequence } = annotatedFixture();
+  add('raw|orf_0-68', sequence, 7, false);
+  add('corrected|orf_0-68', sequence, 7, true, 'corrected');
+  const rows = sequenceRows(data), groups = prepareComparisons(data, rows);
+  const tops = topWindows(data, groups, rows);
+  assert.equal(tops[0].leaders[0].peptides[0].status, 'unconfirmed');
+  assert.equal(tops[1].leaders[0].peptides[0].status, 'contained');
+});
+
+test('RNA evidence tiers follow sample selection and remain independent of method', () => {
+  const { data, variant } = annotatedFixture();
+  variant.rna_support = { category: 'multiple', by_sample: { T1: { category: 'single' }, T2: { category: 'multiple' } } };
+  assert.equal(sequenceRows(data, { rnaTier: 'multiple' }).length, 1);
+  assert.equal(sequenceRows(data, { rnaTier: 'multiple', sample: 'T1' }).length, 0);
+  assert.equal(sequenceRows(data, { rnaTier: 'single', sample: 'T1', method: 'corrected' }).length, 1);
+});

tests/test_web_results.js

@@ -195,3 +195,35 @@ test("an empty later page stops polling that jobs snapshot", async () => {
   assert.equal(node.children.find((child) => child.className === "live-jobs").children.length, 100);
   assert.equal(requests.length, 3);
 });
+
+test('amino-acid failures are named beside the score and show the actual failed residue', () => {
+  const data = result(3, 2, { reads: { outcome: 'proteoform' } });
+  data.variants = [
+    { gene: 'ABI3BP', protein_change: 'p.Arg472Met', rna_support: { category: 'multiple' }, recovery: {
+      outcome: 'proteoform', residue_confirmed: false, expected: { alt_aa: 'M' }, samples: {
+        T2: { arms: { reads: { proteoforms: [{ context: 'AAWBB', context_start: 10,
+          mutant_residue_indices: [11], mutant_residues: 'W' }] } } },
+      },
+    } },
+    { gene: 'CONTROL', rna_support: { category: 'single' }, recovery: { outcome: 'peptide', residue_confirmed: true } },
+    { gene: 'NO_PROTEIN', rna_support: { category: 'zero' }, recovery: { outcome: 'rna_only', residue_confirmed: false } },
+  ];
+  const nodes = render(data, ['renderVerdict', 'renderEvidenceBreakdown']);
+  assert.match(nodes['#verdict'].textContent, /1 amino-acid failures: ABI3BP/);
+  assert.doesNotMatch(nodes['#residue-failures'].textContent, /CONTROL|NO_PROTEIN/);
+  assert.match(nodes['#residue-failures'].textContent, /Expected M; observed W/);
+  const sequence = nodes['#residue-failures'].children[1].children.find(x => x.tag === 'code');
+  assert.equal(sequence.children[2].className, 'sequence-vaccine-difference');
+  assert.equal(sequence.children[2].textContent, 'W');
+  assert.match(nodes['#rna-breakdown'].textContent, /2\+ mutant RNA reads11\/10\/1/);
+  assert.match(nodes['#rna-breakdown'].textContent, /Single-read evidence11\/11\/1/);
+});
+
+test('RNA support groups do not score failed or missing outputs as failed recovery', () => {
+  const data = result(0, 0);
+  data.variants[0].rna_support = { category: 'multiple' };
+  data.variants[0].recovery.outcome = 'no_reads';
+  const nodes = render(data, ['renderEvidenceBreakdown']);
+  assert.match(nodes['#rna-breakdown'].textContent, /2\+ mutant RNA reads1——/);
+  assert.doesNotMatch(nodes['#rna-breakdown'].textContent, /0\/1/);
+});