bioai-platform/backend/app/tools/oligo_qc.py

"""Oligo QC and in-silico PCR helpers for primer design verification.

Implements the checks a wet-lab scientist runs on candidate primers before
ordering them (the IDT OligoAnalyzer + Primer-BLAST + in-silico PCR steps of
the canonical primer design workflow):

  * GC content and salt-adjusted melting temperature
  * Hairpin (intra-molecular stem-loop) prediction
  * Self-dimer / hetero-dimer prediction (complementary annealing)
  * In-silico PCR: where each primer binds in the template, expected amplicon
    length, and specificity (number of binding sites).

These are heuristic approximations suitable for screening. A delta-G is
estimated with a simple nearest-neighbour-style pair score so primers can be
ranked; it is not a substitute for a full NN thermodynamics calculation.
"""

from __future__ import annotations

import math
import re
from typing import Dict, List, Optional, Tuple

_COMP = {ord('A'): 'T', ord('T'): 'A', ord('G'): 'C', ord('C'): 'G',
         ord('N'): 'N', ord('U'): 'A',
         ord('a'): 't', ord('t'): 'a', ord('g'): 'c', ord('c'): 'g',
         ord('n'): 'n', ord('u'): 'a'}

# Per-base-pair stacking energies (kcal/mol) for the complementary stem.
# Only canonical Watson-Crick pairs count towards a stem; wobble pairs are
# ignored so short weak "stems" don't turn every oligo into a dimer/hairpin.
_PAIR_DG = {
    ('A', 'T'): -1.0, ('T', 'A'): -1.0,
    ('G', 'C'): -2.0, ('C', 'G'): -2.0,
}


def clean(seq: str) -> str:
    """Upper-case and strip whitespace (and any non-letter junk)."""
    return re.sub(r"[^A-Za-z]", "", seq).upper().replace("U", "T")


def is_dna(seq: str) -> bool:
    s = clean(seq)
    return bool(s) and all(c in "ATGCN" for c in s)


def reverse_complement(seq: str) -> str:
    return seq.upper().translate(_COMP)[::-1]


def gc_content(seq: str) -> float:
    s = clean(seq)
    if not s:
        return 0.0
    return round((s.count("G") + s.count("C")) / len(s) * 100.0, 1)


def _pair_energy(a: str, b: str) -> float:
    return _PAIR_DG.get((a, b), _PAIR_DG.get((b, a), 0.0))


def salt_adjusted_tm(seq: str, salt_mm: float = 50.0) -> float:
    """IDT-style melting temperature.

    Uses the classic NN-approximation formula Tm = 64.9 + 41*(nGC - 16.4)/N
    (valid at ~50 mM monovalent salt) plus the SantaLucia salt correction.
    """
    s = clean(seq)
    n = len(s)
    if n < 1:
        return 0.0
    n_gc = s.count("G") + s.count("C")
    tm = 64.9 + 41.0 * (n_gc - 16.4) / n
    if salt_mm != 50.0:
        tm += 16.6 * math.log10(max(salt_mm, 0.01) / 50.0)
    return round(tm, 1)


def _risk_from_dg(dg: float) -> Tuple[str, str]:
    if dg <= -7.0:
        return "high", "Strongly stabilizing — high risk of mis-priming."
    if dg <= -4.5:
        return "medium", "Moderately stabilizing — keep an eye on it."
    if dg <= -3.0:
        return "low", "Weakly stabilizing — unlikely to matter in practice."
    return "none", "No significant structure predicted."


def hairpin_analysis(seq: str, min_stem: int = 4, min_loop: int = 3) -> Dict:
    """Predict the strongest intra-molecular hairpin.

    A hairpin forms when the strand folds so that one arm is the reverse
    complement of the other arm, separated by a loop of at least `min_loop`
    bases: seq[i:i+L] must equal reverse_complement(seq[j:j+L]). Returns
    delta-G, stem length, stem sequence, loop sequence and a risk level.
    """
    s = clean(seq)
    n = len(s)
    best: Optional[Tuple[float, int, int, int]] = None  # (dg, stem_len, i, j)
    for i in range(n):
        for j in range(i + min_stem + min_loop, n):
            # loop = j - (i + L) must stay >= min_loop
            lmax = min(j - i - min_loop, n - j, 10)
            if lmax < min_stem:
                continue
            for L in range(min_stem, lmax + 1):
                dg = 0.0
                ok = True
                for k in range(L):
                    e = _pair_energy(s[i + k], s[j + L - 1 - k])
                    if e == 0.0:
                        ok = False
                        break
                    dg += e
                if ok and (best is None or dg < best[0]):
                    best = (dg, L, i, j)

    if best is None:
        return {"dg": 5.0, "stem_length": 0, "stem": "", "loop": "",
                "risk": "none", "note": "No hairpin structure predicted."}
    dg, stem_len, i, j = best
    stem = s[i:i + stem_len]
    loop = s[i + stem_len:j]
    risk, note = _risk_from_dg(dg)
    return {
        "dg": round(dg, 1),
        "stem_length": stem_len,
        "stem": stem,
        "loop": loop,
        "risk": risk,
        "note": note,
    }


