"""Scan a construct end to end and annotate what is not already annotated. Two real problems this solves, neither covered before: 1. **Partially annotated records.** ``plasmid.auto_annotate`` bails the moment a record carries a single feature (``if plasmid.get("features"): return``), so a GenBank file with three curated features got *zero* help with the other twenty kilobases. Here the existing annotations are never touched — they are the ground truth the scan works around, and only the gaps get filled. 2. **Unannotated records.** A raw FASTA or a pasted sequence arrives with no features at all. Everything below applies unchanged; there simply are no gaps to work around. Evidence discipline: every finding says where it came from, and the tiers are kept apart on purpose. ``curated`` is an exact match to a sequence this repo verified against GenBank accessions; ``motif`` is an exact match to a short known element; ``predicted`` is computed from the sequence alone and matches no reference. A predicted ORF is never dressed up as an identified gene — the scan reports "unidentified ORF, 412 aa", because that is all the sequence supports. """ from __future__ import annotations from typing import Any, Dict, List, Optional, Tuple from .plasmid import (COMMON_ENZYMES, _color_for, find_motifs, find_orfs, revcomp) # A finding must sit mostly in unannotated territory to count as a gap fill. # Below this it is treated as commentary on an annotation that already exists. _GAP_FRACTION = 0.5 # Dark regions shorter than this are noise between adjacent features. _MIN_DARK_BP = 150 # An MCS is a run of distinct cutters packed into a short stretch. _MCS_WINDOW = 120 _MCS_MIN_SITES = 4 # An ORF this expected by chance says nothing at all; do not report it. _NOISE_EXPECTATION = 100.0 CURATED, MOTIF, PREDICTED = "curated", "motif", "predicted" _RANK = {CURATED: 0, MOTIF: 1, PREDICTED: 2} # ── interval helpers ───────────────────────────────────────────────────── def _merge(spans: List[Tuple[int, int]]) -> List[Tuple[int, int]]: out: List[Tuple[int, int]] = [] for s, e in sorted(spans): if out and s <= out[-1][1]: out[-1] = (out[-1][0], max(out[-1][1], e)) else: out.append((s, e)) return out def _covered(spans: List[Tuple[int, int]], s: int, e: int) -> int: """Bases of ``[s,e)`` that fall inside ``spans`` (which must be merged).""" n = 0 for a, b in spans: if b <= s: continue if a >= e: break n += min(b, e) - max(a, s) return n def _norm(features: Optional[List[Dict[str, Any]]]) -> List[Dict[str, Any]]: out = [] for f in features or []: try: s, e = int(f.get("start", 0)), int(f.get("end", 0)) except (TypeError, ValueError): continue if e > s: out.append(dict(f, start=s, end=e)) return out # ── detection layers ───────────────────────────────────────────────────── def _curated_promoters() -> List[Tuple[str, str, str]]: """(name, sequence, provenance) across every host in the library. Annotation is host-agnostic: a T7 promoter reads as a T7 promoter whatever organism the user is targeting, so the whole library is searched rather than a host-gated slice of it. """ from .promoters import PROMOTER_LIBRARY seen, out = set(), [] for host_promoters in PROMOTER_LIBRARY.values(): for p in host_promoters: seq = (p.sequence or "").upper() if len(seq) < 12 or seq in seen: continue seen.add(seq) out.append((p.name, seq, p.source or "")) return out def _find_curated(seq: str) -> List[Dict[str, Any]]: up, out = seq.upper(), [] for name, motif, source in _curated_promoters(): for strand, needle in ((1, motif), (-1, revcomp(motif))): i = up.find(needle) while i != -1: out.append({ "name": (name if name.lower().endswith("promoter") else f"{name} promoter"), "type": "promoter", "start": i, "end": i + len(needle), "strand": strand, "color": _color_for("promoter", name), "tier": CURATED, "evidence": f"exact match, {len(needle)} bp. {source}", }) i = up.find(needle, i + 1) return out def _chance_model(seq: str) -> Tuple[float, float]: """(P(ATG), P(stop)) under a zero-order model of this construct's own bases. Deliberately the simplest defensible model: base frequencies measured from the sequence in hand, codons treated as independent draws. It ignores codon bias and dinucleotide structure, so it is an approximation and is labelled as one wherever it surfaces \u2014 but it is measured, not assumed, and it is enough to separate "this ORF is the kind of thing random DNA produces all day" from "a run this long should not be here by chance". """ n = len(seq) or 1 f = {b: seq.count(b) / n for b in "ACGT"} tot = sum(f.values()) or 1.0 a, c, g, t = (f[b] / tot for b in "ACGT") return a * t * g, t * a * a + t * a * g + t * g * a def orf_chance_expectation(length_bp: int, aa: int, p_atg: float, p_stop: float) -> float: """Roughly how many ORFs this long a construct this size would hold if the sequence were random with the same base composition. ``2L`` codon slots across six frames, each needing an ATG followed by ``aa`` consecutive non-stop codons. Calibrated against simulated random sequence at 35-65% GC (36 conditions): observed counts came in at a median 0.83x this figure, spread 0.0-1.8x. So it is an order-of-magnitude estimate, **not** a bound \u2014 individual counts do run above it, and the wording everywhere it surfaces says "on the order of" for that reason. It ignores codon bias and dinucleotide structure, but its inputs are measured from the sequence in hand rather than assumed. Well above 1 means an ORF of this length is unremarkable; well below 1 means chance is a poor explanation for it. """ if p_stop <= 0 or aa <= 0: return float("inf") return 2.0 * max(1, length_bp) * p_atg * ((1.0 - p_stop) ** aa) def _orfs_circular(seq: str, *, min_aa: int, circular: bool) -> List[Dict[str, Any]]: """ORFs, including the ones that cross the origin on a circular plasmid. ``find_orfs`` walks a linear string, so on a circular construct any gene spanning position 0 is invisible to it — a real blind spot, since cloning routinely lands a cassette across the join. Scanning a doubled sequence surfaces those; findings are then folded back into plasmid coordinates and de-duplicated. """ L = len(seq) if not circular or L < 6: return find_orfs(seq, min_aa=min_aa) doubled = find_orfs(seq + seq, min_aa=min_aa) out, seen = [], set() for f in doubled: s, e = f["start"], f["end"] if s >= L or (e - s) > L: # duplicate copy, or longer than the plasmid continue key = (s % L, (e - s), f["strand"]) if key in seen: continue seen.add(key) g = dict(f, start=s, end=e) if e > L: g["wraps"] = True g["notes"] = "auto-detected ORF (crosses the origin)" out.append(g) return out def _find_mcs(seq: str, circular: bool) -> List[Dict[str, Any]]: """A multiple cloning site read structurally: unique cutters packed tight. What makes a polylinker a polylinker is that its enzymes cut *once in the whole plasmid* and sit within a few dozen bases of each other — an enzyme that also cuts three more times elsewhere is useless for directional cloning and is not evidence of an MCS. Nothing here is matched against a reference; the claim is about the arrangement, which is what an MCS is. """ from .plasmid import digest try: d = digest(seq, "circular" if circular else "linear", list(COMMON_ENZYMES)) except Exception: return [] singles = [(e["sites"][0], e["name"]) for e in (d.get("per_enzyme") or []) if len(e.get("sites") or []) == 1] singles.sort() if len(singles) < _MCS_MIN_SITES: return [] best: Optional[Tuple[int, int, List[str]]] = None for i, (pos, _) in enumerate(singles): window = [(p, n) for p, n in singles[i:] if p - pos <= _MCS_WINDOW] if len(window) >= _MCS_MIN_SITES and (best is None or len(window) > len(best[2])): best = (pos, window[-1][0], [n for _, n in window]) if not best: return [] s0, e0, names = best start, end = max(0, s0 - 8), min(len(seq), e0 + 8) return [{ "name": "MCS", "type": "misc_feature", "start": start, "end": end, "strand": 1, "color": _color_for("misc_feature", "MCS"), "tier": PREDICTED, "evidence": (f"{len(names)} enzymes that each cut this construct exactly " f"once, all within {e0 - s0} bp ({', '.join(sorted(names))}) " "\u2014 a polylinker arrangement, not a reference match."), }] def _dedupe_orfs(cands: List[Dict[str, Any]]) -> List[Dict[str, Any]]: """Drop an ORF wholly contained in a longer one (same trick auto_annotate uses).""" cds = sorted([c for c in cands if c["type"] == "CDS"], key=lambda f: f["end"] - f["start"], reverse=True) kept: List[Dict[str, Any]] = [] for f in cds: if not any(k["start"] <= f["start"] and f["end"] <= k["end"] for k in kept): kept.append(f) return [c for c in cands if c["type"] != "CDS"] + kept # ── the scan ───────────────────────────────────────────────────────────── def scan_construct(sequence: str, *, features: Optional[List[Dict[str, Any]]] = None, topology: str = "circular", min_aa: int = 60, include_orfs: bool = True, max_new: int = 200) -> Dict[str, Any]: """Annotate every stretch of ``sequence`` that ``features`` leaves unexplained. Existing features are never modified, reordered or dropped. Returns the findings split by what they mean relative to what was already there: ``new`` — sits in unannotated sequence. These are the gap fills. ``confirms`` — lands on sequence that is already annotated. Named by what it overlaps; overlap is not agreement. ``coverage`` — how much of the construct is accounted for, before and after, plus the stretches still unexplained by anything. """ seq = "".join((sequence or "").split()).upper() L = len(seq) if L == 0: raise ValueError("No sequence to scan.") circular = str(topology or "").lower().startswith("c") existing = _norm(features) prior = _merge([(f["start"], min(f["end"], L)) for f in existing if f["start"] < L]) cands = _find_curated(seq) for m in find_motifs(seq): cands.append(dict(m, tier=MOTIF, evidence=f"exact match to a known {m['type']} element, " f"{m['end'] - m['start']} bp.")) p_atg, p_stop = _chance_model(seq) if include_orfs: for o in _orfs_circular(seq, min_aa=min_aa, circular=circular): aa = (o["end"] - o["start"]) // 3 exp = orf_chance_expectation(L, aa, p_atg, p_stop) if exp > _NOISE_EXPECTATION: continue if exp >= 1.0: verdict = (f"a construct this size would be expected to hold " f"on the order of {exp:.1f} ORFs this long by " "chance alone, so " "the length on its own is not evidence of a gene") elif exp >= 0.05: verdict = (f"expected on the order of {exp:.2f} times by chance in a " "construct this size \u2014 suggestive, not conclusive") else: verdict = (f"expected on the order of {exp:.1e} times by chance in a " "construct this size \u2014 unlikely to be a random " "open frame") cands.append(dict( o, name=f"unidentified ORF ({aa} aa)", tier=PREDICTED, chance=exp, evidence=("open reading frame, ATG to stop, " f"{aa} aa on the {'+' if o['strand'] > 0 else '−'} strand. " "No reference match, so the identity is unknown; " f"{verdict}. (Order-of-magnitude estimate, " "zero-order composition model.)"))) cands += _find_mcs(seq, circular) cands = _dedupe_orfs(cands) # Exact matches nest legitimately: a His tag sits inside its CDS, an # operator inside its promoter, an RBS just before a start codon. So a # curated or motif hit is never suppressed by an overlap and never demoted # to "already known" — it is information either way. Only *predicted* # findings compete, because two overlapping ORF calls really are redundant. # The same element can appear in both the curated library and the motif # list (T7, for one). Same span, same strand — report it once, keeping the # stronger provenance. by_span: Dict[Tuple[int, int, int], Dict[str, Any]] = {} for c in cands: if c.get("tier") == PREDICTED: continue k = (c["start"], c["end"], int(c.get("strand") or 1)) cur = by_span.get(k) if cur is None or _RANK.get(c.get("tier"), 9) < _RANK.get(cur.get("tier"), 9): by_span[k] = c exact = list(by_span.values()) guesses = [c for c in cands if c.get("tier") == PREDICTED] new: List[Dict[str, Any]] = list(exact) overlapping: List[Dict[str, Any]] = [] taken = list(prior) + [(c["start"], min(c["end"], L)) for c in exact] for c in sorted(guesses, key=lambda f: (f.get("chance", 0.0), -(f["end"] - f["start"]))): span = max(1, min(c["end"], L) - c["start"]) if _covered(prior, c["start"], min(c["end"], L)) / span >= _GAP_FRACTION: hit = next((f for f in existing if f["start"] < c["end"] and c["start"] < f["end"]), None) c["overlaps"] = (hit or {}).get("name") or "" overlapping.append(c) continue seen_here = _covered(_merge(taken), c["start"], min(c["end"], L)) if seen_here / span >= _GAP_FRACTION: continue # another prediction already explains it new.append(c) taken.append((c["start"], min(c["end"], L))) if len(new) >= max_new: break new.sort(key=lambda f: (f["start"], _RANK.get(f.get("tier"), 9), f["end"])) kept = new confirms = overlapping after = _merge(taken) dark = [] cur = 0 for a, b in after: if a - cur >= _MIN_DARK_BP: dark.append({"start": cur, "end": a, "length": a - cur}) cur = max(cur, b) if L - cur >= _MIN_DARK_BP: dark.append({"start": cur, "end": L, "length": L - cur}) before_bp = sum(b - a for a, b in prior) after_bp = sum(b - a for a, b in after) return { "ok": True, "length": L, "topology": "circular" if circular else "linear", "existing": existing, "new": kept, "confirms": confirms, "coverage": { "before_bp": before_bp, "after_bp": after_bp, "before_pct": round(100.0 * before_bp / L, 1), "after_pct": round(100.0 * after_bp / L, 1), "dark_regions": dark, "dark_bp": sum(d["length"] for d in dark), }, "counts": { "existing": len(existing), "new": len(kept), "confirms": len(confirms), "curated": sum(1 for c in kept if c.get("tier") == CURATED), "motif": sum(1 for c in kept if c.get("tier") == MOTIF), "predicted": sum(1 for c in kept if c.get("tier") == PREDICTED), }, } def merge_annotations(features: Optional[List[Dict[str, Any]]], found: List[Dict[str, Any]]) -> List[Dict[str, Any]]: """Existing features plus the accepted findings, in coordinate order. Curated annotations keep their exact record; findings are marked so a later reader can always tell a scan result from a human one. """ out = list(_norm(features)) for f in found: g = {k: v for k, v in f.items() if k not in ("overlaps",)} g["notes"] = (g.get("evidence") or g.get("notes") or "") g["auto"] = True out.append(g) out.sort(key=lambda f: (f["start"], f["end"])) return out