Spaces:
Running
Running
| """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 | |