Skip to content

Add allele_set column to faithfully store haplotype-mode predictions #168

Description

@iskandr

Problem

Topiary's long-form schema stores a single allele string per row, which assumes mhc_dependence == "single_allele". MHCflurry's presentation predictor can run in haplotype mode (presentation_allele_mode="haplotype"), where each peptide gets one presentation score against the whole set of alleles, not against each allele individually.

Today's workaround is to store the "best-binding allele within the haplotype" in the allele column for haplotype rows. This is lossy:

  • Two haplotype predictions of the same peptide against the same allele set can land on different "best alleles" depending on which one wins — so identity becomes non-deterministic.
  • A reader of a single row can't tell whether allele=HLA-B*07:02 means "this prediction is against HLA-B07:02 specifically" or "this prediction is against a haplotype, and B07:02 happens to be the best-binding member."
  • Downstream operations (cross-source validation, side-by-side comparison, round-trip through to_tsv/read_tsv) all silently collapse the haplotype meaning.

Proposal

Add an allele_set column to topiary's long-form schema.

Row type allele allele_set
Single-allele binding (mhc_dependence="single_allele") the allele null
Haplotype binding (mhc_dependence="haplotype") best-binding allele within the set (convenience view) canonical set of all alleles in the haplotype
Non-MHC (mhc_dependence="none") null null

Storage format

Comma-joined, sorted allele names: "HLA-A*02:01,HLA-B*07:02,HLA-C*06:02". Pros: round-trips through TSV cleanly, greppable, ordering is canonical (so two haplotype rows over the same set hash equal). Cons: parsing on read, but alleles can't contain commas in any naming convention I know of, so it's unambiguous.

Alternatives considered: frozenset[str] (pythonic but doesn't round-trip), tuple[str, ...] sorted (same problem), JSON-encoded list (verbose, ugly in TSV).

Schema contract

When kind_support[method][kind]["mhc_dependence"] == "haplotype", rows of that (method, kind) MUST populate allele_set. Reader/writer can validate this on the I/O boundary.

Downstream impact

This unlocks several follow-ups that are currently impossible to do cleanly:

  • Identity for cross-source validation. Two haplotype pMHC_presentation predictions of the same peptide against the same allele set can be matched faithfully (via (peptide, allele_set, n_flank, c_flank)) instead of via the lossy (peptide, allele=best, n_flank, c_flank) workaround.
  • Per-mhc_dependence peptide-level projection in the DSL. See companion issue. BestAlleleField exists for mhc_dependence="single_allele". Haplotype and non-MHC rows already-per-peptide-rows fall out naturally once allele_set is in place to disambiguate haplotype rows from single-allele rows that happen to share the allele field.
  • Faithful MHCflurry round-trip. CachedPredictor.from_mhcflurry (and any reader of MHCflurry presentation output) can preserve the haplotype shape end-to-end.
  • Haplotype-aware DSL queries like Presentation.value evaluated on haplotype rows — the row IS the per-peptide value with no aggregation needed.

Ripple effects (in scope for the PR that implements this)

  1. topiary/predictor.py — MHCflurry presentation path needs to populate allele_set when presentation_allele_mode="haplotype".
  2. topiary/cached.pyfrom_mhcflurry and any other reader of haplotype output needs to preserve allele_set.
  3. topiary/result.py / topiary/io.pyto_tsv / read_tsv need to write and parse the column without mangling it; the round-trip test suite needs a haplotype case.
  4. topiary/wide.pyto_wide / from_wide need to treat allele_set as a group-identifying column for haplotype rows (not pivot it).
  5. Existing BestAlleleField and DSL allele-aware nodes need either to skip mhc_dependence="haplotype" rows or to error helpfully (today they implicitly run on allele regardless).
  6. Tests: a haplotype-mode MHCflurry fixture in tests/data/, round-trip tests, validation tests.

Why this is its own PR (not bundled with anything else)

Schema additions ripple through the entire codebase. The pVACseq loader PR (#167) works fine without it because pVACseq's output is always single-allele — but the validator / combine-predictions helpers we sketched on top of #167 can't be done correctly until allele_set exists. This needs to land first; the helpers come after.

See also

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions