fix(blast): poll cap 65min, DNA validation, program/db/max_hits params
cd0c7a9 verified 3 months ago
27.4 kB
from __future__ import annotations

import json
import math
import re
from fastapi import APIRouter, HTTPException, Depends, Request
from pydantic import BaseModel, Field
from typing import Any, Optional

from app.services.supabase import get_client
from app.services.auth import require_user_id
from app.services.ssrf import validate_url
router = APIRouter(prefix="/api/docking", tags=["Docking"])
_TABLE = "docking_jobs"


# ---------------------------------------------------------------------------
# Request / response schemas (match frontend DockingResult type)
# ---------------------------------------------------------------------------

class DockingJobCreate(BaseModel):
    pdb_id: str = ""
    smiles: str
    pdb_url: str = ""
    grid_center: Optional[list[float]] = None
    grid_size: list[float] = Field(default_factory=lambda: [20.0, 20.0, 20.0])
    exhaustiveness: int = 8
    num_modes: int = 9


class DockingJobResponse(BaseModel):
    job_id: str
    status: str
    result: Optional[dict[str, Any]] = None
    error: Optional[str] = None


# ---------------------------------------------------------------------------
# Helpers
# ---------------------------------------------------------------------------

def _prune_old(supabase, max_rows: int = 200):
    try:
        rows = (
            supabase.table(_TABLE)
            .select("id")
            .order("created_at", desc=True)
            .range(max_rows, max_rows + 1000)
            .execute()
            .data
        )
        if rows:
            supabase.table(_TABLE).delete().in_(
                "id", [r["id"] for r in rows]
            ).execute()
    except Exception:
        pass


def _row_to_response(row: dict) -> dict:
    """Convert a Supabase row to the frontend DockingResult shape."""
    result = None

    # Prefer Storage URL (Phase 0c)
    storage_url = row.get("storage_url")
    if storage_url:
        from app.services.artifact_storage import download_json
        result = download_json(storage_url)
    elif row.get("result_sdf"):
        try:
            result = json.loads(row["result_sdf"])
        except Exception:
            pass

    return {
        "job_id": row["id"],
        "status": row["status"],
        "result": result,
        "error": row.get("error"),
    }


def _row_to_list_response(row: dict) -> dict:
    """Lightweight row conversion for list views — skips Storage downloads."""
    return {
        "job_id": row["id"],
        "status": row["status"],
        "result": None,
        "error": row.get("error"),
    }


# ---------------------------------------------------------------------------
# Background worker
# ---------------------------------------------------------------------------

def _run_docking_sync(job_id: str, payload: dict):
    """Run the full docking pipeline synchronously (in a thread)."""
    supabase = get_client()
    try:
        supabase.table(_TABLE).update({"status": "running"}).eq("id", job_id).execute()

        from app.tools.docking import (
            fetch_pdb_from_rcsb,
            compute_grid_center,
            smiles_to_pdbqt,
            pdb_to_pdbqt_receptor,
            run_vina,
        )
        import urllib.request
        from app.services.ssrf import validate_url

        pdb_id = payload.get("pdb_id", "").strip().upper()
        pdb_url = payload.get("pdb_url", "").strip()
        smiles = payload.get("ligand_smiles") or payload.get("smiles")
        if not smiles:
            raise ValueError("Missing ligand_smiles in job payload")

        # 1. Obtain PDB text
        pdb_text: str | None = None
        if pdb_url:
            validate_url(pdb_url)  # SSRF guard even in worker
            try:
                pdb_text = urllib.request.urlopen(pdb_url, timeout=30).read().decode("utf-8", errors="replace")
            except Exception:
                pass
        if not pdb_text and pdb_id:
            pdb_text = fetch_pdb_from_rcsb(pdb_id)

        if not pdb_text:
            raise RuntimeError(
                "Could not obtain a PDB structure. "
                "Provide a valid pdb_id or pdb_url."
            )

        # 2. Strip heteroatoms (keep protein backbone for receptor)
        protein_lines = [
            l for l in pdb_text.splitlines()
            if l.startswith("ATOM") or l.startswith("TER") or l.startswith("END")
        ]
        protein_pdb = "\n".join(protein_lines) if protein_lines else pdb_text

        # 3. Compute grid center if not provided
        grid_center = payload.get("grid_center")
        if not grid_center or all(v == 0 for v in grid_center):
            grid_center = compute_grid_center(protein_pdb)
            # Add a small offset so the center isn't dead on a backbone atom
            grid_center = [round(c + 2.0, 3) for c in grid_center]

        grid_size = payload.get("grid_size", [20.0, 20.0, 20.0])

        # 4. Prepare receptor
        protein_pdbqt = pdb_to_pdbqt_receptor(protein_pdb)

        # 5. Prepare ligand
        lig_pdbqt = smiles_to_pdbqt(smiles)

        # 6. Run AutoDock Vina
        vina_result = run_vina(
            protein_pdbqt=protein_pdbqt,
            ligand_pdbqt=lig_pdbqt,
            grid_center=grid_center,
            grid_size=grid_size,
            exhaustiveness=payload.get("exhaustiveness", 8),
            num_modes=payload.get("num_modes", 9),
        )

        # 7. Compute interaction summary for best pose
        interactions = _compute_interactions(
            protein_pdb, vina_result["ligand_pdb"]
        )
        pose_interactions = _summarize_pose_interactions(
            protein_pdb, vina_result.get("result_sdf", "")
        )

        result_obj = {
            "pdb_id": pdb_id,
            "smiles": smiles,
            "poses": vina_result["poses"],
            "num_poses": vina_result["num_poses"],
            "box_center": {
                "x": grid_center[0],
                "y": grid_center[1],
                "z": grid_center[2],
            },
            "box_size": {
                "x": grid_size[0],
                "y": grid_size[1],
                "z": grid_size[2],
            },
            "vina_log": vina_result.get("vina_log", ""),
            "interactions": interactions,
            "pose_interactions": pose_interactions,
            "ligand_pdb": vina_result.get("ligand_pdb", ""),
        }

        # Offload to Supabase Storage; DB keeps only the URL
        from app.services.artifact_storage import upload_json
        storage_url = upload_json(job_id, "result", result_obj)

        supabase.table(_TABLE).update({
            "status": "complete",
            "storage_url": storage_url,
            "result_sdf": None,  # cleared — data lives in Storage now
        }).eq("id", job_id).execute()

    except Exception as exc:
        import traceback
        tb = traceback.format_exc()
        supabase.table(_TABLE).update({
            "status": "failed",
            "error": f"{exc}\n\n{tb}"[:4000],
        }).eq("id", job_id).execute()
    finally:
        _prune_old(supabase)


# ---------------------------------------------------------------------------
# Geometric interaction detector (H-bonds, hydrophobic, pi-stacking, salt bridges)
# ---------------------------------------------------------------------------

    # Protein atom classification
_HYDROPHOBIC_RES = {"ALA", "VAL", "LEU", "ILE", "MET", "PHE", "TRP", "PRO", "GLY"}
_AROMATIC_RES = {"PHE", "TRP", "TYR", "HIS"}

# Atoms in aromatic rings by residue (PDB atom names)
_AROMATIC_RING_ATOMS = {
    "PHE": ["CG", "CD1", "CD2", "CE1", "CE2", "CZ"],
    "TYR": ["CG", "CD1", "CD2", "CE1", "CE2", "CZ"],
    "HIS": ["CG", "ND1", "CD2", "CE1", "NE2"],
    "TRP": ["CG", "CD1", "CD2", "NE1", "CE2", "CE3", "CZ2", "CZ3", "CH2"],
}

# Two-ring centroids for TRP (5-membered + 6-membered)
_TRP_RING_ATOMS = {
    "five": ["CD1", "NE1", "CE2", "CG", "CD2"],
    "six": ["CE2", "CD2", "CZ2", "CH2", "CZ3", "CE3"],
}

# Polar atoms eligible for H-bonding
_POLAR_ATOMS = {"N", "O", "S"}

# Residue-level charge groups for salt bridges
_ANIONIC_RES = {"ASP", "GLU"}
_CATIONIC_RES = {"LYS", "ARG", "HIS"}

# Atom names that define the charged group center
_ANIONIC_CARBONS = {"ASP": "CG", "GLU": "CD"}
_CATIONIC_NITROGENS = {"LYS": "NZ", "ARG": ["CZ", "NH1", "NH2"]}

_PDB_COORD_RE = re.compile(
    r"^(ATOM|HETATM)\s+\d+\s+(\S+)\s+(\S{3})\s+(\S)\s+(\d+)\s+"
    r"([-\d.]+)\s+([-\d.]+)\s+([-\d.]+)"
)


def _parse_atom_coords(pdb_text: str) -> list[tuple[str, str, str, str, int, float, float, float]]:
    """Parse PDB into (record, atom_name, res_name, chain, res_seq, x, y, z)."""
    atoms = []
    for line in pdb_text.splitlines():
        m = _PDB_COORD_RE.match(line)
        if m:
            atoms.append((
                m.group(1), m.group(2), m.group(3), m.group(4),
                int(m.group(5)),
                float(m.group(6)), float(m.group(7)), float(m.group(8)),
            ))
    return atoms


def _distance(a: tuple[float, float, float], b: tuple[float, float, float]) -> float:
    return math.sqrt(sum((x - y) ** 2 for x, y in zip(a, b)))


def _angle(a: tuple[float, float, float], b: tuple[float, float, float],
           c: tuple[float, float, float]) -> float:
    """Angle at vertex b between segments b→a and b→c, in degrees."""
    ba = tuple(x - y for x, y in zip(a, b))
    bc = tuple(x - y for x, y in zip(c, b))
    dot = sum(x * y for x, y in zip(ba, bc))
    mag_ba = math.sqrt(sum(x * x for x in ba))
    mag_bc = math.sqrt(sum(x * x for x in bc))
    if mag_ba < 1e-9 or mag_bc < 1e-9:
        return 0.0
    cos_angle = max(-1.0, min(1.0, dot / (mag_ba * mag_bc)))
    return math.degrees(math.acos(cos_angle))


def _vec_sub(a: tuple[float, float, float], b: tuple[float, float, float]) -> tuple[float, float, float]:
    return (a[0] - b[0], a[1] - b[1], a[2] - b[2])


def _vec_cross(a: tuple[float, float, float], b: tuple[float, float, float]) -> tuple[float, float, float]:
    return (
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0],
    )


def _vec_norm(v: tuple[float, float, float]) -> float:
    return math.sqrt(sum(x * x for x in v))


def _ring_centroid(coords: list[tuple[float, float, float]]) -> tuple[float, float, float]:
    n = len(coords)
    if n == 0:
        return (0.0, 0.0, 0.0)
    return (
        sum(c[0] for c in coords) / n,
        sum(c[1] for c in coords) / n,
        sum(c[2] for c in coords) / n,
    )


def _ring_normal(coords: list[tuple[float, float, float]]) -> tuple[float, float, float]:
    """Compute the normal vector of a planar ring via cross product of two edges."""
    if len(coords) < 3:
        return (0.0, 0.0, 1.0)
    v1 = _vec_sub(coords[1], coords[0])
    v2 = _vec_sub(coords[2], coords[0])
    cross = _vec_cross(v1, v2)
    n = _vec_norm(cross)
    if n < 1e-9:
        return (0.0, 0.0, 1.0)
    return (cross[0] / n, cross[1] / n, cross[2] / n)


def _build_residue_map(atoms: list[tuple]) -> dict[tuple[str, str, int], list[tuple]]:
    """Group atoms by (chain, res_name, res_seq)."""
    res_map: dict[tuple[str, str, int], list[tuple]] = {}
    for a in atoms:
        key = (a[3], a[2], a[4])  # chain, res_name, res_seq
        res_map.setdefault(key, []).append(a)
    return res_map