def dimer_analysis(a: str, b: str, min_stem: int = 4) -> Dict:
    """Predict the strongest dimer between two primers.

    A dimer forms where primer `a` anneals antiparallel to primer `b`, i.e.
    where a region of `a` is complementary to a region of reversed `b`. Reports
    the most stabilizing complementary stem, whether the stem involves the 3'
    end of either primer (extension-competent dimers are the dangerous ones
    for PCR), and a risk level.
    """
    pa = clean(a)
    rev_b = clean(b)[::-1]
    n, m = len(pa), len(rev_b)

    best: Optional[Tuple[float, int, int, int]] = None  # (dg, stem_len, i, j)
    for i in range(n):
        for j in range(m):
            k = 0
            dg = 0.0
            while i + k < n and j + k < m:
                e = _pair_energy(pa[i + k], rev_b[j + k])
                if e == 0.0:
                    break
                dg += e
                k += 1
                if k >= 12:
                    break
            if k >= min_stem and (best is None or dg < best[0]):
                best = (dg, k, i, j)

    if best is None:
        return {"dg": 5.0, "stem_length": 0, "stem": "", "involves_a3": False,
                "involves_b3": False, "risk": "none",
                "note": "No significant dimer predicted."}
    dg, stem_len, i, j = best
    stem = pa[i:i + stem_len]
    a3 = (i + stem_len == n)  # stem reaches the 3' end of primer a
    b3 = (j == 0)             # rev_b[0] complements b's 3' terminal base
    risk, note = _risk_from_dg(dg)
    if a3 or b3:
        note = ("3' end involved — extension-competent, more likely to "
                "interfere with PCR. ") + note
    return {
        "dg": round(dg, 1),
        "stem_length": stem_len,
        "stem": stem,
        "involves_a3": a3,
        "involves_b3": b3,
        "risk": risk,
        "note": note,
    }


def find_binding_sites(template: str, primer: str) -> List[int]:
    """All 0-based start positions where `primer` matches `template`."""
    t = clean(template)
    p = clean(primer)
    if not t or not p:
        return []
    return [m.start() for m in re.finditer(r"(?=" + re.escape(p) + r")", t)]


def in_silico_pcr(template: str, left: str, right: str,
                  expected_product: Optional[int] = None,
                  left_expected: Optional[int] = None,
                  right_expected: Optional[int] = None) -> Dict:
    """Simulate the PCR reaction against a template sequence.

    The forward primer binds as-is; the reverse primer (given 5'->3' as
    Primer3 returns it) is the reverse complement of the forward-strand region
    it anneals to. Amplicons are [forward_site, reverse_site + len(right)).
    """
    t = clean(template)
    lseq = clean(left)
    rseq = clean(right)

    left_sites = find_binding_sites(t, lseq)
    right_target = reverse_complement(rseq)
    right_sites = find_binding_sites(t, right_target)

    amplicons: List[Tuple[int, int, int]] = []
    for f in left_sites:
        for r in right_sites:
            end = r + len(rseq)
            if end > f:
                amplicons.append((f, end, end - f))
    amplicons.sort(key=lambda x: x[2])

    specific = len(left_sites) == 1 and len(right_sites) == 1
    primer3_consistent = False
    if left_expected is not None and right_expected is not None and right_expected > 0:
        right_bind = right_expected - len(rseq) + 1
        primer3_consistent = (
            left_expected in left_sites and right_bind in right_sites
        )

    result: Dict = {
        "template_length": len(t),
        "forward_binding_sites": len(left_sites),
        "reverse_binding_sites": len(right_sites),
        "forward_positions": left_sites[:20],
        "reverse_positions": right_sites[:20],
        "specific": specific,
        "amplicons": [{"start": f, "end": e, "length": ln}
                      for f, e, ln in amplicons[:10]],
        "primer3_consistent": primer3_consistent,
        "note": None,
    }

    if expected_product is not None:
        # Primer3 product size == right_pos - left_pos + 1; allow 1 bp slack.
        match = any(abs(a[2] - expected_product) <= 1 for a in amplicons)
        result["matches_product_size"] = match
    else:
        result["matches_product_size"] = None

    if not left_sites and not right_sites:
        result["note"] = "Neither primer binds the template."
    elif not left_sites:
        result["note"] = "Forward primer does not bind the template."
    elif not right_sites:
        result["note"] = "Reverse primer does not bind the template."
    elif not specific:
        result["note"] = ("Multiple binding sites detected — this primer pair "
                          "may amplify off-target products.")
    else:
        result["note"] = "Specific: each primer binds exactly once."
    return result


def oligo_report(seq: str) -> Dict:
    """Full QC report for a single primer."""
    s = clean(seq)
    hp = hairpin_analysis(s)
    sd = dimer_analysis(s, s)
    return {
        "sequence": s,
        "length": len(s),
        "gc": gc_content(s),
        "tm_50mM": salt_adjusted_tm(s),
        "hairpin": hp,
        "self_dimer": sd,
    }