Version Compatibility

Reference examples tested with: Ensembl VEP 111+, pVACtools 4.1+, MHCflurry 2.1+, VAtools 5+, pandas 2.2+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Notes specific to this skill: pVACtools is now at 7.x; positional CLI args are stable across 4.x-7.x but defaults and the supported-algorithm list change between releases — always run pvacseq run --help against the installed build. pVACseq requires the Wildtype and Frameshift VEP plugins; the Downstream plugin was replaced by Frameshift in pVACtools 2.0, so 4.x+ pipelines must NOT use Downstream. A local IEDB install (--iedb-install-directory) is strongly preferred over the rate-limited public API for patient data.

Neoantigen Prediction

"Find neoantigens from my tumor mutations" -> Translate somatic variants into mutant peptides, predict patient-HLA presentation, and rank by tumor-specific quality for vaccine/biomarker use.

• CLI: pvacseq run on a VEP-annotated, expression/readcount-annotated somatic VCF + patient HLA (pVACtools)
• CLI: pvacfuse (fusions via AGFusion/Arriba), pvacbind (arbitrary peptides), pvacview (manual re-tiering)
• Python: VAtools annotation, LOHHLA/CCF integration, aggregate-report parsing

The Single Most Important Modern Insight -- binding is the easy part; PPV lives downstream

The visible surface of the field — NetMHCpan, MHCflurry, the IC50 column everyone sorts on — is the binding step, and binding is the one step the field has genuinely cracked. The positive predictive value of a binding-only neoantigen pipeline is single-digit percent: of peptides confidently called strong binders, the large majority are never presented, and of those presented, the large majority never elicit a T-cell response (TESLA; Wells 2020). This is structural, not a bad IC50 cutoff — each step of the presentation-and-recognition cascade multiplies a low conditional probability. The corrective: spend the analysis on the filters and features that govern the predicted->presented->immunogenic attrition (clonality/CCF, HLA LOH, expression, agretopicity, foreignness, processing, validation tiers) and treat the choice of binding algorithm as a near-afterthought with sane defaults. TESLA's five features that actually separated immunogenic peptides from binders: HLA binding affinity, source-gene expression ("tumor abundance"), peptide-HLA binding stability, hydrophobicity, and the two recognition features — agretopicity and foreignness.

The pVACtools Suite

Sub-tool	Input that defines the peptide	Use case
pVACseq	VEP-annotated somatic VCF (SNV + indel/frameshift)	The workhorse: point mutations and frameshifts
pVACfuse	AGFusion / Arriba fusion output	Fusion-junction novel-ORF neoantigens
pVACbind	a plain peptide FASTA	Score arbitrary peptides (MS hits, splice peptides); no WT/agretopicity
pVACvector	chosen epitopes	Order epitopes into a vaccine string, minimizing junctional neo-epitopes
pVACview	*.all_epitopes.aggregated.tsv	Human-in-the-loop review and re-tiering (the decision step)

Upstream Chain (every link can silently poison the output)

Step	Tool(s)	Failure if skipped/wrong
Somatic calling (T/N)	Mutect2, Strelka2 (consensus)	Germline leak -> false neoantigens; indels matter most (frameshifts)
VEP annotation	VEP + Wildtype + Frameshift plugins, --fasta, --tsl, --symbol	Most error-prone step; wrong plugins -> no WT peptide / no frameshift ORF
HLA typing (I and II)	OptiType (class I, WES), arcasHLA (RNA), HLA-HD (II)	Wrong allele = confident garbage; type at 4-digit; reconcile DNA vs RNA
Expression	kallisto/salmon TPM, vcf-expression-annotator	--expn-val passes everything if unannotated -> ships unexpressed "neoantigens"
Read counts	bam-readcount, vcf-readcount-annotator	VAF/coverage filters pass everything if unannotated
Phasing	merge somatic+germline, WhatsHap / GATK ReadBackedPhasing	Proximal in-cis variants -> peptides the patient never makes (neoepiscope, Wood 2020)

Decision Tree by Scenario

Scenario	Recommended	Why
SNV + indel neoantigens	pVACseq, all_class_i	The workhorse; frameshifts via Frameshift plugin
Gene fusions	pVACfuse (AGFusion/Arriba, with STAR-Fusion read support)	Junction novel ORFs; demand junction read support
Proximal germline/somatic variants nearby	pVACseq --phased-proximal-variants-vcf	Otherwise the peptide sequence is wrong
Need quality features (DAI, foreignness, dissimilarity)	NeoFox / antigen.garnish on pVAC candidates	pVAC tiers; NeoFox computes the ~16 published features
Reproducible end-to-end	nextNEOpi (HLA + VEP + pVACseq + NeoFox + LOHHLA)	Wires the whole chain including LOHHLA and purity
Final candidate selection	pVACview manual re-tiering	Tiers say WHY a candidate failed; human triage