def _find_hydrogens(atoms: list[tuple]) -> list[tuple]:
    """Return only hydrogen atoms from parsed PDB."""
    return [a for a in atoms if a[1].startswith("H") or a[1] in ("1H", "2H", "3H")]


def _compute_interactions(protein_pdb: str, ligand_pdb: str) -> dict:
    """
    Compute protein-ligand interactions using proper geometry.
    H-bonds:     donor-H···acceptor angle > 120°, distance < 3.5Å
    Hydrophobic: ligand carbon near protein carbon in hydrophobic residue, < 4.5Å
    Pi-stacking: aromatic ring centroids, distance < 5.5Å, inter-ring angle
    Salt bridges: charged group centroids, distance < 4.0Å
    """
    if not ligand_pdb:
        return {"hbonds": [], "hydrophobic": [], "pi_stacking": [], "salt_bridges": []}

    prot_atoms = _parse_atom_coords(protein_pdb)
    lig_atoms = _parse_atom_coords(ligand_pdb)
    prot_h = _find_hydrogens(prot_atoms)
    lig_h = _find_hydrogens(lig_atoms)
    prot_heavy = [a for a in prot_atoms if not (a[1].startswith("H") or a[1] in ("1H", "2H", "3H"))]
    lig_heavy = [a for a in lig_atoms if not (a[1].startswith("H") or a[1] in ("1H", "2H", "3H"))]

    hbonds: list[dict] = []
    hydrophobic: list[dict] = []
    pi_stacking: list[dict] = []
    salt_bridges: list[dict] = []

    seen_hbonds: set[tuple] = set()
    seen_hydrophobic: set[tuple] = set()
    seen_salt: set[tuple] = set()

    # --- H-bonds with angle check ---
    for la in lig_heavy:
        l_elem = la[1][0] if la[1] else ""
        if l_elem not in _POLAR_ATOMS:
            continue
        lcoord = (la[5], la[6], la[7])

        # Find nearest H on ligand for angle reference
        lig_h_near = None
        min_h_dist = 1.5
        for h in lig_h:
            hd = _distance(lcoord, (h[5], h[6], h[7]))
            if hd < min_h_dist:
                min_h_dist = hd
                lig_h_near = (h[5], h[6], h[7])

        for pa in prot_heavy:
            p_elem = pa[1][0] if pa[1] else ""
            if p_elem not in _POLAR_ATOMS:
                continue
            pcoord = (pa[5], pa[6], pa[7])
            d = _distance(lcoord, pcoord)

            if d > 3.5 or d < 1.0:
                continue

            # Find nearest H on protein donor for angle check
            prot_h_near = None
            min_ph_dist = 1.5
            for h in prot_h:
                hd = _distance(pcoord, (h[5], h[6], h[7]))
                if hd < min_ph_dist:
                    min_ph_dist = hd
                    prot_h_near = (h[5], h[6], h[7])

            # Check angle if we have hydrogen positions
            angle_ok = True
            if lig_h_near and prot_h_near:
                # H-bond angle: ligand-H···protein or protein-H···ligand
                a1 = _angle(lig_h_near, lcoord, pcoord)
                a2 = _angle(prot_h_near, pcoord, lcoord)
                angle_ok = max(a1, a2) > 120.0
            elif lig_h_near:
                a1 = _angle(lig_h_near, lcoord, pcoord)
                angle_ok = a1 > 120.0
            elif prot_h_near:
                a1 = _angle(prot_h_near, pcoord, lcoord)
                angle_ok = a1 > 120.0
            # If no H found at all, accept based on distance + element only

            if not angle_ok:
                continue

            key = (la[4], pa[4])  # (lig_res_seq, prot_res_seq)
            if key in seen_hbonds:
                continue
            seen_hbonds.add(key)

            hbonds.append({
                "type": "hbond",
                "ligand_atom": la[1],
                "ligand_coords": [la[5], la[6], la[7]],
                "protein_residue": pa[2],
                "protein_residue_seq": pa[4],
                "protein_chain": pa[3],
                "protein_atom": pa[1],
                "protein_coords": [pa[5], pa[6], pa[7]],
                "distance": round(d, 2),
                "confidence": "high" if d < 3.0 else "medium",
            })
            if len(hbonds) >= 20:
                break
        if len(hbonds) >= 20:
            break

    # --- Hydrophobic contacts ---
    for la in lig_heavy:
        if la[1][0] != "C":
            continue
        lcoord = (la[5], la[6], la[7])
        for pa in prot_heavy:
            if pa[1][0] != "C":
                continue
            pres = pa[2]
            if pres not in _HYDROPHOBIC_RES:
                continue
            pcoord = (pa[5], pa[6], pa[7])
            d = _distance(lcoord, pcoord)
            if d < 4.5:
                key = (la[4], pa[4])
                if key in seen_hydrophobic:
                    continue
                seen_hydrophobic.add(key)
                hydrophobic.append({
                    "type": "hydrophobic",
                    "ligand_atom": la[1],
                    "ligand_coords": [la[5], la[6], la[7]],
                    "protein_residue": pres,
                    "protein_residue_seq": pa[4],
                    "protein_chain": pa[3],
                    "protein_atom": pa[1],
                    "protein_coords": [pa[5], pa[6], pa[7]],
                    "distance": round(d, 2),
                })
                if len(hydrophobic) >= 20:
                    break
        if len(hydrophobic) >= 20:
            break

    # --- Pi-stacking (aromatic ring centroid geometry) ---
    prot_res_map = _build_residue_map(prot_heavy)

    for res_key, res_atoms in prot_res_map.items():
        chain, res_name, res_seq = res_key
        if res_name not in _AROMATIC_RES:
            continue

        ring_atom_names = _AROMATIC_RING_ATOMS[res_name]
        ring_atoms_by_name = {a[1]: a for a in res_atoms}
        ring_coords = []
        for rn in ring_atom_names:
            if rn in ring_atoms_by_name:
                a = ring_atoms_by_name[rn]
                ring_coords.append((a[5], a[6], a[7]))

        if len(ring_coords) < 3:
            continue

        centroid = _ring_centroid(ring_coords)
        normal = _ring_normal(ring_coords)

        # For TRP, also check the 5-membered ring
        rings_to_check = [(ring_coords, centroid, normal)]
        if res_name == "TRP":
            for ring_name in ("five", "six"):
                ring_atom_names_2 = _TRP_RING_ATOMS[ring_name]
                coords_2 = []
                for rn in ring_atom_names_2:
                    if rn in ring_atoms_by_name:
                        a = ring_atoms_by_name[rn]
                        coords_2.append((a[5], a[6], a[7]))
                if len(coords_2) >= 3:
                    rings_to_check.append((coords_2, _ring_centroid(coords_2), _ring_normal(coords_2)))

        for ring_coords_r, centroid_r, normal_r in rings_to_check:
            # Find aromatic atoms in ligand (heuristic: C/N in a flat region)
            lig_aromatic_coords = []
            for la in lig_heavy:
                if la[1][0] in ("C", "N"):
                    lig_aromatic_coords.append((la[5], la[6], la[7]))

            if len(lig_aromatic_coords) < 3:
                continue

            # Use all ligand heavy atoms as a pseudo-centroid
            lig_centroid = _ring_centroid(lig_aromatic_coords)

            dist = _distance(centroid_r, lig_centroid)
            if dist > 6.5:
                continue

            # Compute angle between ring normal and vector to ligand centroid
            v_to_lig = _vec_sub(lig_centroid, centroid_r)
            v_norm = _vec_norm(v_to_lig)
            if v_norm < 1e-9:
                continue
            cos_angle = abs(sum(x * y for x, y in zip(normal_r, v_to_lig))) / (
                _vec_norm(normal_r) * v_norm
            )
            ring_angle = math.degrees(math.acos(max(0, min(1, cos_angle))))

            # Parallel: ring normal ~parallel to centroid-centroid vector (angle < 30°)
            # T-shaped: ring normal ~perpendicular (angle 60-90°)
            stacking_type = "unknown"
            if ring_angle < 30 and dist < 5.5:
                stacking_type = "parallel"
            elif 60 < ring_angle < 90 and dist < 6.5:
                stacking_type = "perpendicular"

            if stacking_type == "unknown":
                continue

            pi_stacking.append({
                "type": "pi_stacking",
                "protein_residue": res_name,
                "protein_residue_seq": res_seq,
                "protein_chain": chain,
                "ring_centroid": [round(c, 3) for c in centroid_r],
                "ring_normal": [round(c, 3) for c in normal_r],
                "ligand_centroid": [round(c, 3) for c in lig_centroid],
                "distance": round(dist, 2),
                "angle": round(ring_angle, 1),
                "stacking_type": stacking_type,
                "confidence": "high" if dist < 4.5 else "medium",
            })
            if len(pi_stacking) >= 10:
                break
        if len(pi_stacking) >= 10:
            break

    # --- Salt bridges (charged group centroid distance) ---
    for la in lig_heavy:
        l_elem = la[1][0] if la[1] else ""
        if l_elem not in ("N", "O", "S", "C"):
            continue
        lcoord = (la[5], la[6], la[7])

        for pa in prot_heavy:
            pres = pa[2]
            if pres in _ANIONIC_RES and pa[1] in ("OD1", "OD2", "OE1", "OE2"):
                d = _distance(lcoord, pa[1:8] if False else (pa[5], pa[6], pa[7]))
                if d < 4.0 and l_elem in ("N",):
                    key = (la[4], pa[4])
                    if key not in seen_salt:
                        seen_salt.add(key)
                        salt_bridges.append({
                            "type": "salt_bridge",
                            "ligand_atom": la[1],
                            "ligand_coords": [la[5], la[6], la[7]],
                            "protein_residue": pres,
                            "protein_residue_seq": pa[4],
                            "protein_chain": pa[3],
                            "protein_atom": pa[1],
                            "protein_coords": [pa[5], pa[6], pa[7]],
                            "distance": round(d, 2),
                            "charge_pair": "positive-negative",
                        })

            if pres in _CATIONIC_RES:
                cat_atoms = _CATIONIC_NITROGENS.get(pres, [])
                if isinstance(cat_atoms, str):
                    cat_atoms = [cat_atoms]
                if pa[1] in cat_atoms:
                    d = _distance(lcoord, (pa[5], pa[6], pa[7]))
                    if d < 4.0 and l_elem in ("O",):
                        key = (la[4], pa[4])
                        if key not in seen_salt:
                            seen_salt.add(key)
                            salt_bridges.append({
                                "type": "salt_bridge",
                                "ligand_atom": la[1],
                                "ligand_coords": [la[5], la[6], la[7]],
                                "protein_residue": pres,
                                "protein_residue_seq": pa[4],
                                "protein_chain": pa[3],
                                "protein_atom": pa[1],
                                "protein_coords": [pa[5], pa[6], pa[7]],
                                "distance": round(d, 2),
                                "charge_pair": "negative-positive",
                            })

    return {
        "hbonds": hbonds[:20],
        "hydrophobic": hydrophobic[:20],
        "pi_stacking": pi_stacking[:10],
        "salt_bridges": salt_bridges[:10],
    }


