Consumer guide¶

Prediction, evidence and ranking interfaces for downstream consumers, including Vaxrank. For release history and migrations, see the changelog.

For table-only combination, ORF/RNA evidence, additive re-scoring and the Vaxrank handoff, see Combining sources (5.68.0+).

The shortest useful version¶

from topiary import read_lens, evaluate_scores, resolve_default_methods
from topiary.ranking import parse

df = read_lens(path).to_long().df

scores = evaluate_scores(
    df,
    parse("affinity.value.logistic_normalized(350, 150) * (n_rna_alt > 5)"),
    default_methods=resolve_default_methods(df),
)

Every reader produces the same column vocabulary, so that expression runs unchanged against LENS, pVACseq, or a predictor's own output.

RNA and DNA evidence¶

Columns come in three layers. Filter against layer 1.

Canonical cross-source — same meaning from every reader.
Canonical unit-specific — n_rna_alt_reads / n_rna_alt_fragments, present only where a source reports both units and they differ.
Source-prefixed originals — lens_vaf, pvacseq_tumor_dna_vaf: exactly the number the tool printed, never reinterpreted. Find them with source_columns(df), or source_columns(df, "lens") for one tool.

The same vocabulary across readers, not the same columns — a reader emits one only where its source can answer:

Column	Meaning
`n_rna_alt`	Evidence supporting the variant allele
`n_rna_ref`	Evidence supporting the reference allele
`n_rna_other`	Evidence supporting neither — a third allele, error, or nearby indel
`n_rna_overlapping`	Evidence covering the position (total coverage)
`rna_vaf`	Variant allele fraction
`rna_evidence_subject`	`"reads"` or `"fragments"` — what those counts count
`rna_evidence_method`	How they were obtained
`rna_alt_expression`	Abundance attributed to the variant allele
`rna_alt_expression_method`	How that was obtained
`gene_expression`	Gene-level abundance
`transcript_expression`	Transcript-level abundance, where separately stated
`sequence_source`	How the protein sequence came to exist

DNA support is the same shape, so a filter written against RNA depth becomes the DNA one by changing three letters: n_dna_alt, n_dna_ref, n_dna_other, n_dna_overlapping, dna_vaf, dna_evidence_subject, dna_evidence_method. There is no DNA expression pair — abundance is a transcript property.

n_rna_other / n_dna_other are absent unless the source counted the reference independently. Where topiary derived ref as depth - alt, that figure already absorbs the other alleles, so reporting a difference of zero would assert a clean locus nobody checked.

A column is omitted, not nulled, where the source cannot answer it. A pVACseq aggregated report states a DNA VAF but no DNA depth, so it gets dna_vaf and no n_dna_* at all — a column full of nulls would make available_evidence_columns() report a capability the source lacks.

This holds on every path. LENS emits no rna_alt_expression (it has no DNA VAF to scale abundance by, and its vaf never names its assay), and a prediction run where nobody supplied expression emits no gene_expression. This changes filter behaviour: gene_expression > 1 against such a frame now raises instead of evaluating to NaN and silently emptying your candidate list. Check with available_evidence_columns(df) before naming a column in a config that runs against more than one source.

Check which are present before naming one in a config that has to run against more than one source:

from topiary import available_evidence_columns, EVIDENCE_COLUMNS

missing = set(EVIDENCE_COLUMNS) - set(available_evidence_columns(df))

Naming an absent column in an expression raises rather than evaluating to NaN:

transcript_expression > 1  on a pVACseq aggregated frame
  ValueError: Column 'transcript_expression' not found in DataFrame.
              Did you mean: ['gene_expression', 'rna_alt_expression', ...]

That is the intended behaviour — a silent NaN would drop every row in a filter and say nothing. The alternative, emitting all-null columns everywhere, is worse: a column that is present and empty asserts the question was asked and answered as nothing. Absent beats substituted here as everywhere else.

Concretely, a pVACseq aggregated report has gene-level RNA Expr but no separately stated transcript-level abundance, so it has no transcript_expression; a LENS file without a tpm column has no gene_expression either.

n_rna_alt, rna_evidence_subject and rna_evidence_method are present wherever a source can determine variant support, which is what makes a threshold written against n_rna_alt portable — whatever unit that source counts in. Coverage without a usable fraction can still populate n_rna_overlapping and its subject while correctly omitting the unavailable alt count and derivation.

