diff --git a/README.md b/README.md index 7d92c33cd646377fc1c80f5a2376ac24d22c66cb..cbfe3f03f8e1332c0d3c141320281016c3eacaf7 100644 --- a/README.md +++ b/README.md @@ -37,7 +37,7 @@ Recommended Dataset repository files: ```text MuSProt.db -MuSProt-node-lookup.db +MuSProt-index.db musprot_summary.json MuSProt_documentation.md plots/ @@ -47,7 +47,7 @@ Optional Space variables: ```text MUSPROT_DB_PATH=/data/MuSProt.db -MUSPROT_NODE_DB_PATH=/data/MuSProt-node-lookup.db +MUSPROT_INDEX_DB_PATH=/data/MuSProt-index.db MUSPROT_SUMMARY_PATH=/data/musprot_summary.json MUSPROT_DOCS_PATH=/data/MuSProt_documentation.md MUSPROT_PLOTS_DIR=/data/plots @@ -61,13 +61,18 @@ For a private Dataset repository, add `HF_TOKEN` as a Space secret. ## Updating the database -When publishing a new `MuSProt.db`, regenerate the bundled summary sidecar: +When publishing a new `MuSProt.db`, rebuild the website index sidecar +(sequence catalog, search, state pairs, ECOD groups and dataset statistics). +It reads the DB once and takes about 6 minutes: ```bash -cd backend -python scripts/generate_summary.py /path/to/MuSProt.db -o assets/musprot_summary.json +python backend/scripts/build_index.py /path/to/MuSProt.db --out-dir \ + --ecod /path/to/ecod.latest.domains.txt ``` +Upload `MuSProt.db`, `MuSProt-index.db` and `MuSProt_documentation.md` to the +storage bucket mounted at `/data`, then restart the Space. + The deployed application does not create indexes or modify the SQLite file. Create any required indexes before uploading the database to the Dataset repository. @@ -81,6 +86,7 @@ npm run build cd ../backend MUSPROT_DB_PATH=/path/to/MuSProt.db \ +MUSPROT_INDEX_DB_PATH=/path/to/MuSProt-index.db \ MUSPROT_FRONTEND_DIR=../frontend/build \ uvicorn space_app:app --host 0.0.0.0 --port 8501 ``` diff --git a/backend/app/protein/api/routes/explore.py b/backend/app/protein/api/routes/explore.py new file mode 100644 index 0000000000000000000000000000000000000000..f11fd3010a8df2b7d54f47c90938a22dad77a99b --- /dev/null +++ b/backend/app/protein/api/routes/explore.py @@ -0,0 +1,143 @@ +"""Explorer API: overview, search, sequence catalog, sequence and chain detail.""" +from __future__ import annotations + +import logging +from typing import Optional + +from fastapi import APIRouter, HTTPException, Query +from starlette.concurrency import run_in_threadpool + +from app.protein import catalog, records + +logger = logging.getLogger(__name__) +router = APIRouter() + + +def _index_required(value): + if value is None: + raise HTTPException(status_code=503, detail="The MuSProt index is not available yet") + return value + + +@router.get("/overview") +async def get_overview(): + """Dataset-level counts and distributions (precomputed).""" + return _index_required(await run_in_threadpool(catalog.overview)) + + +@router.get("/search") +async def search(q: str = Query(..., min_length=1, max_length=200), limit: int = Query(12, ge=1, le=50)): + """Search by chain (3vgk.A), PDB entry, sequence ID, UniProt accession or function text.""" + return await run_in_threadpool(catalog.search, q, limit) + + +@router.get("/sequences") +async def list_sequences( + sort: str = Query("states"), + order: str = Query("desc", pattern="^(asc|desc)$"), + q: Optional[str] = Query(None, max_length=200), + min_states: Optional[int] = Query(None, ge=1), + min_obs: Optional[int] = Query(None, ge=1), + method: Optional[str] = Query(None, max_length=64), + ligand: Optional[str] = Query(None, pattern="^(both|holo|apo)$"), + cath: Optional[str] = Query(None, max_length=32), + offset: int = Query(0, ge=0), + limit: int = Query(25, ge=1, le=100), + ids: Optional[str] = Query(None, max_length=1000, description="Comma-separated sequence IDs"), + ecod: Optional[str] = Query(None, max_length=200, description="ECOD H-group (e.g. 11.1) or exact H-set (11.1;2003.1)"), + pfam: Optional[str] = Query(None, max_length=20), +): + """Browse the sequence catalog with filters, sorting and pagination.""" + id_list = [i.strip().upper() for i in ids.split(",") if i.strip()][:100] if ids else None + result = await run_in_threadpool( + catalog.list_sequences, sort, order, q, min_states, min_obs, method, ligand, cath, offset, limit, id_list, + ecod, pfam, + ) + return _index_required(result) + + +@router.get("/sequence/{sequence_id}") +async def get_sequence(sequence_id: str): + """One sequence: states, member chains and the state-state similarity matrix.""" + seq = _index_required(await run_in_threadpool(catalog.get_sequence, sequence_id)) + if not seq: + raise HTTPException(status_code=404, detail=f"Sequence {sequence_id} not found") + rep = await run_in_threadpool(records.get_node, seq["rep_pdb"], seq["rep_chain"]) + if rep: + seq["sequence"] = rep["sequence"] + seq["functions"] = rep["functions"] + seq["cath_superfamilies"] = rep["cath_superfamilies"] + return seq + + +@router.get("/chain/{pdb_id}/{chain_id}") +async def get_chain(pdb_id: str, chain_id: str): + """Full record for one chain plus a summary of its sequence family.""" + node = await run_in_threadpool(records.get_node, pdb_id, chain_id) + if node is None: + raise HTTPException(status_code=404, detail=f"{pdb_id.upper()}.{chain_id} is not in MuSProt") + node["family"] = await run_in_threadpool(catalog.sequence_summary, node["sequence_id"]) + node["ecod_names"] = await run_in_threadpool(catalog.ecod_names, None, node["ecod_fids"]) + return node + + +@router.get("/chain/{pdb_id}/{chain_id}/transitions") +async def get_chain_transitions(pdb_id: str, chain_id: str): + """Every comparison from this chain to other observations of the same sequence.""" + node = await run_in_threadpool(records.get_node, pdb_id, chain_id) + if node is None: + raise HTTPException(status_code=404, detail=f"{pdb_id.upper()}.{chain_id} is not in MuSProt") + edges = await run_in_threadpool(records.get_transitions, node["pdb_id"], node["auth_asym_id"]) + members = await run_in_threadpool(catalog.sequence_members, node["sequence_id"]) + for e in edges: + m = members.get((e["pdb_id"], e["auth_asym_id"])) + if m: + e.update( + experimental_method=m["experimental_method"], + resolution=m["resolution"], + binding_status=m["binding_status"], + chain_composition=m["chain_composition"], + binders=m["binders"], + release_date=m["release_date"], + ) + return { + "pdb_id": node["pdb_id"], + "auth_asym_id": node["auth_asym_id"], + "state_id": node["state_id"], + "sequence_id": node["sequence_id"], + "transitions": edges, + } + + +@router.get("/graph/outer") +async def get_outer_graph( + seq: str = Query(..., max_length=400, description="Comma-separated sequence IDs whose ECOD groups to show"), + per_group: int = Query(14, ge=2, le=40), +): + """Outer layer of the knowledge graph around the given sequences (one clique per ECOD H-set).""" + ids = [i.strip().upper() for i in seq.split(",") if i.strip()][:12] + listed = _index_required(await run_in_threadpool(catalog.list_sequences, "sequence", "asc", None, None, None, + None, None, None, 0, len(ids), ids)) + hsets = [] + for s_ in sorted(listed["items"], key=lambda x: ids.index(x["sequence_id"])): + hs = ";".join(s_["ecod_hset"]) + if hs and hs not in hsets: + hsets.append(hs) + return await run_in_threadpool(catalog.outer_graph, hsets, per_group, ids) + + +@router.get("/graph/inner/{sequence_id}") +async def get_inner_graph( + sequence_id: str, + max_nodes: int = Query(60, ge=2, le=150), + max_states: int = Query(7, ge=1, le=40), +): + """Inner layer for one sequence: observations of its largest states and their pairwise similarities.""" + data = _index_required( + await run_in_threadpool(catalog.inner_graph_nodes, sequence_id.upper(), max_nodes, max_states) + ) + if not data: + raise HTTPException(status_code=404, detail=f"Sequence {sequence_id} not found") + chains = [(n["pdb_id"], n["auth_asym_id"]) for n in data["nodes"]] + data["edges"] = await run_in_threadpool(records.pair_similarities, chains) + return data diff --git a/backend/app/protein/catalog.py b/backend/app/protein/catalog.py new file mode 100644 index 0000000000000000000000000000000000000000..a4bdd5c0e0ed27252560f696c5dce7df3337b9ae --- /dev/null +++ b/backend/app/protein/catalog.py @@ -0,0 +1,467 @@ +""" +Queries against the website index sidecar (MuSProt-index.db). + +The sidecar is built by scripts/build_index.py and holds the sequence catalog, +per-chain membership, state sizes, state-pair similarities and an FTS index of +function text. Every function degrades to None / [] when the sidecar is absent, +so chain-level pages keep working off the main DB alone. +""" +from __future__ import annotations + +import json +import re +import sqlite3 +from functools import lru_cache +from typing import Any, Dict, List, Optional, Tuple + +from app.protein.config import get_index_database_path +from app.protein.database import connect_readonly + +MAX_HEATMAP_STATES = 40 + +SORTABLE = { + "states": "n_states", + "observations": "n_obs", + "entries": "n_entries", + "length": "length", + "released": "last_release", + "sequence": "sequence_id", +} + +CHAIN_RE = re.compile(r"^([0-9][A-Za-z0-9]{3})\s*[.:_\s/-]\s*([A-Za-z0-9]{1,4})$") +PDB_RE = re.compile(r"^[0-9][A-Za-z0-9]{3}$") +SEQ_RE = re.compile(r"^SEQ\d{1,6}$", re.IGNORECASE) +UNIPROT_RE = re.compile( + r"^([OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2})(-\d+)?$", + re.IGNORECASE, +) + + +def available() -> bool: + return get_index_database_path() is not None + + +def _connect() -> Optional[sqlite3.Connection]: + path = get_index_database_path() + if path is None: + return None + conn = connect_readonly(path) + conn.row_factory = sqlite3.Row + return conn + + +@lru_cache(maxsize=1) +def overview() -> Optional[Dict[str, Any]]: + conn = _connect() + if conn is None: + return None + try: + row = conn.execute("SELECT value FROM meta WHERE key = 'overview'").fetchone() + return json.loads(row[0]) if row else None + finally: + conn.close() + + +def _sequence_dict(r: sqlite3.Row) -> Dict[str, Any]: + return { + "sequence_id": r["sequence_id"], + "uniprot_id": r["uniprot_id"], + "length": r["length"], + "n_obs": r["n_obs"], + "n_states": r["n_states"], + "n_entries": r["n_entries"], + "n_holo": r["n_holo"], + "n_apo": r["n_apo"], + "n_cross_pairs": r["n_cross_pairs"], + "largest_state": r["largest_state"], + "rep_pdb": r["rep_pdb"], + "rep_chain": r["rep_chain"], + "cath_superfamily": r["cath_superfamily"], + "methods": [m for m in (r["methods"] or "").split(";") if m], + "first_release": r["first_release"], + "last_release": r["last_release"], + "top_function": r["top_function"], + "ecod_hset": [h for h in (r["ecod_hset"] or "").split(";") if h], + "pfam_ids": [h for h in (r["pfam_ids"] or "").split(";") if h], + } + + +def _member_dict(r: sqlite3.Row) -> Dict[str, Any]: + return { + "pdb_id": r["pdb_id"].lower(), + "auth_asym_id": r["auth_asym_id"], + "state_id": r["state_id"], + "uniprot_id": r["uniprot_id"], + "experimental_method": r["experimental_method"], + "resolution": r["resolution"], + "binding_status": r["binding_status"], + "chain_composition": r["chain_composition"], + "binders": [b for b in (r["binders"] or "").split(";") if b], + "release_date": r["release_date"], + "resolved_coverage": r["resolved_coverage"], + "cath_id": r["cath_id"], + "pfam_id": r["pfam_id"], + "ecod_id": r["ecod_id"], + "ecod_fid": r["ecod_fid"], + } + + +def ecod_names(conn: Optional[sqlite3.Connection], f_ids: List[str] = (), h_ids: List[str] = ()) -> Dict[str, Any]: + """Names for ECOD family ids (x.h.t.f) and homology ids (x.h).""" + own = conn is None + conn = conn or _connect() + if conn is None: + return {"families": {}, "homologies": {}} + try: + h_all = set(h_ids) | {".".join(f.split(".")[:2]) for f in f_ids} + fams = {} + for f in f_ids: + r = conn.execute("SELECT t_name, f_name FROM ecod_f WHERE f_id = ?", (f,)).fetchone() + if r: + fams[f] = {"t_name": r[0], "f_name": r[1]} + homs = {} + for h in h_all: + r = conn.execute("SELECT architecture, x_name, h_name FROM ecod_h WHERE h_id = ?", (h,)).fetchone() + if r: + homs[h] = {"architecture": r[0], "x_name": r[1], "h_name": r[2]} + return {"families": fams, "homologies": homs} + finally: + if own: + conn.close() + + +def group_info(conn: sqlite3.Connection, hset: List[str]) -> Optional[Dict[str, Any]]: + if not hset: + return None + r = conn.execute("SELECT * FROM ecod_group WHERE hset = ?", (";".join(hset),)).fetchone() + if r is None: + return None + return { + "hset": r["hset"].split(";"), "label": r["label"], "n_sequences": r["n_sequences"], + "n_obs": r["n_obs"], "n_multistate": r["n_multistate"], "n_states": r["n_states"], + } + + +def _fts_query(q: str) -> str: + """Turn free text into a safe FTS5 prefix query (every term must match).""" + terms = re.findall(r"[A-Za-z0-9]+", q) + return " ".join(f'"{t}"*' for t in terms[:8]) + + +# ─────────────────────────────────────────────────────────── browse +def list_sequences( + sort: str = "states", + order: str = "desc", + q: Optional[str] = None, + min_states: Optional[int] = None, + min_obs: Optional[int] = None, + method: Optional[str] = None, + ligand: Optional[str] = None, + cath: Optional[str] = None, + offset: int = 0, + limit: int = 25, + ids: Optional[List[str]] = None, + ecod: Optional[str] = None, + pfam: Optional[str] = None, +) -> Optional[Dict[str, Any]]: + conn = _connect() + if conn is None: + return None + where: List[str] = [] + params: List[Any] = [] + if ids: + where.append(f"s.sequence_id IN ({','.join('?' * len(ids))})") + params.extend(ids) + if q and _fts_query(q): + where.append("s.sequence_id IN (SELECT sequence_id FROM function_fts WHERE function_fts MATCH ?)") + params.append(_fts_query(q)) + if min_states: + where.append("s.n_states >= ?") + params.append(min_states) + if min_obs: + where.append("s.n_obs >= ?") + params.append(min_obs) + if method: + where.append("(';' || s.methods || ';') LIKE ?") + params.append(f"%;{method};%") + if ligand == "both": + where.append("s.n_holo > 0 AND s.n_apo > 0") + elif ligand == "holo": + where.append("s.n_holo > 0") + elif ligand == "apo": + where.append("s.n_apo > 0") + if cath: + where.append("s.cath_superfamily LIKE ?") + params.append(f"{cath}%") + if ecod: + # an exact H-set (outer-graph group, e.g. "11.1;2003.1") or one H-group contained in it + if ";" in ecod: + where.append("s.ecod_hset = ?") + params.append(ecod) + else: + where.append("(';' || s.ecod_hset || ';') LIKE ?") + params.append(f"%;{ecod};%") + if pfam: + where.append("(';' || s.pfam_ids || ';') LIKE ?") + params.append(f"%;{pfam.upper()};%") + where_sql = f"WHERE {' AND '.join(where)}" if where else "" + col = SORTABLE.get(sort, "n_states") + direction = "ASC" if order == "asc" else "DESC" + try: + total = conn.execute(f"SELECT COUNT(*) FROM sequence s {where_sql}", params).fetchone()[0] + rows = conn.execute( + f"SELECT s.* FROM sequence s {where_sql} " + f"ORDER BY s.{col} {direction} NULLS LAST, s.n_obs DESC, s.sequence_id " + f"LIMIT ? OFFSET ?", + [*params, max(1, min(limit, 100)), max(0, offset)], + ).fetchall() + finally: + conn.close() + return {"total": total, "offset": offset, "items": [_sequence_dict(r) for r in rows]} + + +# ─────────────────────────────────────────────────────────── sequence detail +def get_sequence(sequence_id: str) -> Optional[Dict[str, Any]]: + conn = _connect() + if conn is None: + return None + try: + row = conn.execute( + "SELECT * FROM sequence WHERE sequence_id = ? COLLATE NOCASE", (sequence_id,) + ).fetchone() + if row is None: + return {} + seq = _sequence_dict(row) + sid = seq["sequence_id"] + members = [ + _member_dict(r) for r in conn.execute( + "SELECT * FROM member WHERE sequence_id = ? ORDER BY state_id, resolution IS NULL, resolution", + (sid,), + ) + ] + states = [ + {"state_id": r["state_id"], "n_members": r["n_members"]} + for r in conn.execute( + "SELECT state_id, n_members FROM state WHERE sequence_id = ? " + "ORDER BY n_members DESC, state_id", + (sid,), + ) + ] + shown = {s["state_id"] for s in states[:MAX_HEATMAP_STATES]} + pairs = [ + { + "a": r["state_a"], "b": r["state_b"], "similarity": r["similarity"], + "fidelity": r["fidelity"], "n_pairs": r["n_pairs"], + } + for r in conn.execute( + "SELECT state_a, state_b, similarity, fidelity, n_pairs FROM state_pair " + "WHERE sequence_id = ?", + (sid,), + ) + if r["state_a"] in shown and r["state_b"] in shown + ] + seq["ecod_group"] = group_info(conn, seq["ecod_hset"]) + seq["ecod_names"] = ecod_names(conn, h_ids=seq["ecod_hset"]) + finally: + conn.close() + seq.update(states=states, members=members, state_pairs=pairs, + heatmap_truncated=len(states) > MAX_HEATMAP_STATES) + return seq + + +def sequence_members(sequence_id: str) -> Dict[Tuple[str, str], Dict[str, Any]]: + """Map (pdb_id, chain) -> member metadata for one sequence.""" + conn = _connect() + if conn is None: + return {} + try: + return { + (r["pdb_id"].lower(), r["auth_asym_id"]): _member_dict(r) + for r in conn.execute("SELECT * FROM member WHERE sequence_id = ?", (sequence_id,)) + } + finally: + conn.close() + + +def sequence_summary(sequence_id: str) -> Optional[Dict[str, Any]]: + conn = _connect() + if conn is None: + return None + try: + row = conn.execute("SELECT * FROM sequence WHERE sequence_id = ?", (sequence_id,)).fetchone() + states = [ + {"state_id": r["state_id"], "n_members": r["n_members"]} + for r in conn.execute( + "SELECT state_id, n_members FROM state WHERE sequence_id = ? ORDER BY n_members DESC, state_id", + (sequence_id,), + ) + ] + finally: + conn.close() + if row is None: + return None + out = _sequence_dict(row) + out["states"] = states + return out + + +# ─────────────────────────────────────────────────────────── search +def search(q: str, limit: int = 12) -> Dict[str, Any]: + """Classify a free-text query and return typed hits.""" + q = q.strip() + if not q: + return {"query": q, "kind": "empty", "chains": [], "sequences": []} + + m = CHAIN_RE.match(q) + if m: + from app.protein import records + pdb, chain = m.group(1).lower(), m.group(2) + node_chains = [c for c in records.list_entry_chains(pdb)] + exact = [c for c in node_chains if c["auth_asym_id"] == chain] or [ + c for c in node_chains if c["auth_asym_id"].lower() == chain.lower() + ] + return {"query": q, "kind": "chain", "chains": exact, "sequences": [], "entry_chains": node_chains} + + if PDB_RE.match(q) and not SEQ_RE.match(q): + from app.protein import records + chains = records.list_entry_chains(q.lower()) + if chains: + return {"query": q, "kind": "entry", "chains": chains, "sequences": []} + + conn = _connect() + if conn is None: + return {"query": q, "kind": "unavailable", "chains": [], "sequences": []} + try: + if SEQ_RE.match(q): + digits = re.sub(r"\D", "", q) + rows = conn.execute( + "SELECT * FROM sequence WHERE sequence_id = ?", (f"SEQ{int(digits):06d}",) + ).fetchall() + return {"query": q, "kind": "sequence", "chains": [], "sequences": [_sequence_dict(r) for r in rows]} + + if UNIPROT_RE.match(q): + rows = conn.execute( + "SELECT * FROM sequence WHERE sequence_id IN " + "(SELECT DISTINCT sequence_id FROM member WHERE uniprot_id = ? COLLATE NOCASE) " + "ORDER BY n_obs DESC LIMIT ?", + (q, limit), + ).fetchall() + if not rows: + # multi-accession sequences store e.g. "O43318;Q15750"; the FTS column tokenises them + rows = conn.execute( + "SELECT s.* FROM function_fts f JOIN sequence s ON s.sequence_id = f.sequence_id " + "WHERE function_fts MATCH ? ORDER BY s.n_obs DESC LIMIT ?", + (f'uniprot_id : "{q}"', limit), + ).fetchall() + if rows: + return {"query": q, "kind": "uniprot", "chains": [], "sequences": [_sequence_dict(r) for r in rows]} + + fts = _fts_query(q) + if not fts: + return {"query": q, "kind": "text", "chains": [], "sequences": []} + rows = conn.execute( + "SELECT s.* FROM function_fts f JOIN sequence s ON s.sequence_id = f.sequence_id " + "WHERE function_fts MATCH ? ORDER BY bm25(function_fts), s.n_obs DESC LIMIT ?", + (fts, limit), + ).fetchall() + total = conn.execute( + "SELECT COUNT(*) FROM function_fts WHERE function_fts MATCH ?", (fts,) + ).fetchone()[0] + return { + "query": q, "kind": "text", "chains": [], "total": total, + "sequences": [_sequence_dict(r) for r in rows], + } + finally: + conn.close() + + +# ─────────────────────────────────────────────────────────── knowledge graph +def outer_graph(hsets: List[str], per_group: int = 14, include: List[str] = ()) -> Optional[Dict[str, Any]]: + """Outer layer: sequences grouped by identical ECOD H-set (each group is a clique).""" + conn = _connect() + if conn is None: + return None + try: + groups = [] + for hs in hsets: + info = group_info(conn, hs.split(";")) + if info is None: + continue + must = [ + _sequence_dict(r) for r in conn.execute( + f"SELECT * FROM sequence WHERE ecod_hset = ? AND sequence_id IN ({','.join('?' * len(include))})", + (hs, *include), + ) + ] if include else [] + rest = [ + _sequence_dict(r) for r in conn.execute( + "SELECT * FROM sequence WHERE ecod_hset = ? ORDER BY n_states DESC, n_obs DESC LIMIT ?", + (hs, per_group), + ) + ] + seen = {m["sequence_id"] for m in must} + members = (must + [r for r in rest if r["sequence_id"] not in seen])[:per_group] + groups.append({**info, "sequences": members}) + names = ecod_names(conn, h_ids=[h for g in groups for h in g["hset"]]) + finally: + conn.close() + return {"groups": groups, "names": names["homologies"]} + + +def inner_graph_nodes(sequence_id: str, max_nodes: int = 60, max_states: int = 7) -> Optional[Dict[str, Any]]: + """Pick up to max_nodes observations from the sequence's max_states largest states.""" + conn = _connect() + if conn is None: + return None + try: + row = conn.execute("SELECT * FROM sequence WHERE sequence_id = ?", (sequence_id,)).fetchone() + if row is None: + return {} + members = [_member_dict(r) for r in conn.execute( + "SELECT * FROM member WHERE sequence_id = ? ORDER BY state_id, resolution IS NULL, resolution", + (sequence_id,), + )] + states = [ + {"state_id": r[0], "n_members": r[1]} for r in conn.execute( + "SELECT state_id, n_members FROM state WHERE sequence_id = ? ORDER BY n_members DESC, state_id", + (sequence_id,), + ) + ] + pairs = [ + {"a": r[0], "b": r[1], "similarity": r[2], "fidelity": r[3], "n_pairs": r[4]} + for r in conn.execute( + "SELECT state_a, state_b, similarity, fidelity, n_pairs FROM state_pair WHERE sequence_id = ?", + (sequence_id,), + ) + ] + finally: + conn.close() + + order = [s["state_id"] for s in states[:max_states]] + pool = [m for m in members if m["state_id"] in set(order)] + if len(pool) <= max_nodes: + chosen = pool + else: + # round-robin across the kept states (largest first) so smaller states still appear + by_state: Dict[int, List[Dict[str, Any]]] = {} + for m in pool: + by_state.setdefault(m["state_id"], []).append(m) + chosen, depth = [], 0 + while len(chosen) < max_nodes: + added = False + for st in order: + bucket = by_state.get(st, []) + if depth < len(bucket) and len(chosen) < max_nodes: + chosen.append(bucket[depth]) + added = True + if not added: + break + depth += 1 + shown_states = {s["state_id"] for s in states[:MAX_HEATMAP_STATES]} + return { + "sequence": _sequence_dict(row), + "nodes": chosen, + "n_total": len(members), + "states": states, + "state_pairs": [p for p in pairs if p["a"] in shown_states and p["b"] in shown_states], + } diff --git a/backend/app/protein/config.py b/backend/app/protein/config.py index faaff5b673842f352a223589bb875279df64ec4c..0d7752803b28636d0c206b538a6e29990f3106be 100644 --- a/backend/app/protein/config.py +++ b/backend/app/protein/config.py @@ -66,18 +66,18 @@ def get_summary_path() -> Path | None: ) -def get_node_database_path() -> Path: - """Return an optional indexed node lookup sidecar, falling back to the main DB.""" - explicit_path = os.getenv("MUSPROT_NODE_DB_PATH") +def get_index_database_path() -> Path | None: + """Return the website index sidecar (built by scripts/build_index.py), if present.""" + explicit_path = os.getenv("MUSPROT_INDEX_DB_PATH") if explicit_path: path = Path(explicit_path).expanduser().resolve() - if path.exists(): - return path + return path if path.exists() else None - return ( - _existing_path(DEFAULT_VOLUME_DIR / "MuSProt-node-lookup.db") - or get_database_path() - ).resolve() + return _existing_path( + DEFAULT_VOLUME_DIR / "MuSProt-index.db", + ASSETS_DIR / "MuSProt-index.db", + LOCAL_DATA_DIR / "MuSProt-index.db", + ) def get_docs_path() -> Path | None: diff --git a/backend/app/protein/data_loader.py b/backend/app/protein/data_loader.py index 11bbcf3601e56709b38652ef87daf2afbdf0a6f2..8554fc4ac49dc7ac42f7a8c526176822c6ab4827 100755 --- a/backend/app/protein/data_loader.py +++ b/backend/app/protein/data_loader.py @@ -44,6 +44,13 @@ class DataManager: if summary_path: self._summary = json.loads(summary_path.read_text(encoding="utf-8")) self._total_records = int(self._summary.get("total_records", 0)) + from app.protein import catalog + overview = catalog.overview() + if overview: + # the index is rebuilt with every DB release, so its count wins over a stale summary + self._total_records = int(overview["counts"]["transitions"]) + if self._summary: + self._summary["total_records"] = self._total_records conn = connect_readonly(self.db_path) try: @@ -100,7 +107,7 @@ class DataManager: try: cur = conn.cursor() cur.execute( - "SELECT MIN(CAST(TM1 AS REAL)), MAX(CAST(TM1 AS REAL))," + "SELECT MIN(CAST(tm_aln AS REAL)), MAX(CAST(tm_aln AS REAL))," " MIN(CAST(RMSD AS REAL)), MAX(CAST(RMSD AS REAL))" " FROM edge" ) @@ -161,16 +168,16 @@ class DataManager: where_clauses.append("CAST(RMSD AS REAL) <= ?") params.append(rmsd_max) if tm_min is not None: - where_clauses.append("CAST(TM1 AS REAL) >= ?") + where_clauses.append("CAST(tm_aln AS REAL) >= ?") params.append(tm_min) if tm_max is not None: - where_clauses.append("CAST(TM1 AS REAL) <= ?") + where_clauses.append("CAST(tm_aln AS REAL) <= ?") params.append(tm_max) where_sql = ("WHERE " + " AND ".join(where_clauses)) if where_clauses else "" sql = f""" SELECT e.pdb_id_A, e.auth_asym_id_A, e.pdb_id_B, e.auth_asym_id_B, - CAST(e.TM1 AS REAL) AS TM1, + CAST(e.tm_aln AS REAL) AS TM1, CAST(e.RMSD AS REAL) AS RMSD, CAST(e.pair_similarity AS REAL) AS structure_sim, e."delta_Rosetta", e."delta_FoldX", e."delta_EvoEF2", @@ -249,10 +256,10 @@ class DataManager: where_clauses.append("CAST(RMSD AS REAL) <= ?") params.append(rmsd_max) if tm_min is not None: - where_clauses.append("CAST(TM1 AS REAL) >= ?") + where_clauses.append("CAST(tm_aln AS REAL) >= ?") params.append(tm_min) if tm_max is not None: - where_clauses.append("CAST(TM1 AS REAL) <= ?") + where_clauses.append("CAST(tm_aln AS REAL) <= ?") params.append(tm_max) where_sql = ("WHERE " + " AND ".join(where_clauses)) if where_clauses else "" @@ -261,7 +268,7 @@ class DataManager: try: cur = conn.cursor() cur.execute( - f"SELECT COUNT(*), AVG(CAST(RMSD AS REAL)), AVG(CAST(TM1 AS REAL))" + f"SELECT COUNT(*), AVG(CAST(RMSD AS REAL)), AVG(CAST(tm_aln AS REAL))" f" FROM edge {where_sql}", params, ) @@ -292,7 +299,7 @@ class DataManager: ["0-0.5", "0.5-1.0", "1.0-1.5", "1.5-2.0", "2.0-5.0", ">5.0"], ) tm_dist = _dist( - "TM1", + "tm_aln", [0, 0.5, 0.7, 0.85, 0.95, 1.01], ["0-0.5", "0.5-0.7", "0.7-0.85", "0.85-0.95", "0.95-1.0"], ) diff --git a/backend/app/protein/models.py b/backend/app/protein/models.py index de22c470bbb415f0603fdc71f26f1d2c8cf94725..2c6448af83062371f4fc9316fd33cd0ef99a0584 100755 --- a/backend/app/protein/models.py +++ b/backend/app/protein/models.py @@ -1,144 +1,144 @@ -""" -Data models for the protein visualization API. -""" -from typing import Optional, List, Dict, Any -from pydantic import BaseModel, Field - - -class CrystalGrow(BaseModel): - """Crystal growth conditions.""" - pH: Optional[float] = None - pH_structured: Optional[float] = None - pH_from_details: Optional[float] = None - temp: Optional[float] = None - pdbx_details: Optional[str] = None - - -class EntryInfo(BaseModel): - """Entry-level information.""" - rcsb_id: Optional[str] = None - title: Optional[str] = None - initial_release_date: Optional[str] = None - exptl_method: Optional[str] = None - resolution_combined: Optional[List[float]] = None - crystal_grow: Optional[CrystalGrow] = None - - -class PolymerEntity(BaseModel): - """Polymer entity (sequence-level) information.""" - rcsb_id: Optional[str] = None - entity_id: Optional[str] = None - sequence_length: Optional[int] = None - seq_can: Optional[str] = None - reference_sequences: Optional[List[Dict[str, str]]] = None - - -class PolymerEntityInstance(BaseModel): - """Polymer entity instance (chain-level) information.""" - rcsb_id: Optional[str] = None - asym_id: Optional[str] = None - auth_asym_id: Optional[str] = None - binding_affinity: Optional[Any] = None - - -class NonpolymerEntity(BaseModel): - """Nonpolymer entity (ligand/binder) information.""" - rcsb_id: Optional[str] = None - name: Optional[str] = None - formula_weight: Optional[float] = None - smiles: Optional[str] = None - inchi: Optional[str] = None - - -class ChainMetadata(BaseModel): - """Complete chain metadata from RCSB.""" - pdb_id: str - chain_id: str - entry: EntryInfo - polymer_entity: Optional[PolymerEntity] = None - polymer_entity_instance: Optional[PolymerEntityInstance] = None - nonpolymer_entities: Optional[List[NonpolymerEntity]] = None - binding_status: Optional[str] = None - chain_composition: Optional[str] = None - cath_id: Optional[str] = None - cath_superfamily: Optional[str] = None - state_id: Optional[str] = None - - -class FilterParams(BaseModel): - """Request model for filtering protein data.""" - chain_ids: Optional[List[str]] = Field(default=None, description="List of chain IDs to filter") - pdb_id: Optional[str] = Field(default=None, description="PDB ID for single chain search") - auth_asym_id: Optional[str] = Field(default=None, description="Chain ID for single chain search") - rmsd_min: Optional[float] = Field(default=None, ge=0, description="Minimum RMSD value") - rmsd_max: Optional[float] = Field(default=None, ge=0, description="Maximum RMSD value") - tm_score_min: Optional[float] = Field(default=None, ge=0, le=1, description="Minimum TM-score") - tm_score_max: Optional[float] = Field(default=None, ge=0, le=1, description="Maximum TM-score") - length_min: Optional[int] = Field(default=None, ge=0, description="Minimum sequence length") - length_max: Optional[int] = Field(default=None, ge=0, description="Maximum sequence length") - include_duplicates: Optional[bool] = Field(default=None, description="Include duplicate structures") - cluster_id: Optional[str] = Field(default=None, description="Structural similarity cluster ID") - limit: Optional[int] = Field(default=1000, ge=1, le=10000, description="Maximum results to return") - - -class FilterOptions(BaseModel): - """Available filter options and ranges.""" - chain_ids: List[str] = Field(description="Available chain IDs") - rmsd_range: Dict[str, float] = Field(description="Min and max RMSD values") - tm_score_range: Dict[str, float] = Field(description="Min and max TM-score values") - length_range: Dict[str, int] = Field(description="Min and max sequence lengths") - clusters: List[str] = Field(description="Available similarity clusters") - total_records: int = Field(description="Total number of records") - - -class DataRecord(BaseModel): - """Single protein comparison record.""" - pdb_id_a: str - auth_asym_id_a: str - pdb_id_b: str - auth_asym_id_b: str - tm_score: float - rmsd: float - structure_sim: Optional[float] = None - sequence_identity: Optional[float] = None - length_a: Optional[int] = None - length_b: Optional[int] = None - aligned_length: Optional[int] = None - status: Optional[str] = None - cluster_id: Optional[str] = None - exptl_method: Optional[str] = None - temp: Optional[float] = None - pH: Optional[float] = None - # Matched-chain (B) annotations from the node table - binding_status: Optional[str] = None - chain_composition: Optional[str] = None - # Energy score deltas fetched directly from the edge table - delta_rosetta: Optional[float] = None - delta_foldx: Optional[float] = None - delta_evoef2: Optional[float] = None - delta_rw: Optional[float] = None - delta_rw_plus: Optional[float] = None - # Alternative state similarity fields from edge table - state_id_b: Optional[str] = None - state_fidelity: Optional[str] = None - avg_sim: Optional[str] = None - observation_fidelity: Optional[str] = None - - -class DataResponse(BaseModel): - """Response model for filtered data.""" - data: List[DataRecord] - total: int - filtered: int - - -class SummaryStats(BaseModel): - """Aggregate statistics for the dataset.""" - total_records: int - unique_chains: int - avg_rmsd: float - avg_tm_score: float - avg_sequence_length: float - rmsd_distribution: Dict[str, int] - tm_score_distribution: Dict[str, int] - length_distribution: Dict[str, int] +""" +Data models for the protein visualization API. +""" +from typing import Optional, List, Dict, Any +from pydantic import BaseModel, Field + + +class CrystalGrow(BaseModel): + """Crystal growth conditions.""" + pH: Optional[float] = None + pH_structured: Optional[float] = None + pH_from_details: Optional[float] = None + temp: Optional[float] = None + pdbx_details: Optional[str] = None + + +class EntryInfo(BaseModel): + """Entry-level information.""" + rcsb_id: Optional[str] = None + title: Optional[str] = None + initial_release_date: Optional[str] = None + exptl_method: Optional[str] = None + resolution_combined: Optional[List[float]] = None + crystal_grow: Optional[CrystalGrow] = None + + +class PolymerEntity(BaseModel): + """Polymer entity (sequence-level) information.""" + rcsb_id: Optional[str] = None + entity_id: Optional[str] = None + sequence_length: Optional[int] = None + seq_can: Optional[str] = None + reference_sequences: Optional[List[Dict[str, str]]] = None + + +class PolymerEntityInstance(BaseModel): + """Polymer entity instance (chain-level) information.""" + rcsb_id: Optional[str] = None + asym_id: Optional[str] = None + auth_asym_id: Optional[str] = None + binding_affinity: Optional[Any] = None + + +class NonpolymerEntity(BaseModel): + """Nonpolymer entity (ligand/binder) information.""" + rcsb_id: Optional[str] = None + name: Optional[str] = None + formula_weight: Optional[float] = None + smiles: Optional[str] = None + inchi: Optional[str] = None + + +class ChainMetadata(BaseModel): + """Complete chain metadata from RCSB.""" + pdb_id: str + chain_id: str + entry: EntryInfo + polymer_entity: Optional[PolymerEntity] = None + polymer_entity_instance: Optional[PolymerEntityInstance] = None + nonpolymer_entities: Optional[List[NonpolymerEntity]] = None + binding_status: Optional[str] = None + chain_composition: Optional[str] = None + cath_id: Optional[str] = None + cath_superfamily: Optional[str] = None + state_id: Optional[str] = None + + +class FilterParams(BaseModel): + """Request model for filtering protein data.""" + chain_ids: Optional[List[str]] = Field(default=None, description="List of chain IDs to filter") + pdb_id: Optional[str] = Field(default=None, description="PDB ID for single chain search") + auth_asym_id: Optional[str] = Field(default=None, description="Chain ID for single chain search") + rmsd_min: Optional[float] = Field(default=None, ge=0, description="Minimum RMSD value") + rmsd_max: Optional[float] = Field(default=None, ge=0, description="Maximum RMSD value") + tm_score_min: Optional[float] = Field(default=None, ge=0, le=1, description="Minimum TM-score") + tm_score_max: Optional[float] = Field(default=None, ge=0, le=1, description="Maximum TM-score") + length_min: Optional[int] = Field(default=None, ge=0, description="Minimum sequence length") + length_max: Optional[int] = Field(default=None, ge=0, description="Maximum sequence length") + include_duplicates: Optional[bool] = Field(default=None, description="Include duplicate structures") + cluster_id: Optional[str] = Field(default=None, description="Structural similarity cluster ID") + limit: Optional[int] = Field(default=1000, ge=1, le=10000, description="Maximum results to return") + + +class FilterOptions(BaseModel): + """Available filter options and ranges.""" + chain_ids: List[str] = Field(description="Available chain IDs") + rmsd_range: Dict[str, float] = Field(description="Min and max RMSD values") + tm_score_range: Dict[str, float] = Field(description="Min and max TM-score values") + length_range: Dict[str, int] = Field(description="Min and max sequence lengths") + clusters: List[str] = Field(description="Available similarity clusters") + total_records: int = Field(description="Total number of records") + + +class DataRecord(BaseModel): + """Single protein comparison record.""" + pdb_id_a: str + auth_asym_id_a: str + pdb_id_b: str + auth_asym_id_b: str + tm_score: float + rmsd: float + structure_sim: Optional[float] = None + sequence_identity: Optional[float] = None + length_a: Optional[int] = None + length_b: Optional[int] = None + aligned_length: Optional[int] = None + status: Optional[str] = None + cluster_id: Optional[str] = None + exptl_method: Optional[str] = None + temp: Optional[float] = None + pH: Optional[float] = None + # Matched-chain (B) annotations from the node table + binding_status: Optional[str] = None + chain_composition: Optional[str] = None + # Energy score deltas fetched directly from the edge table + delta_rosetta: Optional[float] = None + delta_foldx: Optional[float] = None + delta_evoef2: Optional[float] = None + delta_rw: Optional[float] = None + delta_rw_plus: Optional[float] = None + # Alternative state similarity fields from edge table + state_id_b: Optional[str] = None + state_fidelity: Optional[str] = None + avg_sim: Optional[str] = None + observation_fidelity: Optional[str] = None + + +class DataResponse(BaseModel): + """Response model for filtered data.""" + data: List[DataRecord] + total: int + filtered: int + + +class SummaryStats(BaseModel): + """Aggregate statistics for the dataset.""" + total_records: int + unique_chains: int + avg_rmsd: float + avg_tm_score: float + avg_sequence_length: float + rmsd_distribution: Dict[str, int] + tm_score_distribution: Dict[str, int] + length_distribution: Dict[str, int] diff --git a/backend/app/protein/records.py b/backend/app/protein/records.py new file mode 100644 index 0000000000000000000000000000000000000000..58f49bffc07ff68f0c5fd842740ec80c7afadcc3 --- /dev/null +++ b/backend/app/protein/records.py @@ -0,0 +1,270 @@ +""" +Per-chain reads from the main MuSProt.db. + +Both tables carry a (pdb_id COLLATE NOCASE, chain) index, so every query here +is a single index seek. pdb_id is case-insensitive; chain ids are case-sensitive +(`d` and `D` are different chains), with a case-insensitive fallback for +convenience when the exact chain does not exist. +""" +from __future__ import annotations + +import json +import math +import sqlite3 +from typing import Any, Dict, List, Optional + +from app.protein.config import get_database_path +from app.protein.database import connect_readonly + +ENERGY_TERMS = ["Rosetta", "FoldX", "EvoEF2", "RW", "RW+"] + +EDGE_COLUMNS = [ + "pdb_id_B", "auth_asym_id_B", "tm_aln", "RMSD", "coverage", "pair_similarity", + "delta_Rosetta", "delta_FoldX", "delta_EvoEF2", "delta_RW", "delta_RW+", + "state_id_B", "state_similarity", "state_fidelity", "pair_fidelity", + "n_shared_resolved", "delta_order_valid", "order_invalid_reason", "n_ordered", "n_disordered", + "ordered_ranges", "disordered_ranges", "longest_ordered_segment", "longest_disordered_segment", + "n_terminal_ordered", "n_internal_ordered", "order_evidence", "ordering_with_ligand_change", +] + + +def _connect() -> sqlite3.Connection: + conn = connect_readonly(get_database_path()) + conn.row_factory = sqlite3.Row + return conn + + +def num(value: Any) -> Optional[float]: + if value is None or value == "": + return None + try: + f = float(value) + except (TypeError, ValueError): + return None + return None if math.isnan(f) or math.isinf(f) else f + + +def integer(value: Any) -> Optional[int]: + f = num(value) + return None if f is None else int(f) + + +def text(value: Any) -> Optional[str]: + if value is None: + return None + s = str(value).strip() + return None if s in ("", "None", "nan", "") else s + + +def _parse_list(raw: Any) -> List[Any]: + raw = (raw or "").strip() if isinstance(raw, str) else raw + if not raw: + return [] + try: + val = json.loads(raw) + except ValueError: + import ast + try: + val = ast.literal_eval(raw) + except (ValueError, SyntaxError): + return [s.strip() for s in str(raw).split(";") if s.strip()] + return val if isinstance(val, list) else [] + + +def _parse_ranges(raw: Optional[str]) -> List[List[int]]: + """'5-173,187-353,400' -> [[5, 173], [187, 353], [400, 400]] (1-based, inclusive).""" + out = [] + for part in (raw or "").split(","): + bounds = part.strip().split("-") + if len(bounds) == 2 and bounds[0].isdigit() and bounds[1].isdigit(): + out.append([int(bounds[0]), int(bounds[1])]) + elif len(bounds) == 1 and bounds[0].isdigit(): + out.append([int(bounds[0]), int(bounds[0])]) + return out + + +def _split(raw: Any) -> List[str]: + return [x for x in (text(raw) or "").split(";") if x] + + +def _flag(raw: Any) -> Optional[bool]: + t = text(raw) + return None if t is None else t == "True" + + +def resolve_chain_key(conn: sqlite3.Connection, pdb_id: str, chain: str) -> Optional[str]: + """Return the stored chain id for (pdb_id, chain), preferring an exact match.""" + row = conn.execute( + "SELECT auth_asym_id FROM node WHERE pdb_id = ? COLLATE NOCASE AND auth_asym_id = ? LIMIT 1", + (pdb_id, chain), + ).fetchone() + if row: + return row[0] + candidates = [ + r[0] for r in conn.execute( + "SELECT auth_asym_id FROM node WHERE pdb_id = ? COLLATE NOCASE", (pdb_id,) + ) + ] + for c in candidates: + if c.lower() == chain.lower(): + return c + return None + + +def get_node(pdb_id: str, chain: str) -> Optional[Dict[str, Any]]: + """Full, parsed node record for one chain, or None when absent.""" + conn = _connect() + try: + stored = resolve_chain_key(conn, pdb_id, chain) + if stored is None: + return None + r = conn.execute( + "SELECT * FROM node WHERE pdb_id = ? COLLATE NOCASE AND auth_asym_id = ? LIMIT 1", + (pdb_id, stored), + ).fetchone() + finally: + conn.close() + + functions = [str(f).strip() for f in _parse_list(r["ranked_functions"])] + scores = _parse_list(r["matching_scores"]) + ranked = [] + for i, name in enumerate(functions): + if len(name) <= 10: + continue + ranked.append({"function": name, "score": num(scores[i]) if i < len(scores) else None}) + + return { + "pdb_id": r["pdb_id"].lower(), + "auth_asym_id": r["auth_asym_id"], + "sequence_id": r["sequence_id"], + "uniprot_id": text(r["uniprot_id"]), + "state_id": integer(r["state_id"]), + "sequence": text(r["sequence"]) or "", + "sequence_length": integer(r["sequence_length"]), + "modeled_sequence_length": integer(r["modeled_sequence_length"]), + "resolved_ranges": _parse_ranges(r["resolved_label_seq"]), + "n_resolved_aa": integer(r["n_resolved_aa"]), + "resolved_coverage": num(r["resolved_coverage"]), + "binders": [b for b in (text(r["binders"]) or "").split(";") if b], + "binding_status": text(r["binding_status"]), + "experimental_method": text(r["experimental_method"]), + "resolution": num(r["resolution"]), + "pH": num(r["pH"]), + "temp_K": num(r["temp_K"]), + "experimental_details": text(r["experimental_details"]), + "chain_composition": text(r["chain_composition"]), + "nucleic_acid_binding": text(r["non_protein_polymer_binding"]), + "release_date": text(r["initial_release_date"]), + "cath_domains": _split(r["cath_id"]), + "cath_superfamilies": _split(r["cath_superfamily"]), + "pfam_ids": _split(r["pfam_id"]), + "ecod_domains": _split(r["ecod_id"]), + "ecod_fids": _split(r["ecod_fid"]), + "energies": {term: num(r[term]) for term in ENERGY_TERMS}, + "scoring": { + "n_seqres": integer(r["n_seqres"]), + "n_scored": integer(r["n_scored"]), + "scored_coverage": num(r["scored_coverage"]), + "foldx_n_scored": integer(r["foldx_n_scored"]), + "foldx_coverage": num(r["foldx_coverage"]), + }, + "functions": ranked, + } + + +def get_transitions(pdb_id: str, chain: str, limit: int = 20000) -> List[Dict[str, Any]]: + """All directed edges leaving one chain (A = the given chain).""" + cols = ", ".join(f'"{c}"' for c in EDGE_COLUMNS) + conn = _connect() + try: + stored = resolve_chain_key(conn, pdb_id, chain) or chain + rows = conn.execute( + f"SELECT {cols} FROM edge WHERE pdb_id_A = ? COLLATE NOCASE AND auth_asym_id_A = ? LIMIT ?", + (pdb_id, stored, limit), + ).fetchall() + finally: + conn.close() + + out = [] + for r in rows: + sim = num(r["pair_similarity"]) + out.append({ + "pdb_id": r["pdb_id_B"].lower(), + "auth_asym_id": r["auth_asym_id_B"], + "state_id": integer(r["state_id_B"]), + "tm_aln": num(r["tm_aln"]), + # no shared residues: RMSD/coverage are stored as 0.0 but are really undefined + "rmsd": num(r["RMSD"]) if sim is not None else None, + "coverage": num(r["coverage"]) if sim is not None else None, + "similarity": sim, + "fidelity": text(r["pair_fidelity"]), + "state_similarity": num(r["state_similarity"]), + "state_fidelity": text(r["state_fidelity"]), + "delta": { + "Rosetta": num(r["delta_Rosetta"]), + "FoldX": num(r["delta_FoldX"]), + "EvoEF2": num(r["delta_EvoEF2"]), + "RW": num(r["delta_RW"]), + "RW+": num(r["delta_RW+"]), + }, + "order": { + "valid": _flag(r["delta_order_valid"]), + "reason": text(r["order_invalid_reason"]), + "n_shared_resolved": integer(r["n_shared_resolved"]), + "n_ordered": integer(r["n_ordered"]), + "n_disordered": integer(r["n_disordered"]), + "ordered_ranges": _parse_ranges(r["ordered_ranges"]), + "disordered_ranges": _parse_ranges(r["disordered_ranges"]), + "longest_ordered": integer(r["longest_ordered_segment"]), + "longest_disordered": integer(r["longest_disordered_segment"]), + "n_terminal_ordered": integer(r["n_terminal_ordered"]), + "n_internal_ordered": integer(r["n_internal_ordered"]), + "evidence": text(r["order_evidence"]), + "with_ligand_change": _flag(r["ordering_with_ligand_change"]), + }, + }) + return out + + +def pair_similarities(chains: List[tuple]) -> List[Dict[str, Any]]: + """Undirected similarity edges among a set of (pdb_id, chain) observations of one sequence.""" + keep = {(p.lower(), c): i for i, (p, c) in enumerate(chains)} + out = [] + conn = _connect() + try: + for i, (pdb, chain) in enumerate(chains): + for r in conn.execute( + "SELECT pdb_id_B, auth_asym_id_B, pair_similarity, pair_fidelity FROM edge " + "WHERE pdb_id_A = ? COLLATE NOCASE AND auth_asym_id_A = ?", + (pdb, chain), + ): + j = keep.get((r[0].lower(), r[1])) + if j is not None and j > i: + out.append({"s": i, "t": j, "sim": num(r[2]), "fid": text(r[3])}) + finally: + conn.close() + return out + + +def list_entry_chains(pdb_id: str) -> List[Dict[str, Any]]: + """Chains of one PDB entry present in MuSProt (lightweight columns only).""" + conn = _connect() + try: + rows = conn.execute( + "SELECT pdb_id, auth_asym_id, sequence_id, uniprot_id, state_id, sequence_length " + "FROM node WHERE pdb_id = ? COLLATE NOCASE", + (pdb_id,), + ).fetchall() + finally: + conn.close() + return [ + { + "pdb_id": r["pdb_id"].lower(), + "auth_asym_id": r["auth_asym_id"], + "sequence_id": r["sequence_id"], + "uniprot_id": text(r["uniprot_id"]), + "state_id": integer(r["state_id"]), + "sequence_length": integer(r["sequence_length"]), + } + for r in rows + ] diff --git a/backend/app/protein/tsv_loader.py b/backend/app/protein/tsv_loader.py index cb763b74dbd5e8bd6c82bbf278790abedbc76a05..1e2a74ff1ada0495d9e531ebb530364a219b432c 100755 --- a/backend/app/protein/tsv_loader.py +++ b/backend/app/protein/tsv_loader.py @@ -8,11 +8,10 @@ key the first occurrence is kept. import ast import sqlite3 import logging -import zlib from functools import lru_cache from typing import Dict, List, Optional, Tuple, Any -from app.protein.config import get_database_path, get_node_database_path +from app.protein.config import get_database_path from app.protein.database import connect_readonly logger = logging.getLogger(__name__) @@ -83,23 +82,25 @@ def reset_index() -> None: @lru_cache(maxsize=4096) def _lookup_node_row(pdb_id: str, auth_asym_id: str) -> Optional[_Row]: - """Read one node row without blocking on a full-table in-memory index.""" - conn = connect_readonly(get_node_database_path()) + """Read one node row via the (pdb_id COLLATE NOCASE, auth_asym_id) index.""" + from app.protein.records import resolve_chain_key + + conn = connect_readonly(get_database_path()) try: conn.row_factory = sqlite3.Row + stored = resolve_chain_key(conn, pdb_id, auth_asym_id) + if stored is None: + return None row = conn.execute( 'SELECT binding_status, sequence, sequence_length, cath_id, cath_superfamily,' ' Rosetta, FoldX, EvoEF2, RW, "RW+", ranked_functions,' ' state_id, experimental_method, pH, temp_K, chain_composition, matching_scores' - ' FROM node WHERE LOWER(pdb_id) = ? AND UPPER(auth_asym_id) = ? LIMIT 1', - (pdb_id.lower(), auth_asym_id.upper()), + ' FROM node WHERE pdb_id = ? COLLATE NOCASE AND auth_asym_id = ? LIMIT 1', + (pdb_id, stored), ).fetchone() if row is None: return None - return { - key: zlib.decompress(value).decode("utf-8") if isinstance(value, bytes) else value or "" - for key, value in ((key, row[key]) for key in row.keys()) - } + return {key: row[key] or "" for key in row.keys()} finally: conn.close() diff --git a/backend/assets/MuSProt_documentation.md b/backend/assets/MuSProt_documentation.md index 99f69157d7902dca00ae58b6cf1d1d5340b52f29..4c693dc660329b79251911f1ed9f401dd14635df 100755 --- a/backend/assets/MuSProt_documentation.md +++ b/backend/assets/MuSProt_documentation.md @@ -1,91 +1,192 @@ # MuSProt Dataset Documentation -**MuSProt** (Multistate Protein Database) is a million-scale multimodal database for multistate proteins, designed to support programmable protein design and AI model development. It links experimentally observed conformational states of identical protein sequences from the PDB, organizes them into state clusters and transition relationships, and enriches each record with structural similarity, experimental context, state-specific function rankings and transition fidelity labels. Users can search, browse and download MuSProt records to study conformational diversity, state-dependent functions and feasible protein state transitions. +**MuSProt** (Multistate Protein Database) is a million-scale multimodal database for multistate proteins, designed to support programmable protein design and AI model development. It links experimentally observed conformational states of identical protein sequences from the PDB, organizes them into state clusters and transition relationships, and enriches each record with structural similarity, experimental context, state-specific function rankings and transition fidelity labels. --- -## Download +## At a glance + +| | Count | +|---|---| +| Conformational observations (`node` rows, one per PDB chain) | 863,643 | +| Sequences (identical SEQRES observed ≥ 2 times) | 102,790 | +| Sequences with ≥ 2 conformational states | 24,673 | +| State clusters | 162,601 | +| Transition triplets (`edge` rows, ordered observation pairs) | 113,663,858 | +| Transitions between distinct states | 25,568,258 | + +--- + +## Concepts + +### Observation + +One protein chain (`pdb_id`, `auth_asym_id`) in one PDB entry. Observations whose full deposited sequence (SEQRES) is 100% identical share a `sequence_id`. Only sequences observed at least twice are included. + +### State cluster + +Within each sequence, two observations are linked when their **pair similarity is ≥ 0.95**. The connected components of that graph are the conformational states, numbered `0, 1, 2, …` per sequence (`state_id`). + +### Transition + +Every ordered pair (A → B) of observations of the same sequence is an `edge` row. Both directions are stored; they share all structural fields and differ only in the sign of the energy differences (B − A). + +### Structural similarity + +Pairs are superposed with US-align on a fixed residue correspondence (identical SEQRES numbering), and similarity is computed over the aligned core so that disordered or unmodelled termini do not penalise otherwise identical conformations: + +``` +tm_aln = (TM_A · L_A + TM_B · L_B) / (2 · aligned_length) +pair_similarity = 0.8 · tm_aln + 0.2 · exp(−(RMSD / 5 Å)²) +coverage = aligned_length / max(L_A, L_B) +``` -The full dataset is distributed as a single SQLite database file (`MuSProt.db`, ~31 GB). Download it from the dataset page using the **Download DB** button. +`coverage` is reported separately as the completeness axis. -For detailed protein structures and atom coordinates, users are expected to download [Protein Data Bank (PDB)](https://rcsb.org) and fetch from `mmCIF` or `.pdb` files based on the PDB entry and chain ID. +### Fidelity labels + +`pair_fidelity` (from `pair_similarity`) and `state_fidelity` (from `state_similarity`) bin how far a transition moves: + +| Label | Similarity | Meaning | +|---|---|---| +| `identical` | ≥ 0.95 | Same conformation | +| `low` | 0.70 – 0.95 | Subtle rearrangement | +| `medium` | 0.50 – 0.70 | Clear conformational change | +| `high` | ≤ 0.50 | Large-scale transition | +| empty | — | Not computable (the two chains share no resolved residues) | --- -## Database Schema +## Download + +The full dataset is a single SQLite database file (`MuSProt.db`, ~37 GB) with two tables, **`node`** and **`edge`**. Download it with the **Download** button at the top of the site. -The database contains two tables: **`node`** and **`edge`**. +Atomic coordinates are not included: fetch mmCIF files from the [RCSB PDB](https://www.rcsb.org) or [PDBe](https://www.ebi.ac.uk/pdbe/) using `pdb_id` and `auth_asym_id`. --- +## Database Schema + ### Table: `node` -Each row represents a single protein chain instance (one PDB entry + chain). +One row per observation (PDB chain). | Column | Type | Description | |---|---|---| -| `sequence_id` | TEXT | Internal sequence cluster identifier (e.g. `SEQ087470`); shared by chains with the same sequence | -| `uniprot_id` | TEXT | UniProt accession number | +| `sequence_id` | TEXT | Sequence identifier (e.g. `SEQ049914`), shared by all chains with an identical SEQRES | +| `uniprot_id` | TEXT | UniProt accession(s) from SIFTS; `;`-joined for chimeric chains; empty if none | | `pdb_id` | TEXT | 4-character PDB entry ID (lowercase) | -| `auth_asym_id` | TEXT | Author chain identifier (e.g. `A`) | -| `sequence` | TEXT | SEQRES amino acid sequence | -| `sequence_length` | INT | Number of residues in the SEQRES sequence | -| `modeled_sequence` | TEXT | Amino acid sequence actually resolved (modeled) in the structure | -| `modeled_sequence_length` | INT | Number of modeled residues | -| `binders` | TEXT | Bound ligand / small-molecule component IDs present in the structure; empty when none | -| `binding_status` | TEXT | Ligand-binding state of the chain: `apo` or `holo` | -| `experimental_method` | TEXT | Structure determination method (e.g. `X-RAY DIFFRACTION`, `ELECTRON MICROSCOPY`, `SOLUTION NMR`) | -| `resolution` | FLOAT | Experimental resolution in Å (lower is sharper; empty for methods without a resolution) | -| `pH` | FLOAT | pH of the experimental / crystallization condition | -| `temp_K` | FLOAT | Temperature of the experiment, in Kelvin | -| `experimental_details` | TEXT | Free-text crystallization / sample-preparation details | -| `chain_composition` | TEXT | Quaternary composition of the deposited assembly: `monomeric`, `homomeric`, or `heteromeric` | -| `non_protein_polymer_binding` | TEXT | Non-protein polymer bound in the structure (`DNA`, `RNA`, or `DNA/RNA`); `None` when absent | -| `initial_release_date` | TEXT | PDB entry initial release date (`YYYY-MM-DD`) | -| `cath_id` | TEXT | CATH domain assignment; empty when unassigned | -| `cath_superfamily` | TEXT | Full CATH superfamily code (e.g. `1.10.10.10`) | -| `state_id` | TEXT | Conformational state cluster this chain is assigned to (`0`, `1`, `2`, …) | -| `Rosetta` | FLOAT | Rosetta total energy score | -| `FoldX` | FLOAT | FoldX total energy score | -| `EvoEF2` | FLOAT | EvoEF2 total energy score | -| `RW` | FLOAT | Random-Walk (RW) energy score | -| `RW+` | FLOAT | Random-Walk+ (RWplus) energy score | -| `ranked_functions` | TEXT | JSON-encoded list of ranked GO/functional annotations | -| `matching_scores` | TEXT | JSON-encoded list of function-matching scores, aligned with `ranked_functions` | - - ---- +| `auth_asym_id` | TEXT | Author chain identifier (case-sensitive: `d` and `D` are different chains) | +| `sequence` | TEXT | Full deposited (canonical) sequence, including unmodelled residues | +| `sequence_length` | INT | Length of `sequence` | +| `modeled_sequence` | TEXT | Residues present in the coordinates, concatenated (gaps closed up) | +| `modeled_sequence_length` | INT | Length of `modeled_sequence` | +| `resolved_label_seq` | TEXT | Resolved SEQRES positions as ranges (`label_seq_id`, 1-based), e.g. `5-173,187-353` | +| `n_resolved_aa` | INT | Number of resolved standard amino acids | +| `resolved_coverage` | FLOAT | `n_resolved_aa / sequence_length` | +| `binders` | TEXT | Ligands, metals and cofactors on this chain (component IDs, `;`-joined; water and common buffers excluded) | +| `binding_status` | TEXT | `holo` if the chain has any binder, otherwise `apo` | +| `experimental_method` | TEXT | e.g. `X-RAY DIFFRACTION`, `ELECTRON MICROSCOPY`, `SOLUTION NMR` | +| `resolution` | FLOAT | Resolution in Å; empty when not applicable (e.g. NMR) | +| `pH` | FLOAT | pH of the crystallisation / sample condition | +| `temp_K` | FLOAT | Temperature in Kelvin | +| `experimental_details` | TEXT | Free-text crystallisation / sample-preparation details | +| `chain_composition` | TEXT | Entry composition: `monomeric`, `homomeric` or `heteromeric` | +| `non_protein_polymer_binding` | TEXT | Bound nucleic acid: `DNA`, `RNA`, `DNA/RNA`, `Other polymers` or `None` | +| `initial_release_date` | TEXT | Initial PDB release date (`YYYY-MM-DD`) | +| `cath_id` | TEXT | CATH domain IDs for the chain (e.g. `3vgkA01;3vgkA02`); empty if unassigned | +| `cath_superfamily` | TEXT | CATH superfamily code(s) (e.g. `3.30.420.40`) | +| `pfam_id` | TEXT | Pfam family accession(s), `;`-joined (e.g. `PF00069`); empty if none | +| `ecod_id` | TEXT | ECOD domain ID(s) for the chain, `;`-joined (e.g. `e3vgkA1;e3vgkA2`) | +| `ecod_fid` | TEXT | ECOD family ID(s) `x.h.t.f` (e.g. `2484.1.1.8`); the first two levels (`x.h`) are the homology group | +| `state_id` | TEXT | Conformational state of this observation within its sequence (`0`, `1`, …) | +| `Rosetta` | FLOAT | Rosetta total score | +| `FoldX` | FLOAT | FoldX total energy | +| `EvoEF2` | FLOAT | EvoEF2 total energy | +| `RW` | FLOAT | RW statistical potential | +| `RW+` | FLOAT | RW+ statistical potential | +| `n_seqres` | INT | SEQRES length used for energy scoring | +| `n_resolved` | INT | Resolved residues used for scoring | +| `n_scored` | INT | Residues in the consensus set scored by Rosetta / EvoEF2 / RW / RW+ | +| `scored_coverage` | FLOAT | `n_scored / n_seqres` | +| `foldx_n_scored` | INT | Residues scored by FoldX | +| `foldx_coverage` | FLOAT | `foldx_n_scored / n_seqres` | +| `ranked_functions` | TEXT | JSON list of functional annotations, ranked for this observation | +| `matching_scores` | TEXT | JSON list of matching scores, index-aligned with `ranked_functions` | + +Energies are computed on a consensus residue set shared by all observations of a sequence, so they are comparable within a sequence. Lower is more favourable. ### Table: `edge` -Each row represents a pairwise structural comparison between two chain instances (A and B) that share the same UniProt identity. +One row per ordered pair of observations (A → B) of the same sequence. | Column | Type | Description | |---|---|---| -| `pdb_id_A` | TEXT | PDB ID of chain A | -| `auth_asym_id_A` | TEXT | Chain identifier of chain A | -| `pdb_id_B` | TEXT | PDB ID of chain B | -| `auth_asym_id_B` | TEXT | Chain identifier of chain B | -| `TM1` | FLOAT | TM-score of the alignment (chain A as reference) | -| `RMSD` | FLOAT | Root-mean-square deviation of Cα atoms (Å) | -| `pair_similarity` | FLOAT | Composite structural similarity score for the A–B pair | -| `delta_Rosetta` | FLOAT | Rosetta energy difference (B − A) | -| `delta_FoldX` | FLOAT | FoldX energy difference (B − A) | -| `delta_EvoEF2` | FLOAT | EvoEF2 energy difference (B − A) | -| `delta_RW` | FLOAT | Random-Walk (RW) energy difference (B − A) | -| `delta_RW+` | FLOAT | Random-Walk+ (RWplus) energy difference (B − A) | -| `state_id_A` | TEXT | Conformational state cluster of chain A | -| `state_id_B` | TEXT | Conformational state cluster of chain B | -| `state_similarity` | TEXT | Average structural similarity within the state cluster | -| `state_fidelity` | TEXT | State-level transition fidelity label: `identical`, `high`, `medium`, or `low` | -| `pair_fidelity` | TEXT | Pairwise (observation-level) transition fidelity label: `identical`, `high`, `medium`, or `low` | -| `sequence_id` | TEXT | Internal sequence cluster identifier shared by chains A and B | +| `pdb_id_A` | TEXT | PDB ID of observation A | +| `auth_asym_id_A` | TEXT | Chain of observation A | +| `pdb_id_B` | TEXT | PDB ID of observation B | +| `auth_asym_id_B` | TEXT | Chain of observation B | +| `tm_aln` | FLOAT | TM-score normalised by the aligned length (conformation only) | +| `RMSD` | FLOAT | Cα RMSD over the aligned residues (Å) | +| `coverage` | FLOAT | `aligned_length / max(L_A, L_B)` | +| `pair_similarity` | FLOAT | `0.8 · tm_aln + 0.2 · exp(−(RMSD/5)²)` | +| `delta_Rosetta` | FLOAT | Rosetta(B) − Rosetta(A) | +| `delta_FoldX` | FLOAT | FoldX(B) − FoldX(A) | +| `delta_EvoEF2` | FLOAT | EvoEF2(B) − EvoEF2(A) | +| `delta_RW` | FLOAT | RW(B) − RW(A) | +| `delta_RW+` | FLOAT | RW+(B) − RW+(A) | +| `state_id_A` | TEXT | State of A | +| `state_id_B` | TEXT | State of B | +| `state_similarity` | FLOAT | Mean `pair_similarity` over all observation pairs between the two states (within-state when `state_id_A = state_id_B`) | +| `state_fidelity` | TEXT | Fidelity label of `state_similarity` | +| `pair_fidelity` | TEXT | Fidelity label of `pair_similarity` | +| `n_shared_resolved` | INT | Residues resolved in both A and B | +| `delta_order_valid` | TEXT | `True` when resolved-residue changes are comparable (same method, not both NMR, overlapping models) | +| `order_invalid_reason` | TEXT | `valid`, `cross_method`, `nmr`, `no_shared_residues` or `missing_data` | +| `n_ordered` | INT | Residues resolved in B but not in A (ordered on going A → B) | +| `n_disordered` | INT | Residues resolved in A but not in B | +| `ordered_ranges` | TEXT | SEQRES positions of the ordered residues, e.g. `115-123,187-189` | +| `disordered_ranges` | TEXT | SEQRES positions of the disordered residues | +| `longest_ordered_segment` | INT | Longest contiguous ordered stretch | +| `longest_disordered_segment` | INT | Longest contiguous disordered stretch | +| `n_terminal_ordered` | INT | Ordered residues in segments touching either terminus | +| `n_internal_ordered` | INT | Ordered residues in internal segments | +| `order_evidence` | TEXT | Evidence grade for the ordering: `high`, `medium`, `low` or `not_applicable` (see below) | +| `ordering_with_ligand_change` | TEXT | `True` when the ordering co-occurs with a change in bound ligands | +| `sequence_id` | TEXT | Sequence shared by A and B | + +When A and B share no resolved residues, `tm_aln`, `pair_similarity` and `pair_fidelity` are empty and `RMSD` / `coverage` are stored as `0.0`. + +### Order / disorder labels + +These columns describe the *observed* resolved state (residues with or without modelled coordinates), not intrinsic disorder: a residue can lack density because it is flexible, or because of crystal contacts, occupancy or modelling choices. They are only computed when both structures use the same experimental method and at least one is not NMR (NMR models build every residue); otherwise `delta_order_valid` is `False` and the counts are empty. + +`order_evidence` grades an ordering event from its length, position and the resolution of both structures: + +| Grade | Rule | +|---|---| +| `high` | an internal ordered segment of ≥ 8 residues, both structures better than 2.5 Å | +| `medium` | otherwise, if not `low` | +| `low` | longest ordered segment ≤ 2 residues, purely terminal, or either structure worse than 3.5 Å | +| `not_applicable` | nothing was ordered | + +### Knowledge graph + +MuSProt is organised as a hierarchical graph. In the **outer layer**, each sequence is a node, and two sequences are linked when their chains carry the same set of ECOD homology groups (`x.h` of `ecod_fid`; 88.1% of sequences have an assignment). In the **inner layer**, every observation of a sequence is connected to every other by an `edge` row; contracting each state cluster to one node gives the quotient graph over states, whose edges carry `state_similarity`. + +### Indexes + +Both tables are indexed for per-chain lookup: `node(pdb_id COLLATE NOCASE, auth_asym_id)` and `edge(pdb_id_A COLLATE NOCASE, auth_asym_id_A)`. Use the same collation in queries to hit the index: + +```sql +SELECT * FROM edge WHERE pdb_id_A = '3vgk' COLLATE NOCASE AND auth_asym_id_A = 'A'; +``` --- ## Usage Examples -### Python (sqlite3) +### Python (sqlite3 + pandas) ```python import sqlite3 @@ -93,59 +194,55 @@ import pandas as pd conn = sqlite3.connect("MuSProt.db") -# Load all chains for a UniProt entry -df_nodes = pd.read_sql( - "SELECT * FROM node WHERE uniprot_id = 'P00533'", - conn +# One observation +node = pd.read_sql( + "SELECT * FROM node WHERE pdb_id = ? COLLATE NOCASE AND auth_asym_id = ?", + conn, params=("3vgk", "A"), ) -# Find all structural neighbours of a given chain -df_edges = pd.read_sql( - "SELECT * FROM edge WHERE pdb_id_A = '1ivo' AND auth_asym_id_A = 'A'", - conn +# Every transition from that observation +edges = pd.read_sql( + "SELECT * FROM edge WHERE pdb_id_A = ? COLLATE NOCASE AND auth_asym_id_A = ?", + conn, params=("3vgk", "A"), ) - -conn.close() ``` -### Filter by structural similarity +### Transitions between distinct states with a large change ```python -# Retrieve pairs with high TM-score and low RMSD df = pd.read_sql(""" - SELECT * + SELECT pdb_id_A, auth_asym_id_A, pdb_id_B, auth_asym_id_B, + CAST(pair_similarity AS REAL) AS similarity, + CAST(RMSD AS REAL) AS rmsd, + CAST(delta_Rosetta AS REAL) AS d_rosetta FROM edge - WHERE CAST(TM1 AS REAL) > 0.8 - AND CAST(RMSD AS REAL) < 2.0 - LIMIT 1000 -""", conn) + WHERE pdb_id_A = ? COLLATE NOCASE AND auth_asym_id_A = ? + AND state_id_A <> state_id_B + AND pair_fidelity IN ('medium', 'high') + ORDER BY similarity +""", conn, params=("1anf", "A")) ``` -### Join nodes and edges +### Apo / holo pairs for one UniProt accession ```python -# Get full info for both chains in each pair df = pd.read_sql(""" - SELECT - e.pdb_id_A, e.auth_asym_id_A, - e.pdb_id_B, e.auth_asym_id_B, - e.TM1, e.RMSD, - nA.sequence_length AS len_A, - nB.sequence_length AS len_B, - nA.cath_superfamily - FROM edge e - JOIN node nA ON e.pdb_id_A = nA.pdb_id AND e.auth_asym_id_A = nA.auth_asym_id - JOIN node nB ON e.pdb_id_B = nB.pdb_id AND e.auth_asym_id_B = nB.auth_asym_id - WHERE nA.uniprot_id = 'P00533' - LIMIT 500 + SELECT e.pdb_id_A, e.auth_asym_id_A, e.pdb_id_B, e.auth_asym_id_B, + e.pair_similarity, e.pair_fidelity, nB.binders + FROM node nA + JOIN edge e ON e.pdb_id_A = nA.pdb_id COLLATE NOCASE AND e.auth_asym_id_A = nA.auth_asym_id + JOIN node nB ON nB.pdb_id = e.pdb_id_B COLLATE NOCASE AND nB.auth_asym_id = e.auth_asym_id_B + WHERE nA.uniprot_id = 'P0AEX9' + AND nA.binding_status = 'apo' AND nB.binding_status = 'holo' """, conn) ``` +Scanning `node` by `uniprot_id` or `sequence_id` reads the whole table; for repeated queries create an index first (`CREATE INDEX idx_node_seq ON node(sequence_id)`). + --- ## Notes -- All numeric fields (energies, TM-score, RMSD, lengths) are stored as `TEXT`; cast them with `CAST(col AS REAL)` or `CAST(col AS INTEGER)` as needed. -- `ranked_functions` in the `node` table is a JSON string. Parse it with `json.loads()`. -- The `edge` table is directional: (A→B) and (B→A) are separate rows and may differ slightly in TM-score. -- Energy delta values represent B − A; a negative delta means chain B is lower energy than chain A. +- All columns are stored as `TEXT`; cast numeric fields with `CAST(col AS REAL)` / `CAST(col AS INTEGER)`. +- `ranked_functions` and `matching_scores` are JSON strings; parse with `json.loads()`. +- `edge` is directional: A → B and B → A are separate rows with identical structural fields and opposite energy deltas. diff --git a/backend/assets/musprot_summary.json b/backend/assets/musprot_summary.json index 4f6d478d73e572bb6bae91a35ffb1c94092c4369..6cce643e135289d54a60ddb03b41212c1564a6d2 100644 --- a/backend/assets/musprot_summary.json +++ b/backend/assets/musprot_summary.json @@ -1,31 +1,31 @@ { - "total_records": 27586348, - "unique_chains": 325320, - "avg_rmsd": 0.6675476253688963, - "avg_tm_score": 0.9602746188440746, - "avg_sequence_length": 235.7717984028132, + "total_records": 113582934, + "unique_chains": 793278, + "avg_rmsd": 1.021488560772695, + "avg_tm_score": 0.9077313277727972, + "avg_sequence_length": 304.6724704328117, "rmsd_distribution": { - "0-0.5": 11970560, - "0.5-1.0": 10422936, - "1.0-1.5": 3551996, - "1.5-2.0": 905068, - "2.0-5.0": 730332, - ">5.0": 5456 + "0-0.5": 34431100, + "0.5-1.0": 39665526, + "1.0-1.5": 17562924, + "1.5-2.0": 6594622, + "2.0-5.0": 15208444, + ">5.0": 120318 }, "tm_score_distribution": { - "0-0.5": 144940, - "0.5-0.7": 319272, - "0.7-0.85": 917410, - "0.85-0.95": 4632994, - "0.95-1.0": 21571732 + "0-0.5": 4580963, + "0.5-0.7": 4431774, + "0.7-0.85": 12281664, + "0.85-0.95": 20890577, + "0.95-1.0": 71397956 }, "length_distribution": { - "300-500": 57390, - "500-1000": 19801, - "100-200": 101177, - "0-100": 76718, - "200-300": 65315, - ">1000": 4921 + "300-500": 193657, + "500-1000": 86502, + "200-300": 162341, + ">1000": 20178, + "0-100": 110685, + "100-200": 219915 }, "tm_score_range": { "min": 0.0, @@ -33,10 +33,10 @@ }, "rmsd_range": { "min": 0.0, - "max": 9.83 + "max": 12.44 }, "length_range": { - "min": 1, - "max": 4542 + "min": 3, + "max": 5247 } } \ No newline at end of file diff --git a/backend/scripts/build_index.py b/backend/scripts/build_index.py new file mode 100644 index 0000000000000000000000000000000000000000..4ba8f8d043cfc6e6093bf2b1a6e004b839faf1ae --- /dev/null +++ b/backend/scripts/build_index.py @@ -0,0 +1,656 @@ +#!/usr/bin/env python3 +""" +Build the website sidecars from MuSProt.db. + +Outputs (into --out-dir): + MuSProt-index.db compact, indexed lookup tables the website queries directly + sequence one row per sequence_id (catalog: counts, representative, function) + member one row per chain observation (light columns only) + state one row per (sequence_id, state_id) + state_pair mean similarity / fidelity between two states of a sequence + function_fts FTS5 over UniProt / CATH / ECOD / Pfam / top ranked functions per sequence + ecod_h, ecod_f ECOD homology / family names (from ecod.latest.domains.txt) + ecod_group outer-graph groups: sequences sharing the same set of ECOD H-groups + overview.json dataset-level counts and histograms for the home page + +ECOD names (homology / family) come from ECOD's own domain list via --ecod. + +The main DB is only read. Everything is derived in a few SQL scans, so the +edge table (100M+ rows) is aggregated inside SQLite rather than in Python. + +Usage: + python backend/scripts/build_index.py /path/MuSProt.db --out-dir /tmp/musprot-bucket \ + --ecod /path/ecod.latest.domains.txt +""" +from __future__ import annotations + +import argparse +import json +import math +import sqlite3 +import time +from collections import Counter, defaultdict +from pathlib import Path + +FUNCTION_TEXT_TOP_N = 3 +FUNCTION_TEXT_MAX_CHARS = 600 +MAX_PAIR_STATES = 40 # state_pair rows are kept only among a sequence's largest states + + +def log(msg: str) -> None: + print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True) + + +def to_float(v): + try: + f = float(v) + return None if math.isnan(f) else f + except (TypeError, ValueError): + return None + + +def to_int(v): + f = to_float(v) + return None if f is None else int(f) + + +def parse_list(raw): + raw = (raw or "").strip() + if not raw: + return [] + try: + val = json.loads(raw) + return val if isinstance(val, list) else [] + except ValueError: + import ast + try: + val = ast.literal_eval(raw) + return val if isinstance(val, list) else [] + except (ValueError, SyntaxError): + return [] + + +def hist(values, edges): + """Counts of values in [edges[i], edges[i+1]); last bin is closed.""" + counts = [0] * (len(edges) - 1) + for v in values: + if v is None: + continue + for i in range(len(edges) - 1): + if v < edges[i + 1] or i == len(edges) - 2: + if v >= edges[i]: + counts[i] += 1 + break + return counts + + +# ─────────────────────────────────────────────────────────────── domain inputs +def split_ids(raw) -> list[str]: + return [x for x in (raw or "").split(";") if x] + + +def ecod_hset(fids: str) -> tuple: + """ECOD F-ids ('292.2.1.1;11.1.1.3') -> sorted homology-level set ('11.1', '292.2').""" + return tuple(sorted({".".join(f.split(".")[:2]) for f in split_ids(fids)})) + + +def load_ecod_names(path: Path | None) -> tuple[dict, dict]: + """f_id -> (t_name, f_name) and h_id -> (architecture, x_name, h_name).""" + fam, hom = {}, {} + if path is None: + return fam, hom + header = None + with open(path, encoding="utf-8") as fh: + for line in fh: + if line.startswith("#"): + cols = line.lstrip("#").rstrip("\n").split("\t") + if "f_id" in cols: + header = cols + continue + parts = line.rstrip("\n").split("\t") + if header is None: + if "f_id" in parts: + header = parts + continue + r = dict(zip(header, parts)) + f_id = r.get("f_id", "") + if not f_id or f_id in fam: + continue + clean = lambda v: (v or "").strip('"') if (v or "").strip('"') not in ("NO_X_NAME", "NO_H_NAME", "NO_T_NAME", "F_UNCLASSIFIED") else "" + fam[f_id] = (clean(r.get("t_name")), clean(r.get("f_name"))) + h_id = ".".join(f_id.split(".")[:2]) + if h_id not in hom: + hom[h_id] = (clean(r.get("architecture_name")), clean(r.get("x_name")), clean(r.get("h_name"))) + return fam, hom + + +# ─────────────────────────────────────────────────────────────── node pass +NODE_SQL = """ +SELECT sequence_id, uniprot_id, pdb_id, auth_asym_id, sequence_length, + n_resolved_aa, resolved_coverage, binders, binding_status, + experimental_method, resolution, pH, temp_K, chain_composition, + non_protein_polymer_binding, initial_release_date, cath_superfamily, + state_id, ranked_functions, cath_id, pfam_id, ecod_id, ecod_fid +FROM node +""" + + +def build_members(src: sqlite3.Connection, dst: sqlite3.Connection, limit: int | None): + dst.executescript( + """ + DROP TABLE IF EXISTS member; + CREATE TABLE member ( + sequence_id TEXT NOT NULL, + state_id INTEGER, + pdb_id TEXT NOT NULL, + auth_asym_id TEXT NOT NULL, + uniprot_id TEXT, + experimental_method TEXT, + resolution REAL, + binding_status TEXT, + chain_composition TEXT, + binders TEXT, + release_date TEXT, + resolved_coverage REAL, + n_resolved_aa INTEGER, + pH REAL, + temp_K REAL, + cath_id TEXT, + pfam_id TEXT, + ecod_id TEXT, + ecod_fid TEXT + ); + """ + ) + seqs: dict[str, dict] = {} + functions: dict[str, list[str]] = {} + ov = { + "method": Counter(), "binding": Counter(), "composition": Counter(), + "npp": Counter(), "year": Counter(), "cath": Counter(), + "resolution": [], "length": {}, "uniprot": set(), "entries": set(), + } + sql = NODE_SQL + (f" LIMIT {int(limit)}" if limit else "") + batch = [] + n = 0 + for row in src.execute(sql): + (seq_id, uniprot, pdb, chain, seq_len, n_res, res_cov, binders, binding, + method, resolution, ph, temp, comp, npp, date, cath_sf, state, funcs, cath_id, + pfam, ecod_ids, ecod_fid) = row + state_i = to_int(state) + res_f = to_float(resolution) + batch.append(( + seq_id, state_i, pdb, chain, uniprot or None, method or None, res_f, + binding or None, comp or None, binders or None, date or None, + to_float(res_cov), to_int(n_res), to_float(ph), to_float(temp), + cath_id or None, pfam or None, ecod_ids or None, ecod_fid or None, + )) + s = seqs.get(seq_id) + if s is None: + s = seqs[seq_id] = { + "uniprot": Counter(), "length": to_int(seq_len), "n_obs": 0, + "states": Counter(), "n_holo": 0, "entries": set(), "methods": Counter(), + "cath": Counter(), "dates": [], "rep": None, + "ecod": Counter(), "pfam": Counter(), + } + s["n_obs"] += 1 + if uniprot: + s["uniprot"][uniprot] += 1 + s["states"][state_i] += 1 + s["n_holo"] += binding == "holo" + s["entries"].add(pdb) + if method: + s["methods"][method] += 1 + if cath_sf: + s["cath"][cath_sf] += 1 + hs = ecod_hset(ecod_fid) + if hs: + s["ecod"][hs] += 1 + pf = tuple(sorted(set(split_ids(pfam)))) + if pf: + s["pfam"][pf] += 1 + if date: + s["dates"].append(date) + # representative: best (lowest) resolution among well-resolved chains + key = ((to_float(res_cov) or 0) < 0.9, res_f if res_f is not None else 99.0) + if s["rep"] is None or key < s["rep"][0]: + s["rep"] = (key, pdb, chain) + if seq_id not in functions and funcs: + fl = [str(f).strip() for f in parse_list(funcs) if len(str(f).strip()) > 10] + if fl: + functions[seq_id] = fl[:FUNCTION_TEXT_TOP_N] + + ov["method"][method or "Unknown"] += 1 + ov["binding"][binding or "Unknown"] += 1 + ov["composition"][comp or "Unknown"] += 1 + ov["npp"][npp or "None"] += 1 + if date: + ov["year"][date[:4]] += 1 + if res_f is not None: + ov["resolution"].append(res_f) + if uniprot: + ov["uniprot"].add(uniprot) + ov["entries"].add(pdb.lower()) + ov["length"][seq_id] = to_int(seq_len) + + n += 1 + if len(batch) >= 50000: + dst.executemany(f"INSERT INTO member VALUES ({','.join('?' * 19)})", batch) + batch.clear() + log(f" node rows: {n:,}") + if batch: + dst.executemany(f"INSERT INTO member VALUES ({','.join('?' * 19)})", batch) + log(f" node rows total: {n:,}; sequences: {len(seqs):,}") + dst.executescript( + """ + CREATE INDEX idx_member_seq ON member(sequence_id, state_id); + CREATE INDEX idx_member_pdb ON member(pdb_id COLLATE NOCASE, auth_asym_id); + CREATE INDEX idx_member_uniprot ON member(uniprot_id COLLATE NOCASE); + """ + ) + return seqs, functions, ov, n + + +def consensus(counter: Counter) -> tuple: + """Most common non-empty assignment set across a sequence's chains (ties: larger set).""" + if not counter: + return () + return max(counter.items(), key=lambda kv: (kv[1], len(kv[0]), kv[0]))[0] + + +def write_sequences(dst: sqlite3.Connection, seqs: dict, functions: dict, ecod_names: tuple): + dst.executescript( + """ + DROP TABLE IF EXISTS sequence; + CREATE TABLE sequence ( + sequence_id TEXT PRIMARY KEY, + uniprot_id TEXT, + length INTEGER, + n_obs INTEGER, + n_states INTEGER, + n_entries INTEGER, + n_holo INTEGER, + n_apo INTEGER, + n_cross_pairs INTEGER, + largest_state INTEGER, + rep_pdb TEXT, + rep_chain TEXT, + cath_superfamily TEXT, + methods TEXT, + first_release TEXT, + last_release TEXT, + top_function TEXT, + ecod_hset TEXT, + pfam_ids TEXT + ); + DROP TABLE IF EXISTS state; + CREATE TABLE state ( + sequence_id TEXT NOT NULL, + state_id INTEGER NOT NULL, + n_members INTEGER, + PRIMARY KEY (sequence_id, state_id) + ); + """ + ) + rows, state_rows = [], [] + for seq_id, s in seqs.items(): + sizes = list(s["states"].values()) + n_obs = s["n_obs"] + cross = n_obs * n_obs - sum(x * x for x in sizes) # directed cross-state pairs + dates = sorted(s["dates"]) + rows.append(( + seq_id, + s["uniprot"].most_common(1)[0][0] if s["uniprot"] else None, + s["length"], n_obs, len(sizes), len(s["entries"]), s["n_holo"], n_obs - s["n_holo"], + cross, max(sizes), + s["rep"][1], s["rep"][2], + s["cath"].most_common(1)[0][0] if s["cath"] else None, + ";".join(m for m, _ in s["methods"].most_common()), + dates[0] if dates else None, dates[-1] if dates else None, + (functions.get(seq_id) or [None])[0], + ";".join(consensus(s["ecod"])) or None, + ";".join(consensus(s["pfam"])) or None, + )) + for st, cnt in s["states"].items(): + state_rows.append((seq_id, st, cnt)) + dst.executemany(f"INSERT INTO sequence VALUES ({','.join('?' * 19)})", rows) + dst.executemany("INSERT INTO state VALUES (?,?,?)", state_rows) + dst.executescript( + """ + CREATE INDEX idx_seq_states ON sequence(n_states DESC, n_obs DESC); + CREATE INDEX idx_seq_obs ON sequence(n_obs DESC); + CREATE INDEX idx_seq_uniprot ON sequence(uniprot_id COLLATE NOCASE); + CREATE INDEX idx_seq_cath ON sequence(cath_superfamily); + CREATE INDEX idx_seq_ecod ON sequence(ecod_hset, n_obs DESC); + """ + ) + + fam, hom = ecod_names + dst.executescript( + """ + DROP TABLE IF EXISTS ecod_h; + CREATE TABLE ecod_h (h_id TEXT PRIMARY KEY, architecture TEXT, x_name TEXT, h_name TEXT); + DROP TABLE IF EXISTS ecod_f; + CREATE TABLE ecod_f (f_id TEXT PRIMARY KEY, t_name TEXT, f_name TEXT); + DROP TABLE IF EXISTS ecod_group; + CREATE TABLE ecod_group ( + hset TEXT PRIMARY KEY, + label TEXT, + n_sequences INTEGER, + n_obs INTEGER, + n_multistate INTEGER, + n_states INTEGER + ); + """ + ) + dst.executemany("INSERT INTO ecod_h VALUES (?,?,?,?)", [(k, *v) for k, v in hom.items()]) + dst.executemany("INSERT INTO ecod_f VALUES (?,?,?)", [(k, *v) for k, v in fam.items()]) + groups: dict[tuple, list] = {} + for s in seqs.values(): + hs = consensus(s["ecod"]) + if hs: + g = groups.setdefault(hs, [0, 0, 0, 0]) + g[0] += 1 + g[1] += s["n_obs"] + g[2] += len(s["states"]) > 1 + g[3] += len(s["states"]) + dst.executemany( + "INSERT INTO ecod_group VALUES (?,?,?,?,?,?)", + [ + (";".join(hs), " + ".join(hom.get(h, ("", "", ""))[2] or h for h in hs), *v) + for hs, v in groups.items() + ], + ) + dst.execute("CREATE INDEX idx_group_size ON ecod_group(n_sequences DESC)") + + dst.executescript( + """ + DROP TABLE IF EXISTS function_fts; + CREATE VIRTUAL TABLE function_fts USING fts5( + sequence_id UNINDEXED, uniprot_id, cath_superfamily, domains, text, + tokenize = 'porter unicode61' + ); + """ + ) + fts_rows = [] + for seq_id, s in seqs.items(): + text = " ".join(functions.get(seq_id, []))[:FUNCTION_TEXT_MAX_CHARS * FUNCTION_TEXT_TOP_N] + uni = " ".join(s["uniprot"].keys()) + cath = " ".join(s["cath"].keys()).replace(";", " ") + hs = consensus(s["ecod"]) + dom = " ".join( + [*(h.replace(".", "_") for h in hs), *(hom.get(h, ("", "", ""))[2] for h in hs), + *consensus(s["pfam"])] + ) + fts_rows.append((seq_id, uni, cath, dom, text)) + dst.executemany("INSERT INTO function_fts VALUES (?,?,?,?,?)", fts_rows) + + +# ─────────────────────────────────────────────────────────────── edge passes +def build_state_pairs(src: sqlite3.Connection, dst: sqlite3.Connection, seqs: dict, + limit: int | None) -> dict: + """Store state-pair similarities (top states only); return fidelity counts over all pairs.""" + keep = { + seq_id: {st for st, _ in s["states"].most_common(MAX_PAIR_STATES)} + for seq_id, s in seqs.items() if len(s["states"]) > MAX_PAIR_STATES + } + fid_counts: dict = defaultdict(lambda: {"state_pairs": 0, "observation_pairs": 0}) + dst.executescript( + """ + DROP TABLE IF EXISTS state_pair; + CREATE TABLE state_pair ( + sequence_id TEXT NOT NULL, + state_a INTEGER NOT NULL, + state_b INTEGER NOT NULL, + similarity REAL, + fidelity TEXT, + n_pairs INTEGER, + PRIMARY KEY (sequence_id, state_a, state_b) + ); + """ + ) + src_tbl = f"(SELECT * FROM edge LIMIT {int(limit)})" if limit else "edge" + # directed edges -> keep one direction (A <= B) and count undirected pairs + sql = f""" + SELECT sequence_id, CAST(state_id_A AS INTEGER) a, CAST(state_id_B AS INTEGER) b, + MIN(CAST(state_similarity AS REAL)), MIN(state_fidelity), COUNT(*) + FROM {src_tbl} + WHERE CAST(state_id_A AS INTEGER) <= CAST(state_id_B AS INTEGER) + GROUP BY sequence_id, a, b + """ + cur = src.execute(sql) + n = 0 + while True: + rows = cur.fetchmany(100000) + if not rows: + break + out = [] + for s, a, b, sim, fid, c in rows: + # within-state rows appear in both directions with A==B -> halve the count + c = c // 2 if a == b else c + if a != b: + fc = fid_counts[fid or "NA"] + fc["state_pairs"] += 1 + fc["observation_pairs"] += c + allowed = keep.get(s) + if allowed is None or (a in allowed and b in allowed): + out.append((s, a, b, sim, fid, c)) + dst.executemany("INSERT INTO state_pair VALUES (?,?,?,?,?,?)", out) + n += len(rows) + log(f" state pairs: {n:,}") + return dict(fid_counts) + + +def edge_histograms(src: sqlite3.Connection, limit: int | None) -> dict: + src_tbl = f"(SELECT * FROM edge LIMIT {int(limit)})" if limit else "edge" + sql = f""" + SELECT state_id_A = state_id_B AS same_state, + COALESCE(NULLIF(pair_fidelity, ''), 'NA') AS fid, + -- pairs without shared residues store '' similarity and 0.0 RMSD/coverage: keep them + -- out of every histogram (they are counted under fidelity 'NA') + CASE WHEN pair_similarity = '' THEN NULL + ELSE MIN(CAST(CAST(pair_similarity AS REAL) * 50 AS INTEGER), 49) END AS sim_bin, + CASE WHEN pair_similarity = '' THEN NULL + ELSE MIN(CAST(CAST(RMSD AS REAL) * 2 AS INTEGER), 40) END AS rmsd_bin, + CASE WHEN pair_similarity = '' THEN NULL + ELSE MIN(CAST(CAST(coverage AS REAL) * 20 AS INTEGER), 19) END AS cov_bin, + CASE WHEN tm_aln = '' THEN NULL + ELSE MIN(CAST(CAST(tm_aln AS REAL) * 50 AS INTEGER), 49) END AS tm_bin, + COUNT(*) + FROM {src_tbl} + GROUP BY 1, 2, 3, 4, 5, 6 + """ + out = { + "total": 0, "cross_state": 0, + "pair_fidelity": Counter(), "pair_fidelity_cross": Counter(), + "similarity": [0] * 50, "similarity_cross": [0] * 50, + "tm_aln": [0] * 50, "rmsd": [0] * 41, "coverage": [0] * 20, + } + for same, fid, sim_b, rmsd_b, cov_b, tm_b, c in src.execute(sql): + out["total"] += c + out["pair_fidelity"][fid] += c + if not same: + out["cross_state"] += c + out["pair_fidelity_cross"][fid] += c + if sim_b is not None and sim_b >= 0: + out["similarity"][sim_b] += c + if not same: + out["similarity_cross"][sim_b] += c + if tm_b is not None and tm_b >= 0: + out["tm_aln"][tm_b] += c + if rmsd_b is not None and rmsd_b >= 0: + out["rmsd"][rmsd_b] += c + if cov_b is not None and cov_b >= 0: + out["coverage"][cov_b] += c + return out + + +def order_stats(src: sqlite3.Connection, limit: int | None) -> dict: + """Counts of the observed order/disorder labels over all directed transitions.""" + src_tbl = f"(SELECT * FROM edge LIMIT {int(limit)})" if limit else "edge" + rows = src.execute(f""" + SELECT order_invalid_reason, order_evidence, ordering_with_ligand_change, + state_id_A = state_id_B, COUNT(*) + FROM {src_tbl} GROUP BY 1, 2, 3, 4 + """).fetchall() + reason, evidence, evidence_cross, ligand = Counter(), Counter(), Counter(), Counter() + for r, ev, lig, same, c in rows: + reason[r or "NA"] += c + if r == "valid": + evidence[ev or "NA"] += c + if not same: + evidence_cross[ev or "NA"] += c + if ev in ("low", "medium", "high"): + ligand[(ev, lig == "True")] += c + order = ["high", "medium", "low", "not_applicable"] + return { + "validity": [{"label": k, "value": v} for k, v in reason.most_common()], + "evidence": [{"label": k, "value": evidence.get(k, 0), "cross_state": evidence_cross.get(k, 0)} for k in order], + "with_ligand_change": [ + {"label": k, "with": ligand.get((k, True), 0), "without": ligand.get((k, False), 0)} + for k in order[:3] + ], + } + + +# ─────────────────────────────────────────────────────────────── overview +def outer_graph_stats(dst: sqlite3.Connection, n_sequences: int) -> dict: + rows = dst.execute( + "SELECT hset, label, n_sequences, n_obs, n_multistate FROM ecod_group ORDER BY n_sequences DESC" + ).fetchall() + covered = sum(r[2] for r in rows) + size_bins = Counter() + for r in rows: + n = r[2] + size_bins["1" if n == 1 else "2–5" if n <= 5 else "6–20" if n <= 20 else "21–100" if n <= 100 else ">100"] += 1 + return { + "level": "ECOD homology (H)", + "groups": len(rows), + "covered_sequences": covered, + "coverage": covered / n_sequences if n_sequences else 0, + "edges": sum(r[2] * (r[2] - 1) // 2 for r in rows), + "group_sizes": [{"label": k, "value": size_bins[k]} for k in ["1", "2–5", "6–20", "21–100", ">100"]], + "largest": [ + {"hset": r[0], "label": r[1], "n_sequences": r[2], "n_obs": r[3], "n_multistate": r[4]} + for r in rows[:12] + ], + } + + +def build_overview(seqs, ov, n_nodes, edges, state_fid) -> dict: + states_per_seq = Counter(len(s["states"]) for s in seqs.values()) + obs_per_seq = [s["n_obs"] for s in seqs.values()] + n_states_total = sum(len(s["states"]) for s in seqs.values()) + multi_state = sum(1 for s in seqs.values() if len(s["states"]) > 1) + lengths = [v for v in ov["length"].values() if v] + + def capped(counter, cap): + out = Counter() + for k, v in counter.items(): + out[k if k < cap else cap] += v + return [{"label": (f"{k}+" if k == cap else str(k)), "value": out[k]} for k in sorted(out)] + + obs_edges = [2, 3, 4, 5, 6, 11, 21, 51, 101, 1_000_000] + obs_labels = ["2", "3", "4", "5", "6–10", "11–20", "21–50", "51–100", ">100"] + len_edges = [0, 100, 200, 300, 400, 500, 750, 1000, 1_000_000] + len_labels = ["<100", "100–199", "200–299", "300–399", "400–499", "500–749", "750–999", "≥1000"] + res_edges = [0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 1000] + res_labels = ["<1.5", "1.5–2.0", "2.0–2.5", "2.5–3.0", "3.0–3.5", "3.5–4.0", "4.0–5.0", "≥5.0"] + fid_order = ["identical", "low", "medium", "high", "no_shared_residues", "NA"] + + return { + "generated": time.strftime("%Y-%m-%d"), + "counts": { + "observations": n_nodes, + "sequences": len(seqs), + "multi_state_sequences": multi_state, + "state_clusters": n_states_total, + "pdb_entries": len(ov["entries"]), + "uniprot_accessions": len(ov["uniprot"]), + "transitions": edges["total"], + "cross_state_transitions": edges["cross_state"], + }, + "states_per_sequence": capped(states_per_seq, 10), + "observations_per_sequence": [ + {"label": l, "value": v} for l, v in zip(obs_labels, hist(obs_per_seq, obs_edges)) + ], + "sequence_length": [ + {"label": l, "value": v} for l, v in zip(len_labels, hist(lengths, len_edges)) + ], + "resolution": [ + {"label": l, "value": v} for l, v in zip(res_labels, hist(ov["resolution"], res_edges)) + ], + "experimental_method": [{"label": k, "value": v} for k, v in ov["method"].most_common()], + "binding_status": [{"label": k, "value": v} for k, v in ov["binding"].most_common()], + "chain_composition": [{"label": k, "value": v} for k, v in ov["composition"].most_common()], + "nucleic_acid_binding": [{"label": k, "value": v} for k, v in ov["npp"].most_common()], + "release_year": [{"label": k, "value": ov["year"][k]} for k in sorted(ov["year"])], + "pair_fidelity": [ + {"label": k, "value": edges["pair_fidelity"].get(k, 0), + "cross_state": edges["pair_fidelity_cross"].get(k, 0)} + for k in fid_order if edges["pair_fidelity"].get(k) + ], + "state_fidelity": [ + {"label": k, **state_fid[k]} for k in fid_order if k in state_fid + ], + "pair_similarity": { + "bin_width": 0.02, "all": edges["similarity"], "cross_state": edges["similarity_cross"], + }, + "tm_aln": {"bin_width": 0.02, "all": edges["tm_aln"]}, + "rmsd": {"bin_width": 0.5, "all": edges["rmsd"], "last_bin_open": True}, + "coverage": {"bin_width": 0.05, "all": edges["coverage"]}, + } + + +def main(): + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("db", type=Path) + ap.add_argument("--out-dir", type=Path, required=True) + ap.add_argument("--limit", type=int, default=None, help="debug: only read the first N rows") + ap.add_argument("--ecod", type=Path, default=None, help="ecod.latest.domains.txt (for ECOD names)") + args = ap.parse_args() + + args.out_dir.mkdir(parents=True, exist_ok=True) + index_path = args.out_dir / "MuSProt-index.db" + tmp_path = index_path.with_suffix(".db.tmp") + tmp_path.unlink(missing_ok=True) + + src = sqlite3.connect(f"file:{args.db}?mode=ro", uri=True) + src.execute("PRAGMA temp_store = FILE") + src.execute("PRAGMA cache_size = -2000000") + dst = sqlite3.connect(tmp_path) + dst.execute("PRAGMA journal_mode = OFF") + dst.execute("PRAGMA synchronous = OFF") + + ecod_names = load_ecod_names(args.ecod) + log(f"ECOD names: {len(ecod_names[0]):,} families, {len(ecod_names[1]):,} homology groups") + + log("pass 1/4: node table") + seqs, functions, ov, n_nodes = build_members(src, dst, args.limit) + write_sequences(dst, seqs, functions, ecod_names) + dst.commit() + + log("pass 2/4: state pairs (edge GROUP BY)") + state_fid = build_state_pairs(src, dst, seqs, args.limit) + dst.commit() + + log("pass 3/4: edge histograms") + edges = edge_histograms(src, args.limit) + log("pass 4/4: order / disorder labels") + order = order_stats(src, args.limit) + + overview = build_overview(seqs, ov, n_nodes, edges, state_fid) + overview["outer_graph"] = outer_graph_stats(dst, len(seqs)) + overview["order"] = order + dst.execute("CREATE TABLE meta (key TEXT PRIMARY KEY, value TEXT)") + dst.execute("INSERT INTO meta VALUES ('overview', ?)", (json.dumps(overview),)) + dst.commit() + dst.execute("VACUUM") + dst.close() + tmp_path.replace(index_path) + (args.out_dir / "overview.json").write_text(json.dumps(overview, indent=1)) + log(f"done → {index_path} ({index_path.stat().st_size / 1e6:.1f} MB)") + log(json.dumps(overview["counts"], indent=1)) + + +if __name__ == "__main__": + main() diff --git a/backend/scripts/generate_node_lookup.py b/backend/scripts/generate_node_lookup.py deleted file mode 100644 index 0ed46e3fdc112c45c84bb254b89b7ca49be11621..0000000000000000000000000000000000000000 --- a/backend/scripts/generate_node_lookup.py +++ /dev/null @@ -1,68 +0,0 @@ -"""Generate an indexed read-only node lookup sidecar from MuSProt.db.""" -from __future__ import annotations - -import argparse -import sqlite3 -import zlib -from pathlib import Path - - -NODE_COLUMNS = ( - "pdb_id", "auth_asym_id", "binding_status", "sequence", "sequence_length", - "cath_id", "cath_superfamily", "Rosetta", "FoldX", "EvoEF2", "RW", "RW+", - "ranked_functions", "state_id", "experimental_method", "pH", "temp_K", - "chain_composition", "matching_scores", -) - - -def main() -> None: - parser = argparse.ArgumentParser() - parser.add_argument("source", type=Path) - parser.add_argument("-o", "--output", type=Path, default=Path("MuSProt-node-lookup.db")) - args = parser.parse_args() - - output = args.output.resolve() - output.unlink(missing_ok=True) - - conn = sqlite3.connect(output) - try: - conn.execute("PRAGMA journal_mode = OFF") - conn.execute("PRAGMA synchronous = OFF") - conn.execute( - 'CREATE TABLE node (pdb_id, auth_asym_id, binding_status, sequence BLOB,' - ' sequence_length, cath_id, cath_superfamily, Rosetta, FoldX, EvoEF2,' - ' RW, "RW+", ranked_functions BLOB, state_id, experimental_method, pH, temp_K,' - ' chain_composition, matching_scores BLOB)' - ) - source = sqlite3.connect(f"file:{args.source.resolve()}?mode=ro", uri=True) - try: - select_columns = ", ".join(f'"{column}"' for column in NODE_COLUMNS) - rows = source.execute(f"SELECT {select_columns} FROM node") - placeholders = ", ".join("?" for _ in NODE_COLUMNS) - batch = [] - for row in rows: - row = list(row) - for index in (3, 12, 18): # sequence, ranked_functions, matching_scores - row[index] = zlib.compress((row[index] or "").encode("utf-8"), 1) - batch.append(row) - if len(batch) == 1000: - conn.executemany(f"INSERT INTO node VALUES ({placeholders})", batch) - batch.clear() - if batch: - conn.executemany(f"INSERT INTO node VALUES ({placeholders})", batch) - finally: - source.close() - conn.execute( - "CREATE INDEX idx_node_chain " - "ON node(LOWER(pdb_id), UPPER(auth_asym_id))" - ) - conn.execute("ANALYZE") - conn.commit() - finally: - conn.close() - - print(f"Created {output} ({output.stat().st_size:,} bytes)") - - -if __name__ == "__main__": - main() diff --git a/backend/space_app.py b/backend/space_app.py index a9c0732694e88adf36a96db43e908cd4556aed5d..622789f69380c279cc426a27233a468a93ab006c 100644 --- a/backend/space_app.py +++ b/backend/space_app.py @@ -7,11 +7,13 @@ from contextlib import asynccontextmanager from pathlib import Path from fastapi import FastAPI, HTTPException +from fastapi.middleware.gzip import GZipMiddleware from fastapi.responses import FileResponse, PlainTextResponse, RedirectResponse from fastapi.staticfiles import StaticFiles from app.protein.api.routes.protein import router as protein_router from app.protein.api.routes.protein import set_data_manager +from app.protein.api.routes.explore import router as explore_router from app.protein.config import get_database_path, get_docs_path, get_plots_dir from app.protein.data_loader import DataManager from app.protein import tsv_loader @@ -52,7 +54,9 @@ app = FastAPI( redoc_url="/api/redoc", openapi_url="/api/openapi.json", ) +app.add_middleware(GZipMiddleware, minimum_size=2048) app.include_router(protein_router, prefix="/api/protein", tags=["MuSProt"]) +app.include_router(explore_router, prefix="/api/protein", tags=["MuSProt explorer"]) @app.get("/health") diff --git a/frontend/index.html b/frontend/index.html index 479ea02472fe14f7301a34b3e48db98990f277af..f1700741388593e1d4b684be91d98c96b826def7 100755 --- a/frontend/index.html +++ b/frontend/index.html @@ -1,16 +1,21 @@ - - - - - - - MuSProt: Multistate Protein Structure Database - - - - -
- - - - + + + + + + MuSProt — Multistate Protein Database + + + + + + + + +
+ + + diff --git a/frontend/package-lock.json b/frontend/package-lock.json index 97c7e07b261fca0ebb85fa7eb1d52a7c58c4a6e9..3516fdcbc6df5a86cb944e91efa03516511d6bf2 100755 --- a/frontend/package-lock.json +++ b/frontend/package-lock.json @@ -10,6 +10,7 @@ "license": "MIT", "dependencies": { "axios": "^1.13.4", + "d3-force": "^3.0.0", "lucide-react": "^0.487.0", "react": "^18.3.1", "react-dom": "^18.3.1", @@ -19,6 +20,7 @@ "remark-gfm": "^4.0.1" }, "devDependencies": { + "@types/d3-force": "^3.0.10", "@types/node": "^20.10.0", "@types/react": "^19.2.9", "@types/react-dom": "^19.2.3", @@ -1636,6 +1638,13 @@ "@swc/counter": "^0.1.3" } }, + "node_modules/@types/d3-force": { + "version": "3.0.10", + "resolved": "https://registry.npmjs.org/@types/d3-force/-/d3-force-3.0.10.tgz", + "integrity": "sha512-ZYeSaCF3p73RdOKcjj+swRlZfnYpK1EbaDiYICEEp5Q6sUiqFaFQ9qgoshp5CzIyyb/yD09kD9o2zEltCexlgw==", + "dev": true, + "license": "MIT" + }, "node_modules/@types/debug": { "version": "4.1.13", "resolved": "https://registry.npmjs.org/@types/debug/-/debug-4.1.13.tgz", @@ -2385,6 +2394,47 @@ "integrity": "sha512-z1HGKcYy2xA8AGQfwrn0PAy+PB7X/GSj3UVJW9qKyn43xWa+gl5nXmU4qqLMRzWVLFC8KusUX8T/0kCiOYpAIQ==", "license": "MIT" }, + "node_modules/d3-dispatch": { + "version": "3.0.1", + "resolved": "https://registry.npmjs.org/d3-dispatch/-/d3-dispatch-3.0.1.tgz", + "integrity": "sha512-rzUyPU/S7rwUflMyLc1ETDeBj0NRuHKKAcvukozwhshr6g6c5d8zh4c2gQjY2bZ0dXeGLWc1PF174P2tVvKhfg==", + "license": "ISC", + "engines": { + "node": ">=12" + } + }, + "node_modules/d3-force": { + "version": "3.0.0", + "resolved": "https://registry.npmjs.org/d3-force/-/d3-force-3.0.0.tgz", + "integrity": "sha512-zxV/SsA+U4yte8051P4ECydjD/S+qeYtnaIyAs9tgHCqfguma/aAQDjo85A9Z6EKhBirHRJHXIgJUlffT4wdLg==", + "license": "ISC", + "dependencies": { + "d3-dispatch": "1 - 3", + "d3-quadtree": "1 - 3", + "d3-timer": "1 - 3" + }, + "engines": { + "node": ">=12" + } + }, + "node_modules/d3-quadtree": { + "version": "3.0.1", + "resolved": "https://registry.npmjs.org/d3-quadtree/-/d3-quadtree-3.0.1.tgz", + "integrity": "sha512-04xDrxQTDTCFwP5H6hRhsRcb9xxv2RzkcsygFzmkSIOJy3PeRJP7sNk3VRIbKXcog561P9oU0/rVH6vDROAgUw==", + "license": "ISC", + "engines": { + "node": ">=12" + } + }, + "node_modules/d3-timer": { + "version": "3.0.1", + "resolved": "https://registry.npmjs.org/d3-timer/-/d3-timer-3.0.1.tgz", + "integrity": "sha512-ndfJ/JxxMd3nw31uyKoY2naivF+r29V+Lc0svZxe1JvvIRmi8hUsrMvdOwgS1o6uBHmiz91geQ0ylPP0aj1VUA==", + "license": "ISC", + "engines": { + "node": ">=12" + } + }, "node_modules/debug": { "version": "4.4.3", "resolved": "https://registry.npmjs.org/debug/-/debug-4.4.3.tgz", diff --git a/frontend/package.json b/frontend/package.json index ae0209ab43231e3756ba67f983d2ef74d543bebf..07c86aa550db61e62da28f21c423161b89cd0150 100755 --- a/frontend/package.json +++ b/frontend/package.json @@ -7,6 +7,7 @@ "homepage": "https://huggingface.co/spaces/omaib/MuSProt", "dependencies": { "axios": "^1.13.4", + "d3-force": "^3.0.0", "lucide-react": "^0.487.0", "react": "^18.3.1", "react-dom": "^18.3.1", @@ -16,6 +17,7 @@ "remark-gfm": "^4.0.1" }, "devDependencies": { + "@types/d3-force": "^3.0.10", "@types/node": "^20.10.0", "@types/react": "^19.2.9", "@types/react-dom": "^19.2.3", diff --git a/frontend/public/favicon.svg b/frontend/public/favicon.svg new file mode 100644 index 0000000000000000000000000000000000000000..4a31ecdae7c9d6b49486c08997a2b7d327b2b087 --- /dev/null +++ b/frontend/public/favicon.svg @@ -0,0 +1 @@ + diff --git a/frontend/public/sitemap.xml b/frontend/public/sitemap.xml index c18bd805cfeae34c86cf64a7248e57a7d66c2a7b..6ed322b104d7e77ee42890fe41476237abfc87cb 100644 --- a/frontend/public/sitemap.xml +++ b/frontend/public/sitemap.xml @@ -5,6 +5,16 @@ monthly 1.0 + + https://omaib-musprot.hf.space/browse + monthly + 0.9 + + + https://omaib-musprot.hf.space/statistics + monthly + 0.7 + https://omaib-musprot.hf.space/docs monthly diff --git a/frontend/scripts/prerender.mjs b/frontend/scripts/prerender.mjs index 8e81b07586f6edc924d691a8265da6b24bdd7559..ed8924ba14e7c10644d9d11b342086a32fb47707 100644 --- a/frontend/scripts/prerender.mjs +++ b/frontend/scripts/prerender.mjs @@ -8,14 +8,24 @@ const template = readFileSync(join(BUILD_DIR, 'index.html'), 'utf-8'); const routes = [ { path: '', - title: 'MuSProt: Multistate Protein Structure Database', + title: 'MuSProt — Multistate Protein Database', description: - 'Explore multistate protein structures, pairwise comparisons, energy scores, and functional annotations.', + 'Search and download 863,643 experimentally observed conformations of 102,790 protein sequences, clustered into 162,601 states with 113.7 million annotated state transitions.', + }, + { + path: 'browse', + title: 'Browse sequences — MuSProt', + description: 'Browse protein sequences with multiple experimentally observed conformations, sorted by number of conformational states.', + }, + { + path: 'statistics', + title: 'Statistics — MuSProt', + description: 'Distributions of conformational states, transitions, structural similarity, fidelity and experimental metadata across MuSProt.', }, { path: 'docs', - title: 'MuSProt Documentation', - description: 'MuSProt dataset format, SQLite schema, and usage examples.', + title: 'Documentation — MuSProt', + description: 'MuSProt SQLite schema, field definitions and usage examples.', }, ]; diff --git a/frontend/src/App.tsx b/frontend/src/App.tsx index 30bda2a200a13f9654b323c2fbbcd1e5f0b3c470..7a0b8408b06ea23d0294728ddbaa5904ea198e91 100644 --- a/frontend/src/App.tsx +++ b/frontend/src/App.tsx @@ -1,20 +1,62 @@ -import { Navigate, Route, Routes } from 'react-router-dom'; -import { Header } from './components/layout/Header'; -import { Footer } from './components/layout/Footer'; -import { ProteinDashboard } from './components/datasets/protein/ProteinDashboard'; -import { ProteinDocsPage } from './pages/datasets/ProteinDocsPage'; +import { useEffect } from 'react'; +import { Link, Navigate, Route, Routes, useLocation, useParams } from 'react-router-dom'; +import { Footer, Header } from './components/Layout'; +import { DatasetJsonLd, PageMeta } from './components/PageMeta'; +import { HomePage } from './pages/HomePage'; +import { BrowsePage } from './pages/BrowsePage'; +import { SearchPage } from './pages/SearchPage'; +import { SequencePage } from './pages/SequencePage'; +import { ChainPage } from './pages/ChainPage'; +import { StatisticsPage } from './pages/StatisticsPage'; +import { DocsPage } from './pages/DocsPage'; + +function ScrollToTop() { + const { pathname } = useLocation(); + useEffect(() => { + window.scrollTo(0, 0); + }, [pathname]); + return null; +} + +/** /chain/3vgk.A → /chain/3vgk/A */ +function ChainShorthand() { + const { id = '' } = useParams(); + const m = /^([0-9][A-Za-z0-9]{3})[.:_-](.+)$/.exec(id); + return m ? : ; +} + +function NotFound() { + return ( +
+ +

404

+

This page is in an unobserved state.

+

+ Try the sequence browser or search from the header. +

+
+ ); +} export default function App() { return ( -
+
+ +
-
+
- } /> - } /> + } /> + } /> + } /> + } /> + } /> + } /> + } /> + } /> } /> } /> - } /> + } />