def _summarize_pose_interactions(protein_pdb: str, output_pdbqt: str) -> list[dict]:
    """Per-pose interaction summary."""
    if not output_pdbqt:
        return []
    models: dict[int, list[str]] = {}
    current: int | None = None
    for line in output_pdbqt.splitlines():
        if line.startswith("MODEL"):
            parts = line.split()
            if len(parts) >= 2:
                current = int(parts[1])
                models[current] = []
        elif line.startswith("ENDMDL"):
            current = None
        elif current is not None:
            models.setdefault(current, []).append(line)

    summaries = []
    for mid in sorted(models.keys()):
        lig_pdb = "\n".join(l for l in models[mid] if l.startswith("HETATM")) + "\nEND"
        inter = _compute_interactions(protein_pdb, lig_pdb)
        summaries.append({
            "model": mid,
            "hbonds": len(inter.get("hbonds", [])),
            "hydrophobic": len(inter.get("hydrophobic", [])),
            "pi_stacking": len(inter.get("pi_stacking", [])),
            "salt_bridges": len(inter.get("salt_bridges", [])),
        })
    return summaries


# ---------------------------------------------------------------------------
# API endpoints
# ---------------------------------------------------------------------------

@router.post("/run", response_model=DockingJobResponse)
async def create_docking_job(request: Request, body: DockingJobCreate, user_id: str = Depends(require_user_id)):
    supabase = get_client()
    _prune_old(supabase)

    # SSRF validation on user-supplied URL
    if body.pdb_url:
        validate_url(body.pdb_url)

    import uuid, datetime
    job_id = str(uuid.uuid4())
    now = datetime.datetime.utcnow().isoformat()

    insert_row = {
        "id": job_id,
        "status": "queued",
        "ligand_smiles": body.smiles,
        "user_id": user_id,
        "payload": {
            "pdb_id": body.pdb_id,
            "pdb_url": body.pdb_url,
            "grid_center": body.grid_center or [0, 0, 0],
            "grid_size": body.grid_size,
            "exhaustiveness": body.exhaustiveness,
            "num_modes": body.num_modes,
            "smiles": body.smiles,
            "ligand_smiles": body.smiles,
        },
    }
    try:
        supabase.table(_TABLE).insert(insert_row).execute()
    except Exception as e:
        if "ligand_smiles" in str(e):
            supabase.table(_TABLE).insert({
                "id": job_id, "status": "queued", "user_id": user_id,
                "payload": insert_row["payload"],
            }).execute()
        else:
            raise

    return DockingJobResponse(job_id=job_id, status="queued", result=None)