Read rna_evidence_subject before a number leaves the run. Within one run the unit is consistent and rankings are unaffected. It matters for things that travel — a documented threshold, a config copied between projects, a number in a paper. Five fragments and five reads are different bars.

Where each number came from¶

rna_evidence_method is one of:

Method	Meaning	Provenance
`rna_alignment`	Counted from an RNA alignment	`measured`
`rna_depth_x_vaf`	RNA depth × RNA VAF, rounded	`approximated`
`rna_depth_x_source_vaf`	RNA depth × a VAF whose assay the source did not state	`approximated`
`cds_overlap_reads`	Counted, but of reads overlapping the peptide's CDS	`approximated`
`tpm_x_dna_vaf`	Transcript abundance × DNA VAF	`approximated`
`source_reported`	The source supplied the number without saying how	`approximated`

topiary.provenance_for_method(method) maps a method to its provenance using the definitions above; derived estimates remain approximated.

describe_read_evidence(df) summarises a whole frame without walking rows.

Per source¶

	`n_rna_alt`	subject	method	DNA columns
isovar	30	`fragments`	`rna_alignment`	none — isovar reads RNA
pVACseq (all_epitopes)	429	`reads`	`rna_depth_x_vaf`	yes, from DNA depth × DNA VAF
pVACseq (aggregated)	from `RNA VAF`	`reads`	`rna_depth_x_vaf`	`dna_vaf` only — no DNA depth stated
LENS	from `lens_vaf`	`reads`	`rna_depth_x_source_vaf`	none — see below

LENS's vaf carries no assay qualifier, and lands as lens_vaf. LENS names its read columns rna_* explicitly and leaves vaf bare, so topiary cannot tell whether the fraction is from RNA or DNA. It is used to split the RNA depth — under a method that says so — and not to scale expression or to populate any n_dna_* column, either of which would assert an assay nobody stated.

Multiple samples¶

Evidence is per row, and sample_name is a first-class column and a DSL group key, so a stacked multi-sample frame already carries each sample's own counts and each stays attributable:

merged = stack_results([per_sample_result(name) for name in samples])
merged.df[["sample_name", "n_rna_alt", "n_rna_overlapping", "rna_vaf"]]

Build a separate pooled view with the same canonical count names:

pooled = aggregate_evidence_across_samples(merged.df)
pooled[["n_samples", "n_rna_alt", "n_rna_overlapping", "rna_vaf"]]

The original frame is unchanged, so both the per-sample and pooled answers stay available. Repeated prediction rows within one sample count once. Counts are summed only when every represented sample states them; an unmeasured sample is not silently converted to zero. VAF is recomputed as pooled alternate count / pooled overlapping count, never averaged across samples.

Pooling also requires the samples to agree on evidence subject. Allele-support counts must agree on method as well: a read count cannot be added to a fragment count, and a measured count cannot be flattened together with a derived estimate. Coverage-only depths need no invented allele-derivation method. Expression remains per-sample because it has no generally valid cross-sample sum. n_samples records how many samples actually contain rows for each candidate; a candidate absent from a sample is not treated as a zero row.

Fragments¶

Four sources, one shape:

fragments_from_variants(variants, alignment_file=bam)   # assembled from RNA
fragments_from_variants(variants)                        # translated from reference
fragments_from_effects(effects, padding_around_mutation) # reference arm alone
fragments_from_dataframe(read_lens(path).df)             # from a reader frame
fragment_from_isovar_result(result)                      # an IsovarResult you hold

SEMANTIC_CORE names the fields every source speaks to. They differ only in which they can fill, so one consumer function reads all of them:

def support(fragment):
    return fragment.n_rna_alt          # None where the source has no RNA

isovar is optional in the strong sense — it is not a base requirement, and import topiary does not import it. Install it with pip install 'topiary[isovar]'; only fragments_from_variants with an alignment_file needs it.

The extra requires isovar>=1.39.5,<2, compatible with Osteosarc 0.14 and the current RNA support record and export schemas. This range is enforced at run time as well as at install time: an older release can import cleanly and still return wrong evidence (1.17.x miscounts reads at insertion boundaries), so Topiary refuses it with an upgrade instruction. Upgrade with pip install --upgrade 'topiary[isovar]'.

RNA context for ligands versus vaccine peptides¶

Tell Topiary how long the peptide of interest is. For ligand-only workflows, the default objective is the longest epitope_lengths entry: 11 aa with the default 8–11mer lengths, retaining a 21-aa RNA context target. For long vaccine peptides, set the vaccine size separately from MHC prediction lengths:

fragments = fragments_from_variants(
    variants,
    alignment_file=bam,
    protein_context_peptide_length=25,       # vaccine peptide size
    protein_sequence_preference="balanced",
    min_protein_sequence_support_fraction=0.85,
    min_variant_sequence_coverage=2,
)

Isovar derives a context target of 2*K - 1 residues: 29 for a 15mer, 49 for a 25mer, 59 for a 30mer. This covers all placements around a centered single-residue substitution, not a guarantee for every wider mutation or deletion junction. Available RNA and stop codons can shorten the result. protein_sequence_length=35 explicitly overrides the target, not the peptide size used to evaluate it. protein_sequence_creator=my_creator remains supported unchanged; combining a custom creator with explicit creator options raises an error instead of silently choosing one.

balanced favors useful context while retaining the specified fraction of the best mutant candidate's compatible read-name support. support prioritizes support; context prioritizes context without that relative budget. All three obey the separate per-base RNA read-object floor. Neither fraction nor floor is a confidence probability, total-alt VAF or a claim that every supporting name spans the whole peptide. No RNA sequence is padded from reference and no floor is lowered to reach the requested length.

The new RNA controls require alignment_file. Reference translation and explicit fallback retain their separate padding_around_mutation rule; requesting a longer RNA peptide objective does not change reference padding. Fallback fragments remain reference-derived and carry no invented RNA counts.

RNA fragments record resolved settings in isovar_* annotations, including isovar_version, isovar_protein_context_peptide_length, isovar_protein_sequence_length, the preference and both support controls. These survive fragment TSV/JSON save/reload and become prediction columns. Custom-creator settings are recorded only where the creator exposes them; Topiary does not fabricate missing values. These settings describe the run, not a validation of historical vaccine selection or clinical efficacy.

Reads and fragments¶

isovar reports both units. Both are carried, each under a name that says what it holds:

fragment.n_rna_alt_reads        # 58
fragment.n_rna_alt_fragments    # 30
fragment.n_rna_alt          # 30 — fragments preferred
fragment.rna_evidence_subject()      # "fragments"

Fragments are preferred where a source has them: a paired-end fragment is one molecule read twice, so it is one piece of evidence and two reads. Reads are used where that is all there is. There is no conversion between them — that needs library information no source carries.

Knowing which fields are real¶

fragment.is_known("n_rna_alt_reads")             # populated at all?
fragment.provenance_of("n_rna_alt_reads")        # measured | approximated | synthesized
fragment.is_approximate("n_rna_alt_reads")
fragment.is_usable_as_biology("n_rna_alt_reads") # False for absent *and* synthesized

None means the source could not answer. It is not zero, and the distinction survives a TSV round trip. is_usable_as_biology is the one that matters for correctness: a placeholder the loader invented has a value that means nothing, and anything doing variant-effect annotation on it must refuse rather than compute.

The ranking DSL¶

Resolving ambiguity¶

A frame with two models producing one kind, or one model at two versions, raises rather than guessing. Both have a configured answer:

evaluate_scores(df, node,
                default_methods=resolve_default_methods(df),
                default_versions=resolve_default_versions(df))

describe_default_versions(df) returns the candidates each choice was made between, so a run can report "netMHCpan at 4.1b and 4.2, scoring with 4.2" without re-deriving anything. validate_default_* catches an entry naming a model or version that never ran.

Version ordering is PEP 440 — 4.10 beats 4.9 — and a version that is not PEP 440 sorts before everything that parses, so "newest" means a real release.

Alleles¶

alleles= takes a sequence, a mapping, or a callable:

EvalContext(df, alleles=["HLA-A*02:01", "HLA-B*07:02"])   # every peptide
EvalContext(df, alleles={"SIINFEKLA": ["HLA-A*02:01"]})    # per peptide
EvalContext(df, alleles=lambda keys: genotype_for(keys["peptide"]))

Use a per-peptide form when peptides were not each reported against the whole genotype — a reader emitting one row per (peptide, allele) passing its own threshold produces exactly that. Declaring the union invents a group for every pairing that was never scored, and an expression reading only peptide-level evidence gives each invented group a real number.

Peptide-level evidence¶

A row whose kind is MHC-independent broadcasts across a peptide's allele groups — unless the row names an allele, in which case it lands there and nowhere else. That is the mechanism for allele attribution: writing a row onto chosen alleles is how a policy says "credit this evidence here".

