pirl-unc / DoesExactoWorkYet

Does Exacto work yet?

Every vaccine neoantigen from Sid Sijbrandij's osteosarcoma is a mutant protein someone already believed in enough to manufacture. This runs Exacto over every long-read RNA-seq sample from his tumour — ONT at three biopsies and PacBio Iso-Seq at the first — and asks a single question: does the mutant protein sequence come back out?

Recovery by sample

One column per sequencing sample, not per biopsy: T1 was sequenced on ONT and again on PacBio, and the two can disagree. Each cell is the best outcome Exacto reached for that mutation in that sample, across both arms of the pipeline.

Somatic variants used in a vaccine

somatic mutations from Sid's tumour that were selected into at least one of his personalised vaccines. Click a row for the per-arm Exacto detail, including any mutant protein sequence recovered.

The VAF columns pool every sample the portal sequenced for that assay and timepoint; each cell shows the variant allele fraction over the supporting and total read counts. The table scrolls sideways — gene, change and verdict stay pinned.

Paths to a contig

Exacto does not read FASTQs — it reads a BAM of transcripts, and how those are produced is decided before Exacto ever runs. Two routes are taken here. The funnels show how much sequence survives each one; the ladder below shows how many of the 37 vaccine mutations each route actually reaches. They do not agree, and the disagreement is the most consequential thing this test has measured.

What fraction of the variants each path reaches

Counted across all samples: a mutation is credited to a route if that route reached the rung in any sample. Both routes start from the same reads, so the first rung is identical by construction — everything below it is the route's own doing.

Where mutations drop out, per sample

The funnels above count sequence. This counts mutations: how many of the 37 are still in play at each stage, for every sample and method. The gap between allele present in the RNA and Exacto called it is the cost of everything in between — the filters, the preparation, and the caller's own thresholds — and the gap after that is translation.

Methods benchmark: sensitivity against specificity

Every way of turning reads into transcripts, scored on the same mutations. Sensitivity is measured against the mutations whose allele is actually in that sample's RNA, not against all 37 — scoring a method for mutations that are not expressed would punish it for the tumour's biology. The other three columns are what that sensitivity costs. Methods and their parameter sets are defined in pipeline/methods.py; a threshold sweep is a row here, not a code path.

Why assembly is cleaner and still misses more

The tables below are measured. The reading of them that follows is conjecture, labelled as such, and is the kind of claim this harness could be extended to test rather than argue.

What the numbers support

Assembly does produce better reading frames, and the gradient is exactly what sequencing error predicts. Raw ONT reads frameshift most; PacBio Iso-Seq reads, which are already circular-consensus corrected, frameshift far less; RNA-Bloom2 contigs, which are a consensus over many reads, do not frameshift at all in this run. An indel miscalled in a homopolymer shifts the frame and everything downstream of the mutation is wrong, and collapsing many reads into one contig is precisely the operation that removes an error appearing in a minority of them.

What it does not buy is recovery. Every mutation the assembly route finds, the reads route also finds, and the frameshifts on the reads route largely do not matter for this question: the mutated residue itself still comes out correct in the great majority of proteoforms, and wherever a published vaccine epitope existed alongside a proteoform, that epitope was found verbatim. A frameshift downstream of the mutation ruins the tail of the protein without moving the codon under test.

Conjecture: why it misses so much

Assembly is a consensus, and a subclonal mutation is a minority signal. The VAF table is the tell — the mutations both routes recover sit near a third of reads, the ones only the reads route recovers sit near a seventh, and the ones nobody recovers are down near one in a hundred. A contig is built from many reads at a locus; if the mutant allele is a minority of them, the consensus is the wild-type sequence and the mutation is averaged away before Exacto ever sees it. Handing reads over individually cannot lose a variant that way, because each read carries its own allele.

The support filter compounds it. nexus_filter_rnabloom2_transcripts requires three reads behind a contig. For a clonal variant that is trivial; for one at 5% VAF it demands depth the sample may not have at that locus, and the mutant contig is dropped while the wild-type one passes. The filter is not wrong — it exists to remove assembly artefacts — but its cost falls specifically on the low-VAF variants a neoantigen test most wants.

And the splice filter is the largest single loss. About three quarters of what survives assembly and alignment is then removed by remove-unspliced-rnas. A plausible reading is that subsampled input produces fragmentary contigs which no longer span an intron and so no longer look spliced — in which case the loss is partly an artefact of this harness's read caps, not of the method.

Where assembly should still win

