Protein Fragments¶

Identity and RNA outcomes¶

A record is one observation of a candidate: the pair (sample_name, fragment_id). fragment_id names the candidate, and sample_name names which sample (or library, or analysis policy) observed it. The same candidate seen in two samples is two records sharing one ID. They may differ in evidence, but not in the candidate itself (CANDIDATE_FIELDS: sequence, target intervals, reference and germline sequences).

fragments_from_dataframe derives the ID from everything it groups rows by: source, variant, transcript (transcript_id or pVACseq's transcript), reported peptide and sequence. A frame's sample_name column fills the field. A LENS context reported for several peptides becomes one fragment per peptide, each with that peptide's evidence and annotations["reported_peptide"]. Rows describing the same fragment that disagree about its evidence raise.
fragments_from_variants(..., sample_name="T1") and fragments_for_sample(fragments, "T1") label observations from any source. The label is kept exactly and is checked before any alignment is read.
Prediction scans each candidate once and writes one row per observation. sample_name comes from the fragment, never from a model or cache (blank when unlabelled), so aggregate_evidence_across_samples pools a candidate across samples with its default keys.

unique_fragments(records) keeps the first record per (sample_name, fragment_id) and coalesces repeats that hold the same content. It raises, naming the fields that differ, on conflicting content or on one ID naming different candidates. Content is compared by passing it through fragment IO itself, so a record and its own write_fragments round trip always agree (5 and 5.0, NaN, blank or "None" text and None, tuples and lists, numbers in text fields). A value IO cannot store cannot be compared, and fails explicitly. Prediction runs this check before any model, preventing silent last-record overwrites.

describe_isovar_result(result) returns a JSON-compatible diagnostic record for every completed upstream result, including empty or filtered reconstructions. It preserves separate read/template counts, unknowns as null, sequence/edit interval, transcript IDs and named failed filters. passing means the supplied result has a protein and passes its recorded filters, not clinical eligibility or proof that default filters were used. Fragment adapters remain diagnostic: do not promote a reconstructed but filtered sequence into accepted inputs. Acquisition errors must be reported separately, never as biological negatives.

ProteinFragment is a universal record for a protein / peptide sequence with source-type, target-region, and comparator metadata. It's the substrate that lets Topiary handle antigens from any origin — somatic variants, structural variants, ERVs, CTAs, viral proteins, allergens, autoantigens, synthetic constructs — through one pipeline, and threads identity through predictions so downstream tools (vaxrank, vaccine-window selection) can group peptides back to their source.

The dataclass¶

from topiary import ProteinFragment

f = ProteinFragment.from_variant(
    sequence="MAAVTDVGMAVATGSWDSFLKIWN",
    reference_sequence="MAAVTDVGMAAATGSWDSFLKIWN",   # WT protein
    mutation_start=10, mutation_end=11, inframe=True,
    variant="chr7:140753336",
    effect="p.Val600Glu",
    gene="BRAF",
    annotations={"vaf": 0.42, "ccf": 0.9},
)

Every field except fragment_id and sequence is optional. target_intervals carries the half-open regions within sequence considered targetable — its meaning depends on source type (see vocabulary below).

Running predictions¶

from topiary import TopiaryPredictor

predictor = TopiaryPredictor(models=[...], alleles=[...])
df = predictor.predict_from_fragments(fragments)

Output DataFrame columns, beyond the standard prediction fields:

Column	Meaning
`fragment_id`	Source identity — threads back to `ProteinFragment.fragment_id`
`source_type`, `variant`, `effect`, `effect_type`, `gene`, `gene_id`, `transcript_id`, `transcript_name`	Propagated from the fragment
`gene_expression`, `transcript_expression`	Propagated from the fragment
`overlaps_target`	`True` / `False` / NaN — whether the peptide overlaps any of the fragment's target intervals
`contains_mutant_residues`	Backwards-compat alias — `True` iff `source_type.startswith("variant")` AND `overlaps_target` is True
`wt_peptide`, `wt_peptide_length`	Derived by slicing `effective_baseline` at the peptide's offset. Only populated for substitution-compatible fragments (baseline and `sequence` the same length); `None` otherwise or when no baseline exists.
`wt_value`, `wt_score`, `wt_affinity`, `wt_percentile_rank`, `wt_prediction_method_name`, `wt_predictor_version`	Populated when `TopiaryPredictor(predict_wt=True)` scores non-null `wt_peptide` values with the configured MHC model(s). Rows without a length-compatible WT peptide keep NaN values.
(each annotation key)	Flattened from every fragment's `annotations` dict. Underscore-prefixed keys (e.g. `_subsequence_offset`) are reserved for internal plumbing and never surface as columns.

Building a fragment from a variant effect¶

fragment_from_effect(effect, padding_around_mutation) is the varcode arm of the multi-source story — the same ProteinFragment a LENS or pVACseq reader produces, built from a translated variant effect instead:

from topiary import fragment_from_effect

fragment = fragment_from_effect(effect, padding_around_mutation=8)

Two rules worth knowing before you rely on it:

The window is clipped at the protein's first stop codon, so a mutation near the end yields a shorter fragment rather than a padded one. If the stop falls before the reported mutation, there is nothing to present and the function returns None.
reference_sequence is populated only when the pre- and post-mutation proteins are the same length. An indel or frameshift shifts every downstream residue, so the window sliced at the same offsets would be a different piece of protein presented as a comparator — the same restriction wt_peptide applies.

It returns None when the effect has no mutant protein sequence at all (silent, non-coding, untranslatable). That is an absence, not an error.

Building fragments from variants, with or without RNA¶

from topiary import fragments_from_variants

# RNA available: the context around each mutation is assembled from reads
fragments = fragments_from_variants(variants, alignment_file=bam)

# No RNA: the same variants translated from the reference
fragments = fragments_from_variants(variants)

The two arms are interchangeable. Both return ProteinFragments with the same core, so the rest of a pipeline does not change shape when the RNA does or does not exist. What changes is what the fragments can tell you: an assembled sequence carries the patient's other variants and whatever phasing the reads support, and comes with counted read support; a translated one carries the reference everywhere except the variant itself, and no read counts at all. annotations["sequence_source"] says which.

protein_context_peptide_length sets the RNA reconstruction objective; it defaults to the longest requested epitope_lengths (11 aa, targeting 21 aa of context). Long vaccine workflows can request a separate peptide size. protein_sequence_length overrides the context target directly, while padding_around_mutation controls reference translation independently. See the consumer guide for selection and support controls.

allow_reference_fallback=True translates variants isovar could not support rather than dropping them. The fragments stay distinguishable by sequence_source, which is the reason to record it rather than blend an RNA-backed candidate with an inferred one.

isovar is needed only when alignment_file is given. fragments_from_effects is the reference arm on its own, public because a caller with variants and no alignment file wants exactly that. Reference-only and configuration tests run without Isovar; marked integration tests exercise real RNA reconstruction with both the minimum supported and latest Isovar releases.

Every source reaches a fragment¶

Four paths, one shape. They differ only in which fields they can fill:

Source	Builder	Sequence is	Read counts
isovar	`fragment_from_isovar_result(result)`	assembled from RNA reads	counted (`measured`)
varcode	`fragment_from_effect(effect, padding)`	translated from the reference	none
LENS	`fragments_from_dataframe(read_lens(p).df)`	each reported peptide's context	overlapping counted; support is CDS-overlap (`approximated`)
pVACseq	`fragments_from_dataframe(read_pvacseq(p).df)`	the peptide itself	depth × VAF (`approximated`)

def rna_support(fragment):
    subject = fragment.rna_evidence_subject()
    if subject is None:
        return None
    field = f"n_rna_alt_{subject}"
    if not fragment.is_usable_as_biology(field):
        return None
    return fragment.n_rna_alt, fragment.is_approximate(field)

That function reads all four without knowing which it has: isovar returns (30, False), pVACseq returns a count with True, varcode returns None because it has no RNA data at all. No branching on source_type — which is the whole point of the abstraction, and why source_type stays biological.

SEMANTIC_CORE names the fields every source is expected to speak to, whether or not it can populate them.

isovar is optional in the strong sense¶

Not imported at module scope and not in Topiary's base requirements, so import topiary does not import it. Install the optional integration with pip install 'topiary[isovar]'. Only fragments_from_variants with an alignment_file imports it. fragment_from_isovar_result needs nothing from isovar at all — it reads an already-built result by attribute. A consumer that only reads LENS reports should not pay for a package it never calls.

isovar is also the only source that counts the reads supporting an assembled sequence. Everything else derives them or counts something adjacent, which is what the derivation names record.

One vocabulary across readers¶

Both readers describe expression the same way, so a filter does not have to know which produced the frame:

Column	Meaning
`gene_expression`	Gene-level abundance. On LENS also available as `gene_tpm` (its native spelling), with the original string in `gene_tpm_raw`
`transcript_expression`	Transcript-level abundance, where the source has it
`variant_allele_expression`	Abundance attributed to the variant allele — an estimate; see below
`*_method`	How each of those was obtained

RNA evidence, and knowing which fields are real¶

Beyond gene_expression / transcript_expression, a fragment can carry read-level evidence:

Field	Meaning
`n_rna_overlapping_reads` / `n_rna_overlapping_fragments`	Reads or fragments spanning the variant position
`n_rna_alt_reads` / `n_rna_alt_fragments`	Reads or fragments supporting the variant allele
`n_rna_ref_reads` / `n_rna_ref_fragments`	Reads or fragments supporting the reference allele
`n_rna_other_reads` / `n_rna_other_fragments`	Reads or fragments supporting neither reference nor alt
`n_rna_alt_reads_supporting_protein_sequence` / `n_rna_alt_fragments_supporting_protein_sequence`	Reads or fragments supporting this assembled protein sequence, not merely the allele

These are not derivable from a TPM, which is why they are fields rather than annotations.

None is unknown; it is not 0. A source with no read data leaves them None. A source that looked and found no support sets 0. A consumer must be able to tell those apart, so ask fragment.is_known("n_rna_alt_reads") rather than testing truthiness — and the distinction survives writing to and reading from a TSV.

Different sources populate different subsets, and some of them estimate. So a fragment can also say how real each populated field is:

from topiary import ProteinFragment, APPROXIMATED, SYNTHESIZED

ProteinFragment(
    fragment_id="...",
    sequence="...",
    variant="chr1:100:N>N",       # invented by the loader; not real alleles
    n_rna_alt_reads=12,           # reconstructed as depth x VAF, not counted
    field_provenance={
        "variant": SYNTHESIZED,
        "n_rna_alt_reads": APPROXIMATED,
    },
)

Provenance	Meaning
`"measured"`	Observed directly from data
`"approximated"`	Derived or estimated — e.g. read counts as depth × VAF
`"synthesized"`	A placeholder the loader invented because the source supplied none

A field not named in the mapping is unqualified: it means what it says.

Read it through the accessors rather than the dict:

Call	Answers
`is_known(name)`	Does the field carry a value at all?
`provenance_of(name)`	How real is it, or `None` if unqualified
`is_approximate(name)`	Was it estimated rather than observed?
`is_usable_as_biology(name)`	May it be interpreted as a fact about the sample?

is_usable_as_biology is the one that matters for correctness: it is False for an absent field and for a synthesized one. Anything that annotates variant effects, reannotates transcripts, or otherwise treats a value as biology must check it and refuse rather than compute on a placeholder.

The point of all this is that a consumer never branches on source_type. Every source produces the same shape; they differ only in which fields are populated and how real those fields are:

def rna_support(fragment):
    subject = fragment.rna_evidence_subject()
    if subject is None:
        return None
    field = f"n_rna_alt_{subject}"
    if not fragment.is_usable_as_biology(field):
        return None
    return fragment.n_rna_alt

source_type vocabulary (recommended, not enforced)¶

Free-form string. Topiary never interprets it; used for display and DSL filtering. Colon subtyping is convention.

Category	Values
Variant, small	`variant:snv`, `variant:indel`, `variant:frameshift`, `variant:stop_gain`, `variant:stop_loss`, `variant:start_loss`, `variant:exon_loss`, `variant:alternate_start`
Structural variant	`sv:fusion`, `sv:tandem_duplication`, `sv:inversion`, `sv:translocation`, `sv:cryptic_exon`, `sv:large_insertion`, `sv:large_deletion`
Aberrant expression	`erv`, `cta`, `tumor_overexpressed`, `intron_retention`, `utr`, `novel_orf`
Pathogen	`viral`, `viral:hpv16`, `viral:hiv`, `bacterial`, `parasitic`
Environmental	`allergen`, `allergen:plant`, `allergen:food`, `allergen:mold`, `allergen:dander`
Self / autoimmunity	`self`, `autoantigen`, `autoantigen:myelin`
Synthetic	`synthetic`, `designed`

Producers are free to invent new subtypes.

target_intervals — geometry per source type¶

Reader frames retain recoverable target geometry through fragments_from_dataframe. pVACseq's reported Pos / Mutation Position is normalized into JSON-encoded, zero-based half-open mutation_intervals_in_peptide; disjoint changed positions stay disjoint. The legacy start/end columns are populated only for a contiguous interval. Missing or unsupported positions remain unknown. A historical representative position does not imply the complete altered tail of a frameshift.

When a fragment uses longer context, peptide-relative intervals shift through an exact, unique peptide occurrence. Repeated or absent matches remain unknown; peptide_offset is not assumed to be an offset in that context. A negative reported peptide does not mark the surrounding context negative. Inspect annotations["target_interval_status"] for the mapping outcome.

LENS mut_aa_pos remains unresolved because its coordinate meaning varies by antigen category. Its raw value is preserved as reported_mut_aa_pos, with pipeline/file source and biological antigen_source / source_type retained separately. Unknown coordinates never turn into whole-sequence targets.

only_novel_epitopes=True keeps windows overlapping known target intervals, including junctions, and excludes unknown geometry. overlaps_target is the source-independent selection field; contains_mutant_residues remains a variant-only compatibility field. Neither asserts tumor specificity. For unresolved reports, ordinary rescanning and exact reported-peptide rescoring remain available; no context scan is required for table-only ranking.

Reported WT peptides fill reference_sequence only when the fragment is that reported peptide. They are retained as reported_wt_peptide for longer contexts; Topiary does not invent WT flanks or claim a patient-specific germline baseline. Explicit reference_sequence and germline_sequence columns are preserved.

The public helpers mutation_intervals_from_positions and map_peptide_intervals provide the same normalization and mapping for other report adapters. Synthetic report builders for the regression tests are kept in tests/report_geometry_helpers.py.

The producer computes target_intervals; Topiary never interprets. Meaning varies by source type:

source_type	`target_intervals`
`variant:snv` at position k	`[(k, k+1)]`
`variant:indel` (in-frame insertion) at k, length L	`[(k, k+L)]`
`variant:indel` (in-frame deletion) at k	`[(k, k)]` — the junction where formerly-distant residues now sit together
`variant:frameshift` at k	`[(k, len(sequence))]` — everything downstream is novel (sequence should be truncated at the new stop)
`sv:fusion` (in-frame, coding-coding) with junction at k	`[(k-1, k+1)]` — junction residues only; internal partner residues are self
`sv:fusion` onto non-coding partner with junction at k	`[(k, len(sequence))]` — readthrough translation is all novel
`sv:tandem_duplication` with breakpoints at k1, k2	`[(k1, k1+1), (k2, k2+1)]` — breakpoints only; duplicated bulk is self
`sv:inversion` within coding region [a, b]	`[(a, b)]` — reversed translation is entirely novel
`sv:cryptic_exon` (in-frame inclusion) at [a, b]	`[(a, b)]`
`sv:cryptic_exon` (frameshift inclusion) at [a, b]	`[(a, len(sequence))]`
`erv`, `cta`	Producer-computed non-self regions (based on the producer's definition of "self" — healthy-tissue expression, homology to non-CTA proteins). `None` when the producer can't decide.
`viral`, `allergen`	Immunodominant / IgE-reactive hotspots if known; `None` otherwise

None means unspecified — downstream tools decide whether to treat as "whole sequence." Empty list [] explicitly means "nothing targetable."

Reference vs germline¶

reference_sequence: str | None      # canonical (Ensembl, RefSeq, reference strain)
germline_sequence: str | None       # patient / strain-specific baseline

The DSL's wt.* scope reads effective_baseline:

@property
def effective_baseline(self) -> str | None:
    return self.germline_sequence if self.germline_sequence is not None else self.reference_sequence

Germline takes precedence when populated; reference is the fallback. Both None → wt.* returns NaN.

source_type	typical `reference_sequence`	typical `germline_sequence`
`variant:*`	Canonical WT from varcode/Ensembl	Patient's non-tumor protein if available
`sv:*`	Usually `None`	`None`
`viral[:strain]`	Reference-strain protein	`None` (patient has no germline virus)
`erv`	`None`	`None`
`cta`	Canonical protein (equals `sequence`)	Same as reference (CTAs are non-neoantigens)
`autoantigen`	Canonical (UniProt MBP etc.)	Patient-specific with SNPs — can matter for TCR specificity
`allergen`	Canonical isoform	`None` (patient doesn't have it)
`synthetic`	Natural parent if any, else `None`	`None`

Reserved DSL scope: `self_nearest`¶

For cross-reactivity filtering — "what's the closest peptide in essential healthy tissues, and does it also bind this MHC?"

Topiary computes these when you give it a reference proteome: TopiaryPredictor(self_proteome=..., predict_self_nearest=True) fills the similarity columns from SelfProteome.nearest() and then scores each self_nearest_peptide at its row's own allele, filling self_nearest_value / _score / _percentile_rank. That second half is the one a cross-reactivity judgement turns on — a near-identical self peptide the patient's MHC never presents is not the same risk as one it does.

The self peptide is scored without flanking context: it comes from the reference proteome, and nearest() reports its gene, transcript and offset but not the residues either side. Kinds that read flanks — antigen processing, and presentation where its model uses them — are therefore scored on the peptide alone; affinity and stability are unaffected.

Without predict_self_nearest, or for producers populating externally (via BLAST / edit distance against a healthy-tissue proteome with their own definition of "self"). The DSL scope just reads self_nearest_* columns. When columns are absent, self_nearest.* returns NaN.

Reserved column namespace:

Column	Meaning
`self_nearest_peptide`	Closest healthy-tissue peptide at the same length
`self_nearest_peptide_length`	(Trivially same as the mutant)
`self_nearest_edit_distance`	Producer-chosen distance metric (Hamming / Levenshtein / BLAST score)
`self_nearest_gene`	Source gene of the nearest-self hit
`self_nearest_gene_id`, `self_nearest_transcript_id`	Provenance
`self_nearest_tissues`	Which healthy tissues the source gene is expressed in
`self_nearest_value`, `self_nearest_score`, `self_nearest_percentile_rank`	MHC binding of the nearest-self peptide, paired to the same allele

from topiary import Affinity, Column, apply_filter, self_nearest

# Drop neoepitopes too similar to healthy-tissue self
df = apply_filter(
    df,
    (Affinity.score >= 0.5) & (Column("self_nearest_edit_distance") >= 3),
)

# Ranking that penalizes cross-reactivity
ranking = Affinity.score - 0.5 * self_nearest.Affinity.score

IO¶

from topiary import read_fragments, write_fragments

write_fragments(fragments, "fragments.tsv")
loaded = read_fragments("fragments.tsv")

TSV format: one row per fragment. Scalar fields map to same-named columns. target_intervals, field_provenance, and annotations are JSON-encoded in their own columns. Missing columns on read fall back to field defaults; unknown columns raise. Evidence names from 5.47 and earlier are migrated through RENAMED_COLUMNS, including direct constructor keywords, compatibility attribute reads, and keys inside field_provenance. New dictionaries and TSVs always emit only the assay-scoped names.

For single-fragment / API use: fragment.to_dict(), fragment.to_json(), and the from_dict / from_json classmethods.

NumPy scalars work like their Python equivalents: np.bool_(True) becomes True, integer scalars become int, and ordinary floating scalars become float. String-backed enums preserve their underlying text: a member whose value is "variant:snv" stays "variant:snv", not an enum display name. Construction and JSON/TSV serialization apply the same conversion, including inside nested annotations and to annotations added later. Callers' input data is not modified. Strings such as "False" stay strings, and unsupported objects still raise serialization errors. The shared conversion is also available as topiary.normalize_python_types(value).

Normalization is not a universal encoder. Dates, durations, decimals, fractions, extended-precision floats, complex numbers and arrays remain intact rather than losing units or precision. Use fragment.to_json(default=your_encoder) to choose an explicit representation for these values. Non-JSON prediction-file metadata retains its existing text fallback; reading that fallback returns text, not an automatically reconstructed date or duration.

Nested dataclass annotations still serialize as field dictionaries. This conversion shares the same traversal (dataclasses_as_dict=True) and does not invoke custom deep-copy hooks that could alter values before encoding.

Identity¶

fragment_id names the candidate; (sample_name, fragment_id) names the record. Two fragments with the same pair are equal and hash-equal, regardless of other content (see Identity and RNA outcomes). Use make_fragment_id(prefix, sequence, variant=..., qualifiers=...) for a deterministic content-derived id with a readable prefix:

BRAF_p.Val600Glu__a1b2c3d4
EWSR1-FLI1_fusion__3c8e4b91
erv_Hsap38.chr7.64991215__7f2e89a1
HPV16_E6__5f6a1c23
__4f9c2a8e                  # no metadata → hash-only fallback

Prefix is sanitized to [A-Za-z0-9._:-]; runs of other characters collapse to _. Hash is 8 hex chars of SHA-1 over sequence, the optional variant and any qualifiers — the other values a producer groups records by, such as a reported peptide. Without qualifiers the id is unchanged from earlier releases.

What's not in this release¶

Indel / frameshift wt_peptide — wt_peptide is only populated when the baseline is the same length as the mutant sequence (substitution-compatible). Length-changing edits yield None. Whether to define a remapped comparator remains an open design question (#411).
Binding-aware cross-reactivity axes (self_mimic_*, self_strongest_nearby_*, self_nearest_candidates) — the sequence-level self_nearest_* columns are computed by SelfProteome / predict_self_nearest; the axes that also require MHC prediction on each candidate are not implemented (#412).
Dedicated Isovar fragment-file loader (read_isovar_fragments) — separate work. Existing in-memory results compose through fragment_from_isovar_result. Native Exacto now has read_exacto_fragments; pVACseq composes through fragments_from_dataframe(read_pvacseq(...).df).