KIND_MHC_DEPENDENCE and mhc_dependence() say which kinds are peptide-level. There are five, not one.

ctx = EvalContext(df, group_keys=gk)
a = evaluate_scores(df, node_a, context=ctx)
b = evaluate_scores(df, node_b, context=ctx)

apply_filter and apply_sort return new frames, so a context cannot be threaded down a pipeline — reuse applies to several operations on one unchanged frame. Passing a stale one raises.

Shared helpers, so you do not reimplement them¶

is_stated(value)        # did the source say anything here?
stated_values(series)   # the same rule over a column
is_named_version(value) # the same rule, named for predictor_version

None, NaN, blank, and the literal strings "nan" / "None" / "<NA>" all mean not stated. The obvious version of this test — if str(v).strip() — excludes only the blank spellings, since str(None) is "None" and str(float("nan")) is "nan", both truthy. That mistake shipped in topiary and independently in a consumer that had reimplemented it.

Also public rather than reimplemented: fragment_from_effect, resolve_default_methods, CANONICAL_METHOD_PREFERENCE, KIND_MHC_DEPENDENCE, derive_mhc_class.

Migrating renamed columns¶

Every column renamed since 5.46.0 is in RENAMED_COLUMNS, and renamed_column(name) looks one up:

from topiary import RENAMED_COLUMNS, renamed_column

renamed_column("vaf")              # 'lens_vaf'
renamed_column("tumor_rna_depth")  # 'pvacseq_tumor_rna_depth'
renamed_column("gene_expression")  # None — not renamed

The DSL consults it before fuzzy matching, so a stale expression says what to do:

Column 'vaf' not found in DataFrame. It was renamed to 'lens_vaf'.

Do not fuzzy-match these yourself. vaf is the trap: the closest surviving name is rna_vaf, and that is the wrong answer — rna_vaf is the canonical cross-source fraction, while vaf became lens_vaf, LENS's own fraction whose assay the file never states.

If you read reader-frame columns with row.get(...) or df[...] rather than through the DSL, topiary cannot warn you — a missed .get() returns None and becomes a silent zero. Check your column names against RENAMED_COLUMNS once at startup. Serialized ProteinFragment JSON and TSV are the exception: their old unit-specific evidence names are migrated on load, including field_provenance keys.

Reader frames have no compatibility aliases. Two output columns for one quantity is the ambiguity the renames existed to remove. ProteinFragment accepts its old unit-specific names when loading, constructing, and reading attributes so the 5.x API remains compatible, but serialization emits only the new names.

Reader escape hatches¶

Topiary automatically recognizes common prediction vocabulary in new pVACtools and LENS columns: EL means presentation, BA / Aff / Affinity mean affinity, IM means immunogenicity, and explicit quantity words such as Processing take precedence over a model suffix. MT/WT score and percentile modifiers are recognized for affinity, processing, presentation, and immunogenicity. Consumers that inspect external headers directly can call parse_prediction_metric(model_name, metric_name) to apply the same rule. For an otherwise unknown pVACtools predictor, an unqualified percentile can inherit the kind of its sole explicit companion, such as an IC50 column. Columns that remain ambiguous are warned about and preserved under their original names rather than discarded.

binding_metrics remains the escape hatch for a LENS column whose meaning cannot be inferred safely:

read_lens(
    path,
    binding_metrics={("newtool", "opaque_metric"): ("affinity", "value")},
)

Merged over the built-in table, keyed on (tool, metric) — the pair the unmapped-column warning prints, and version-free so one entry covers a tool however a file spells its release. None as a value declares a column a non-prediction, silencing the warning without remapping it.

LENS splits <tool>_<version>.<metric> at the final underscore before the version, so tool names such as foo_2 remain intact in foo_2_1.0.MT_Presentation_Score and may be used directly in an override key. Inferred WT metrics normalize to the corresponding _wt_value, _wt_score, or _wt_rank wide column rather than being relabeled as MT. Wide-form round-trips also preserve distinct WT predictor metadata in _wt_method and _wt_version columns; older frames without those columns continue to infer that metadata from the MT model.

Cache compatibility¶

CachedPredictor normalizes older stores on load. Identical repeated keys coalesce, including rows that differ only in source/sample labels. Conflicting values for one prediction key raise on construction and concatenation; conflicting_predictions(df) returns the offending rows. The key includes flanks and allele_set, so context-dependent predictions and distinct genotypes remain separate. See Cached predictions for the complete contract.