EOCRC Exercise 5: From Splicing Events to Neoantigen Candidates

1. What you will practice

The output of EOCRC Exercise 4 contains genomic coordinates and PSI differences. This exercise converts those coordinates into short peptide candidates.

Step	Input	Output
Extend coordinates	Tumour-enriched splice region	Strand-aware BED12
Extract sequence	BED12 and hg38 FASTA	Spliced nucleotide FASTA
Translate	Nucleotide sequence	Amino-acid sequence by reading frame
Create peptides	Amino-acid sequence	Overlapping 9-mers
Predict MHC-I binding	9-mer and HLA allele	Binding affinity and percentile rank

Correct coordinates with an incorrect reading frame still produce the wrong peptide. Predicted MHC-I binding also does not establish that the peptide is produced, presented on a cell, or recognized by T cells.

2. Audit how far the public input can proceed

The paper started with 82 significant EOCRC Tumor-versus-Normal splice events. It added 27 flanking nucleotides on both sides of tumour-enriched alternatively spliced regions, then used strand and split options in bedtools getfasta to exclude introns.

Every Location in Supplementary Table S6 can be parsed, but the table does not publish the adjacent exon blocks that form each BED12 record. Those blocks are required by bedtools getfasta -split. Guessing them from one start and end coordinate could create nucleotide and peptide sequences different from the paper, so the public-input run stops before extraction.

Execution prompt: audit neoantigen input completeness

Audit whether the public Exercise 4 event table is sufficient for sequence extraction.

1. Read Gene name, Location, delta PSI, FDR, and Event type from outputs/paper-table-s6-82-splicing-events.csv.
2. Parse Location into chromosome, strand, start, and end, recording successes and failures.
3. Count tumour-increased events with positive delta PSI and retain selection_direction.
4. Check each event for BED12 blockCount, blockSizes, and blockStarts or equivalent adjacent-exon structure.
5. If exact exon blocks are absent, do not treat the whole coordinate interval as one exon. Record exact_bed12_records=0.
6. Save total events, tumour-increased events, and exact BED12 records able to proceed to outputs/neoantigen-input-audit.csv and make a bar chart.
7. If exact BED12 records are zero, skip hg38 download, getfasta, translation, 9-mer generation, and NetMHCpan. Explain the missing input and stop.
8. Continue with the conditional steps below only after obtaining author-provided rMATS or Whippet tables or BED12 records.

Actual run: all 82 coordinates parsed and 31 events had increased tumour PSI. Zero records contained enough exon-block information for an exact BED12 reconstruction, so the workflow stopped before extraction, translation, 9-mer generation, and NetMHCpan.

Without -s, a reverse-strand gene is translated in the wrong direction. Without -split, introns between BED blocks may be concatenated into the candidate. These options determine the peptide sequence.

3. Translate three frames and create 9-mers only when BED12 is available

A short genomic sequence around a splice junction does not always identify the true coding frame. The paper used EMBOSS transeq to translate three forward reading frames and removed amino-acid sequences containing stop codons.

Step 2 prompt: translation and 9-mer generation

Convert the nucleotide FASTA into amino-acid and 9-mer candidates.

1. Check whether EMBOSS transeq is installed. If absent, do not install it automatically. State that a separate Biopython function implements the same three forward-frame translation.
2. Translate frames 1, 2, and 3 for every event and audit nucleotide remainder and amino-acid length.
3. Exclude amino-acid sequences containing stop codons, but record excluded event-frame pairs in outputs/translation-audit.tsv.
4. Move a one-residue window to create every 9-mer. Save event ID, gene, frame, peptide_start, peptide, and junction_overlap.
5. Mark flank-only 9-mers that do not cross the splice junction. Do not call every flank peptide tumour-specific merely because it came from a tumour-specific event.
6. Do not discard duplicate peptide sequences before preserving all event provenance.
7. Save outputs/eocrc-splice-9mers.tsv and FASTA-format outputs/eocrc-splice-9mers.fasta.
8. Print event count, valid-frame count, all 9-mers, junction-overlap 9-mers, and duplicate-sequence count, then stop.

When an exact transcript CDS frame is available, using that frame is preferable to retaining all three. This exercise follows the paper’s broad candidate procedure while preserving frame and junction overlap as audit fields.

4. Read NetMHCpan output only after peptide FASTA exists

The paper submitted 9-mers to NetMHCpan 4.0 for four HLA alleles.

HLA-A*01:01
HLA-A*02:01
HLA-A*03:01
HLA-A*24:02

It defined peptides in the top 0.5% for an HLA as potential strong binders. This is the paper’s candidate-selection rule, not a clinical threshold that automatically applies to every version and analysis.

This illustrative output only demonstrates the columns.

peptide	HLA	`EL_Rank`	`BA_Rank`	Classification in this exercise
`SLDVLKQMG`	`HLA-A*02:01`	0.23	0.41	paper-rule strong binder
`YLKQMGDVL`	`HLA-A*02:01`	1.80	2.10	does not pass

Step 3 prompt: audit NetMHCpan output

Audit NetMHCpan MHC-I binding results for the 9-mer candidates.

1. Check for a local NetMHCpan executable and its terms of use. If absent, do not download it or accept a licence automatically. Point to the official service and academic download, then stop after preparing the input FASTA.
2. If it can run, record version, model mode, peptide length, four alleles, and complete command in outputs/netmhcpan-run.txt.
3. Treat NetMHCpan 4.0 as the paper-reproduction setting. If current 4.1 is used, state the EL and BA output differences and version difference.
4. Run outputs/eocrc-splice-9mers.fasta for the four alleles and preserve unmodified raw output.
5. Parse peptide, allele, EL score, EL rank, BA score, BA rank, predicted affinity, and binder level when available to outputs/netmhcpan-parsed.tsv.
6. Only for comparison with the paper, mark the top-0.5-percent rule in a paper_strong_binder column. Keep it separate from the tool's current SB/WB labels.
7. Join peptide provenance back to retain event ID, gene, frame, and junction_overlap.
8. Count genes, events, and strong-binder peptides. Compare with the paper's 18 genes and seven EOCRC-specific events without hiding differences.
9. Do not infer expression, proteasomal processing, cell-surface presentation, or T-cell response from MHC-I binding prediction alone, then stop.

5. Preserve provenance through the final row

Every final candidate row needs enough information to return to its source.

Field	What to verify
`event_id`	Which splicing event started the path?
`coordinates`, `strand`	Which genomic blocks were extracted?
`frame`	Which reading frame was translated?
`peptide_start`	Where does the 9-mer start in the amino-acid sequence?
`junction_overlap`	Does the 9-mer cross the new junction?
`HLA`	Which allele was evaluated?
`rank`, `tool_version`	Which model and selection rule were used?

If provenance is broken, a strong binder cannot be traced back to the event for verification.

6. Questions to answer

Why do BED strand and block information change the peptide?
Why can peptides from all three reading frames not all be called real proteins?
Why should flank-only 9-mers be marked separately?
Are percentile rank and predicted affinity the same value?
What else must be tested before calling a strong MHC-I binder a confirmed neoantigen?

Rosie PoC provides a second starting point that creates 9-mers from DNA variants. Rosie starts from a variant and this exercise starts from an RNA splicing event, but both require preserved peptide provenance.

7. Sources and reproducibility

Record reference FASTA, BED coordinate convention, flank length, strand and split options, translation implementation, peptide length, HLA alleles, NetMHCpan version, mode, command, and rank threshold.