bio-immunoinformatics-mhc-binding-prediction

What this skill does

## Version Compatibility

Reference examples tested with: MHCflurry 2.1+, 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: MHCflurry 2.2.0+ switched its backend from TensorFlow to PyTorch (Python 3.10+) — confirm `mhcflurry-downloads fetch` succeeded and the backend imports before scoring. NetMHCpan-4.1 and MixMHCpred are standalone academic binaries, not pip-installable; the IEDB REST API wraps NetMHCpan if a local install is unavailable. Tool versions move fast — re-verify supported-allele lists and default %Rank thresholds against current docs.

# MHC Binding Prediction

**"Predict which peptides bind/are presented by MHC class I"** -> Score peptide-HLA class I binding affinity and natural-presentation likelihood to nominate candidate CD8 epitopes.
- Python: `mhcflurry.Class1PresentationPredictor.load().predict()` (pip-installable, forgiving allele parser)
- CLI: `netMHCpan` (field default; EL score by default, `-BA` adds affinity) or `MixMHCpred` (MS-deconvolution, EL-only)

## The Single Most Important Modern Insight -- a strong predicted binder is a candidate for the next experiment, not an epitope

Binding to MHC is necessary but nowhere near sufficient for immunogenicity. The real path is a funnel: expression -> proteasomal processing -> TAP transport and loading -> stable surface display -> a cognate T cell that survived thymic selection and activates. Binding prediction addresses essentially one stage. Each downstream stage discards a large fraction of binders, so the precision of "predicted binder -> validated epitope" is low even when the binding model itself is excellent. Two operational corollaries follow. First, never report a presentation score to a collaborator as an "immunogenicity" or "epitope" probability — that is a different, far weaker prediction (immunoinformatics/immunogenicity-scoring). Second, the modern EL/MS models that now define the field learned natural presentation from mass-spec immunopeptidomes, which over-represent peptides from highly expressed proteins; the model therefore partly learns "comes from an abundant protein" as a proxy for "is presented." That bias is exactly backwards for neoantigen discovery, where the targets are mutated and often lowly expressed, living in the under-detected tail the model systematically under-ranks.

## Tool Taxonomy (Class I)

| Tool | Citation | Score type | Form | Use when |
|------|----------|-----------|------|----------|
| NetMHCpan-4.1 | Reynisson 2020 | EL (default) + BA (`-BA`) | standalone/web/IEDB | Field default; broadest allele coverage; presentation discovery |
| MHCflurry 2.0 | O'Donnell 2020 | BA + processing + presentation | pip Python | Scripting, messy allele strings, integrated presentation score |
| MixMHCpred 3.0 | Tadros 2025 | EL only (MS motifs) | standalone | MS-grounded presentation; cross-allele/species extrapolation study |
| NetMHC-4.0 | Andreatta 2016 | BA only | standalone/web | Legacy reproducibility; allele-specific, data-rich common alleles |
| MHCnuggets | Shao 2020 | BA (IC50) | pip Python | High-throughput TCGA-scale screens; rare-allele transfer learning |

## BA vs EL -- the conceptual axis that determines which score to read

BA (binding affinity) models train on in-vitro competitive-binding IC50 assays and measure only whether the groove can hold the peptide thermodynamically. EL (eluted-ligand / presentation) models train on mass-spec immunopeptidomics — peptides actually eluted from MHC on real cells — so the label implicitly folds in processing, transport, editing, and surface stability. Read BA for "could this bind the groove if delivered there"; read EL/presentation for "is this likely naturally presented," which is the default and recommended output of NetMHCpan-4.1, MHCflurry's presentation predictor, and MixMHCpred. IEDB codifies the split: `netmhcpan_ba` = recommended-binding, `netmhcpan_el` = recommended-epitope. Pick by intent.

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Quick scriptable scan, messy allele names | MHCflurry presentation predictor | pip-installable, normalizes `A*02:01`/`A0201`/`HLA-A0201` |
| Maximal accuracy / broadest alleles | NetMHCpan-4.1 EL | Field default; concurrent MS motif deconvolution |
| Neoantigen candidate scoring | EL/presentation + expression check | EL under-ranks low-expression mutants (abundance bias) |
| "Can it physically bind" (engineered/delivered peptide) | BA mode (`-BA`, MHCflurry affinity) | Question is thermodynamic, not presentation |
| Rare / non-European allele | NetMHCpan-4.1, but verify training support | Pan-models extrapolate; confidence drops off the manifold |
| Cross-allele ranking in a multi-HLA patient | %Rank, never raw nM | nM scales differ per allele; nM cutoffs are allele-biased |

## Predict Presentation with MHCflurry

**Goal:** Score peptides against a patient genotype and report the best-presenting allele per peptide.

**Approach:** Load the presentation predictor; pass `alleles` as a sample->genotype dict so the model reports `best_allele`, `affinity` (nM), `affinity_percentile` (%Rank), and `presentation_score` (0-1, higher = more likely presented). Supply real `n_flanks`/`c_flanks` only if the genomic context is known.

```python
from mhcflurry import Class1PresentationPredictor

predictor = Class1PresentationPredictor.load()
df = predictor.predict(
    peptides=['SIINFEKL', 'GILGFVFTL', 'NLVPMVATV'],
    alleles={'patient1': ['HLA-A*02:01', 'HLA-A*24:02', 'HLA-B*07:02']},
    include_affinity_percentile=True,   # required: %Rank column is off by default
    verbose=0,
)
# columns: peptide, sample_name, affinity, best_allele, processing_score,
#          presentation_score, and affinity_percentile (only with the flag above)
# affinity nM: LOWER is stronger. presentation_score: HIGHER is more likely presented.
```

## Interpret with %Rank, Not Raw nM

**Goal:** Classify binding strength in a way that is comparable across alleles.

**Approach:** Threshold on %Rank (percentile of the score against random peptides for that same allele), not on absolute IC50. The 500 nM convention is allele-biased — it over-calls permissive alleles and under-calls restrictive ones, skewing a multi-HLA patient's epitope list toward a subset of the genotype.

```python
def classify_by_percentile(affinity_percentile):
    '''Class I %Rank cutoffs (NetMHCpan convention). LOWER percentile = stronger.
    Strong binder <= 0.5%; weak binder <= 2.0%. Use %Rank for any cross-allele
    comparison; raw nM is only meaningful within a single allele.'''
    if affinity_percentile <= 0.5:
        return 'strong'
    elif affinity_percentile <= 2.0:
        return 'weak'
    return 'non-binder'
```

## Scan a Protein for Class I Epitopes

**Goal:** Enumerate candidate epitopes across a protein for a patient genotype.

**Approach:** Tile 8-11mers (9mers dominate real ligands), score all windows in one batched call, keep windows under the 2% weak-binder cutoff. See examples/mhc_binding.py for the full tiling-and-rank script.

```python
def scan_protein(protein_seq, genotype, lengths=(8, 9, 10, 11)):
    from mhcflurry import Class1PresentationPredictor
    predictor = Class1PresentationPredictor.load()
    peptides = [protein_seq[i:i + k] for k in lengths for i in range(len(protein_seq) - k + 1)]
    df = predictor.predict(peptides=peptides, alleles={'patient': list(genotype)},
                           include_affinity_percentile=True, verbose=0)
    return df[df['affinity_percentile'] <= 2.0].sort_values('affinity_percentile')
```

## Per-Method Failure Modes

### Pa