Run VEP, Then pVACseq

Goal: Produce the VEP annotation pVACseq actually consumes, then call neoantigens.

Approach: Run VEP with the Wildtype + Frameshift plugins and a protein FASTA; annotate expression and read counts with VAtools; supply a phased proximal-variants VCF; then pvacseq run with the patient HLA and sane filters.

pvacseq install_vep_plugin $VEP_PLUGINS          # installs Wildtype + Frameshift
vep --input_file somatic.vcf --output_file somatic.vep.vcf --format vcf --vcf \
    --symbol --terms SO --tsl --hgvs --fasta GRCh38.fa --offline --cache --dir_cache $VEP_CACHE \
    --plugin Frameshift --plugin Wildtype --pick

vcf-expression-annotator somatic.vep.vcf kallisto.tsv custom transcript -s TUMOR \
    --id-column target_id --expression-column tpm -o somatic.vep.expn.vcf

pvacseq run somatic.vep.expn.vcf TUMOR \
    "HLA-A*02:01,HLA-A*24:02,HLA-B*07:02,HLA-B*44:02,HLA-C*07:02,DRB1*01:01" \
    all_class_i pvac_out/ \
    -e1 8,9,10,11 --iedb-install-directory $IEDB \
    --phased-proximal-variants-vcf phased.vcf.gz \
    --normal-vaf 0.02 --tdna-vaf 0.25 --trna-vaf 0.25 --expn-val 1.0 -t 8

Key flags: -e1/-e2 epitope lengths; -b/--binding-threshold (default 500 nM); --percentile-threshold (recommend 2); -m/--top-score-metric median|lowest; --allele-specific-binding-thresholds (preferred over flat 500 nM); --net-chop-method/--netmhc-stab (processing + stability features).

Compute Agretopicity (DAI) Correctly

Goal: Quantify how much more foreign the mutant looks than its wild-type counterpart.

Approach: Agretopicity (the fitness-model amplitude; Łuksza 2017) is the WT/MT binding ratio. A high value means the mutant binds while the WT does not — the surface is new to the immune system, so reactive T cells were not deleted in the thymus. The original differential agretopicity index (DAI; Duan 2014) is the difference form; both forms share the traps below. Requires the matched WT peptide (the Wildtype plugin), so pVACbind cannot compute it.

import pandas as pd

def add_agretopicity(df, wt='Median WT IC50 Score', mt='Median MT IC50 Score'):
    '''Agretopicity (amplitude) = IC50_WT / IC50_MT (ratio > 1 = mutant binds better -> favorable).
    Anchor-position mutations inflate DAI without changing the TCR-facing surface, so
    pair DAI with anchor evaluation rather than trusting it alone.'''
    out = df.copy()
    out['agretopicity'] = out[wt] / out[mt]
    out['dai_favorable'] = out['agretopicity'] > 1
    return out

Drop Candidates on Lost HLA Alleles (LOHHLA)

Goal: Remove neoantigens predicted to be presented by an HLA allele the tumor has deleted.

Approach: HLA LOH is an immune-escape mechanism in ~40% of NSCLC (McGranahan 2017) and is invisible to binding/expression/clonality filters. Run LOHHLA (or a subclonal-sensitive equivalent like DASH) with the HLA type and tumor purity/ploidy, then filter the aggregate report. This step sits outside pVACtools and errors silently if skipped.

def drop_lost_allele_candidates(df, lost_alleles, allele_col='HLA Allele'):
    '''lost_alleles: set of alleles called as LOH-lost by LOHHLA. A peptide assigned
    to a lost allele is not weakly presented - it is not presented at all.'''
    return df[~df[allele_col].isin(set(lost_alleles))].copy()

Per-Method Failure Modes

HLA LOH silent invalidation

Trigger: ranking candidates without running LOHHLA. Mechanism: tumor deletes the haplotype that would present its neoantigens; upstream signals all look fine. Symptom: beautiful candidates on an absent allele. Fix: mandatory separate LOHHLA step; drop lost-allele candidates.

Subclonal mis-tiering from raw VAF

Trigger: using VAF as clonality without purity/CN correction. Mechanism: clonality needs cancer cell fraction (CCF = f(VAF, purity, local CN)). Symptom: clonal mutation in low-purity sample read as subclonal (and vice versa in amplified regions). Fix: estimate purity (ASCAT/Sequenza/PURPLE) and CCF (PyClone) before tiering.

Unphased proximal variants

Trigger: running pVACseq with only the somatic VCF when nearby in-cis variants exist. Mechanism: the translated peptide depends on both variants on the haplotype. Symptom: predicted/synthesized peptides the tumor never makes. Fix: supply --phased-proximal-variants-vcf (merge somatic+germline, phase with WhatsHap/GATK).

Silent filter pass-through

