bio-immunoinformatics-neoantigen-prediction

What this skill does

## 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.

```bash
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.

```python
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