@router.get("/status/{job_id}", response_model=DockingJobResponse)
async def get_docking_job(job_id: str, user_id: str = Depends(require_user_id)):
    supabase = get_client()
    result = supabase.table(_TABLE).select("*").eq("id", job_id).eq("user_id", user_id).single().execute()
    if not result.data:
        raise HTTPException(status_code=404, detail="Docking job not found")
    return DockingJobResponse(**_row_to_response(result.data))


@router.get("/result/{job_id}/pdb")
async def get_docking_pdb(job_id: str, user_id: str = Depends(require_user_id)):
    supabase = get_client()
    row = supabase.table(_TABLE).select("result_sdf,storage_url").eq("id", job_id).eq("user_id", user_id).single().execute()
    if not row.data:
        raise HTTPException(status_code=404, detail="Docking result not found")

    data = None
    if row.data.get("storage_url"):
        from app.services.artifact_storage import download_json
        data = download_json(row.data["storage_url"])
    elif row.data.get("result_sdf"):
        try:
            data = json.loads(row.data["result_sdf"])
        except Exception:
            pass

    if not data:
        raise HTTPException(status_code=404, detail="Docking result not found")
    ligand_pdb = data.get("ligand_pdb", "")
    if not ligand_pdb:
        raise HTTPException(status_code=404, detail="No ligand PDB available")
    from fastapi.responses import PlainTextResponse
    return PlainTextResponse(ligand_pdb, media_type="text/plain")


@router.get("")
async def list_docking_jobs(limit: int = 50, user_id: str = Depends(require_user_id)):
    supabase = get_client()
    rows = (
        supabase.table(_TABLE)
        .select("*")
        .eq("user_id", user_id)
        .order("created_at", desc=True)
        .limit(limit)
        .execute()
        .data
    )
    return {"jobs": [_row_to_list_response(r) for r in rows]}