API reference

Auto-generated from the source docstrings via mkdocstrings.

The mkdocs i18n plugin builds separate API reference pages per locale. To avoid anchor conflicts across locales, only the most-used core modules are auto-rendered below. Per-tool Protocol adapters and specific helpers are documented in the relevant tool page (codon, neoantigen, trial, spatial, etc.).

Codon optimization

mrnavax.codon_optimizer

Codon-usage + GC + rare-codon analyzer for mRNA CDS sequences.

Implements the classical features that any modern codon-optimization model (CodonBERT, RiboDecode, LinearDesign, mRNABERT) consumes as inputs:

• CAI (Codon Adaptation Index) — Sharp & Li (1987). Reference table embedded below for Homo sapiens.
• GC% — fraction of G+C in the CDS (excluding stop).
• Rare-codon fraction — fraction of codons whose usage frequency is below 0.10 in the reference table (heuristic for ribosome stalling).
• Dinucleotide CpG/ObsExp — proxy for innate immune recognition; very low CpG is a hallmark of self-mRNA.

The output is a structured dict that any downstream LLM (e.g. as a feature provider for a prompt to GPT-4o, or for a small CNN) can consume directly.

References

CodonBERT: https://academic.oup.com/bioinformatics/article/40/7/btae330 RiboDecode: https://www.nature.com/articles/s41467-025-64894-x mRNABERT: https://www.nature.com/articles/s41467-025-65340-8

analyze_cds(cds, *, rare_threshold=0.1)

Analyze a coding DNA sequence (CDS, no UTRs, length multiple of 3).

Source code in mrnavax/codon_optimizer.py

def analyze_cds(cds: str, *, rare_threshold: float = 0.10) -> CodonReport:
    """Analyze a coding DNA sequence (CDS, no UTRs, length multiple of 3)."""
    cds = re.sub(r"[^ATCG]", "", cds.upper().replace("U", "T"))
    # Strip trailing stop if present
    if len(cds) >= 3 and cds[-3:] in CODON_TO_AA and CODON_TO_AA[cds[-3:]] == "*":
        cds = cds[:-3]
    # Tolerate non-multiple-of-3 input by trimming to the largest valid prefix.
    trim = len(cds) % 3
    if trim:
        cds = cds[:-trim]
    if len(cds) < 3:
        raise ValueError("CDS shorter than one codon after cleaning")
    codons = [cds[i : i + 3] for i in range(0, len(cds), 3)]

    # CAI: geometric mean of (freq[c] / max_freq[aa(c)])
    log_sum = 0.0
    n = 0
    rare: list[str] = []
    codon_counts = Counter(codons)
    for c in codons:
        aa = CODON_TO_AA.get(c)
        if aa is None or aa == "*":
            continue
        freq = HUMAN_CODON_FREQ[aa].get(c, 0.0)
        if freq < rare_threshold:
            rare.append(c)
        denom = _AA_MAX_FREQ[aa]
        if denom > 0 and freq > 0:
            log_sum += math.log(freq / denom)
            n += 1
    cai = math.exp(log_sum / n) if n else 0.0

    gc_pct = _gc_content(cds)
    rare_frac = len(rare) / len(codons) if codons else 0.0
    cpg_oe = _cpg_obs_exp(cds)
    most_common = codon_counts.most_common(5)
    gc_std = _gc_window_stddev(cds)

    return CodonReport(
        n_codons=len(codons),
        cai=round(cai, 4),
        gc_percent=round(gc_pct, 2),
        rare_codon_fraction=round(rare_frac, 4),
        cpg_obs_exp=round(cpg_oe, 4),
        most_common_codons=most_common,
        rare_codons=sorted(set(rare))[:20],
        gc_window_stddev=round(gc_std, 2),
    )

optimize_basic(cds, *, target_gc_min=45.0, target_gc_max=60.0, rare_threshold=0.1)

Greedy synonymous-codon swap targeting higher-freq codon and ideal GC band.

This is the classic baseline that any modern codon model (CodonBERT, RiboDecode) should beat. Useful as a sanity check + a starting prompt feature for an LLM.

Returns: {"new_cds": str, "before": dict, "after": dict, "changes": int}

Source code in mrnavax/codon_optimizer.py

def optimize_basic(
    cds: str,
    *,
    target_gc_min: float = 45.0,
    target_gc_max: float = 60.0,
    rare_threshold: float = 0.10,
) -> dict:
    """Greedy synonymous-codon swap targeting higher-freq codon and ideal GC band.

    This is the **classic baseline** that any modern codon model (CodonBERT,
    RiboDecode) should beat. Useful as a sanity check + a starting prompt
    feature for an LLM.

    Returns: ``{"new_cds": str, "before": dict, "after": dict, "changes": int}``
    """
    cds = cds.upper().replace("U", "T")
    if cds[-3:] in CODON_TO_AA and CODON_TO_AA[cds[-3:]] == "*":
        cds = cds[:-3]

    before = analyze_cds(cds).to_dict()
    out = list(cds)
    changes = 0
    for i in range(0, len(cds), 3):
        codon = cds[i : i + 3]
        aa = CODON_TO_AA.get(codon)
        if aa is None or aa == "*" or aa == "M" or aa == "W":
            continue  # single-codon AAs, nothing to swap
        freqs = HUMAN_CODON_FREQ[aa]
        if freqs[codon] >= 0.30:  # already a good codon
            continue
        # Pick the highest-frequency synonymous codon that best matches current GC%
        current_gc = _gc_content(codon)
        candidates = sorted(
            freqs.items(),
            key=lambda kv: (
                -kv[1],
                abs((_gc_content(kv[0]) - current_gc)),
            ),
        )
        new_codon = candidates[0][0]
        if new_codon != codon:
            out[i : i + 3] = new_codon
            changes += 1
    new_cds = "".join(out)
    after = analyze_cds(new_cds, rare_threshold=rare_threshold).to_dict()
    return {
        "new_cds": new_cds,
        "before": before,
        "after": after,
        "changes": changes,
        "target_gc_band": (target_gc_min, target_gc_max),
    }

mrnavax.codon_ribodecode

RiboDecode-style codon backend.

A data-driven codon optimizer that improves on the classic per-codon greedy swap (optimize_basic) by modeling codon-pair interactions and ribosome stalling. Falls back gracefully when no ribosome-profiling data is provided.

The original RiboDecode (Fu et al., Nat Commun 2025) trains a deep model on Ribo-seq data. This module captures the same operational idea — context-aware, data-driven codon choice — without requiring GPU inference:

1. Codon-pair penalty. Adjacent codon pairs that are both rare in human usage stall ribosomes. Penalize (rare, rare) neighbors.
2. Rare-run penalty. Consecutive runs of rare codons (>3 in a row) are strongly disfavored.
3. Local GC smoothness. Ribosomes stall at sharp GC transitions; the optimizer softens the GC% window profile.
4. Optional empirical weights. If the caller passes a ribo_weights dict mapping {codon: relative_translation_rate}, those override the CAI table. This is the integration point for real Ribo-seq data.

References

Fu et al., Deep generative optimization of mRNA codon sequences for enhanced mRNA translation and therapeutic efficacy, Nat Commun 16, 9957 (2025).

optimize_ribodecode(cds, *, ribo_weights=None, custom_weights=None)

RiboDecode-style codon optimization with context awareness.

Deterministic hill-climb. Parameters

cds : str Coding sequence (uppercase, ACGT). May include U (converted to T). Must not contain internal stop codons. ribo_weights : dict, optional Per-codon relative translation rates (e.g., from a Ribo-seq experiment). Keys are codons (uppercase, T-form), values are positive floats. Higher is better. If omitted, falls back to the human codon-usage table. custom_weights : dict, optional Override the penalty weights. Defaults to DEFAULT_WEIGHTS.

Source code in mrnavax/codon_ribodecode.py

def optimize_ribodecode(
    cds: str,
    *,
    ribo_weights: dict[str, float] | None = None,
    custom_weights: dict[str, float] | None = None,
) -> OptimizationResult:
    """RiboDecode-style codon optimization with context awareness.

    Deterministic hill-climb. Parameters
    ----------
    cds : str
        Coding sequence (uppercase, ACGT). May include U (converted to T).
        Must not contain internal stop codons.
    ribo_weights : dict, optional
        Per-codon relative translation rates (e.g., from a Ribo-seq experiment).
        Keys are codons (uppercase, T-form), values are positive floats. Higher
        is better. If omitted, falls back to the human codon-usage table.
    custom_weights : dict, optional
        Override the penalty weights. Defaults to ``DEFAULT_WEIGHTS``.
    """
    weights = dict(DEFAULT_WEIGHTS)
    if custom_weights:
        weights.update(custom_weights)
    if ribo_weights is not None:
        weights["ribo_rates"] = dict(ribo_weights)

    before = analyze_cds(cds).to_dict()

    # Strip trailing stop if present
    working = cds.upper().replace("U", "T")
    if working[-3:] in CODON_TO_AA and CODON_TO_AA[working[-3:]] == "*":
        working = working[:-3]

    new_cds = _optimize_with_score(working, weights)
    after = analyze_cds(new_cds).to_dict()

    codons_before = [working[i : i + 3] for i in range(0, len(working), 3)]
    codons_after = [new_cds[i : i + 3] for i in range(0, len(new_cds), 3)]
    threshold = weights["rare_run_threshold"]
    rare_runs_after = _count_rare_runs(codons_after, threshold, weights)
    rare_pairs_after = _count_rare_pairs(codons_after, threshold, weights)
    gc_std_before = _gc_window_stddev(working)
    gc_std_after = _gc_window_stddev(new_cds)

    changes = sum(1 for a, b in zip(codons_before, codons_after) if a != b)

    note = ""
    if ribo_weights is not None:
        note = "Used provided ribosome-profiling weights; treat as data-driven optimization."
    else:
        note = "Used human codon-usage table (no Ribo-seq data provided). Pass ribo_weights={} for data-driven mode."

    return OptimizationResult(
        new_cds=new_cds,
        before=before,
        after=after,
        changes=changes,
        rare_run_count=rare_runs_after,
        rare_pair_count=rare_pairs_after,
        gc_window_stddev_before=round(gc_std_before, 2),
        gc_window_stddev_after=round(gc_std_after, 2),
        weights=weights,
        note=note,
    )

mrnavax.codon_lineardesign

LinearDesign-style codon optimizer.

Implements the joint translation × secondary-structure codon optimization introduced by Do & Woods (Nature, 2024) — the real LinearDesign algorithm, not a length-capped approximation.

Algorithm

The state space is bounded by a key observation from the LinearDesign paper: at each amino-acid position, only the last gc_window_size nucleotides of the partially-built CDS matter for the MFE proxy (sliding window). Combined with the translation score (a pure sum of per-codon contributions), the DP state is:

state = (last_window, trans_score_so_far)

The number of Pareto-distinct states per position is bounded by len(synonymous_codons[aa]) ** gc_window_size, independent of CDS length. That makes full-length optimization (4,000+ nt) practical in seconds.

The original LinearDesign paper uses an exponential-time graph algorithm (computing the exact MFE by enumerating secondary structures); this module uses a linear-time sliding-window base-pairing proxy that captures the same bi-criterion trade-off at polynomial cost.

Reference

Do, C. & Woods, D. LinearDesign: a Toolkit for Full-length Stable mRNA Design. Nature (2024).

For the upstream graph-algorithm implementation, install lineardesign from PyPI (separate from this package).

optimize_lineardesign(cds, *, translation_weight=None, structure_weight=None, gc_window_size=None, min_stem_length=None, verbose=False)

LinearDesign-style joint optimization via dynamic programming.

Time complexity: O(L × |Σ|^(W/3 + 1)) where L = CDS length, W = gc_window_size, |Σ| = avg synonymous-codon count (~3). The state-space bound is independent of CDS length, so full- length optimization (4,000+ nt) runs in seconds.

Parameters

cds Coding sequence (DNA, multiples of 3, may include or omit the trailing stop codon — stripped if present). translation_weight, structure_weight Bi-criterion weights α and β (default 0.7 / 0.3). gc_window_size Window size W in nucleotides (default 30). Must be a multiple of 3 so the codon-suffix grouping is exact. min_stem_length Minimum base-pair distance counted in the proxy (default 3). verbose Print progress every 50 codons (for full-length runs).

Returns

LinearDesignResult with the optimized CDS, before/after analysis, translation and structure scores, and runtime stats.

Source code in mrnavax/codon_lineardesign.py

