Extend NMD detection and escape-rule determination to all variant sources (currently fusion-only)

[riasc]

Motivation

NMD (nonsense-mediated decay) detection is currently wired only for fusion events. For SNV / indel (and exitron / alt-splicing) variants the NMD and NMD_escape_rule output columns come out empty.

NMD escape determines whether a PTC-introducing variant's truncated/frameshifted protein is actually translated and presentable. Indels and frameshifts are a large share of ScanNeo2's neoantigen yield, so not evaluating NMD for them systematically over- or under-calls those neoantigens. This issue: extend NMD detection and escape-rule determination to all variant sources.

Current state

The machinery substantially exists already — this is "extend + wire up", not greenfield:

• workflow/scripts/prioritization/effects.py has find_stop_codon, annotate_stop_codon, and check_escape — check_escape already returns escape-rule codes (e.g. 2/3/4, or -1 for no escape) — but the PTC/escape block is exercised only on the fusion path.
• The neoepitope output table already carries NMD, PTC_dist_ejc, PTC_exon_number, NMD_escape_rule columns; for non-fusion variants they are empty.
• effects.py carries a # TODO: check for NMD out of VEP.

Approach

Two routes for the non-fusion variants:

• (a) VEP NMD annotation — what the existing # TODO points at. Low effort, but VEP tends to give a yes/no rather than the specific escape rule.
• (b) Reuse effects.py's own escape logic — feed the PTC genomic position and the transcript exon model (already loaded as self.exome) into find_stop_codon / check_escape for SNV / indel / exitron / alt-splicing variants, exactly as the fusion path already does.

Recommended: (b) — one uniform NMD path across all variant sources: transparent, consistent with the fusion path, and it produces the explicit escape rule (not just a yes/no).

Escape rules to cover

check_escape already encodes some of these; the standard set is:

1. PTC in the last exon (no downstream exon-exon junction).
2. PTC within ~50–55 nt of the last exon–exon junction.
3. Start-proximal PTC (close to the start codon → translation re-initiation).
4. Single-exon / intronless transcript.
5. (optional) long last exon.

The NMD_escape_rule codes are currently bare numbers — document which code maps to which rule.

Scope / acceptance

• NMD / NMD_escape_rule populated for SNV, indel, exitron and alt-splicing variants — not just fusions.
• A documented escape-rule code set.
• Verify on a real run that the NMD columns are filled for the non-fusion variant types.

[riasc]

Engineering plan (drafted during attempted v0.3.15 work, deferred to v0.4.x)

While drafting the wiring fix (remove the VEP NMD pre-filter so all frameshift variants reach determine_NMD), an investigation surfaced a deeper correctness bug that has to be fixed first.

The deeper bug

For non-fusion variants, variants.py:80-90 constructs the "transcript" by slicing the genome:

transcript = annotation.ref[chrom][bnds[0]+1:start+2]   # genomic span
transcript += field["Allele"]
transcript += annotation.ref[chrom][stop+1:bnds[1]+1]   # genomic span

bnds[0] and bnds[1] come from Annotation.parse_transcriptome, which reads the GTF "transcript" feature row — i.e., genomic start/end including introns. The constructed transcript is therefore pre-mRNA, not spliced mRNA.

effects.py:find_stop_codon then walks this pre-mRNA in 3-nt steps. Because intron sequence is effectively random, the algorithm typically finds a spurious in-frame TAA/TAG/TGA triplet within the first intron downstream of the variant — long before reaching a real exonic PTC. The genomic coordinate it returns (stop_coord) is then intronic, and annotate_stop_codon finds no exon containing that coord, so it returns (-1, -1) and determine_NMD leaves the NMD/PTC columns as None.

So the empirical "NMD is empty for non-fusion variants" is not (only) a wiring issue — it's the algorithm walking unspliced sequence and silently bailing out. The fusion path is correct because fusions.py constructs the transcript from arriba's column 27, which is the already-spliced fusion peptide.

What the proper fix needs to do

1. Track strand per transcript. parse_transcriptome stores [start, end]; needs to also store strand from GTF column 7. Bumps the data structure to [start, end, strand] — consumers need to be audited and updated.

2. Add Annotation.splice_transcript(transcript_id, chrom) that builds spliced mRNA: iterate self.exome[tid] by exon number, take ref[chrom][s:e+1] for each exon, concatenate, and reverse-complement if strand is -.

3. Replace transcript construction in variants.py:80-90: get the spliced mRNA via splice_transcript, map the variant's genomic position to its mRNA-relative position by walking exons (strand-aware), splice the variant allele in at that mRNA position, set transcript_bp to the mRNA position.

4. Adapt find_stop_codon / annotate_stop_codon: stop position is now mRNA-relative, so the round-trip back to a genomic coordinate has to walk the exons again. annotate_stop_codon's dist_ejc = exon_end - stop_coord calculation is +-strand-specific; for - strand it should be stop_coord - exon_start.

5. Leave the fusion path alone — arriba already provides spliced fusion peptide sequence.

6. Then the wiring fix (remove the early continue at variants.py:74-75 so all frameshift variants reach determine_NMD) lands cleanly on top.

Scope estimate

~150-300 lines across reference.py, variants.py, effects.py. Real-data SLURM test mandatory — output table will gain rows (NMD-target variants previously dropped); NMD column should be populated meaningfully across all vartypes; spot-check 5-10 variants by hand against transcript exon structure.

Open uncertainties to think through during implementation

• Variant in 5' UTR: PTC concept doesn't really apply; current code's re.search(r'ATG', transcript) for start codon detection can land on a UTR ATG. Defensive handling needed.
• Variant at exon-intron boundary: genomic-to-mRNA mapping needs a clear side-choice. Rare but worth a defensive branch.
• VEP transcript-ID matching: field['Feature'] must match a key in self.exome and self.transcriptome. If GTF has versioned IDs (ENST00000123.5) but VEP gives unversioned (ENST00000123), we strip versions; needs verification on the project's specific VEP setup.

Why we deferred

This is genuine multi-file engineering with non-trivial test surface. v0.3.15 is otherwise small/safe; doing this properly in v0.3.15 would balloon scope. v0.4.x will have the SLURM-testing bandwidth this kind of behavioural change needs.

The escape-rule docstring update I drafted for effects.py:check_escape (mapping return codes 1-4 / -1 to rule descriptions) is also deferred — would land naturally alongside the splicing fix.