neoantigen-pipeline

name	neoantigen-pipeline
description	Build a bioinformatics pipeline to identify neoantigen candidates from tumor whole-genome sequencing data and produce IVT plasmid sequences for mRNA vaccine production. Use when working with VCF files, somatic variant filtering, protein annotation, MHC binding prediction, epitope generation, or GenBank IVT plasmid design for cancer immunotherapy.

Neoantigen Pipeline — Agent Skill

Guide for building a pipeline that goes from a multi-sample Mutect2 VCF to ranked neoantigen candidates to IVT plasmid GenBank files for mRNA vaccine production. Based on hard-won lessons from the canine mammary tumor pipeline.

Pipeline Steps

0. Download / validate input data (VCF + GFF3 genome annotation)
1. Filter somatic variants (germline removal from raw Mutect2)
2. Annotate variants genome-wide (GFF3 gene/CDS mapping + NCBI CDS fetch)
3. Predict MHC binding (MHCflurry epitope scoring for all coding variants)
4. Rank neoantigen candidates (best candidate per gene, genome-wide)
5. Cross-reference target genes (filter candidates against curated gene list)
6. Design IVT plasmid (select specific genes → pUC19 + T7 cassette → GenBank)

Target gene filtering is applied after genome-wide ranking (step 5), so the full neoantigen landscape is available before any filtering. Vaccine gene selection (step 6) picks specific genes from the genome-wide list by name.

For detailed technical reference on each step, see reference.md.

Critical Rules

VCF Validation (Step 0–1)

• Check filtering status first. If the FILTER column is . everywhere, the VCF is unfiltered Mutect2. You MUST apply your own germline filter — do not assume variants are somatic.
• Check the genome assembly. Read ##reference= and ##contig= from the VCF header. All gene coordinates, NCBI fetches, and annotation must use the same assembly. GSD data uses UU_Cfam_GSD_1.0 (GCF_011100685.1), not CanFam3.1.
• Normal AF alone is unreliable. In raw Mutect2, the normal sample often has 0 alt reads (AD = XX,0), making normal AF = 0 even for germline variants. Use cross-tumor AF patterns instead.

Germline Filtering Heuristic (Step 1)

When FilterMutectCalls was not run, classify by tumor AF pattern:

Class	Rule
GERMLINE_HOM	Mean AF > 0.85, range < 0.15, present in >= 5/7 tumors
GERMLINE_HET	Mean AF 0.30–0.65, range < 0.35, present in >= 5/7 tumors
SOMATIC	Everything else (variable AF or limited tumor presence)

GFF3 Annotation & NCBI CDS Fetch (Step 2)

• Use the assembly GFF3 for gene/CDS mapping. Parse gene, mRNA, and CDS features from the GFF3 to build a genome-wide gene/CDS interval index. This is faster and more reliable than per-gene NCBI feature table queries.
• GFF3 attribute parsing: The transcript_id for CDS records is embedded in the Parent attribute (e.g., Parent=rna-XM_...). Strip the rna- prefix to get the actual transcript accession.
• Cache the GFF3 parse result (e.g., as JSON) — re-parsing 700K+ CDS records is expensive.
• Use protein accession for CDS sequence fetch. efetch?db=protein&id={NP_acc}&rettype=fasta_cds_na is more reliable than fetching via mRNA accession. Cache responses to disk.
• Minus-strand genes: CDS exons must be reverse-complemented. Invert CDS offset: cds_offset = total_cds_len - 1 - forward_offset. Apply reverse_complement() to REF/ALT bases.
• NCBI rate limits: Add 0.35s delay between API requests. Handle transient HTTP errors (502, 429) gracefully — log and continue.

NCBI Feature Table Parsing (legacy fallback)

If GFF3 is unavailable and you must query NCBI feature tables per gene:

• Continuation vs feature lines: Continuation lines have exactly 2 columns (start, end). Feature lines have 3+ columns (start, end, feature_type). Never treat a 3-column line as an exon continuation — this causes overlapping gene corruption (e.g., WRAP53 exons appended to TP53).
• Partial coordinates: Strip < and > before int(): int(s.lstrip("<>")).
• Exclude non-coding transcripts. Filter out XR_ accessions. Prefer NM_ (curated) over XM_ (predicted).

MHCflurry (Step 3)

• Check predictor.supported_alleles before running. DLA alleles may not all be supported.
• The output DataFrame may use best_allele instead of allele — check column names dynamically.
• Filter out peptides containing X (non-standard amino acids) or * (stop codons) before prediction.
• Affinity tiers: < 50 nM strong, 50–500 nM moderate, 500–5000 nM weak.

IVT Plasmid Design (Step 6)

• Poly(A) tail: Use segmented A(30)-GCATATGACT-A(70) (BNT162b2 design), not a homopolymeric run. Pure poly(A) > 50 nt is unstable during bacterial propagation.
• Linearization site: Verify the chosen restriction enzyme cuts exactly once in the final plasmid. Test both forward and reverse complement. SpeI (ACTAGT) is absent from pUC19 and most coding sequences.
• pUC19 features: Map backbone features (bla, ori, lacZ-alpha) by marker sequence search, not hardcoded positions — positions shift after cassette insertion.
• GenBank format: LOCUS line has strict column formatting. Feature locations are 1-based. Wrap long qualifier values. Use complement() for minus-strand features.

Common Mistakes to Avoid

1. Skipping germline filtering because the VCF has a FILTER column — check if it's actually populated.
2. Using CanFam3.1 coordinates with a GSD-aligned VCF (or vice versa) — all coordinates silently fail.
3. Treating all feature table continuation lines as exons — overlapping genes corrupt the exon list.
4. Fetching CDS by mRNA accession when transcript_id is on the parent mRNA feature, not the CDS — use protein_id from the CDS feature instead.
5. Misreading GFF3 Parent attribute — Parent=rna-XM_... needs the rna- prefix stripped to get the actual transcript accession.
6. Using a homopolymeric poly(A) tail in a plasmid — it will shorten during bacterial propagation.
7. Assuming MHCflurry column names are stable — allele vs best_allele changes between versions.
8. Applying target gene filtering too early — annotate and rank genome-wide first; filter to targets at the end so the full landscape is available for review.

Architecture Advice

• Annotate genome-wide, filter late. Run somatic filtering, GFF3 annotation, MHC prediction, and ranking across the entire genome. Apply target gene selection as a separate, late step. This preserves the full neoantigen landscape for review and lets you pick vaccine genes independently.
• Separate VCF scanning from annotation. The VCF scan is expensive (~5 min). Write somatic variants to a TSV once; run all downstream analyses on the TSV.
• Use JSON for inter-step data. Gene annotations (protein sequences, exon coordinates) serialize cleanly to JSON for consumption by downstream steps.
• One script per logical step. Each script should accept input/output paths via argparse and be runnable independently.
• Composite scoring: Combine tumor penetrance (AF × n_tumors) with MHC binding affinity: rank_score / log10(affinity + 1). This balances clonal abundance against immunogenicity.
• Vaccine gene selection should be explicit. Step 6 selects genes by name from the genome-wide candidates, not by taking the top N from a target-filtered list. This gives full control over which genes enter the final vaccine constructs.