Spaces:
Running
Running
Download bioai-platform/backend/app/routers/phylo.py from Samad14/bio-nexus-api: direct link, hf CLI and curl.
- Browser
- Download file 15.5 kB
-
https://huggingface.co/spaces/Samad14/bio-nexus-api/resolve/main/bioai-platform/backend/app/routers/phylo.py
- Command line
-
hf download hf://spaces/Samad14/bio-nexus-api/bioai-platform/backend/app/routers/phylo.py
-
curl -L -o phylo.py https://huggingface.co/spaces/Samad14/bio-nexus-api/resolve/main/bioai-platform/backend/app/routers/phylo.py
15.5 kB
| """ | |
| Phylogenetic tree router β three methods: | |
| nj : Neighbor-Joining via Clustal Omega guide tree (~60s) | |
| upgma : UPGMA computed locally from Clustal alignment (adds ~0s after MSA) | |
| ml : Maximum Likelihood via local PhyML binary (~3-5 min, includes bootstrap) | |
| PhyML binary must be installed at build time (see Dockerfile). | |
| Download source from: https://github.com/stephaneguindon/phyml | |
| Compile: ./configure --enable-phyml && make && make install | |
| Job lifecycle (in-memory, thread-safe): | |
| queued -> msa_running -> msa_done -> tree_running -> complete | error | |
| """ | |
| from __future__ import annotations | |
| import asyncio | |
| import logging | |
| import threading | |
| import time | |
| import uuid | |
| from typing import Literal, Optional | |
| import httpx | |
| from fastapi import APIRouter, BackgroundTasks, Depends, HTTPException | |
| from pydantic import BaseModel, Field | |
| logger = logging.getLogger(__name__) | |
| router = APIRouter(prefix="/phylo", tags=["phylo"]) | |
| # ββ EBI base URLs ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| _EBI_CLUSTALO = "https://www.ebi.ac.uk/Tools/services/rest/clustalo" | |
| _EMAIL = "bionexus@example.com" | |
| # ββ PhyML protein models (most-used first) ββββββββββββββββββββββββββββββββββββ | |
| PROTEIN_MODELS = ["LG", "WAG", "JTT", "Blosum62", "MtREV", "Dayhoff"] | |
| DNA_MODELS = ["GTR", "HKY85", "K80", "F81", "TN93", "SYM"] | |
| # βββ Models βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| Method = Literal["nj", "ml", "upgma"] | |
| SeqType = Literal["protein", "dna"] | |
| JobPhase = Literal["queued", "msa_running", "msa_done", "tree_running", "complete", "error"] | |
| class PhyloRequest(BaseModel): | |
| sequences: list[dict] | |
| method: Method = "nj" | |
| seq_type: SeqType = "protein" | |
| # ML-only options | |
| model: str = "LG" | |
| bootstrap: int = Field(100, ge=0, le=1000) | |
| class PhyloJob(BaseModel): | |
| job_id: str | |
| method: Method | |
| seq_type: SeqType | |
| model: Optional[str] | |
| bootstrap: Optional[int] | |
| phase: JobPhase | |
| aln_fasta: Optional[str] = None | |
| newick: Optional[str] = None | |
| stats: Optional[str] = None | |
| error: Optional[str] = None | |
| created_at: float = 0.0 | |
| msa_done_at: Optional[float] = None | |
| done_at: Optional[float] = None | |
| class RunResponse(BaseModel): | |
| job_id: str | |
| status: str | |
| # βββ In-memory store ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| _jobs: dict[str, dict] = {} | |
| _lock = threading.Lock() | |
| def _init(job_id: str, req: PhyloRequest) -> None: | |
| with _lock: | |
| _jobs[job_id] = { | |
| "job_id": job_id, | |
| "method": req.method, | |
| "seq_type": req.seq_type, | |
| "model": req.model, | |
| "bootstrap": req.bootstrap, | |
| "phase": "queued", | |
| "aln_fasta": None, | |
| "newick": None, | |
| "stats": None, | |
| "error": None, | |
| "created_at": time.time(), | |
| "msa_done_at": None, | |
| "done_at": None, | |
| "_req": req.model_dump(), | |
| } | |
| def _patch(job_id: str, **kw) -> None: | |
| with _lock: | |
| if job_id in _jobs: | |
| _jobs[job_id].update(kw) | |
| def _read(job_id: str) -> dict | None: | |
| with _lock: | |
| return dict(_jobs[job_id]) if job_id in _jobs else None | |
| # βββ EBI helpers (same pattern as alignment.py) βββββββββββββββββββββββββββββββ | |
| async def _ebi_submit(url: str, data: dict) -> str: | |
| async with httpx.AsyncClient(timeout=30) as client: | |
| r = await client.post(url + "/run", data=data) | |
| r.raise_for_status() | |
| return r.text.strip() | |
| async def _ebi_poll(url: str, ebi_job: str, interval: float = 2.0, max_polls: int = 150) -> str: | |
| async with httpx.AsyncClient(timeout=10) as client: | |
| for _ in range(max_polls): | |
| await asyncio.sleep(interval) | |
| r = await client.get(f"{url}/status/{ebi_job}") | |
| status = r.text.strip() | |
| if status in ("FINISHED", "FAILED", "ERROR", "NOT_FOUND"): | |
| return status | |
| return "TIMEOUT" | |
| async def _ebi_result(url: str, ebi_job: str, result_type: str, retries: int = 3) -> str: | |
| async with httpx.AsyncClient(timeout=60) as client: | |
| for attempt in range(retries): | |
| try: | |
| r = await client.get(f"{url}/result/{ebi_job}/{result_type}") | |
| if r.status_code == 200: | |
| return r.text | |
| await asyncio.sleep(2 ** attempt) | |
| except httpx.RequestError: | |
| await asyncio.sleep(2 ** attempt) | |
| raise RuntimeError(f"Failed to fetch EBI result {result_type} after {retries} attempts") | |
| # βββ MSA via Clustal Omega ββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| def _to_fasta(sequences: list[dict]) -> str: | |
| return "\n".join(f">{s['id']}\n{s['sequence']}" for s in sequences) | |
| async def _run_clustalo(job_id: str, sequences: list[dict], stype: str) -> tuple[str, str] | None: | |
| _patch(job_id, phase="msa_running") | |
| fasta = _to_fasta(sequences) | |
| try: | |
| ebi_job = await _ebi_submit(_EBI_CLUSTALO, { | |
| "email": _EMAIL, | |
| "sequence": fasta, | |
| "outfmt": "fa", | |
| "stype": "protein" if stype == "protein" else "dna", | |
| }) | |
| logger.info(f"[{job_id}] Clustal Omega job: {ebi_job}") | |
| except Exception as e: | |
| _patch(job_id, phase="error", error=f"Clustal Omega submission failed: {e}") | |
| return None | |
| status = await _ebi_poll(_EBI_CLUSTALO, ebi_job) | |
| if status != "FINISHED": | |
| _patch(job_id, phase="error", error=f"Clustal Omega ended with status: {status}") | |
| return None | |
| try: | |
| aln_fasta = await _ebi_result(_EBI_CLUSTALO, ebi_job, "fa") | |
| nj_newick = await _ebi_result(_EBI_CLUSTALO, ebi_job, "phylotree") | |
| except Exception as e: | |
| _patch(job_id, phase="error", error=f"Clustal Omega result fetch failed: {e}") | |
| return None | |
| _patch(job_id, phase="msa_done", aln_fasta=aln_fasta, msa_done_at=time.time()) | |
| return aln_fasta, nj_newick | |
| # βββ UPGMA (pure Python) ββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| def _parse_aligned_fasta(fasta: str) -> dict[str, str]: | |
| seqs: dict[str, str] = {} | |
| cur = None | |
| for line in fasta.strip().splitlines(): | |
| stripped = line.strip() | |
| if stripped.startswith(">"): | |
| cur = stripped[1:].split()[0] | |
| seqs[cur] = "" | |
| elif cur: | |
| seqs[cur] += stripped | |
| return seqs | |
| def _p_distance(s1: str, s2: str) -> float: | |
| pairs = [(a, b) for a, b in zip(s1, s2) if a != "-" and b != "-"] | |
| if not pairs: | |
| return 1.0 | |
| return sum(1 for a, b in pairs if a != b) / len(pairs) | |
| def _upgma_newick(aln_fasta: str) -> str: | |
| seqs = _parse_aligned_fasta(aln_fasta) | |
| names = list(seqs.keys()) | |
| n = len(names) | |
| if n < 2: | |
| return f"({names[0]}:0.0);" if names else "();" | |
| dist: dict[tuple[str, str], float] = {} | |
| for i in range(n): | |
| for j in range(i + 1, n): | |
| d = _p_distance(seqs[names[i]], seqs[names[j]]) | |
| dist[(names[i], names[j])] = d | |
| dist[(names[j], names[i])] = d | |
| clusters: dict[str, dict] = { | |
| nm: {"newick": nm, "height": 0.0, "size": 1} for nm in names | |
| } | |
| counter = 0 | |
| while len(clusters) > 1: | |
| ckeys = list(clusters.keys()) | |
| min_d = float("inf") | |
| best = ("", "") | |
| for i in range(len(ckeys)): | |
| for j in range(i + 1, len(ckeys)): | |
| a, b = ckeys[i], ckeys[j] | |
| d = dist.get((a, b), float("inf")) | |
| if d < min_d: | |
| min_d = d | |
| best = (a, b) | |
| a, b = best | |
| new_h = min_d / 2.0 | |
| bl_a = max(0.0, new_h - clusters[a]["height"]) | |
| bl_b = max(0.0, new_h - clusters[b]["height"]) | |
| new_nw = f"({clusters[a]['newick']}:{bl_a:.6f},{clusters[b]['newick']}:{bl_b:.6f})" | |
| new_sz = clusters[a]["size"] + clusters[b]["size"] | |
| counter += 1 | |
| new_id = f"__c{counter}" | |
| for c in ckeys: | |
| if c in (a, b): | |
| continue | |
| da = dist.get((a, c), dist.get((c, a), 0.0)) | |
| db = dist.get((b, c), dist.get((c, b), 0.0)) | |
| nd = (da * clusters[a]["size"] + db * clusters[b]["size"]) / new_sz | |
| dist[(new_id, c)] = nd | |
| dist[(c, new_id)] = nd | |
| del clusters[a] | |
| del clusters[b] | |
| clusters[new_id] = {"newick": new_nw, "height": new_h, "size": new_sz} | |
| root = next(iter(clusters.values())) | |
| return root["newick"] + ";" | |
| # βββ PhyML local (subprocess) ββββββββββββββββββββββββββββββββββββββββββββββββ | |
| def fasta_to_phylip(fasta: str) -> str: | |
| """Convert aligned FASTA to relaxed PHYLIP (names up to 100 chars).""" | |
| seqs: dict[str, str] = {} | |
| cur: str | None = None | |
| for line in fasta.strip().splitlines(): | |
| t = line.strip() | |
| if t.startswith(">"): | |
| cur = t[1:].split()[0][:100] | |
| seqs[cur] = "" | |
| elif cur: | |
| seqs[cur] += t.upper() | |
| if not seqs: | |
| return "" | |
| n = len(seqs) | |
| L = len(next(iter(seqs.values()))) | |
| lines = [f"{n} {L}"] | |
| for name, s in seqs.items(): | |
| lines.append(f"{name:<100}{s}") | |
| return "\n".join(lines) | |
| async def _run_phyml_local(job_id: str, aln_fasta: str, req: PhyloRequest) -> None: | |
| """Run PhyML as a subprocess on a temp PHYLIP file.""" | |
| _patch(job_id, phase="tree_running") | |
| import os | |
| import shutil | |
| import tempfile | |
| phyml_path = shutil.which("phyml") | |
| if not phyml_path: | |
| _patch(job_id, phase="error", | |
| error="PhyML binary not found. ML method requires PhyML compiled from " | |
| "https://github.com/stephaneguindon/phyml. Try NJ or UPGMA instead.") | |
| return | |
| fd, phy_path = tempfile.mkstemp(suffix=".phy") | |
| os.close(fd) | |
| try: | |
| with open(phy_path, "w") as f: | |
| f.write(fasta_to_phylip(aln_fasta)) | |
| datatype = "aa" if req.seq_type == "protein" else "nt" | |
| model = req.model if req.model else ("LG" if req.seq_type == "protein" else "GTR") | |
| proc = await asyncio.create_subprocess_exec( | |
| "phyml", | |
| "-i", phy_path, | |
| "-d", datatype, | |
| "-m", model, | |
| "-b", str(req.bootstrap if req.bootstrap else 0), | |
| "-o", "tlr", | |
| "--no_memory_check", | |
| stdout=asyncio.subprocess.PIPE, | |
| stderr=asyncio.subprocess.PIPE, | |
| ) | |
| try: | |
| _, stderr = await asyncio.wait_for(proc.communicate(), timeout=900) | |
| except asyncio.TimeoutError: | |
| proc.kill() | |
| await proc.communicate() | |
| _patch(job_id, phase="error", error="PhyML timed out after 15 minutes") | |
| return | |
| if proc.returncode != 0: | |
| err_text = stderr.decode("utf-8", errors="replace")[:500] if stderr else "" | |
| _patch(job_id, phase="error", | |
| error=f"PhyML failed (exit {proc.returncode}): {err_text}") | |
| return | |
| tree_path = phy_path + "_phyml_tree.txt" | |
| stats_path = phy_path + "_phyml_stats.txt" | |
| newick = None | |
| stats = None | |
| if os.path.exists(tree_path): | |
| with open(tree_path) as f: | |
| newick = f.read().strip() | |
| if os.path.exists(stats_path): | |
| with open(stats_path) as f: | |
| stats = f.read().strip() | |
| if not newick: | |
| _patch(job_id, phase="error", error="PhyML produced no output tree") | |
| return | |
| _patch(job_id, phase="complete", newick=newick, stats=stats, done_at=time.time()) | |
| except Exception as e: | |
| _patch(job_id, phase="error", error=f"PhyML error: {e}") | |
| finally: | |
| for suffix in ["", "_phyml_tree.txt", "_phyml_stats.txt", "_phyml_boot_trees.txt"]: | |
| p = phy_path + suffix | |
| if os.path.exists(p): | |
| os.remove(p) | |
| # βββ Main pipeline worker βββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| async def _worker(job_id: str) -> None: | |
| job = _read(job_id) | |
| if not job: | |
| return | |
| req = PhyloRequest(**job["_req"]) | |
| result = await _run_clustalo(job_id, req.sequences, req.seq_type) | |
| if result is None: | |
| return | |
| aln_fasta, nj_newick = result | |
| if req.method == "nj": | |
| _patch(job_id, phase="complete", | |
| newick=nj_newick.strip(), done_at=time.time()) | |
| elif req.method == "upgma": | |
| _patch(job_id, phase="tree_running") | |
| try: | |
| newick = _upgma_newick(aln_fasta) | |
| _patch(job_id, phase="complete", | |
| newick=newick, done_at=time.time()) | |
| except Exception as e: | |
| _patch(job_id, phase="error", error=f"UPGMA computation failed: {e}") | |
| elif req.method == "ml": | |
| await _run_phyml_local(job_id, aln_fasta, req) | |
| else: | |
| _patch(job_id, phase="error", error=f"Unknown method: {req.method}") | |
| # βββ API endpoints βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ | |
| async def run_phylo( | |
| req: PhyloRequest, | |
| background_tasks: BackgroundTasks, | |
| ): | |
| if len(req.sequences) < 2: | |
| raise HTTPException(400, detail="At least 2 sequences are required") | |
| if len(req.sequences) > 50: | |
| raise HTTPException(400, detail="Maximum 50 sequences per run") | |
| valid_models = PROTEIN_MODELS if req.seq_type == "protein" else DNA_MODELS | |
| if req.method == "ml" and req.model not in valid_models: | |
| raise HTTPException( | |
| 400, | |
| detail=f"Model '{req.model}' not valid for {req.seq_type}. " | |
| f"Choose from: {', '.join(valid_models)}" | |
| ) | |
| job_id = str(uuid.uuid4()) | |
| _init(job_id, req) | |
| background_tasks.add_task(_worker, job_id) | |
| return RunResponse(job_id=job_id, status="queued") | |
| async def get_status(job_id: str): | |
| job = _read(job_id) | |
| if not job: | |
| raise HTTPException(404, detail=f"Job {job_id} not found") | |
| return {k: v for k, v in job.items() if not k.startswith("_")} | |
| async def list_models(seq_type: SeqType = "protein"): | |
| return { | |
| "seq_type": seq_type, | |
| "models": PROTEIN_MODELS if seq_type == "protein" else DNA_MODELS, | |
| "default": "LG" if seq_type == "protein" else "GTR", | |
| } | |