Spaces:
Running
Running
| """Sequence utilities toolkit. | |
| Computes the everyday metrics a biologist needs on a single sequence: | |
| * GC content (nucleotides) | |
| * Reverse complement (nucleotides, IUPAC-aware) | |
| * Molecular weight (ssDNA / ssRNA / protein average-residue) | |
| * Six-frame-free translation of the forward three frames with the longest | |
| open reading frame (ORF) flagged | |
| * Amino-acid composition (proteins / translated CDS) | |
| * Restriction-enzyme site scan against a curated set of common, unambiguous | |
| (palindromic) recognition sequences | |
| Pure local computation — no network calls. Mirrors the "Sequence Stats" | |
| workflow every bench tool (ExPASy, BioPython scripts) implements. | |
| """ | |
| from __future__ import annotations | |
| import re | |
| from Bio.SeqUtils.ProtParam import ProteinAnalysis | |
| from app.services.sequence_utils import detect_sequence_type | |
| class SequenceUtilitiesError(ValueError): | |
| pass | |
| # Standard genetic code — stop codons map to '*', ambiguous codons to 'X'. | |
| _CODON_TABLE = { | |
| "TTT": "F", "TTC": "F", "TTA": "L", "TTG": "L", | |
| "TCT": "S", "TCC": "S", "TCA": "S", "TCG": "S", | |
| "TAT": "Y", "TAC": "Y", "TAA": "*", "TAG": "*", | |
| "TGT": "C", "TGC": "C", "TGA": "*", "TGG": "W", | |
| "CTT": "L", "CTC": "L", "CTA": "L", "CTG": "L", | |
| "CCT": "P", "CCC": "P", "CCA": "P", "CCG": "P", | |
| "CAT": "H", "CAC": "H", "CAA": "Q", "CAG": "Q", | |
| "CGT": "R", "CGC": "R", "CGA": "R", "CGG": "R", | |
| "ATT": "I", "ATC": "I", "ATA": "I", "ATG": "M", | |
| "ACT": "T", "ACC": "T", "ACA": "T", "ACG": "T", | |
| "AAT": "N", "AAC": "N", "AAA": "K", "AAG": "K", | |
| "AGT": "S", "AGC": "S", "AGA": "R", "AGG": "R", | |
| "GTT": "V", "GTC": "V", "GTA": "V", "GTG": "V", | |
| "GCT": "A", "GCC": "A", "GCA": "A", "GCG": "A", | |
| "GAT": "D", "GAC": "D", "GAA": "E", "GAG": "E", | |
| "GGT": "G", "GGC": "G", "GGA": "G", "GGG": "G", | |
| } | |
| # Monoisotopic-ish single-strand base weights (g/mol). | |
| _SSDNA_MW = {"A": 313.21, "C": 289.18, "G": 329.21, "T": 304.20} | |
| _SSRNA_MW = {"A": 329.21, "C": 305.18, "G": 345.21, "U": 306.17} | |
| # Curated common restriction enzymes — all palindromic so a forward-strand | |
| # scan finds every cut site. Recognition site is given 5' -> 3'. | |
| _RESTRICTION_ENZYMES = [ | |
| ("EcoRI", "GAATTC"), | |
| ("BamHI", "GGATCC"), | |
| ("HindIII", "AAGCTT"), | |
| ("SalI", "GTCGAC"), | |
| ("XbaI", "TCTAGA"), | |
| ("XhoI", "CTCGAG"), | |
| ("NotI", "GCGGCCGC"), | |
| ("KpnI", "GGTACC"), | |
| ("SmaI", "CCCGGG"), | |
| ("PstI", "CTGCAG"), | |
| ("SacI", "GAGCTC"), | |
| ] | |
| _RNA_TO_DNA = str.maketrans("Uu", "Tt") | |
| _AA_ALPHABET = "ACDEFGHIKLMNPQRSTVWY" | |
| def _strip_fasta(seq: str) -> str: | |
| """Return just the sequence body of a raw or FASTA-formatted input.""" | |
| lines = (seq or "").strip().splitlines() | |
| lines = [ln.strip() for ln in lines if not ln.strip().startswith(">")] | |
| return "".join(lines) | |
| def clean_sequence(seq: str, seq_type: str) -> str: | |
| """Normalize to uppercase alpha-only. For nucleotides, drop ambiguous IUPAC | |
| codes so downstream math (GC%, MW) only counts real bases.""" | |
| body = _strip_fasta(seq) | |
| letters = "".join(re.findall(r"[A-Za-z]", body)).upper() | |
| if not letters: | |
| raise SequenceUtilitiesError("Sequence is empty") | |
| if seq_type in ("dna", "rna"): | |
| allowed = set("ACGTRYSWKMBDHVN") if seq_type == "dna" else set("ACGURSYKMWBDHVN") | |
| kept = "".join(c for c in letters if c in allowed) | |
| if not kept: | |
| raise SequenceUtilitiesError(f"Sequence contains no valid {seq_type.upper()} bases") | |
| return kept.translate(_RNA_TO_DNA) if seq_type == "rna" else kept | |
| valid = set("ACDEFGHIKLMNPQRSTVWY") | |
| kept = "".join(c for c in letters if c in valid) | |
| if not kept: | |
| raise SequenceUtilitiesError("Sequence contains no valid amino acids") | |
| return kept | |
| def _translate_frame(seq: str, frame: int) -> str: | |
| codons = [seq[i:i + 3] for i in range(frame, len(seq) - 2, 3)] | |
| return "".join(_CODON_TABLE.get(c, "X") for c in codons) | |
| def _best_orf(seq: str, frame: int, translated: str) -> dict | None: | |
| """Longest ORF (M -> stop/end) in a translated frame, with 1-based start.""" | |
| best = None | |
| for m in re.finditer("M[^*]*", translated): | |
| length = len(m.group(0)) | |
| start = frame + m.start() * 3 + 1 | |
| if best is None or length > best["length"]: | |
| best = { | |
| "frame": frame + 1, | |
| "protein": m.group(0), | |
| "start": start, | |
| "length": length, | |
| "has_stop": len(translated) > m.end() and translated[m.end()] == "*", | |
| "starts_with_m": True, | |
| } | |
| return best | |
| def _protein_mw(seq: str) -> float: | |
| try: | |
| return round(ProteinAnalysis(seq).molecular_weight(), 2) | |
| except Exception: | |
| avg = 110.0 | |
| return round(sum(avg for _ in seq), 2) | |
| def _nucleotide_mw(seq: str, seq_type: str) -> float: | |
| table = _SSDNA_MW if seq_type == "dna" else _SSRNA_MW | |
| n = len(seq) | |
| if n == 0: | |
| return 0.0 | |
| total = sum(table.get(c, table["A"]) for c in seq) | |
| return round(total - 61.96 * (n - 1) + 18.02, 2) | |
| def _aa_composition(seq: str) -> list[dict]: | |
| counts: dict[str, int] = {} | |
| for c in seq.upper(): | |
| if c in _AA_ALPHABET: | |
| counts[c] = counts.get(c, 0) + 1 | |
| total = sum(counts.values()) | |
| comp = [ | |
| {"aa": aa, "count": count, "pct": round(count / total * 100, 1) if total else 0.0} | |
| for aa, count in sorted(counts.items(), key=lambda kv: (-kv[1], kv[0])) | |
| ] | |
| return comp | |
| def _restriction_scan(seq: str) -> list[dict]: | |
| sites = [] | |
| for name, recognition in _RESTRICTION_ENZYMES: | |
| positions = [m.start() + 1 for m in re.finditer(recognition, seq)] | |
| if positions: | |
| sites.append({"name": name, "recognition": recognition, "count": len(positions), "positions": positions}) | |
| return sites | |
| def analyze_sequence(sequence: str, seq_type: str = "auto") -> dict: | |
| """Analyze a single sequence and return a flat report dict. | |
| ``seq_type`` may be 'auto', 'dna', 'rna' or 'protein'. Auto-detection uses | |
| the shared alphabet classifier (a pure ACGTN/U string is treated as DNA). | |
| """ | |
| seq_type = (seq_type or "auto").lower() | |
| if seq_type not in ("auto", "dna", "rna", "protein"): | |
| raise SequenceUtilitiesError("seq_type must be auto, dna, rna or protein") | |
| raw = _strip_fasta(sequence) | |
| detected = detect_sequence_type(raw) if raw else "unknown" | |
| effective = seq_type if seq_type != "auto" else detected | |
| if effective == "unknown": | |
| raise SequenceUtilitiesError( | |
| "Could not detect sequence type — expected nucleotide or protein characters" | |
| ) | |
| seq = clean_sequence(raw, effective) | |
| issues: list[str] = [] | |
| report: dict = { | |
| "sequence_type": effective, | |
| "detected_type": detected, | |
| "length": len(seq), | |
| "gc_content": None, | |
| "molecular_weight": None, | |
| "reverse_complement": None, | |
| "translation": None, | |
| "aa_composition": None, | |
| "restriction_sites": None, | |
| "issues": issues, | |
| } | |
| if effective in ("dna", "rna"): | |
| report["gc_content"] = round((seq.count("G") + seq.count("C")) / len(seq) * 100.0, 1) | |
| report["molecular_weight"] = _nucleotide_mw(seq, effective) | |
| if effective == "dna": | |
| comp = {"A": "T", "T": "A", "G": "C", "C": "G", "N": "N"} | |
| else: | |
| comp = {"A": "U", "U": "A", "G": "C", "C": "G", "N": "N"} | |
| report["reverse_complement"] = "".join(comp.get(c, "N") for c in reversed(seq)) | |
| dna_seq = seq.translate(_RNA_TO_DNA) | |
| if len(dna_seq) < 3: | |
| issues.append("Sequence too short for translation (<3 nt)") | |
| else: | |
| frames = {} | |
| best = None | |
| for frame in (0, 1, 2): | |
| translated = _translate_frame(dna_seq, frame) | |
| frames[str(frame + 1)] = translated | |
| orf = _best_orf(dna_seq, frame, translated) | |
| if orf and (best is None or orf["length"] > best["length"]): | |
| best = orf | |
| report["translation"] = {"frames": frames, "best": best} | |
| if best is None: | |
| issues.append("No in-frame methionine (ATG) found — no ORF to report") | |
| else: | |
| report["aa_composition"] = _aa_composition(best["protein"]) | |
| if effective == "dna": | |
| report["restriction_sites"] = _restriction_scan(seq) | |
| else: | |
| issues.append("Restriction-site scan is DNA-only") | |
| else: | |
| report["molecular_weight"] = _protein_mw(seq) | |
| report["aa_composition"] = _aa_composition(seq) | |
| if len(seq) < 2: | |
| issues.append("Protein sequence very short — composition may be uninformative") | |
| return report | |