Trigger: expression/VAF/coverage filters set but the values never annotated into the VCF. Mechanism: the filter passes everything when the field is absent. Symptom: unexpressed/low-coverage candidates in the output. Fix: annotate with VAtools first; confirm the FORMAT/INFO fields exist.

Frameshift/fusion over-trust and MS-gap

Trigger: treating frameshift/fusion presentation scores like canonical SNV scores. Mechanism: EL/MS training is dominated by canonical 8-11mers from point mutations. Symptom: narrow-looking CIs on a poorly-supported class; no MS evidence misread as absence. Fix: widen confidence on these high-value classes; build a personalized MS search DB before claiming MS absence.

Quantitative Thresholds

Threshold	Source	Rationale
Binding threshold 500 nM (default)	pVACseq -b default	Entry gate only; prefer --allele-specific-binding-thresholds / %Rank
--normal-vaf 0.02	pVACseq default	Germline-leak guard (esp. tumor-only-ish setups)
--tdna-vaf / --trna-vaf 0.25	pVACseq default	Min tumor DNA/RNA VAF to keep
--expn-val 1.0 TPM	pVACseq default	Unexpressed mutation is not a neoantigen
Coverage normal/tDNA/tRNA 5/10/10	pVACseq defaults	Below this, VAF/clonality calls are noise
Clonal CCF ~1 (clonal >> subclonal)	McGranahan 2016	Subclonal targets select for resistant majority
Agretopicity/DAI > 1 favorable	Łuksza 2017; TESLA	WT binds poorly -> surface not tolerized
Validate beyond tier 1	Wells 2020; Ott/Sahin 2017	Predicted != presented != immunogenic

Common Errors

Error / symptom	Cause	Solution
pVACseq misses frameshifts / no WT peptide	Used Downstream plugin or omitted Wildtype	Install Wildtype + Frameshift (Downstream dropped in pVACtools 2.0)
Everything passes the expression filter	TPM never annotated	vcf-expression-annotator before run
Non-overlapping neoantigen lists across labs	Different HLA typers/resolution	Type at 4-digit; reconcile DNA vs RNA; WES preferred
Candidates on a deleted allele	LOHHLA skipped	Run LOHHLA; drop lost-allele candidates
Wrong mutant peptide sequence	Proximal variants unphased	--phased-proximal-variants-vcf
Subclonal target promoted	Ranked by VAF/IC50, no CCF	Estimate purity + CCF; respect the Subclonal tier

References

• Wells DK, van Buuren MM, Dang KK, et al. 2020. Key parameters of tumor epitope immunogenicity revealed through a consortium approach improve neoantigen prediction (TESLA). Cell 183(3):818-834.
• Hundal J, Kiwala S, McMichael J, et al. 2020. pVACtools: a computational toolkit to identify and visualize cancer neoantigens. Cancer Immunology Research 8(3):409-420.
• McGranahan N, Furness AJS, Rosenthal R, et al. 2016. Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade. Science 351(6280):1463-1469.
• McGranahan N, Rosenthal R, Hiley CT, et al. 2017. Allele-specific HLA loss and immune escape in lung cancer evolution (LOHHLA). Cell 171(6):1259-1271.
• Łuksza M, Riaz N, Makarov V, et al. 2017. A neoantigen fitness model predicts tumour response to checkpoint blockade immunotherapy. Nature 551:517-520.
• Balachandran VP, Łuksza M, Zhao JN, et al. 2017. Identification of unique neoantigen qualities in long-term survivors of pancreatic cancer. Nature 551:512-516.
• Richman LP, Vonderheide RH, Rech AJ. 2019. Neoantigen dissimilarity to the self-proteome predicts immunogenicity and response to immune checkpoint blockade. Cell Systems 9(4):375-382.
• Wood MA, Nguyen A, Struck AJ, et al. 2020. neoepiscope improves neoepitope prediction with multivariant phasing. Bioinformatics 36(3):713-720.
• Lang F, Riesgo-Ferreiro P, Löwer M, Sahin U, Schrörs B. 2021. NeoFox: annotating neoantigen candidates with neoantigen features. Bioinformatics 37(22):4246-4247.
• Ott PA, Hu Z, Keskin DB, et al. 2017. An immunogenic personal neoantigen vaccine for patients with melanoma. Nature 547:217-221.

Related Skills

• immunoinformatics/mhc-binding-prediction - the binding step (the solved, low-leverage part); EL abundance bias bites here
• immunoinformatics/mhc-class-ii-prediction - class II neoantigens for CD4 help (compounded uncertainty)
• immunoinformatics/immunogenicity-scoring - quality ranking (DAI, foreignness, dissimilarity) of the candidate list
• clinical-databases/hla-typing - the genotype substrate; wrong calls poison everything
• clinical-databases/somatic-signatures - clonal neoantigen burden predicts ICI response (McGranahan 2016)
• variant-calling/variant-calling - upstream somatic SNV/indel calls
• workflows/neoantigen-pipeline - the end-to-end orchestration