None of this says assembly is useless. It is the only route here that reconstructs a full-length transcript, so it is the one that can say which isoform carries a mutation and what the complete protein looks like — a question this test does not currently ask. Its proteoforms are frameshift-free, which matters for any downstream use that depends on sequence far from the mutation. And its advantages should grow with depth: the failure mode conjectured above is a coverage problem, so a clonal variant in a deeply sequenced sample is where assembly should be at its strongest, which is consistent with the three it did recover being the highest-VAF ones in the set.

Testable next: raise the per-variant read caps and re-run the assembly arm to separate genuine method loss from harness-induced loss; add a second assembler to distinguish RNA-Bloom2's behaviour from assembly in general; and relax the Nexus support threshold to see whether low-VAF mutant contigs reappear.

How this run is configured

Every option this pipeline sets, measured against the canonical pipeline — Andy Lee's PEPTIDE_PREDICTION_EXACTO subworkflow in Nexus, or Exacto's own defaults where Nexus passes nothing. Anything not marked stock is a deliberate choice with a reason attached.

What broke along the way

Issues in Exacto that this harness has run into, and what the pipeline does about each. Full write-ups, with the failing commands and their output, are on the bug reports page.

Pipeline runs

One run per sequencing sample per arm. Every Exacto invocation's exit status and wall time is recorded, so a crash shows up as a crash rather than a zero.

Engineering log

Every change to this harness, newest first, straight from the commit log — what broke, what the evidence was, and what was done about it. Published so a number on this page can be traced back to the reasoning that produced it. Distinct from bug reports, which are about Exacto; these are about the test.

How this works

  1. Variants. The vaccine neoantigen list and the per-assay read counts are pulled live from osteosarc.com, merged on locus, and written to results/vaccine_variants.json.
  2. Reference. Each mutation's GENCODE gene body (±10 kb) is fetched from hg38 over HTTP byte ranges and written back at its true coordinate inside otherwise all-N chromosomes. Coordinates stay identical to hg38; only a few megabases have to be indexed.
  3. Reads. The three dedup ONT BAMs and the PacBio Iso-Seq BAM total ~159 GB and stay on Backblaze — htslib range-reads only the blocks over our loci. Coverage spans five orders of magnitude across these genes, so reads whose alignment covers a variant are capped per variant and everything else is capped per region, both sampled with a fixed seed. The caps sit far above what any caller needs to decide a locus.
  4. Alignment. The portal's BAMs were aligned without minimap2's --cs tag, which Exacto needs, so everything is realigned.
  5. Exacto. Two arms. assembly follows the documented pipeline — RNA-Bloom2 assembles spanning and context reads into transcripts, remove-unspliced-rnas filters, then call-rna-vars. reads skips the assembler and hands the variant-spanning reads straight in as transcripts. Both then run annotate-vars → integrate-vars → translate-structs → call-peptide-vars.
  6. Verdict. A mutation counts as recovered when a translated primary structure carries an amino acid Exacto marked mutant, traced back to an RNA variant Exacto called at that exact locus with that exact allele — never via integrate-vars, which links a DNA variant to anything within 10 kb of a transcript edge and would inflate the count several-fold. For missense mutations the residue is additionally checked against the portal's HGVS annotation — a change at the right codon but the wrong amino acid is reported, not counted as a win.
  7. Does it match the vaccine? The portal also publishes the curated pVACtools run the vaccine designs were picked from. For the ten mutations it covers, its MT Epitope Seq values are searched for inside Exacto's translated proteoform — the strictest form of the question available, since those are the peptides that were actually made.

What this test does not do

  • It does not re-discover the somatic variants from DNA. Sid's WGS is short-read; Exacto's DNA callers want long reads. The portal's curated calls are supplied as the somatic callset, and Exacto is asked whether it can find them in the RNA and translate them.
  • Reads are only pulled from the vaccine genes, so this measures sensitivity at known loci, not genome-wide precision. Masking the rest of the genome also removes paralogues that would otherwise compete for alignments.
  • call-peptide-vars is run against the reference proteins of the tested genes, so "novel" means absent from the gene's own wild-type isoforms.
  • Read counts are capped, so the assembly arm sees less depth than a whole-sample run would at the most highly expressed loci. Both the kept and the available counts are shown in each variant's detail.
  • The exact-peptide check only covers the ten mutations the published pVACtools run includes. For the other twenty-seven the strongest available check is the amino acid itself.
  • Exacto 0.4.6a1 cannot finish its own documented pipeline on this data without the workarounds listed above. Each is applied deliberately and recorded in the run's JSON — none of them is silent.
  • A recovered proteoform is not automatically the neoantigen. longest_orf has no reference CDS to anchor on and often picks the wrong frame on a truncated transcript, which is why the residue is checked separately and mismatches are flagged rather than counted.