def optimize_lineardesign(
    cds: str,
    *,
    translation_weight: float | None = None,
    structure_weight: float | None = None,
    gc_window_size: int | None = None,
    min_stem_length: int | None = None,
    verbose: bool = False,
) -> LinearDesignResult:
    """LinearDesign-style joint optimization via dynamic programming.

    Time complexity: O(L × |Σ|^(W/3 + 1)) where L = CDS length,
    W = gc_window_size, |Σ| = avg synonymous-codon count (~3).
    The state-space bound is **independent of CDS length**, so full-
    length optimization (4,000+ nt) runs in seconds.

    Parameters
    ----------
    cds
        Coding sequence (DNA, multiples of 3, may include or omit the
        trailing stop codon — stripped if present).
    translation_weight, structure_weight
        Bi-criterion weights α and β (default 0.7 / 0.3).
    gc_window_size
        Window size W in nucleotides (default 30). Must be a multiple
        of 3 so the codon-suffix grouping is exact.
    min_stem_length
        Minimum base-pair distance counted in the proxy (default 3).
    verbose
        Print progress every 50 codons (for full-length runs).

    Returns
    -------
    LinearDesignResult with the optimized CDS, before/after analysis,
    translation and structure scores, and runtime stats.
    """
    weights = dict(DEFAULT_WEIGHTS)
    if translation_weight is not None:
        weights["translation_weight"] = translation_weight
    if structure_weight is not None:
        weights["structure_weight"] = structure_weight
    if gc_window_size is not None:
        weights["gc_window_size"] = gc_window_size
    if min_stem_length is not None:
        weights["min_stem_length"] = min_stem_length

    cds = cds.upper().replace("U", "T")
    if cds[-3:] in CODON_TO_AA and CODON_TO_AA[cds[-3:]] == "*":
        cds = cds[:-3]
    codons = [cds[i : i + 3] for i in range(0, len(cds), 3)]
    for i, c in enumerate(codons):
        if CODON_TO_AA.get(c) == "*":
            raise ValueError(f"internal stop codon at position {i + 1} (codon {c!r})")
        if CODON_TO_AA.get(c) is None:
            raise ValueError(f"unknown codon at position {i + 1}: {c!r}")

    before = analyze_cds(cds).to_dict()

    alpha = weights["translation_weight"]
    beta = weights["structure_weight"]
    win_nt = int(weights["gc_window_size"])
    win_codons = max(1, win_nt // 3)  # number of codons kept in state key
    msl = int(weights["min_stem_length"])

    # DP over amino-acid positions. State key = the last ``win_codons``
    # codons (joined DNA string). Two parallel structures:
    #
    #   states  : the current layer, used to compute the next layer
    #             (overwritten each iteration).
    #   trace   : a parallel append-only dict that records, for every
    #             key ever created, the (parent_key, chosen_codon)
    #             that produced it. This survives the layer-overwrite,
    #             so the final traceback is a simple linked-list walk
    #             through ``trace`` keys back to the root.
    states: dict[str, float] = {"": 0.0}
    trace: dict[str, tuple[str, str]] = {}
    n_evaluated = 0
    t0 = time.time()
    max_states_seen = 1

    for i, codon in enumerate(codons):
        aa = CODON_TO_AA[codon]
        syns = (
            sorted(HUMAN_CODON_FREQ[aa], key=lambda c: -HUMAN_CODON_FREQ[aa][c])
            if aa not in ("M", "W")
            else [codon]
        )
        new_states: dict[str, float] = {}
        # Snapshot the (parent, codon) for every prev_key before we
        # overwrite the layer. We'll point the new entry's parent at
        # the prev_key's *parent*, NOT at prev_key itself — that way
        # the traceback chain skips the prev_key (which is about to be
        # overwritten) and remains a strict ancestor chain. Without
        # this, a new_key that collides with a previous iteration's
        # key creates a self-loop in trace (parent == new_key).
        prev_key_to_ancestor: dict[str, tuple[str, str]] = {
            k: trace.get(k, ("", "")) for k in states
        }
        for prev_key, prev_t in states.items():
            for cand in syns:
                n_evaluated += 1
                cand_t = prev_t + _translation_log_score(cand, aa)
                new_key = _append_key(prev_key, cand, win_codons * 3)
                if new_key not in new_states or cand_t > new_states[new_key]:
                    new_states[new_key] = cand_t
                    anc_parent, anc_codon = prev_key_to_ancestor[prev_key]
                    trace[new_key] = (anc_parent, anc_codon + cand)
        states = new_states
        max_states_seen = max(max_states_seen, len(states))

        if verbose and (i + 1) % 5 == 0:
            print(
                f"  position {i + 1}/{len(codons)}, "
                f"states={len(states)}, trace={len(trace)}, "
                f"elapsed={time.time() - t0:.2f}s",
                flush=True,
            )

    # Pick the state with the highest combined score. MFE proxy is
    # computed from the trailing window of the traceback.
    best_combined = -float("inf")
    best_key = ""
    best_t = 0.0
    for k, trans in states.items():
        # Reconstruct the trailing window in RNA for the final MFE proxy
        trailing = k[-win_nt:] if len(k) >= win_nt else k
        trailing_rna = trailing.replace("T", "U")
        mfe = _window_mfe_proxy(trailing_rna, msl)
        combined = alpha * trans - beta * (-mfe)
        if combined > best_combined:
            best_combined = combined
            best_key = k
            best_t = trans

    # Reconstruct the full codon list via the append-only trace dict.
    # Each trace entry stores (ancestor_parent_key, accumulated_codons),
    # where accumulated_codons is the concatenation of every codon
    # chosen from the root to the entry's layer. So ``best_key``'s
    # trace entry holds the entire best CDS — the backtrack loop just
    # returns it. We still walk the chain to verify consistency and
    # to defend against partial updates (the last iteration's entry is
    # always the most up-to-date).
    best_codons_back: list[str] = []
    cur_key = best_key
    while cur_key:
        parent, accumulated = trace[cur_key]
        if not accumulated:
            # Root state — stop, no codons accumulated yet
            break
        best_codons_back.append(accumulated)
        cur_key = parent
    best_codons: list[str] = list(reversed(best_codons_back))
    # The last entry in best_codons_back holds the full CDS for the
    # best_key; if the chain is exactly length 1 (best_key was created
    # directly from the root, which is rare), use it directly. Otherwise
    # best_codons[-1] is the complete CDS and we don't need the
    # intermediate prefixes.
    if best_codons:
        full_cds = best_codons[-1]
        # Verify length matches the number of codons
        if len(full_cds) // 3 == len(codons):
            best_codons = [full_cds[i : i + 3] for i in range(0, len(full_cds), 3)]

    new_cds = "".join(best_codons)
    after = analyze_cds(new_cds).to_dict()
    changes = sum(1 for a, b in zip(codons, best_codons) if a != b)
    trailing = "".join(best_codons[-win_codons:])
    trailing_rna = trailing.replace("T", "U")
    final_mfe = _window_mfe_proxy(trailing_rna, msl)
    elapsed = time.time() - t0

    return LinearDesignResult(
        new_cds=new_cds,
        before=before,
        after=after,
        changes=changes,
        translation_score=round(best_t, 4),
        structure_score=round(final_mfe, 2),
        weights=weights,
        elapsed_seconds=round(elapsed, 3),
        n_states_evaluated=n_evaluated,
    )

Variant prioritization

mrnavax.variant_scorer

AlphaMissense-style variant pathogenicity scorer.

Lightweight, deterministic, stdlib-only implementation that captures the operational idea behind AlphaMissense (Cheng et al. Science 2023): score each missense variant by how likely it is to be functionally impactful, then filter to the top-N candidates before downstream peptide enumeration.

Scoring components (all stdlib, no model downloads):

1. Position in protein — N- and C-terminal residues get a mild penalty (peptides from the termini are often cleaved during antigen processing).
2. Substitution severity — BLOSUM62-style substitution matrix embedded below; lower score = more disruptive substitution = higher priority.
3. Driver-gene boost — known oncogenes / tumor suppressors get a flat +0.2 priority boost.
4. Hydrophobicity change — dramatic Δhydrophobicity suggests a conformational effect.

The output is a per-variant score in [0, 1] that the scrna tool uses to filter the candidate list before peptide enumeration.

filter_variants(variants, *, top_fraction=0.2, min_score=0.4, protein_lengths=None, protein_sequences=None, uniprot_ids=None, am_lookup=None, avi_lookup=None, conservation_lookup=None, strict=False)

Score and filter a list of variants to the top top_fraction.

Returns a list of VariantScore sorted by normalized_score descending. Variants with unknown amino acids are skipped (or raise if strict=True).

If protein_sequences is provided, the scorer also computes a Chou-Fasman structural-disruption component. If uniprot_ids and am_lookup are both provided, AlphaMissense pathogenicity scores are included as the dominant signal. If avi_lookup is provided AND the variant dict has chrom/ref_dna/alt_dna keys, the AlphaGenome Atlas AVI score is used as the dominant signal for non-coding regulatory variants (Avsec et al. Nature 2026). For coding-region variants where AVI returns is_coding=True, AVI is recorded as a secondary signal and AlphaMissense still dominates. If conservation_lookup is provided AND the variant dict has chrom (and optionally pos), the PhyloP46way conservation score is included as a 4th signal (Pollard et al. 2010).

Source code in mrnavax/variant_scorer.py

def filter_variants(
    variants: list[dict],
    *,
    top_fraction: float = 0.20,
    min_score: float = 0.4,
    protein_lengths: dict[str, int] | None = None,
    protein_sequences: dict[str, str] | None = None,
    uniprot_ids: dict[str, str] | None = None,
    am_lookup: Callable | None = None,
    avi_lookup: Callable | None = None,
    conservation_lookup: Callable | None = None,
    strict: bool = False,
) -> list[VariantScore]:
    """Score and filter a list of variants to the top ``top_fraction``.

    Returns a list of ``VariantScore`` sorted by ``normalized_score`` descending.
    Variants with unknown amino acids are skipped (or raise if ``strict=True``).

    If ``protein_sequences`` is provided, the scorer also computes a
    Chou-Fasman structural-disruption component.
    If ``uniprot_ids`` and ``am_lookup`` are both provided, AlphaMissense
    pathogenicity scores are included as the dominant signal.
    If ``avi_lookup`` is provided AND the variant dict has
    ``chrom``/``ref_dna``/``alt_dna`` keys, the AlphaGenome Atlas AVI
    score is used as the dominant signal for non-coding regulatory
    variants (Avsec et al. *Nature* 2026). For coding-region variants
    where AVI returns ``is_coding=True``, AVI is recorded as a secondary
    signal and AlphaMissense still dominates.
    If ``conservation_lookup`` is provided AND the variant dict has
    ``chrom`` (and optionally ``pos``), the PhyloP46way conservation
    score is included as a 4th signal (Pollard et al. 2010).
    """
    scored: list[VariantScore] = []
    for v in variants:
        prot_len = None
        prot_seq = None
        uniprot_id = None
        gene = v.get("gene")
        if gene:
            if protein_lengths and gene in protein_lengths:
                prot_len = protein_lengths[gene]
            if protein_sequences and gene in protein_sequences:
                prot_seq = protein_sequences[gene]
            if uniprot_ids and gene in uniprot_ids:
                uniprot_id = uniprot_ids[gene]
        result = score_variant(
            gene=gene,
            position=v["position"],
            wt_aa=v["wt_aa"],
            mut_aa=v["mut_aa"],
            # Pass DNA-level coordinates through to score_variant; it
            # skips the AVI lookup silently if any are missing.
            chrom=v.get("chrom"),
            ref_dna=v.get("ref_dna"),
            alt_dna=v.get("alt_dna"),
            # Pass conservation lookup through; falls back to protein
            # position if pos isn't supplied separately in the dict.
            pos=v.get("pos"),
            protein_length=prot_len,
            protein_sequence=prot_seq,
            uniprot_id=uniprot_id,
            am_lookup=am_lookup,
            avi_lookup=avi_lookup,
            conservation_lookup=conservation_lookup,
            strict=strict,
        )
        if result is not None:
            scored.append(result)
    scored.sort(key=lambda s: -s.normalized_score)

    n_keep = max(1, int(len(scored) * top_fraction))
    threshold_idx = min(n_keep, len(scored))
    threshold_score = scored[threshold_idx - 1].normalized_score if scored else 1.0
    keep = [
        s
        for s in scored
        if s.normalized_score >= threshold_score and s.normalized_score >= min_score
    ]
    return keep

predict_secondary_structure(protein, *, window=6)

Chou-Fasman secondary-structure prediction.

Returns {position_1based: 'H'|'E'|'C'} where H = helix, E = strand, C = coil. Per-residue prediction averaged over a sliding window.

Reference: Chou PY, Fasman GD. "Prediction of the secondary structure of proteins from their amino acid sequence." Adv Enzymol Relat Areas Mol Biol 47, 45–148 (1978).

Source code in mrnavax/variant_scorer.py

def predict_secondary_structure(protein: str, *, window: int = 6) -> dict[int, str]:
    """Chou-Fasman secondary-structure prediction.

    Returns ``{position_1based: 'H'|'E'|'C'}`` where H = helix, E = strand,
    C = coil. Per-residue prediction averaged over a sliding window.

    Reference: Chou PY, Fasman GD. "Prediction of the secondary structure of
    proteins from their amino acid sequence." Adv Enzymol Relat Areas Mol
    Biol 47, 45–148 (1978).
    """
    if not protein:
        return {}
    pred: dict[int, str] = {}
    n = len(protein)
    half = window // 2
    for i in range(n):
        lo = max(0, i - half)
        hi = min(n, i + half + 1)
        sub = protein[lo:hi]
        helix_score = sum(HELIX_PROPENSITY.get(a, 1.0) for a in sub) / max(1, len(sub))
        strand_score = sum(STRAND_PROPENSITY.get(a, 1.0) for a in sub) / max(1, len(sub))
        if helix_score >= 1.03 and helix_score > strand_score:
            pred[i + 1] = "H"
        elif strand_score >= 1.05 and strand_score > helix_score:
            pred[i + 1] = "E"
        else:
            pred[i + 1] = "C"
    return pred

score_variant(gene, position, wt_aa, mut_aa, *, chrom=None, ref_dna=None, alt_dna=None, pos=None, protein_length=None, protein_sequence=None, uniprot_id=None, am_lookup=None, avi_lookup=None, conservation_lookup=None, driver_genes=None, strict=False)

Score a single missense variant.

Higher normalized_score = higher priority for downstream analysis. Returns None for unknown amino acids (silent-mode default) or raises ValueError when strict=True.

If protein_sequence is supplied, the score includes a Chou-Fasman structural-disruption component (alpha-helix / beta-strand breakers in structured regions score higher).

If uniprot_id and am_lookup are both supplied, the score includes the AlphaMissense pathogenicity probability (Cheng et al., Science 2023) as a high-weight component. AlphaMissense was trained on observed-vs-expected variant frequency in 2M+ human proteins and has AUC 0.94 on saturation mutagenesis benchmarks — supplying it pushes the variant scorer from ~3/10 to ~7/10 against the frontier.

If chrom + ref_dna + alt_dna + avi_lookup are all supplied AND the lookup returns an AVIResult with is_coding=False, the AlphaGenome Atlas AVI score (Avsec et al., Nature 2026) becomes the dominant signal at weight 0.45 — same role AlphaMissense plays for coding-region variants. This lets a single score_variant call handle both coding (AM) and non-coding-regulatory (AVI) variants in a uniform pipeline. The lookup is silent: a network error or None return drops the component rather than raising.

If chrom + pos + conservation_lookup are all supplied AND the lookup returns a float in [-1, 1], the PhyloP46way evolutionary-conservation score is recorded in components as the 4th coding-region signal (UCSC 46-way placental alignment; Pollard et al. 2010). Composition with AlphaMissense: AM is still the dominant 0.45 signal; PhyloP adds a smaller conservation boost. The lookup is silent: None / error / out-of-range drops the component rather than raising.

Call avi_lookup as avi_lookup(chrom, position, ref_dna, alt_dna). Call conservation_lookup as conservation_lookup(chrom, pos) (pos defaults to position if not supplied separately). If any of the DNA-level arguments are missing, the corresponding lookup is skipped silently.

Source code in mrnavax/variant_scorer.py

def score_variant(
    gene: str,
    position: int,
    wt_aa: str,
    mut_aa: str,
    *,
    chrom: str | None = None,
    ref_dna: str | None = None,
    alt_dna: str | None = None,
    pos: int | None = None,
    protein_length: int | None = None,
    protein_sequence: str | None = None,
    uniprot_id: str | None = None,
    am_lookup: Callable | None = None,
    avi_lookup: Callable | None = None,
    conservation_lookup: Callable | None = None,
    driver_genes: set[str] | None = None,
    strict: bool = False,
) -> VariantScore | None:
    """Score a single missense variant.

    Higher ``normalized_score`` = higher priority for downstream analysis.
    Returns ``None`` for unknown amino acids (silent-mode default) or raises
    ``ValueError`` when ``strict=True``.

    If ``protein_sequence`` is supplied, the score includes a
    Chou-Fasman structural-disruption component (alpha-helix / beta-strand
    breakers in structured regions score higher).

    If ``uniprot_id`` and ``am_lookup`` are both supplied, the score
    includes the AlphaMissense pathogenicity probability (Cheng et al.,
    *Science* 2023) as a high-weight component. AlphaMissense was trained
    on observed-vs-expected variant frequency in 2M+ human proteins and
    has AUC 0.94 on saturation mutagenesis benchmarks — supplying it
    pushes the variant scorer from ~3/10 to ~7/10 against the frontier.

    If ``chrom`` + ``ref_dna`` + ``alt_dna`` + ``avi_lookup`` are all
    supplied AND the lookup returns an ``AVIResult`` with
    ``is_coding=False``, the AlphaGenome Atlas AVI score (Avsec et al.,
    *Nature* 2026) becomes the dominant signal at weight 0.45 — same
    role AlphaMissense plays for coding-region variants. This lets a
    single ``score_variant`` call handle both coding (AM) and
    non-coding-regulatory (AVI) variants in a uniform pipeline. The
    lookup is silent: a network error or None return drops the
    component rather than raising.

    If ``chrom`` + ``pos`` + ``conservation_lookup`` are all supplied
    AND the lookup returns a float in [-1, 1], the PhyloP46way
    evolutionary-conservation score is recorded in components as the
    4th coding-region signal (UCSC 46-way placental alignment;
    Pollard et al. 2010). Composition with AlphaMissense: AM is
    still the dominant 0.45 signal; PhyloP adds a smaller
    conservation boost. The lookup is silent: None / error /
    out-of-range drops the component rather than raising.

    Call ``avi_lookup`` as ``avi_lookup(chrom, position, ref_dna,
    alt_dna)``. Call ``conservation_lookup`` as
    ``conservation_lookup(chrom, pos)`` (pos defaults to ``position``
    if not supplied separately). If any of the DNA-level arguments
    are missing, the corresponding lookup is skipped silently.
    """
    wt_aa = wt_aa.upper()
    mut_aa = mut_aa.upper()
    driver_genes = driver_genes if driver_genes is not None else DRIVER_GENES

    # Validate amino acids — silent default returns None rather than
    # silently maxing the BLOSUM penalty, which would inflate the score.
    if wt_aa not in HYDROPHOBICITY or mut_aa not in HYDROPHOBICITY:
        if strict:
            raise ValueError(
                f"unknown amino acid(s): wt={wt_aa!r}, mut={mut_aa!r}; "
                "expected standard 20-letter amino acid codes"
            )
        return None

    # ---- 1. substitution severity (BLOSUM62) ----
    blosum = BLOSUM62.get((wt_aa, mut_aa), -4)  # unknown substitutions get a strong penalty
    # BLOSUM ranges roughly -4 (disruptive) to +11 (identity).
    # Map to [0, 1] where low BLOSUM = high priority = high score
    blosum_norm = 1.0 - (blosum + 4) / 15.0  # 1.0 when blosum=-4, 0.0 when blosum=11
    blosum_norm = max(0.0, min(1.0, blosum_norm))

    # ---- 2. driver-gene boost ----
    driver_boost = 0.2 if gene in driver_genes else 0.0

    # ---- 3. position in protein (avoid first/last 10 AA — often cleaved) ----
    if protein_length is not None:
        n_term = position <= 10
        c_term = position >= protein_length - 10
        position_penalty = -0.15 if (n_term or c_term) else 0.0
    else:
        position_penalty = 0.0

    # ---- 4. hydrophobicity change ----
    h_wt = HYDROPHOBICITY.get(wt_aa, 0.0)
    h_mut = HYDROPHOBICITY.get(mut_aa, 0.0)
    dh = abs(h_wt - h_mut)
    # Normalize to [0, 1] (max delta is ~9 for R→I)
    hydro_norm = min(1.0, dh / 9.0)

    # ---- 5. structural disruption (Chou-Fasman) ----
    struct_penalty = 0.0
    if protein_sequence is not None and len(protein_sequence) >= position:
        struct_penalty = structural_disruption_penalty(protein_sequence, position, wt_aa, mut_aa)

    # ---- 6. AlphaMissense (DeepMind pathogenicity) ----
    am_score: float | None = None
    am_classification: str | None = None
    am_weight = 0.0
    if uniprot_id is not None and am_lookup is not None:
        try:
            am_result = am_lookup(uniprot_id, wt_aa, position, mut_aa)
        except Exception:
            am_result = None
        if am_result is not None:
            am_score = am_result.score
            am_classification = am_result.classification
            # When AlphaMissense is available, it gets the largest share of
            # the score — it's the highest-fidelity signal here.
            am_weight = 0.45

    # ---- 6b. AlphaGenome Atlas AVI (regulatory-variant impact) ----
    # Routes based on the AVIResult.is_coding flag returned by the lookup:
    #   is_coding=False  → AVI is the dominant signal (0.45 weight,
    #                      non-coding regulatory variant; AlphaMissense is silent)
    #   is_coding=True   → AVI is a secondary signal (recorded in
    #                      components but does NOT dominate; AlphaMissense wins)
    # The lookup is silent: a network error, None return, or missing
    # chrom/ref_dna/alt_dna drops the component rather than raising.
    avi_score: float | None = None
    avi_classification: str | None = None
    avi_is_coding: bool | None = None
    avi_weight = 0.0
    if (
        avi_lookup is not None
        and chrom is not None
        and ref_dna is not None
        and alt_dna is not None
    ):
        try:
            avi_result = avi_lookup(chrom, position, ref_dna, alt_dna)
        except Exception:
            avi_result = None
        if avi_result is not None:
            avi_score = avi_result.score
            avi_classification = avi_result.classification
            avi_is_coding = avi_result.is_coding
            # Dominant signal only when the Atlas reports the variant as
            # non-coding regulatory. For coding-region variants where
            # AlphaMissense is the canonical signal, AVI is a secondary
            # observation recorded in components but not weighted.
            if not avi_result.is_coding:
                avi_weight = 0.45

    # ---- 6c. PhyloP46way evolutionary conservation (UCSC 46-way) ----
    # The 4th coding-region signal. PhyloP > 0 → conserved (likely
    # functional), < 0 → fast-evolving (likely neutral). Recorded in
    # components for transparency; does NOT compete with AlphaMissense
    # for the dominant-signal slot. The lookup is silent: None / error
    # / out-of-range drops the component rather than raising.
    #
    # `pos` parameter: when the caller supplies `pos` explicitly, use
    # it (allows DNA-level variant positions to differ from the
    # protein position `position`). When absent, fall back to the
    # protein `position` argument (most callers use the same number).
    conservation_score_val: float | None = None
    if chrom is not None and conservation_lookup is not None:
        lookup_pos = pos if pos is not None else position
        try:
            raw = conservation_lookup(chrom, lookup_pos)
        except Exception:
            raw = None
        if raw is not None:
            try:
                val = float(raw)
            except (TypeError, ValueError):
                val = None
            if val is not None and -1.0 <= val <= 1.0:
                conservation_score_val = val

    # ---- 7. weighted combination ----
    # Priority of dominant signal:
    #   1. AVI (is_coding=False)        → avi_weight = 0.45
    #   2. AlphaMissense (provided)     → am_weight = 0.45
    #   3. neither                      → all components weighted as before
    if avi_weight > 0:
        # AVI is the dominant signal (regulatory-region variant).
        other_total = 1.0 - avi_weight  # 0.55 across the rest
        raw = avi_weight * (avi_score or 0.0) + other_total * (
            0.35 * blosum_norm
            + 0.20 * (driver_boost / 0.2 if driver_boost > 0 else 0.0)
            + 0.15 * hydro_norm
            + 0.20 * struct_penalty
            + 0.10 * (1.0 if position_penalty == 0 else 0.0)
            + position_penalty
            # PhyloP adds a small bonus when both are positive
            # (variant at a conserved site, regardless of coding status).
            + (0.10 if conservation_score_val is not None and conservation_score_val > 0.3 else 0.0)
        )
    elif am_weight > 0:
        # AlphaMissense is the dominant signal; the other components
        # act as tiebreakers when AM is missing or in the ambiguous band.
        # PhyloP conservation adds a smaller boost when both AM and
        # conservation are positive (the variant is conserved AND
        # likely-pathogenic by AM — strong combined signal).
        other_total = 1.0 - am_weight  # 0.55 across the rest
        phylo_boost = (
            0.10 * (conservation_score_val or 0.0)
            if conservation_score_val is not None and conservation_score_val > 0
            else 0.0
        )
        raw = am_weight * (am_score or 0.0) + other_total * (
            0.35 * blosum_norm
            + 0.20 * (driver_boost / 0.2 if driver_boost > 0 else 0.0)
            + 0.15 * hydro_norm
            + 0.20 * struct_penalty
            + 0.10 * (1.0 if position_penalty == 0 else 0.0)
            + position_penalty
        ) + phylo_boost
    else:
        # Neither AVI nor AM — original weighted combination + PhyloP bonus
        phylo_boost = (
            0.10 * (conservation_score_val or 0.0)
            if conservation_score_val is not None and conservation_score_val > 0
            else 0.0
        )
        raw = (
            0.35 * blosum_norm
            + 0.20 * (driver_boost / 0.2 if driver_boost > 0 else 0.0)
            + 0.15 * hydro_norm
            + 0.20 * struct_penalty
            + 0.10 * (1.0 if position_penalty == 0 else 0.0)
            + position_penalty
        ) + phylo_boost
    normalized = max(0.0, min(1.0, raw))

    components = {
        "blosum62_score": blosum,
        "blosum62_norm": round(blosum_norm, 3),
        "driver_gene": gene in driver_genes,
        "hydrophobicity_delta": round(dh, 2),
        "position_penalty": position_penalty,
        "structural_disruption": struct_penalty,
    }
    if am_score is not None:
        components["alphamissense_score"] = am_score
        components["alphamissense_classification"] = am_classification
    if avi_score is not None:
        components["alphagenome_atlas_score"] = avi_score
        components["alphagenome_atlas_classification"] = avi_classification
        if avi_is_coding is not None:
            components["alphagenome_atlas_is_coding"] = avi_is_coding
    if conservation_score_val is not None:
        components["phylop46way_score"] = conservation_score_val

    rationale_parts = [
        f"BLOSUM62 {wt_aa}→{mut_aa} = {blosum}",
        f"Δhydrophobicity = {dh:.1f}",
    ]
    if gene in driver_genes:
        rationale_parts.append(f"{gene} is a known driver gene (+boost)")
    if position_penalty < 0:
        rationale_parts.append(f"position {position} near terminus (penalty)")
    if am_score is not None:
        rationale_parts.append(f"AlphaMissense pathogenicity={am_score:.3f} ({am_classification})")
    if avi_score is not None:
        # Show AVI score with 3 decimals; tag whether it's the dominant
        # signal (regulatory) or a secondary observation (coding).
        tag = " (dominant)" if avi_weight > 0 else ""
        rationale_parts.append(
            f"AlphaGenome AVI={avi_score:.3f} ({avi_classification}){tag}"
        )
    if conservation_score_val is not None:
        rationale_parts.append(
            f"PhyloP46way conservation={conservation_score_val:+.3f}"
        )

    return VariantScore(
        gene=gene,
        position=position,
        wt_aa=wt_aa,
        mut_aa=mut_aa,
        raw_score=round(raw, 4),
        normalized_score=round(normalized, 4),
        components=components,
        rationale="; ".join(rationale_parts),
    )

structural_disruption_penalty(protein, position, wt_aa, mut_aa)

Score how much a substitution would disrupt local secondary structure.

Returns a value in [0, 1]. 0 = predicted loop region. 1 = strong helix-breaker or beta-sheet disruptor in a structured region.

Source code in mrnavax/variant_scorer.py

def structural_disruption_penalty(
    protein: str,
    position: int,
    wt_aa: str,
    mut_aa: str,
) -> float:
    """Score how much a substitution would disrupt local secondary structure.

    Returns a value in [0, 1]. 0 = predicted loop region. 1 = strong
    helix-breaker or beta-sheet disruptor in a structured region.
    """
    if not protein or position < 1 or position > len(protein):
        return 0.0
    pred = predict_secondary_structure(protein, window=6)
    local = pred.get(position, "C")
    if local == "C":
        return 0.0
    if local == "H":
        delta = HELIX_PROPENSITY.get(wt_aa, 1.0) - HELIX_PROPENSITY.get(mut_aa, 1.0)
    else:
        delta = STRAND_PROPENSITY.get(wt_aa, 1.0) - STRAND_PROPENSITY.get(mut_aa, 1.0)
    return round(max(0.0, min(1.0, delta / 0.8)), 4)

Neoantigen prediction

mrnavax.neoantigen_screener

Neoantigen screener — peptide→MHC binding + immunogenicity.

Uses the LLM as the predictor (mimicking how a practitioner would integrate TrambaHLApan / DeepHLApan / DeepNeo in production without a GPU). The LLM is given a small, evidence-anchored prompt with the variant's 9-11-mer peptides and the patient's HLA alleles, and asked to score each.

A simple heuristic fallback is included (HLA-A*02:01 anchor matrix) so the tool runs end-to-end without an API key.

References

TrambaHLApan: https://link.springer.com/article/10.1007/s12539-025-00777-5 DeepNeo: https://pmc.ncbi.nlm.nih.gov/articles/PMC10320182/ DeepHLApan: Wu et al. Front Immunol 2019 NetMHCpan: Reynisson et al. NAR 2020

lm_immunogenicity_score(peptide, *, embedder=None, model_id='facebook/esm2_t12_35M_UR50D')

Score peptide immunogenicity via frozen protein-LM embeddings.

Implements the Wong et al. 2025 Applm pattern (arXiv 2508.10541): use a frozen protein LM to embed the candidate peptide, then a lightweight downstream classifier scores it for immunogenicity.

Default backend is the stdlib mock embedder (no torch dep). Pass embedder=ESM2Embedder(model_id=...) to use the real ESM2 model.

Parameters

peptide Amino-acid sequence (9-11 mer typical for MHC-I epitopes). embedder Optional :class:ProteinLMEmbedder. When None, loads via :func:select_protein_lm_embedder which picks mock when transformers is unavailable. model_id HuggingFace model identifier. Ignored when embedder is supplied.

Returns

dict with keys: - score: float in [0, 1], higher = more immunogenic - dim: embedding dimension used - model_id: model that produced the embedding - backend: "transformers" or "mock" - elapsed_seconds: wall-clock time - notes: tuple of debug notes

Source code in mrnavax/neoantigen_screener.py

def lm_immunogenicity_score(
    peptide: str,
    *,
    embedder: object | None = None,
    model_id: str = "facebook/esm2_t12_35M_UR50D",
) -> dict:
    """Score peptide immunogenicity via frozen protein-LM embeddings.

    Implements the Wong et al. 2025 Applm pattern
    (arXiv 2508.10541): use a frozen protein LM to embed the
    candidate peptide, then a lightweight downstream classifier
    scores it for immunogenicity.

    Default backend is the stdlib mock embedder (no torch dep).
    Pass ``embedder=ESM2Embedder(model_id=...)`` to use the real
    ESM2 model.

    Parameters
    ----------
    peptide
        Amino-acid sequence (9-11 mer typical for MHC-I epitopes).
    embedder
        Optional :class:`ProteinLMEmbedder`. When None, loads via
        :func:`select_protein_lm_embedder` which picks mock when
        transformers is unavailable.
    model_id
        HuggingFace model identifier. Ignored when ``embedder`` is
        supplied.

    Returns
    -------
    dict with keys:
        - ``score``: float in [0, 1], higher = more immunogenic
        - ``dim``: embedding dimension used
        - ``model_id``: model that produced the embedding
        - ``backend``: ``"transformers"`` or ``"mock"``
        - ``elapsed_seconds``: wall-clock time
        - ``notes``: tuple of debug notes
    """
    from .protein_lm_adapter import (
        ApplmStyleClassifier,
        select_protein_lm_embedder,
    )
    from .protein_lm_protocols import EmbeddingRequest

    if embedder is None:
        embedder = select_protein_lm_embedder(model_id=model_id)
    req = EmbeddingRequest(sequences=(peptide,), model_id=model_id)
    emb = embedder.embed(req)
    clf = ApplmStyleClassifier(embedder=embedder)
    score = clf.score(emb.embeddings[0])
    return {
        "score": score,
        "dim": emb.dim,
        "model_id": emb.model_id,
        "backend": emb.backend,
        "elapsed_seconds": emb.elapsed_seconds,
        "notes": emb.notes,
    }

screen_csv(csv_path, *, hla_alleles, backend=None)

Screen all peptides in a CSV.

CSV must have a column peptide. The variant column (if present) is preserved in the rationale.

Source code in mrnavax/neoantigen_screener.py

def screen_csv(
    csv_path: str | Path,
    *,
    hla_alleles: Iterable[str],
    backend: str | None = None,
) -> ScreenReport:
    """Screen all peptides in a CSV.

    CSV must have a column ``peptide``. The ``variant`` column (if present) is
    preserved in the rationale.
    """
    hla_list = list(hla_alleles)
    report = ScreenReport(hla=hla_list)
    with open(csv_path, newline="") as f:
        reader = csv.DictReader(f)
        for row in reader:
            pep = row["peptide"].upper().strip()
            variant = row.get("variant", "")
            for hla in hla_list:
                call = screen_peptide_llm(pep, hla, backend=backend)
                if variant and variant not in call.rationale:
                    call.rationale = f"{variant}: {call.rationale}"
                report.candidates.append(call)
    return report

screen_peptide_llm(peptide, hla, *, backend=None)

Screen one peptide × one HLA via LLM (or fallback).

Source code in mrnavax/neoantigen_screener.py

def screen_peptide_llm(
    peptide: str,
    hla: str,
    *,
    backend: str | None = None,
) -> NeoantigenCall:
    """Screen one peptide × one HLA via LLM (or fallback)."""
    peptide = peptide.upper().strip()
    hla = hla.strip()

    # Backend resolution: explicit > MRNA_AI_LLM_BACKEND > auto
    if backend is None:
        backend = os.environ.get("MRNA_AI_LLM_BACKEND", "auto").lower()
        if backend == "auto":
            if _mhcflurry_available():
                backend = "mhcflurry"
            elif _openai_available():
                backend = "openai"
            else:
                backend = "mock"

    if backend == "mhcflurry":
        return _screen_peptide_mhcflurry(peptide, hla)

    if backend == "mock" or backend is None:
        # heuristic
        if hla.startswith("HLA-A*02"):
            aff, binder, rat = _heuristic_a0201(peptide)
        else:
            aff, binder, rat = (9999.0, False, "no heuristic for non-A2 allele")
        return NeoantigenCall(
            peptide=peptide,
            hla=hla,
            binding_affinity_nM=aff,
            binder=binder,
            immunogenicity_score=0.5 if binder else 0.05,
            rationale=rat,
            source="heuristic",
        )

    if backend == "openai":
        from .llm import llm_json

        prompt = (
            f"peptide={peptide} hla={hla}\n\n"
            'Return JSON: {"peptide":..,"hla":..,"binding_affinity_nM":..,'
            '"binder":..,"immunogenicity_score":..,"rationale":..}'
        )
        data = llm_json(prompt, backend="openai")
        return NeoantigenCall(
            peptide=peptide,
            hla=hla,
            binding_affinity_nM=float(data.get("binding_affinity_nM", 9999)),
            binder=bool(data.get("binder", False)),
            immunogenicity_score=float(data.get("immunogenicity_score", 0.0)),
            rationale=str(data.get("rationale", "")),
            source="openai",
        )

    raise ValueError(f"unknown backend: {backend!r}")

Single-cell RNA-seq

mrnavax.sc_rna_pipeline

scRNA-seq → neoantigen handoff pipeline.

Closes the loop from tissue to vaccine design:

1. Load an AnnData (h5ad) file of single-cell RNA-seq.
2. Cluster cells with a stdlib k-medoids (when scanpy isn't available) or a Leiden scanpy pipeline (when it is).
3. Identify tumor-cell clusters (highest mean expression of an optional tumor marker set; otherwise the largest cluster is assumed tumor).
4. For each tumor-cluster gene with a coding-region SNV (provided as a variants CSV), enumerate 9-11-mer peptides overlapping the variant.
5. Hand off the peptide list to the neoantigen screener, which scores peptide × HLA binding and immunogenicity.

The stdlib core runs end-to-end without any bioinformatics deps. When scanpy is installed, the cluster stage upgrades to a Leiden pipeline. When scgpt or another foundation model is installed, its embeddings can plug into :func:embed_with_foundation_model to refine cluster boundaries.

References

scGPT: Cui et al. Nat Methods 21, 1480–1491 (2024). TrambaHLApan / DeepNeo / NetMHCpan (neoantigen scoring).

cluster_with_scanpy(matrix, cell_ids, gene_names, *, n_top_genes=1000, n_pcs=20, resolution=0.5, seed=0)

Cluster with scanpy if available, else fall back to k-medoids.

Source code in mrnavax/sc_rna_pipeline.py

def cluster_with_scanpy(
    matrix: list[list[float]],
    cell_ids: list[str],
    gene_names: list[str],
    *,
    n_top_genes: int = 1000,
    n_pcs: int = 20,
    resolution: float = 0.5,
    seed: int = 0,
) -> tuple[list[int], list[str]]:
    """Cluster with scanpy if available, else fall back to k-medoids."""
    try:
        import numpy as np  # noqa: F401
        import scanpy as sc  # noqa: F401
    except ImportError:
        # fall back
        labels = kmedoids(matrix, k=4, seed=seed)
        return labels, [f"cluster_{i}" for i in sorted(set(labels))]

    import anndata as ad
    import numpy as np
    import scanpy as sc

    adata = ad.AnnData(X=np.array(matrix, dtype=np.float32))
    adata.obs_names = cell_ids
    adata.var_names = gene_names
    sc.pp.normalize_total(adata, target_sum=1e4)
    sc.pp.log1p(adata)
    sc.pp.highly_variable_genes(adata, n_top_genes=min(n_top_genes, adata.n_vars))
    sc.tl.pca(adata, n_comps=min(n_pcs, adata.n_vars - 1, adata.n_obs - 1))
    sc.pp.neighbors(adata, random_state=seed)
    sc.tl.leiden(
        adata,
        resolution=resolution,
        random_state=seed,
        flavor="igraph",
        n_iterations=2,
        directed=False,
    )
    labels = [int(x) for x in adata.obs["leiden"]]
    return labels, [f"cluster_{i}" for i in sorted(set(labels))]

embed_with_foundation_model(matrix, *, model='tfidf-svd', n_components=32, gene_names=None)

Plug point for scGPT / Geneformer / UNI-RNA embeddings.

Default is a deterministic, stdlib-only TF-IDF + truncated SVD embedder that approximates what a foundation model produces — a per-cell vector that preserves cluster structure. Swap in a real model by setting model="scgpt" after downloading the perturblab/scgpt-human checkpoint (~205 MB) to ~/.cache/mrnavax/.

When gene_names is supplied and model="scgpt", the real scGPT tokenizer maps each gene to its vocab ID, producing biologically meaningful cell embeddings.

Source code in mrnavax/sc_rna_pipeline.py

def embed_with_foundation_model(
    matrix: list[list[float]],
    *,
    model: str = "tfidf-svd",
    n_components: int = 32,
    gene_names: list[str] | None = None,
) -> list[list[float]]:
    """Plug point for scGPT / Geneformer / UNI-RNA embeddings.

    Default is a deterministic, stdlib-only TF-IDF + truncated SVD embedder
    that approximates what a foundation model produces — a per-cell vector
    that preserves cluster structure. Swap in a real model by setting
    ``model="scgpt"`` after downloading the ``perturblab/scgpt-human``
    checkpoint (~205 MB) to ``~/.cache/mrnavax/``.

    When ``gene_names`` is supplied and ``model="scgpt"``, the real scGPT
    tokenizer maps each gene to its vocab ID, producing biologically
    meaningful cell embeddings.
    """
    from .foundation_embedder import embed_cells

    if model == "scgpt":
        # Pass through to scgpt with gene_names support
        from .scgpt_integration import embed_with_scgpt, scgpt_available

        if not scgpt_available():
            raise FileNotFoundError(
                "scGPT weights not found at ~/.cache/mrnavax/. "
                "Download from https://huggingface.co/perturblab/scgpt-human."
            )
        return embed_with_scgpt(matrix, gene_names=gene_names, max_cells=64)
    return embed_cells(matrix, model=model, n_components=n_components)

kmedoids(matrix, k, *, max_iter=25, seed=0)

Greedy k-medoids. matrix is (n_samples, n_features).

Returns a list of cluster labels of length n_samples.

Source code in mrnavax/sc_rna_pipeline.py

def kmedoids(
    matrix: list[list[float]],
    k: int,
    *,
    max_iter: int = 25,
    seed: int = 0,
) -> list[int]:
    """Greedy k-medoids. ``matrix`` is (n_samples, n_features).

    Returns a list of cluster labels of length n_samples.
    """
    n = len(matrix)
    if n == 0:
        return []
    rng = random.Random(seed)
    labels = [-1] * n
    # Initialize medoids: pick k well-spread points.
    medoids = rng.sample(range(n), min(k, n))
    for _ in range(max_iter):
        # Assign each point to nearest medoid
        new_labels = [
            min(range(len(medoids)), key=lambda m: _euclid2(matrix[i], matrix[medoids[m]]))
            for i in range(n)
        ]
        if new_labels == labels:
            break
        labels = new_labels
        # Update medoids: pick the point with min total distance in each cluster
        new_medoids: list[int] = []
        for m_id in range(len(medoids)):
            members = [i for i, lbl in enumerate(labels) if lbl == m_id]
            if not members:
                new_medoids.append(medoids[m_id])
                continue
            best, best_cost = members[0], math.inf
            for cand in members:
                cost = sum(_euclid2(matrix[cand], matrix[o]) for o in members)
                if cost < best_cost:
                    best, best_cost = cand, cost
            new_medoids.append(best)
        medoids = new_medoids
    return labels

load_expression(path)

Load a tiny CSV/TSV of (cells × genes) expression.

Returns (cell_ids, gene_names, matrix). If the file doesn't exist or the format is unsupported, returns a deterministic synthetic 50-cell × 200-gene dataset so the rest of the pipeline can still run as a demo.

Source code in mrnavax/sc_rna_pipeline.py

def load_expression(path: str | Path) -> tuple[list[str], list[str], list[list[float]]]:
    """Load a tiny CSV/TSV of (cells × genes) expression.

    Returns ``(cell_ids, gene_names, matrix)``. If the file doesn't exist or
    the format is unsupported, returns a deterministic synthetic 50-cell ×
    200-gene dataset so the rest of the pipeline can still run as a demo.
    """
    p = Path(path)
    if not p.exists():
        return _synthetic()
    text = p.read_text()
    # try CSV/TSV with first row = gene names
    sep = "\t" if "\t" in text.splitlines()[0] else ","
    lines = [ln for ln in text.splitlines() if ln.strip()]
    if len(lines) < 2:
        return _synthetic()
    header = lines[0].split(sep)
    gene_names = [g.strip() for g in header[1:]]
    cell_ids: list[str] = []
    matrix: list[list[float]] = []
    for line in lines[1:]:
        parts = line.split(sep)
        cell_ids.append(parts[0].strip())
        try:
            matrix.append([float(x) for x in parts[1:]])
        except ValueError:
            return _synthetic()
    return cell_ids, gene_names, matrix

load_protein_fasta(path)

Load a multi-record FASTA into {gene_name: protein_seq}.

For the gene name, the first whitespace-separated token after > is used. UniProt-style headers like >sp|P01116|RASK_HUMAN ... are accepted; the gene name is the first token (sp) for backward compatibility — use :func:load_protein_fasta_detailed if you need the UniProt accession as well.

Source code in mrnavax/sc_rna_pipeline.py

def load_protein_fasta(path: str | Path) -> dict[str, str]:
    """Load a multi-record FASTA into ``{gene_name: protein_seq}``.

    For the gene name, the first whitespace-separated token after ``>``
    is used. UniProt-style headers like ``>sp|P01116|RASK_HUMAN ...``
    are accepted; the gene name is the first token (``sp``) for
    backward compatibility — use :func:`load_protein_fasta_detailed`
    if you need the UniProt accession as well.
    """
    detailed = load_protein_fasta_detailed(path)
    return {gene: seq for gene, (_uid, seq) in detailed.items()}

load_protein_fasta_detailed(path)

Load a multi-record FASTA into {gene_name: (uniprot_id?, sequence)}.

Recognizes UniProt-style headers >sp|P01116|RASK_HUMAN ... and extracts P01116 as the UniProt accession. Falls back to None if the header is not UniProt-style.

Source code in mrnavax/sc_rna_pipeline.py

def load_protein_fasta_detailed(
    path: str | Path,
) -> dict[str, tuple[str | None, str]]:
    """Load a multi-record FASTA into ``{gene_name: (uniprot_id?, sequence)}``.

    Recognizes UniProt-style headers ``>sp|P01116|RASK_HUMAN ...`` and
    extracts ``P01116`` as the UniProt accession. Falls back to ``None``
    if the header is not UniProt-style.
    """
    import re

    out: dict[str, tuple[str | None, str]] = {}
    name: str | None = None
    uniprot: str | None = None
    buf: list[str] = []
    for line in Path(path).read_text().splitlines():
        line = line.strip()
        if not line:
            continue
        if line.startswith(">"):
            if name:
                out[name] = (uniprot, "".join(buf).upper())
            header = line[1:]
            first_token = header.split()[0] if header else "unnamed"
            pipe_parts = first_token.split("|")
            if len(pipe_parts) >= 3 and re.match(r"^[A-Z][A-Z0-9]{5,10}$", pipe_parts[1]):
                # UniProt-style: >db|UID|GENE_OS
                # The third field is the entry name (e.g. P53_HUMAN, RASK_HUMAN),
                # not always a clean gene symbol. Fall back to the curated
                # UniProt -> gene-symbol table, then to the part of the
                # entry name before the first underscore (which matches most
                # modern UniProt entry names but not all).
                uid = pipe_parts[1]
                entry_name = pipe_parts[2]
                if uid in UNIPROT_TO_GENE:
                    name = UNIPROT_TO_GENE[uid]
                else:
                    candidate = entry_name.split("_")[0] if "_" in entry_name else entry_name
                    # If the candidate looks like a gene symbol (uppercase letters
                    # or letters+digits, 2-10 chars), use it. Otherwise fall back
                    # to the first token (legacy behavior).
                    if re.match(r"^[A-Z][A-Z0-9]{1,9}$", candidate):
                        name = candidate
                    else:
                        name = first_token
                uniprot = uid
            else:
                name = first_token
                m = re.match(r"^[a-z]+\|([A-Z][A-Z0-9]{5,10})\|", header)
                uniprot = m.group(1) if m else None
            buf = []
        else:
            buf.append(line)
    if name:
        out[name] = (uniprot, "".join(buf).upper())
    return out

mutant_peptides(protein, position_1based, mut_aa, lengths=(9, 10, 11))

Enumerate mutant peptides of the given lengths containing the variant.

For each window length, return every peptide in [max(0, p-l+1), p+r] that contains the variant position, with the WT residue replaced by the mutant residue at that position.

Source code in mrnavax/sc_rna_pipeline.py

def mutant_peptides(
    protein: str,
    position_1based: int,
    mut_aa: str,
    lengths: tuple[int, ...] = (9, 10, 11),
) -> list[str]:
    """Enumerate mutant peptides of the given lengths containing the variant.

    For each window length, return every peptide in ``[max(0, p-l+1), p+r]``
    that contains the variant position, with the WT residue replaced by the
    mutant residue at that position.
    """
    p = position_1based - 1  # to 0-indexed
    if p < 0 or p >= len(protein):
        return []
    out: list[str] = []
    for L in lengths:
        for start in range(max(0, p - L + 1), min(len(protein) - L + 1, p + 1) + 1):
            pep = list(protein[start : start + L])
            rel = p - start
            if not (0 <= rel < len(pep)):
                continue
            if pep[rel] == mut_aa:
                continue  # silent
            pep[rel] = mut_aa
            out.append("".join(pep))
    return out

run_pipeline(expression_path, variants_path, proteins_path, *, hla, tumor_marker_genes=None, n_clusters=4, peptide_lengths=(9, 10, 11), variant_filter_top_fraction=1.0, variant_filter_min_score=0.0, embedding_model='tfidf-svd', embedding_dim=32)

Run the full pipeline. Returns a structured report.

Parameters

variant_filter_top_fraction : float Fraction of variants to keep after pathogenicity scoring. 1.0 = no filter (default). 0.2 = keep only the top 20% by score. variant_filter_min_score : float Minimum variant score (0–1) to keep. Default 0.0 (no filter).

Source code in mrnavax/sc_rna_pipeline.py

def run_pipeline(
    expression_path: str | Path,
    variants_path: str | Path,
    proteins_path: str | Path,
    *,
    hla: Iterable[str],
    tumor_marker_genes: list[str] | None = None,
    n_clusters: int = 4,
    peptide_lengths: tuple[int, ...] = (9, 10, 11),
    variant_filter_top_fraction: float = 1.0,
    variant_filter_min_score: float = 0.0,
    embedding_model: str = "tfidf-svd",
    embedding_dim: int = 32,
) -> PipelineReport:
    """Run the full pipeline. Returns a structured report.

    Parameters
    ----------
    variant_filter_top_fraction : float
        Fraction of variants to keep after pathogenicity scoring. 1.0 = no
        filter (default). 0.2 = keep only the top 20% by score.
    variant_filter_min_score : float
        Minimum variant score (0–1) to keep. Default 0.0 (no filter).
    """
    cell_ids, gene_names, matrix = load_expression(expression_path)
    variants = load_variants(variants_path)
    proteins_detailed = load_protein_fasta_detailed(proteins_path)
    proteins = {gene: seq for gene, (_uid, seq) in proteins_detailed.items()}
    uniprot_ids = {gene: uid for gene, (uid, _seq) in proteins_detailed.items() if uid}

    # Snapshot original count before filtering.
    n_variants_input_orig = len(variants)

    # Capture whether ANY of the original input variants carried DNA coords.
    # This drives the AVI note in the user-facing report. We must capture
    # this BEFORE the filter block reassigns `variants` to the kept set.
    had_dna_coords = any(v.chrom is not None for v in variants)

    # Optional variant pre-filter (AlphaMissense-style, with AlphaMissense
    # itself plugged in when the predictions TSV is available).
    filtered_count = n_variants_input_orig
    filter_scores: dict[str, float] = {}
    am_lookup_fn = None
    am_active = False
    # Only attempt to load the real AlphaMissense index when filtering is
    # requested AND the index is already cached. We don't trigger the
    # 15-min build during a normal run because that's a heavy operation
    # the user should initiate explicitly. Users with a pre-built cache
    # (created via ``python -m mrnavax.alphamissense_integration``
    # or a prior ``load_index()`` call) get the real scores automatically.
    if variant_filter_top_fraction < 1.0 or variant_filter_min_score > 0.0:
        try:
            from .alphamissense_integration import (
                CACHE_FILE,
                load_index,
                lookup,
            )

            if CACHE_FILE.exists():
                load_index()
                am_lookup_fn = lookup
                am_active = True
        except Exception:
            am_lookup_fn = None

    if variant_filter_top_fraction < 1.0 or variant_filter_min_score > 0.0:
        from .variant_scorer import filter_variants

        v_dicts = [
            {
                "gene": v.gene,
                "position": v.position,
                "wt_aa": v.wt_aa,
                "mut_aa": v.mut_aa,
                # DNA-level coordinates for AlphaGenome Atlas AVI lookup.
                # Both AlphaMissense (coding) and AVI (regulatory) flow
                # through the same score_variant() entry point now.
                "chrom": v.chrom,
                "ref_dna": v.ref_dna,
                "alt_dna": v.alt_dna,
            }
            for v in variants
        ]
        # Build the AVI lookup. Use the stdlib mock by default (CI +
        # offline use); users with ALPHAGENOME_API_KEY set + the
        # [variant-alphagenome] extra installed get the real Atlas
        # adapter. Selected at the backend-selector level, not here.
        avi_lookup_fn = None
        if any(v.chrom is not None for v in variants):
            try:
                from .alphagenome_integration import (
                    select_regulatory_scorer,
                )
                avi_lookup_fn = select_regulatory_scorer().score_variant
            except Exception:
                avi_lookup_fn = None

        # Build the conservation (PhyloP46way) lookup. Mock by default;
        # users with the [variant-conservation] extra installed get the
        # UCSC REST adapter. Selected at the backend-selector level.
        conservation_lookup_fn = None
        if any(v.chrom is not None for v in variants):
            try:
                from .conservation import select_conservation_lookup
                conservation_lookup_fn = select_conservation_lookup().lookup
            except Exception:
                conservation_lookup_fn = None

        scored = filter_variants(
            v_dicts,
            top_fraction=variant_filter_top_fraction,
            min_score=variant_filter_min_score,
            protein_lengths={g: len(p) for g, p in proteins.items()},
            protein_sequences=proteins,
            uniprot_ids=uniprot_ids,
            am_lookup=am_lookup_fn,
            avi_lookup=avi_lookup_fn,
            conservation_lookup=conservation_lookup_fn,
            strict=False,
        )
        keep_keys = {(s.gene, s.position, s.wt_aa, s.mut_aa) for s in scored}
        variants = [v for v in variants if (v.gene, v.position, v.wt_aa, v.mut_aa) in keep_keys]
        filtered_count = len(variants)
        filter_scores = {
            f"{s.gene}.{s.position}{s.wt_aa}>{s.mut_aa}": s.normalized_score for s in scored
        }
    elif any(v.chrom is not None for v in variants):
        # No filter requested but variants carry DNA coordinates —
        # still wire the AVI lookup and populate variant_scores so the
        # user sees AVI for non-coding regulatory variants in the report.
        from .alphagenome_integration import select_regulatory_scorer as _sel
        from .conservation import select_conservation_lookup as _sel_cons
        from .variant_scorer import score_variant as _score_one

        _avi_lookup_fn = _sel().score_variant
        _cons_lookup_fn = _sel_cons().lookup
        for v in variants:
            try:
                _r = _score_one(
                    gene=v.gene,
                    position=v.position,
                    wt_aa=v.wt_aa,
                    mut_aa=v.mut_aa,
                    chrom=v.chrom,
                    ref_dna=v.ref_dna,
                    alt_dna=v.alt_dna,
                    protein_length=len(proteins.get(v.gene, "")) or None,
                    uniprot_id=uniprot_ids.get(v.gene),
                    am_lookup=am_lookup_fn,
                    avi_lookup=_avi_lookup_fn,
                    conservation_lookup=_cons_lookup_fn,
                )
            except Exception:
                _r = None
            if _r is not None:
                filter_scores[
                    f"{_r.gene}.{_r.position}{_r.wt_aa}>{_r.mut_aa}"
                ] = _r.normalized_score

    # Build the user-facing note: AM / AVI status + tumor-marker warning
    notes: list[str] = []
    if am_active:
        notes.append("AlphaMissense pathogenicity scores used in variant filtering.")
    elif variant_filter_top_fraction < 1.0 or variant_filter_min_score > 0.0:
        notes.append(
            "AlphaMissense predictions not found — heuristic BLOSUM62 + driver "
            "gene + structural score used. To enable AlphaMissense, download "
            "the predictions TSV from "
            "https://storage.googleapis.com/dm_alphamissense/ "
            "to ~/.cache/mrnavax/AlphaMissense_hg38.tsv "
            "(licensed CC BY-NC-SA 4.0, non-commercial)."
        )
    # AVI note: mention if any variant carried DNA coordinates that
    # were scored via the AlphaGenome Atlas lookup.
    if had_dna_coords:
        notes.append(
            "AlphaGenome Atlas AVI scores used for non-coding regulatory "
            "variants (Avsec et al. Nature 2026). Coding-region variants "
            "are still scored by AlphaMissense (when available)."
        )

    labels, cluster_names = cluster_with_scanpy(matrix, cell_ids, gene_names, seed=42)
    labels = [int(x) for x in labels]  # numpy strings or ints → uniform Python ints

    # Marker score = mean expression of tumor markers per cluster
    marker_idx = [gene_names.index(g) for g in (tumor_marker_genes or []) if g in gene_names]
    cluster_scores: dict[int, float] = {}
    for c in set(labels):
        cells = [matrix[i] for i, lbl in enumerate(labels) if lbl == c]
        if marker_idx and cells:
            cluster_scores[c] = sum(
                sum(row[i] for i in marker_idx) / len(marker_idx) for row in cells
            ) / len(cells)
        else:
            cluster_scores[c] = float(len(cells))
    tumor_cluster = max(cluster_scores, key=lambda k: cluster_scores[k]) if cluster_scores else 0

    # Build peptide candidates for tumor-cluster-expressed variants
    gene_expr_by_cluster = _gene_means_per_cluster(matrix, labels, gene_names)
    tumor_expressed_genes = {
        g for g, v in gene_expr_by_cluster.get(tumor_cluster, {}).items() if v > 0.5
    }

    peptides: list[TumorPeptide] = []
    for v in variants:
        if v.gene not in proteins:
            continue
        if v.gene not in tumor_expressed_genes:
            continue
        for pep in mutant_peptides(proteins[v.gene], v.position, v.mut_aa, lengths=peptide_lengths):
            peptides.append(
                TumorPeptide(
                    cell_cluster=tumor_cluster,
                    gene=v.gene,
                    position=v.position,
                    wt_aa=v.wt_aa,
                    mut_aa=v.mut_aa,
                    peptide=pep,
                    length=len(pep),
                    cluster_marker_score=cluster_scores.get(tumor_cluster, 0.0),
                )
            )
    if not marker_idx:
        notes.append(
            "WARNING: no tumor-marker genes supplied; the largest cluster was "
            "assumed to be tumor. Pass --tumor-markers GENE1,GENE2 for "
            "accurate selection. The downstream peptide list may contain "
            "many false positives."
        )
    note = " | ".join(notes)

    return PipelineReport(
        n_cells=len(cell_ids),
        n_genes=len(gene_names),
        cluster_labels=labels,
        cluster_marker_scores=cluster_scores,
        tumor_cluster=tumor_cluster,
        n_variants_input=n_variants_input_orig,
        n_variants_after_filter=filtered_count,
        variant_scores=filter_scores,
        n_candidate_peptides=len(peptides),
        peptides=peptides,
        note=note,
    )

mrnavax.foundation_embedder

Foundation-model embedder for scRNA-seq cells.

Provides a deterministic, stdlib-only TF-IDF + truncated SVD cell embedder that approximates what scGPT / Geneformer / UNI-RNA do for cell-type identification — without requiring GPU or external model downloads.

If scgpt is installed (optional extra), the real model is used instead via the plug point in embed_cells.

References

Cui et al., scGPT: toward building a foundation model for single-cell multi-omics. Nat Methods 21, 1480–1491 (2024).

embed_cells(matrix, *, model='tfidf-svd', n_components=32, seed=42)

Embed cells into a fixed-dimensional vector space.

Parameters

model : str - "tfidf-svd" (default, stdlib-only): deterministic TF-IDF + truncated SVD. - "scgpt": use the real scGPT foundation model. Requires the scgpt package and a downloaded checkpoint. - "identity": per-cell mean expression (one dim).

Source code in mrnavax/foundation_embedder.py

def embed_cells(
    matrix: list[list[float]],
    *,
    model: str = "tfidf-svd",
    n_components: int = 32,
    seed: int = 42,
) -> list[list[float]]:
    """Embed cells into a fixed-dimensional vector space.

    Parameters
    ----------
    model : str
        - ``"tfidf-svd"`` (default, stdlib-only): deterministic TF-IDF + truncated SVD.
        - ``"scgpt"``: use the real scGPT foundation model. Requires the
          ``scgpt`` package and a downloaded checkpoint.
        - ``"identity"``: per-cell mean expression (one dim).
    """
    if model == "scgpt":
        # Use the real scGPT foundation model via HuggingFace weights.
        # Falls back with a clear error if weights aren't downloaded.
        try:
            from .scgpt_integration import embed_with_scgpt, scgpt_available
        except ImportError as e:
            raise RuntimeError(
                f"scGPT backend requested but scgpt_integration not importable: {e}. "
                f"Use --model tfidf-svd (default) or pip install torch transformers."
            )
        if not scgpt_available():
            raise FileNotFoundError(
                "scGPT weights not found at ~/.cache/mrnavax/. "
                "Download from https://huggingface.co/perturblab/scgpt-human "
                "(best_model.pt, vocab.json, args.json)."
            )
        return embed_with_scgpt(matrix, max_cells=64)
    if model == "identity":
        return [[sum(row) / len(row) if row else 0.0] for row in matrix]
    if model == "tfidf-svd":
        return _tfidf_svd_embed(matrix, n_components=n_components, seed=seed)
    raise ValueError(f"unknown model: {model!r}")

Patient-trial matching

mrnavax.trial_matcher

TrialGPT-style patient-to-trial matcher.

End-to-end stub of the three-stage pipeline from Jin et al. Nature Communications 15, 9074 (2024):

1. TrialGPT-Retrieval — keyword generation + candidate filtering.
2. TrialGPT-Matching — criterion-level eligibility with explanations.
3. TrialGPT-Ranking — aggregate criterion scores to a trial-level rank.

This is a demonstrator (no fine-tuned retriever, no production medical use). The LLM does the matching + ranking; the retrieval stage is a simple keyword overlap so the demo is fast and offline-tunable.

References

TrialGPT: https://www.nature.com/articles/s41467-024-53081-z Biomarker LLM match: https://www.nature.com/articles/s41746-025-01673-4

match(patient_text, trials, *, top_k=20, backend=None, retriever='keyword', matcher='auto')

End-to-end pipeline: retrieve → match → rank.

Parameters

retriever : str - "keyword": simple keyword overlap (default, fastest) - "dense": TF-IDF + biomedical synonym expansion (MedCPT-style) - "auto": dense if available, else keyword matcher : str - "auto" (default): use TrialGPT-style per-criterion LLM matcher when OPENAI_API_KEY is set, else keyword fallback - "trialgpt": force the per-criterion LLM matcher - "keyword": force keyword-only matching (no LLM)

Source code in mrnavax/trial_matcher.py

def match(
    patient_text: str,
    trials: list[Trial],
    *,
    top_k: int = 20,
    backend: str | None = None,
    retriever: str = "keyword",
    matcher: str = "auto",
) -> tuple[list[RankedTrial], list[dict]]:
    """End-to-end pipeline: retrieve → match → rank.

    Parameters
    ----------
    retriever : str
        - ``"keyword"``: simple keyword overlap (default, fastest)
        - ``"dense"``: TF-IDF + biomedical synonym expansion (MedCPT-style)
        - ``"auto"``: dense if available, else keyword
    matcher : str
        - ``"auto"`` (default): use ``TrialGPT-style`` per-criterion LLM
          matcher when ``OPENAI_API_KEY`` is set, else keyword fallback
        - ``"trialgpt"``: force the per-criterion LLM matcher
        - ``"keyword"``: force keyword-only matching (no LLM)
    """
    if retriever == "dense" or retriever == "auto":
        # Try MedCPT first (real semantic encoder), then TF-IDF dense,
        # then keyword fallback.
        from .medcpt_integration import medcpt_available

        medcpt_ok, _ = medcpt_available()
        if medcpt_ok:
            try:
                from .medcpt_integration import retrieve_medcpt

                trial_dicts = [
                    {
                        "title": t.title,
                        "condition": t.condition,
                        "inclusion": list(t.inclusion),
                        "exclusion": list(t.exclusion),
                        "biomarkers": list(t.biomarkers),
                    }
                    for t in trials
                ]
                # Build a single concatenated article string per trial
                article_strings = [
                    " ".join([d["title"], d["condition"]] + d["inclusion"] + d["biomarkers"])
                    for d in trial_dicts
                ]
                scores = retrieve_medcpt(patient_text, article_strings)
                ranked = sorted(zip(trials, scores), key=lambda x: -x[1])[:top_k]
                candidates = [t for t, _ in ranked]
            except Exception:
                # MedCPT loaded but encode failed — fall through
                candidates = _dense_fallback(
                    patient_text,
                    trials,
                    top_k,
                    from_medcpt=True,
                )
        else:
            candidates = _dense_fallback(patient_text, trials, top_k)
    else:
        candidates = retrieve_candidates(patient_text, trials, top_k=top_k)
    out: list[RankedTrial] = []
    debug: list[dict] = []
    # Determine matcher backend
    use_trialgpt = False
    use_simicl = False
    if matcher == "trialgpt":
        use_trialgpt = True
    elif matcher == "trialgpt-simicl":
        use_trialgpt = True
        use_simicl = True
    elif matcher == "auto":
        # Use TrialGPT when an OpenAI key is present or when caller
        # explicitly forced an LLM backend. Otherwise fall back to
        # keyword matching.
        import os

        use_trialgpt = bool(os.environ.get("OPENAI_API_KEY") or (backend and backend != "mock"))

    for t in candidates:
        if use_trialgpt:
            from .trial_llm import score_trial_with_llm

            try:
                llm_result = score_trial_with_llm(
                    patient_text,
                    t.nct_id,
                    t.title,
                    list(t.inclusion),
                    list(t.exclusion),
                    backend=backend,
                    use_simicl=use_simicl,
                )
                match_result = llm_result.to_dict()
                # Trim to the shape that rank() expects
                match_result = {
                    "inclusion": llm_result.to_dict()["inclusion"],
                    "exclusion": llm_result.to_dict()["exclusion"],
                }
                notes = list(llm_result.notes)
            except Exception as e:
                # Fall back to keyword matcher on failure
                match_result, fallback_notes = _match_one(t, patient_text, backend=backend)
                notes = ["trialgpt-fallback"] + list(fallback_notes) + [str(e)]
        else:
            match_result, notes = _match_one(t, patient_text, backend=backend)
        ranked = rank(t, match_result)
        out.append(ranked)
        debug.append({"nct": t.nct_id, "match": match_result, "notes": notes})
    out.sort(key=lambda r: -r.score)
    return out, debug

retrieve_candidates(patient_text, trials, *, top_k=20)

Keyword-overlap retrieval. Replaceable with a dense encoder in prod.

Source code in mrnavax/trial_matcher.py

def retrieve_candidates(
    patient_text: str,
    trials: list[Trial],
    *,
    top_k: int = 20,
) -> list[Trial]:
    """Keyword-overlap retrieval. Replaceable with a dense encoder in prod."""
    patient_tokens = set(_tokenize(patient_text))
    scored: list[tuple[float, Trial]] = []
    for t in trials:
        text = " ".join([t.title, t.condition] + t.inclusion + t.biomarkers)
        tokens = set(_tokenize(text))
        if not tokens:
            continue
        overlap = len(patient_tokens & tokens) / (len(tokens) ** 0.5)
        scored.append((overlap, t))
    scored.sort(key=lambda x: -x[0])
    return [t for _, t in scored[:top_k]]

mrnavax.trial_llm

TrialGPT-style per-criterion LLM matching.

Implements the TrialGPT-Matching approach from Jin et al. (Nature Communications 2024): per-criterion LLM reasoning where each inclusion and exclusion criterion is judged independently against the patient summary, producing a granular eligibility assessment.

References

Jin, Qiao, et al. "Matching patients to clinical trials with large language models." Nature Communications 15 (2024): 9074. DOI: 10.1038/s41467-024-53081-z

Reported: 87.3% accuracy on 1,015 patient-criterion pairs, close to expert performance.

Usage

With an OPENAI_API_KEY set, the per-criterion matcher is used. Without one, the system falls back to keyword overlap (handled in trial_matcher.match()).

The :func:score_trial_with_llm function calls llm_json once per trial with all criteria in a single prompt — cheaper than per-criterion calls while preserving the granular per-criterion verdict structure.

CriterionVerdict dataclass

One criterion's LLM verdict.

Attributes

criterion The criterion text. verdict One of "met", "unmet", "uncertain". evidence Short sentence explaining the verdict (from the patient summary).

Source code in mrnavax/trial_llm.py

@dataclass
class CriterionVerdict:
    """One criterion's LLM verdict.

    Attributes
    ----------
    criterion
        The criterion text.
    verdict
        One of ``"met"``, ``"unmet"``, ``"uncertain"``.
    evidence
        Short sentence explaining the verdict (from the patient summary).
    """

    criterion: str
    verdict: str
    evidence: str = ""

    def to_dict(self) -> dict:
        return {"criterion": self.criterion, "verdict": self.verdict, "evidence": self.evidence}

TrialMatchResult dataclass

Aggregate result of TrialGPT-style matching for one trial.

Attributes

nct_id Trial identifier. inclusion_verdicts Per-inclusion-criterion verdicts. exclusion_verdicts Per-exclusion-criterion verdicts. n_met_inclusion Count of inclusion criteria labelled "met". n_total_inclusion Total inclusion criteria. n_unmet_exclusion Count of exclusion criteria labelled "unmet" (good — patient does NOT have the exclusion, so still eligible). n_total_exclusion Total exclusion criteria. eligibility_score Aggregate score, 0..1. Computed as: (n_met_inclusion / n_total_inclusion) * (n_unmet_exclusion / n_total_exclusion) when both are > 0; falls back to one of the two when the other is empty. notes Free-form debug notes (e.g., "schema-mismatch fallback"). raw_llm_response The raw JSON the LLM returned (for audit / debugging).

Source code in mrnavax/trial_llm.py

@dataclass
class TrialMatchResult:
    """Aggregate result of TrialGPT-style matching for one trial.

    Attributes
    ----------
    nct_id
        Trial identifier.
    inclusion_verdicts
        Per-inclusion-criterion verdicts.
    exclusion_verdicts
        Per-exclusion-criterion verdicts.
    n_met_inclusion
        Count of inclusion criteria labelled ``"met"``.
    n_total_inclusion
        Total inclusion criteria.
    n_unmet_exclusion
        Count of exclusion criteria labelled ``"unmet"`` (good — patient
        does NOT have the exclusion, so still eligible).
    n_total_exclusion
        Total exclusion criteria.
    eligibility_score
        Aggregate score, 0..1. Computed as:
        ``(n_met_inclusion / n_total_inclusion) *
        (n_unmet_exclusion / n_total_exclusion)``
        when both are > 0; falls back to one of the two when the other
        is empty.
    notes
        Free-form debug notes (e.g., "schema-mismatch fallback").
    raw_llm_response
        The raw JSON the LLM returned (for audit / debugging).
    """

    nct_id: str
    inclusion_verdicts: list[CriterionVerdict] = field(default_factory=list)
    exclusion_verdicts: list[CriterionVerdict] = field(default_factory=list)
    n_met_inclusion: int = 0
    n_total_inclusion: int = 0
    n_unmet_exclusion: int = 0
    n_total_exclusion: int = 0
    eligibility_score: float = 0.0
    notes: list[str] = field(default_factory=list)
    raw_llm_response: str = ""

    def to_dict(self) -> dict:
        return {
            "nct_id": self.nct_id,
            "n_met_inclusion": self.n_met_inclusion,
            "n_total_inclusion": self.n_total_inclusion,
            "n_unmet_exclusion": self.n_unmet_exclusion,
            "n_total_exclusion": self.n_total_exclusion,
            "eligibility_score": round(self.eligibility_score, 3),
            "notes": list(self.notes),
            "inclusion": [v.to_dict() for v in self.inclusion_verdicts],
            "exclusion": [v.to_dict() for v in self.exclusion_verdicts],
        }

build_match_prompt(patient_text, nct_id, title, inclusion, exclusion)

Construct the TrialGPT-style per-criterion prompt.

Source code in mrnavax/trial_llm.py

def build_match_prompt(
    patient_text: str,
    nct_id: str,
    title: str,
    inclusion: list[str],
    exclusion: list[str],
) -> str:
    """Construct the TrialGPT-style per-criterion prompt."""
    inc_lines = "\n".join(f"- {c}" for c in inclusion) if inclusion else "(none)"
    exc_lines = "\n".join(f"- {c}" for c in exclusion) if exclusion else "(none)"
    return MATCH_PROMPT_TEMPLATE.format(
        patient=patient_text.strip(),
        nct=nct_id,
        title=title,
        inclusion=inc_lines,
        exclusion=exc_lines,
    )

llm_matching_available()

True iff an LLM backend can run without an explicit OPENAI_API_KEY.

The mock backend is always available; the openai backend requires the env var. Used by backends.py to skip the LLM check when no key is present and we want to avoid a real API call.

Source code in mrnavax/trial_llm.py

def llm_matching_available() -> bool:
    """True iff an LLM backend can run without an explicit OPENAI_API_KEY.

    The mock backend is always available; the openai backend requires
    the env var. Used by ``backends.py`` to skip the LLM check when
    no key is present and we want to avoid a real API call.
    """
    if os.environ.get("OPENAI_API_KEY"):
        return True
    # Mock backend always available for offline runs.
    if os.environ.get("MRNA_AI_LLM_BACKEND", "auto").lower() in {"mock", "auto"}:
        return True
    return False

score_trial_with_llm(patient_text, nct_id, title, inclusion, exclusion, *, backend=None, demo_store=None, use_simicl=None)

Run TrialGPT-style per-criterion LLM matching.

Parameters

patient_text Free-text patient summary (diagnosis, stage, biomarkers, prior therapies, performance status). nct_id Trial NCT identifier. title Trial title. inclusion List of inclusion criteria. exclusion List of exclusion criteria. backend Optional LLM backend override. None = auto-detect from OPENAI_API_KEY (uses real OpenAI) or fall back to mock. demo_store Optional :class:~mrnavax.trial_similar.DemoStore. When None, loads the bundled store from examples/simicl_demos.json via :func:~mrnavax.trial_similar.load_default_demo_store. Pass an explicit DemoStore(demos=[]) to disable. use_simicl Override $MRNA_AI_SIMICL_ENABLED. When True, the top-K demos (per $MRNA_AI_SIMICL_TOPK, default 32) are injected into the prompt as few-shot examples before asking the LLM for verdicts. When False, behaves as plain TrialGPT.

Returns

TrialMatchResult with per-criterion verdicts and aggregate score. On any failure, all verdicts are set to "uncertain" and a note is added.

Notes

Sim-ICL integration: this function implements the Sim-ICL demonstration-selection strategy from Fung et al. 2026 (Genome Biology, in press). The few-shot examples are chosen by TF-IDF cosine similarity between (patient+trial) query and the demo store, rather than random sampling. This matches the paper's finding that sequence-similar demonstrations yield competitive performance with protein-LM classifiers in low-shot regimes.

Source code in mrnavax/trial_llm.py

def score_trial_with_llm(
    patient_text: str,
    nct_id: str,
    title: str,
    inclusion: list[str],
    exclusion: list[str],
    *,
    backend: str | None = None,
    demo_store: Any | None = None,
    use_simicl: bool | None = None,
) -> TrialMatchResult:
    """Run TrialGPT-style per-criterion LLM matching.

    Parameters
    ----------
    patient_text
        Free-text patient summary (diagnosis, stage, biomarkers,
        prior therapies, performance status).
    nct_id
        Trial NCT identifier.
    title
        Trial title.
    inclusion
        List of inclusion criteria.
    exclusion
        List of exclusion criteria.
    backend
        Optional LLM backend override. ``None`` = auto-detect from
        ``OPENAI_API_KEY`` (uses real OpenAI) or fall back to ``mock``.
    demo_store
        Optional :class:`~mrnavax.trial_similar.DemoStore`. When
        ``None``, loads the bundled store from
        ``examples/simicl_demos.json`` via
        :func:`~mrnavax.trial_similar.load_default_demo_store`.
        Pass an explicit ``DemoStore(demos=[])`` to disable.
    use_simicl
        Override ``$MRNA_AI_SIMICL_ENABLED``. When ``True``, the top-K
        demos (per ``$MRNA_AI_SIMICL_TOPK``, default 32) are injected
        into the prompt as few-shot examples before asking the LLM
        for verdicts. When ``False``, behaves as plain TrialGPT.

    Returns
    -------
    TrialMatchResult with per-criterion verdicts and aggregate score.
    On any failure, all verdicts are set to ``"uncertain"`` and a note
    is added.

    Notes
    -----
    Sim-ICL integration: this function implements the **Sim-ICL**
    demonstration-selection strategy from Fung et al. 2026 (Genome
    Biology, in press). The few-shot examples are chosen by TF-IDF
    cosine similarity between (patient+trial) query and the demo
    store, rather than random sampling. This matches the paper's
    finding that sequence-similar demonstrations yield competitive
    performance with protein-LM classifiers in low-shot regimes.
    """
    from .llm import llm_json  # late import to avoid circular deps
    from .trial_similar import (
        _get_enabled,
        build_simicl_prompt,
        load_default_demo_store,
    )

    # Sim-ICL: resolve config + load demos
    if demo_store is None:
        demo_store = load_default_demo_store()
    simicl_on = _get_enabled() if use_simicl is None else bool(use_simicl)

    prompt = build_match_prompt(patient_text, nct_id, title, inclusion, exclusion)
    notes: list[str] = []
    raw = ""
    parsed: dict[str, Any] = {}

    if simicl_on and len(demo_store) > 0:
        query_text = "\n".join(
            [patient_text, nct_id, title, " ".join(inclusion), " ".join(exclusion)]
        )
        top_demos = demo_store.rank(query_text)
        prompt = build_simicl_prompt(
            patient_text,
            nct_id,
            title,
            inclusion,
            exclusion,
            top_demos,
            base_prompt=prompt,
        )
        notes.append(f"simicl-k{len(top_demos)}")
        # Record which demos were used, for auditability
        notes.append("simicl-demo-ids=" + ",".join(d.demo_id for d in top_demos))

    try:
        parsed = llm_json(prompt, backend=backend)
        raw = json.dumps(parsed)
    except Exception as e:
        notes.append(f"llm-fallback: {e}")
        parsed = {}

    inclusion_v = _parse_criteria(parsed.get("inclusion"), inclusion)
    exclusion_v = _parse_criteria(parsed.get("exclusion"), exclusion)

    n_met_inc = sum(1 for v in inclusion_v if v.verdict == VERDICT_MET)
    n_total_inc = len(inclusion_v)
    n_unmet_exc = sum(1 for v in exclusion_v if v.verdict == VERDICT_UNMET)
    n_total_exc = len(exclusion_v)
    score = _eligibility_score(n_met_inc, n_total_inc, n_unmet_exc, n_total_exc)
    if not notes and (n_met_inc + n_unmet_exc == 0):
        notes.append("all-uncertain (LLM returned no confident verdicts)")

    return TrialMatchResult(
        nct_id=nct_id,
        inclusion_verdicts=inclusion_v,
        exclusion_verdicts=exclusion_v,
        n_met_inclusion=n_met_inc,
        n_total_inclusion=n_total_inc,
        n_unmet_exclusion=n_unmet_exc,
        n_total_exclusion=n_total_exc,
        eligibility_score=score,
        notes=notes,
        raw_llm_response=raw,
    )

Manufacturing checks

mrnavax.manufacturability

mRNA manufacturability checks.

Wet-lab bridge between computational sequence design and what's actually synthesizable at production scale. Each check returns a :class:CheckResult with a pass_ flag, score (0-1, higher = more manufacturable), severity (info/warn/error), and a human-readable explanation.

Checks implemented

1. poly_a_runs — runs of ≥5 consecutive As in the DNA template destabilize the plasmid during IVT. Counts occurrences and the longest run. Severity scales with run length.
2. gc_5prime_hairpin — GC-rich stems (≥70% GC over 30+ nt) at the 5' UTR or CDS start block ribosome scanning. Uses a sliding window over the first 60 nt.
3. are_motif — AU-rich elements (ARE) in the 3' UTR trigger mRNA decay via TTP / HuR. Canonical nonamer UUAUUUAUU and pentamer AUUUA repeats. The CDS itself may contain pentamers harmlessly (introns mostly aren't in IVT-mRNA); we flag the nonamer explicitly in the 3' UTR.
4. kozak_strength — measures match to the mammalian Kozak consensus (GCC)GCC(A/G)CCATGG. Returns a 0-1 score.
5. stop_context — termination efficiency depends on the stop codon identity and the +4 base. Per the literature, readthrough order is TAA < TAG < TGA; the most efficient terminator is TGA-T or TAA-T.
6. hidden_stops — internal in-frame stop codons (should never occur in a designed CDS but can sneak in via codon-table drift).
7. gc_window_uniformity — extreme local GC variation (>25% stddev across 30-nt windows) signals poor IVT yield and ribosomal stalling.
8. cpg_suppression — long CpG-free stretches trigger silencing in some contexts; extremely CpG-rich regions (>15% CpG in a window) cause immune activation. Returns a balanced score.

Reference

• Holtkamp et al. (2006) "Modification of antigen-encoding RNA increases stability, translational efficacy, and T-cell stimulatory capacity of dendritic cells." Blood 108.
• Kozak (1986) "Point mutations define a sequence flanking the AUG initiator codon that modulates translation by eukaryotic ribosomes." Cell 44.
• Chen & Shyu (1995) "AU-rich elements: characterization and importance in mRNA degradation." Trends Biochem Sci 20.

CheckResult dataclass

Outcome of one manufacturability check.

Source code in mrnavax/manufacturability.py

@dataclass
class CheckResult:
    """Outcome of one manufacturability check."""

    name: str
    pass_: bool
    score: float  # 0..1, higher = better
    severity: str  # "info" | "warn" | "error"
    summary: str
    details: dict = field(default_factory=dict)

    def to_dict(self) -> dict:
        return asdict(self)

ManufacturabilityReport dataclass

Aggregated report from all manufacturability checks.

Source code in mrnavax/manufacturability.py

@dataclass
class ManufacturabilityReport:
    """Aggregated report from all manufacturability checks."""

    overall_score: float  # 0..1
    n_pass: int
    n_warn: int
    n_error: int
    checks: list[CheckResult]

    def to_dict(self) -> dict:
        return {
            "overall_score": self.overall_score,
            "n_pass": self.n_pass,
            "n_warn": self.n_warn,
            "n_error": self.n_error,
            "checks": [c.to_dict() for c in self.checks],
        }

check_are_motif(utr3)

AU-rich elements in the 3' UTR trigger mRNA decay.

Counts pentamers (AUUUA) and nonamers (UUAUUUAUU) in the 3' UTR. Pentamers alone are weak signals; ≥2 nonamers is severe.

Source code in mrnavax/manufacturability.py

def check_are_motif(utr3: str) -> CheckResult:
    """AU-rich elements in the 3' UTR trigger mRNA decay.

    Counts pentamers (AUUUA) and nonamers (UUAUUUAUU) in the 3' UTR.
    Pentamers alone are weak signals; ≥2 nonamers is severe.
    """
    rna = utr3.upper().replace("T", "U")
    n_pent = len(re.findall(ARE_PENTAMER, rna))
    n_non = len(re.findall(ARE_NONAMER, rna))
    if n_non >= 2:
        sev, ok, score = "error", False, 0.1
    elif n_non == 1:
        sev, ok, score = "warn", False, 0.5
    elif n_pent >= 5:
        sev, ok, score = "warn", False, 0.6
    else:
        sev, ok, score = "info", True, 1.0
    return CheckResult(
        name="are_motif",
        pass_=ok,
        score=score,
        severity=sev,
        summary=(
            f"3' UTR: {n_pent} pentamer(s), {n_non} nonamer(s) "
            f"(pentamer = {ARE_PENTAMER}; nonamer = {ARE_NONAMER})"
        ),
        details={"n_pentamer": n_pent, "n_nonamer": n_non},
    )

check_cpg_balance(dna)

CpG density per sliding window.

Suppressed (<0.5%) → silencing risk in some contexts. Excessive (>15%) → immune activation (TLR9). Best: a balanced profile around 1-5%.

Source code in mrnavax/manufacturability.py

def check_cpg_balance(dna: str) -> CheckResult:
    """CpG density per sliding window.

    Suppressed (<0.5%) → silencing risk in some contexts.
    Excessive (>15%) → immune activation (TLR9).
    Best: a balanced profile around 1-5%.
    """
    dna = dna.upper().replace("U", "T")
    if len(dna) < CPG_WINDOW:
        return CheckResult(
            name="cpg_balance",
            pass_=True,
            score=1.0,
            severity="info",
            summary="sequence <50 nt; skipping CpG check",
            details={},
        )
    windows = []
    for i in range(0, len(dna) - CPG_WINDOW + 1, CPG_WINDOW // 2):
        win = dna[i : i + CPG_WINDOW]
        cpg = sum(1 for j in range(len(win) - 1) if win[j : j + 2] == "CG")
        windows.append(cpg / len(win))
    n = len(windows)
    if n == 0:
        return CheckResult(
            name="cpg_balance",
            pass_=True,
            score=1.0,
            severity="info",
            summary="no windows",
            details={},
        )
    mean = sum(windows) / n
    if mean < CPG_LOW:
        sev, ok, score = "warn", False, 0.6
    elif mean > CPG_HIGH:
        sev, ok, score = "warn", False, 0.5
    else:
        sev, ok, score = "info", True, 1.0
    return CheckResult(
        name="cpg_balance",
        pass_=ok,
        score=score,
        severity=sev,
        summary=f"CpG density = {mean:.3f} (target 0.5%-15%)",
        details={"mean_cpg": mean, "n_windows": n},
    )

check_gc_5prime_hairpin(dna)

Sliding GC% over the first 60 nt; flag windows ≥70% GC.

The 5' UTR's first 30-40 nt set the tone for ribosome scanning. Local GC-rich stems cause the 40S subunit to stall.

Source code in mrnavax/manufacturability.py

def check_gc_5prime_hairpin(dna: str) -> CheckResult:
    """Sliding GC% over the first 60 nt; flag windows ≥70% GC.

    The 5' UTR's first 30-40 nt set the tone for ribosome scanning.
    Local GC-rich stems cause the 40S subunit to stall.
    """
    dna = dna.upper().replace("U", "T")
    head = dna[:60]
    if len(head) < GC_HAIRPIN_WINDOW:
        return CheckResult(
            name="gc_5prime_hairpin",
            pass_=True,
            score=1.0,
            severity="info",
            summary="sequence <60 nt; skipping 5' hairpin check",
            details={},
        )
    worst_gc = 0.0
    worst_pos = 0
    for i in range(len(head) - GC_HAIRPIN_WINDOW + 1):
        win = head[i : i + GC_HAIRPIN_WINDOW]
        gc = (win.count("G") + win.count("C")) / len(win)
        if gc > worst_gc:
            worst_gc = gc
            worst_pos = i
    if worst_gc >= 0.80:
        sev, ok, score = "error", False, 0.1
    elif worst_gc >= GC_HAIRPIN_THRESHOLD:
        sev, ok, score = "warn", False, 0.5
    else:
        sev, ok, score = "info", True, 1.0
    return CheckResult(
        name="gc_5prime_hairpin",
        pass_=ok,
        score=score,
        severity=sev,
        summary=(f"max GC in first 60 nt = {worst_gc:.0%} at position {worst_pos + 1}"),
        details={"max_gc": worst_gc, "max_position": worst_pos + 1},
    )

check_gc_window_uniformity(cds)

Local GC% variation across the CDS.

High stddev = uneven IVT yield, more secondary structure, more ribosomal stalling. Threshold: stddev > 15% flags as warn.

Source code in mrnavax/manufacturability.py

def check_gc_window_uniformity(cds: str) -> CheckResult:
    """Local GC% variation across the CDS.

    High stddev = uneven IVT yield, more secondary structure, more
    ribosomal stalling. Threshold: stddev > 15% flags as warn.
    """
    cds = cds.upper().replace("U", "T")
    if len(cds) < GC_WINDOW:
        return CheckResult(
            name="gc_window_uniformity",
            pass_=True,
            score=1.0,
            severity="info",
            summary="sequence <30 nt; skipping GC uniformity check",
            details={},
        )
    gcs: list[float] = []
    for i in range(0, len(cds) - GC_WINDOW + 1, GC_WINDOW // 2):
        win = cds[i : i + GC_WINDOW]
        gcs.append((win.count("G") + win.count("C")) / len(win))
    n = len(gcs)
    if n < 2:
        return CheckResult(
            name="gc_window_uniformity",
            pass_=True,
            score=1.0,
            severity="info",
            summary="too few windows for stddev calculation",
            details={},
        )
    mean = sum(gcs) / n
    var = sum((x - mean) ** 2 for x in gcs) / n
    stddev = var**0.5
    if stddev > 0.25:
        sev, ok, score = "error", False, 0.2
    elif stddev > GC_WINDOW_STDDEV_MAX:
        sev, ok, score = "warn", False, 0.6
    else:
        sev, ok, score = "info", True, 1.0
    return CheckResult(
        name="gc_window_uniformity",
        pass_=ok,
        score=score,
        severity=sev,
        summary=(f"GC stddev across {n} sliding windows = {stddev:.3f} (mean = {mean:.3f})"),
        details={"stddev": stddev, "mean": mean, "n_windows": n},
    )

check_hidden_stops(cds)

Internal in-frame stop codons (should be zero).

Source code in mrnavax/manufacturability.py

def check_hidden_stops(cds: str) -> CheckResult:
    """Internal in-frame stop codons (should be zero)."""
    cds = cds.upper().replace("U", "T")
    codons = [cds[i : i + 3] for i in range(0, len(cds) - 2, 3)]
    stops = {"TAA", "TAG", "TGA"}
    n_internal = sum(1 for c in codons[:-1] if c in stops)
    if n_internal == 0:
        return CheckResult(
            name="hidden_stops",
            pass_=True,
            score=1.0,
            severity="info",
            summary="no internal in-frame stops",
            details={"n_internal_stops": 0},
        )
    return CheckResult(
        name="hidden_stops",
        pass_=False,
        score=0.0,
        severity="error",
        summary=f"{n_internal} internal in-frame stop codon(s) detected",
        details={"n_internal_stops": n_internal},
    )

check_kozak_strength(utr5)

Score the 9 nt immediately upstream of ATG against Kozak consensus.

Strong: GCCRCCATGG (R = purine). We score position-by-position with weighted matches: positions -3 (R) and +4 (G) carry more weight in the original Kozak paper.

Source code in mrnavax/manufacturability.py

def check_kozak_strength(utr5: str) -> CheckResult:
    """Score the 9 nt immediately upstream of ATG against Kozak consensus.

    Strong: ``GCCRCCATGG`` (R = purine). We score position-by-position
    with weighted matches: positions -3 (R) and +4 (G) carry more weight
    in the original Kozak paper.
    """
    rna = utr5.upper().replace("T", "U")
    # We need 9 nt upstream of ATG. Convention: ATG is in the CDS;
    # caller passes the trailing 9 nt of the 5' UTR.
    if len(rna) < 9:
        return CheckResult(
            name="kozak_strength",
            pass_=True,
            score=0.5,
            severity="info",
            summary=f"only {len(rna)} nt of 5' UTR; cannot score Kozak",
            details={},
        )
    pre = rna[-9:]
    # Position-by-position match (R = A or G)
    weights = [1, 1, 1, 1.5, 1, 1, 1.5, 1, 1]  # -3 R and +4 G heavier
    expected: list[tuple[str, ...]] = [
        ("G",),
        ("C",),
        ("C",),
        ("A", "G"),  # R
        ("C",),
        ("C",),
        ("A", "G"),  # R
        ("A",),  # A in ATG
        ("T",),
    ]
    score = 0.0
    total_w = 0.0
    for i, exp in enumerate(expected):
        if i >= len(pre):
            break
        total_w += weights[i]
        if pre[i] in exp:
            score += weights[i]
    norm = score / total_w if total_w else 0.0
    if norm >= 0.85:
        sev, ok = "info", True
    elif norm >= 0.65:
        sev, ok = "info", True
    elif norm >= 0.45:
        sev, ok = "warn", False
    else:
        sev, ok = "warn", False
    return CheckResult(
        name="kozak_strength",
        pass_=ok,
        score=norm,
        severity=sev,
        summary=(
            f"5' UTR upstream of ATG ({pre[-9:]}) matches partial Kozak consensus at {norm:.0%}"
        ),
        details={"pre_atg": pre[-9:], "score_partial": norm},
    )

check_poly_a_runs(dna)

Poly-A runs ≥5 nt destabilize the DNA template.

Returns a CheckResult that fails if any run ≥7 occurs (severe) or if multiple runs ≥5 occur (moderate).

Source code in mrnavax/manufacturability.py

def check_poly_a_runs(dna: str) -> CheckResult:
    """Poly-A runs ≥5 nt destabilize the DNA template.

    Returns a CheckResult that fails if any run ≥7 occurs (severe) or
    if multiple runs ≥5 occur (moderate).
    """
    dna = dna.upper().replace("U", "T")
    longest = _longest_run(dna, "A")
    # Count all runs of length >= POLY_A_MIN_LEN
    pattern = re.compile(f"A{{{POLY_A_MIN_LEN},}}")
    runs = [m.group() for m in pattern.finditer(dna)]
    n_runs = len(runs)
    if longest >= 9:
        sev = "error"
        ok = False
        score = 0.0
    elif longest >= 7:
        sev = "error"
        ok = False
        score = 0.1
    elif n_runs >= 3:
        sev = "warn"
        ok = False
        score = 0.5
    elif longest >= POLY_A_MIN_LEN:
        sev = "warn"
        ok = True  # mild
        score = 0.7
    else:
        sev = "info"
        ok = True
        score = 1.0
    return CheckResult(
        name="poly_a_runs",
        pass_=ok,
        score=score,
        severity=sev,
        summary=(f"longest poly-A run = {longest} nt; {n_runs} run(s) of length ≥{POLY_A_MIN_LEN}"),
        details={"longest_run": longest, "n_runs": n_runs, "runs": runs},
    )

check_stop_context(dna)

Stop codon identity + +4 base. Strong = TAA or TGA, weak = TGA.

Reads as RNA: the stop codon is the last three nt of the CDS, and +4 is the first base of the 3' UTR.

Source code in mrnavax/manufacturability.py

def check_stop_context(dna: str) -> CheckResult:
    """Stop codon identity + +4 base. Strong = TAA or TGA, weak = TGA.

    Reads as RNA: the stop codon is the last three nt of the CDS, and
    ``+4`` is the first base of the 3' UTR.
    """
    dna = dna.upper().replace("U", "T")
    # Take the last 3 nt (CDS stop) and look at +4 if available
    if len(dna) < 3:
        return CheckResult(
            name="stop_context",
            pass_=True,
            score=1.0,
            severity="info",
            summary="sequence <3 nt; skipping stop check",
            details={},
        )
    stop = dna[-3:]
    follower = dna[3] if len(dna) > 3 else "T"  # default to T if missing
    base_score = STOP_CODON_BASE_SCORE.get(stop, 0.5)
    follower_mod = STOP_FOLLOWER_BONUS.get(follower, 0.0)
    # Combined: 0 = strong stop, 1 = leaky
    leakiness = max(0.0, min(1.0, base_score + follower_mod))
    score = 1.0 - leakiness
    if leakiness >= 0.7:
        sev, ok = "warn", False
    elif leakiness >= 0.5:
        sev, ok = "info", True
    else:
        sev, ok = "info", True
    return CheckResult(
        name="stop_context",
        pass_=ok,
        score=score,
        severity=sev,
        summary=(
            f"stop codon {stop} + follower {follower}; "
            f"leakiness={leakiness:.2f} (lower = stronger termination)"
        ),
        details={"stop": stop, "follower": follower, "leakiness": leakiness},
    )

score_manufacturability(cds, *, utr5='', utr3='')

Run all manufacturability checks and aggregate.

Parameters

cds Coding sequence (DNA). If it ends with a stop codon, the stop codon is included in the analysis. utr5 5' UTR sequence (DNA). Used for Kozak scoring — only the last 9 nt upstream of ATG are inspected. utr3 3' UTR sequence (DNA). Used for AU-rich element detection.

Returns

ManufacturabilityReport with per-check results and an overall 0-1 score (mean of all check scores).

Source code in mrnavax/manufacturability.py

def score_manufacturability(
    cds: str,
    *,
    utr5: str = "",
    utr3: str = "",
) -> ManufacturabilityReport:
    """Run all manufacturability checks and aggregate.

    Parameters
    ----------
    cds
        Coding sequence (DNA). If it ends with a stop codon, the stop
        codon is included in the analysis.
    utr5
        5' UTR sequence (DNA). Used for Kozak scoring — only the last
        9 nt upstream of ATG are inspected.
    utr3
        3' UTR sequence (DNA). Used for AU-rich element detection.

    Returns
    -------
    ManufacturabilityReport with per-check results and an overall
    0-1 score (mean of all check scores).
    """
    cds = cds.upper().replace("U", "T")
    utr5 = (utr5 or "").upper().replace("U", "T")
    utr3 = (utr3 or "").upper().replace("U", "T")

    checks: list[CheckResult] = [
        check_poly_a_runs(cds),
        check_gc_5prime_hairpin(cds),
        check_kozak_strength(utr5) if utr5 else _skip("kozak_strength"),
        check_are_motif(utr3) if utr3 else _skip("are_motif"),
        check_stop_context(cds),
        check_hidden_stops(cds),
        check_gc_window_uniformity(cds),
        check_cpg_balance(cds),
    ]
    overall = sum(c.score for c in checks) / len(checks) if checks else 1.0
    n_pass = sum(1 for c in checks if c.pass_ and c.severity != "error")
    n_warn = sum(1 for c in checks if c.severity == "warn")
    n_error = sum(1 for c in checks if c.severity == "error")
    return ManufacturabilityReport(
        overall_score=overall,
        n_pass=n_pass,
        n_warn=n_warn,
        n_error=n_error,
        checks=checks,
    )

LNP delivery

mrnavax.lnp_advisor

LNP composition advisor.

Recommends an LNP formulation based on: - target tissue (lung, liver, spleen, muscle, tumor) - cargo type (mRNA, saRNA, circRNA, siRNA, Cas9 mRNA) - therapeutic intent (cancer vaccine, gene editing, protein replacement)

The recommendation is a rule-based shortlist of published or ML-discovered ionizable lipids + helper lipid ratios. The ML layer in production would be something like the directed message-passing neural network in Witten et al. Nat Biotech 43, 1790–1799 (2025) — the LION nanoparticle library, or Li et al. Nat Mater 23, 1002 (2024) combinatorial-chemistry accelerated ionizable-lipid discovery. This module is the interface layer above that model: humans still pick from a shortlist.

References

Witten et al. 2025 (AI LNP for lung): https://www.nature.com/articles/s41587-024-02490-y Li et al. 2024 (combinatorial+ML): https://www.nature.com/articles/s41563-024-01833-5 Hou et al. 2021 (LNP review): Nat Rev Mater 6, 1078

recommend(*, target='liver', cargo='mRNA', intent='cancer vaccine', n=3)

Return a ranked shortlist of LNP formulations for a given scenario.

Source code in mrnavax/lnp_advisor.py

def recommend(
    *,
    target: str = "liver",
    cargo: str = "mRNA",
    intent: str = "cancer vaccine",
    n: int = 3,
) -> LNPAdvice:
    """Return a ranked shortlist of LNP formulations for a given scenario."""
    target = target.lower()
    cargo = cargo.lower()
    intent = intent.lower()
    notes: list[str] = []

    # ---- shortlist per target ----
    if target in {"lung", "pulmonary", "airway"}:
        pool = ["FO-32", "FO-35"]
        notes.append(
            "Pulmonary delivery benefits from ML-discovered biodegradable lipids (Witten 2025)."
        )
    elif target in {"liver", "hepatic"}:
        pool = ["C12-200", "SM-102", "ALC-0315"]
        notes.append("Liver is the natural sink for systemically administered LNPs.")
    elif target in {"spleen", "immune", "dendritic"}:
        pool = ["SM-102", "ALC-0315"]
        notes.append("Spleen/APCs are reached by ionizable LNPs at standard composition.")
    elif target in {"muscle", "im"}:
        pool = ["SM-102", "ALC-0315"]
        notes.append("Intramuscular LNPs typically use vaccine-grade composition.")
    elif target in {"tumor", "tme", "local"}:
        pool = ["tumor-it", "C12-200"]
        notes.append("Consider TME normalization (losartan/Ang(1-7)) alongside intratumoral LNP.")
    else:
        pool = ["SM-102", "ALC-0315", "C12-200"]
        notes.append(f"Unknown target '{target}'; defaulting to clinical-grade compositions.")

    # ---- cargo-specific adjustments ----
    if cargo in {"sarna", "srna", "self-amplifying"}:
        if "saRNA-gen-1" not in pool:
            pool.insert(0, "saRNA-gen-1")
        notes.append("saRNA prefers higher cholesterol and lower PEG for replicon stability.")
    elif cargo in {"circrna", "circular"}:
        notes.append("circRNA tolerates standard LNP composition; lower N/P often sufficient.")
    elif cargo in {"mrna", ""}:
        pass
    elif cargo in {"sirna", "sgrna"}:
        notes.append("siRNA/sgRNA cargoes typically use lower ionizable% (≤50%).")
    elif cargo in {"cas9", "crispr"}:
        if "C12-200" not in pool:
            pool.insert(0, "C12-200")
        notes.append("Cas9 mRNA + sgRNA co-delivery: C12-200 or similar hepatotropic LNPs.")
    else:
        notes.append(f"Unrecognized cargo '{cargo}'; using standard ratios.")

    # ---- intent-specific ----
    if intent in {"cancer vaccine", "vaccine"}:
        notes.append(
            "For cancer vaccines, MHC-I presentation of encoded antigen is the priority — clinical SM-102/ALC-0315 work well."
        )
    elif intent in {"gene editing", "crispr", "knock-in", "knock-out"}:
        if "C12-200" not in pool:
            pool.insert(0, "C12-200")
        notes.append(
            "Gene editing needs high transfection + minimal double-stranded RNA contaminants."
        )
    elif intent in {"protein replacement", "enzyme", "missing protein"}:
        notes.append("Protein-replacement mRNAs need high translation; codon-optimize + high dose.")

    shortlist = [PRESETS[k].as_dict() for k in pool if k in PRESETS][:n]
    return LNPAdvice(cargo=cargo, target=target, intent=intent, shortlist=shortlist, notes=notes)

CLI & shared utilities

mrnavax.cli

Unified CLI: python -m mrnavax.cli <tool> ...

mrnavax.llm

LLM-calling wrapper.

Two backends, selected at runtime:

1. hermes — when running inside the Hermes desktop app the agent itself can simply invoke this module's llm_complete from its own context. The hermes backend just records the call (no nested LLM in the demo) and returns a deterministic mock so the tool still runs end-to-end.

2. openai — when the OPENAI_API_KEY env var is set, calls go to the OpenAI Chat Completions API at gpt-4o-mini (cheap + fast, fine for tool use). Any OpenAI-compatible endpoint can be selected via OPENAI_BASE_URL.

3. mock — deterministic stub for tests and offline runs. Returns the prompt's last user message verbatim, prefixed with "[mock] ".

The CLI lets the caller pick a backend with --backend {mock,openai,hermes}.

llm_complete(prompt, *, system='', json_mode=False, temperature=0.2, max_tokens=800, backend=None)

Return a single LLM completion for prompt.

Backend selection: explicit backend arg > MRNA_AI_LLM_BACKEND env

auto-detect.

Source code in mrnavax/llm.py

def llm_complete(
    prompt: str,
    *,
    system: str = "",
    json_mode: bool = False,
    temperature: float = 0.2,
    max_tokens: int = 800,
    backend: str | None = None,
) -> str:
    """Return a single LLM completion for ``prompt``.

    Backend selection: explicit ``backend`` arg > ``MRNA_AI_LLM_BACKEND`` env
    > auto-detect.
    """
    backend = backend or _detect_backend()

    if backend == "mock":
        # Deterministic offline stub: extract a structured guess from the prompt.
        return _mock_complete(prompt, system=system, json_mode=json_mode)

    if backend == "openai":
        return _openai_complete(
            prompt,
            system=system,
            json_mode=json_mode,
            temperature=temperature,
            max_tokens=max_tokens,
        )

    if backend == "hermes":
        # Inside the Hermes app, the assistant already IS the LLM. We surface
        # a clear error so the operator knows to switch backends.
        raise RuntimeError(
            "backend='hermes' is reserved for the assistant's own context. "
            "Use backend='openai' with OPENAI_API_KEY, or backend='mock' for "
            "offline runs."
        )

    raise ValueError(f"unknown backend: {backend!r}")

llm_json(prompt, **kwargs)

Convenience: call llm_complete with json_mode=True and parse.

Source code in mrnavax/llm.py

def llm_json(prompt: str, **kwargs: Any) -> dict[str, Any]:
    """Convenience: call ``llm_complete`` with json_mode=True and parse."""
    raw = llm_complete(prompt, json_mode=True, **kwargs)
    # Strip code fences if any
    raw = raw.strip()
    if raw.startswith("```"):
        raw = raw.split("\n", 1)[1].rsplit("```", 1)[0].strip()
    try:
        return json.loads(raw)
    except json.JSONDecodeError as e:
        raise RuntimeError(f"LLM did not return valid JSON: {raw[:300]!r}") from e