Blame for engines/molgate_module4_engineA.py

"""
MolGate — Module 4: Pocket Definition (Engine A)
Stage 21: HETATM scan & filter

輸入：protein/cleaned.pdb（Module 3 Engine A）、prep_decisions.json（須 Module 3 Engine B 完成並
      ready_for_engine_e=true，與 Module 3 Engine C 前置條件一致）
輸出：stage_21_hetatm.json（候選配體列表；含 ligand_role=target/candidate）
規則：Master Index 之 het_id／cocrystal_het 白名單優先；糖基化／輔因子黑名單（輔因子可於 Index
      設 keep_cofactors 保留）；其餘套用緩衝離子黑名單與 MW 門檻。
依賴：biopython, rdkit；讀 Index 時與 molgate_module3_engineA 共用載入邏輯。
"""

from __future__ import annotations

import argparse
import json
import sys
from datetime import datetime
from pathlib import Path
from typing import Any

from Bio.PDB import PDBParser

_PKG_ROOT = Path(__file__).resolve().parent
_UI_ROOT = _PKG_ROOT.parent
if str(_PKG_ROOT) not in sys.path:
    sys.path.insert(0, str(_PKG_ROOT))
if str(_UI_ROOT) not in sys.path:
    sys.path.insert(0, str(_UI_ROOT))

from molgate_ui.console_utf8 import configure_console_utf8

configure_console_utf8()

SCHEMA_VERSION = "0.2"

# ── 黑名單：緩衝劑、離子、溶劑 ───────────────────────
BLACKLIST = {
    "HOH","WAT","GOL","EDO","PEG","MPD","DMS","FMT",
    "ACT","TRS","MES","EPE","BME","DTT","TCE","IMD",
    "SO4","PO4","CIT","ACY","EOH","IPA","MOH","BOG",
    "CL","NA","MG","ZN","CA","K","FE","MN","CU","NI",
    "CO","CD","HG","PT","AU","AG","SR","BA","CS","RB",
    "IOD","BR","F","LI","AL","SE",
}

# ── 糖基化／寡糖常見殘基（非口袋藥物配體）────────────────
GLYCAN_BLACKLIST = {
    "NAG", "NGA", "BMA", "MAN", "GAL", "GLC", "GCS", "MAL", "MAB", "MBG", "BGC",
    "FUC", "FUL", "RAM", "LFR", "GLO", "XYP", "XYS", "XYF", "XUL",
    "SIA", "SLB", "A2G", "M6D", "MMA", "MGL", "NGK", "IDR",
}

# ── 常見輔因子／核苷酸／金屬簇（非口袋配體；除非 Index 標記 keep_cofactors）──
COFACTOR_BLACKLIST = {
    "HEM", "HEA", "HEB", "HEC", "HAS", "CHL", "BCL", "CLA", "FAD", "FMN",
    "NAD", "NAP", "NAI", "NAO", "NAM", "NDP", "NDG", "NT2", "NAJ", "NAC",
    "ATP", "ADP", "AMP", "ANP", "GTP", "GDP", "GMP", "GNP",
    "UTP", "UDP", "UMP", "CTP", "CDP", "CMP", "TTP", "TDP", "TMP",
    "COA", "SAM", "SAH", "ACO", "SCA",
    "SF4", "FES", "F3S", "CLF", "CFM", "FSO", "FS4",
    "PLP", "PLH", "PMP", "PSP", "MGD", "GTN", "Q10", "U10", "3CO", "CUA", "CU1",
}

MW_CUTOFF = 150.0  # Da 以下排除（Master Index 標記之目標配體不受此限）

# ── 工具函式 ──────────────────────────────────────────
def save_stage(session_dir: Path, stage: int, name: str, data: dict[str, Any]) -> dict[str, Any]:
    data = dict(data)
    data["_stage"] = stage
    data["_name"] = name
    data["_timestamp"] = datetime.now().isoformat()
    data["_schema_version"] = SCHEMA_VERSION
    path = session_dir / f"stage_{stage:02d}_{name}.json"
    path.write_text(json.dumps(data, indent=2, ensure_ascii=False), encoding="utf-8")
    return data

def load_prep_decisions(session_dir: Path) -> dict[str, Any]:
    path = session_dir / "prep_decisions.json"
    if not path.is_file():
        raise FileNotFoundError(
            "找不到 prep_decisions.json，請先完成 Module 3 Engine B（人機決策 → ready_for_engine_e）"
        )
    return json.loads(path.read_text(encoding="utf-8"))


def validate_module3_prerequisites(
    session_dir: Path,
    target: str | None,
    *,
    skip_ready_check: bool = False,
    skip_manifest_target_check: bool = False,
) -> dict[str, Any]:
    """
    與 molgate_module3_engineC 對齊：須 prep_decisions、ready_for_engine_e；
    若提供 --target 則與 manifest / prep_decisions.target 交叉驗證。
    """
    if target and not skip_manifest_target_check:
        try:
            from molgate_manifest import validate_target_matches_manifest

            validate_target_matches_manifest(session_dir, target)
        except FileNotFoundError:
            pass
        except ValueError as e:
            print(f"❌ {e}")
            sys.exit(1)

    try:
        prep = load_prep_decisions(session_dir)
    except FileNotFoundError as e:
        print(f"❌ {e}")
        sys.exit(1)

    if not skip_ready_check and not prep.get("ready_for_engine_e"):
        print(
            "❌ prep_decisions.json 尚未標記 ready_for_engine_e=true，請先完成 Module 3 Engine B"
        )
        sys.exit(1)

    if target:
        pt = (prep.get("target") or "").strip().upper()
        if pt and pt != target.strip().upper():
            print(
                f"❌ prep_decisions.target={prep.get('target')!r} 與 --target={target!r} 不一致"
            )
            sys.exit(1)

    pdb_id = (prep.get("pdb_id") or "").strip().upper()
    if not pdb_id:
        print("❌ prep_decisions 缺少 pdb_id（請先完成 Module 3 Engine B）")
        sys.exit(1)

    return prep


def load_cleaned_pdb(session_dir: Path) -> Path:
    p = session_dir / "protein" / "cleaned.pdb"
    if not p.is_file():
        raise FileNotFoundError(
            "找不到 protein/cleaned.pdb，請先完成 Module 3 Engine A（產出 cleaned.pdb）"
        )
    return p


def resolve_index_path(session_dir: Path, cli_index: Path | None) -> Path:
    if cli_index is not None:
        return cli_index
    meta = session_dir / "session_meta.json"
    if meta.is_file():
        try:
            data = json.loads(meta.read_text(encoding="utf-8"))
            p = data.get("index_path")
            if isinstance(p, str) and p.strip():
                return Path(p)
        except (json.JSONDecodeError, OSError):
            pass
    mp = session_dir / "molgate_run_manifest.json"
    if mp.is_file():
        try:
            data = json.loads(mp.read_text(encoding="utf-8"))
            p = data.get("index_path")
            if isinstance(p, str) and p.strip():
                return Path(p)
        except (json.JSONDecodeError, OSError):
            pass
    return _PKG_ROOT / "master_index.json"


def load_target_index_entry(
    session_dir: Path,
    prep: dict[str, Any],
    index_path: Path | None,
) -> dict[str, Any]:
    """Master Index 條目（與 Module 3 對齊）；失敗時以 prep 的 cocrystal_het 等為後備。"""
    target = (prep.get("target") or "").strip()
    if not target:
        return {}
    resolved = resolve_index_path(session_dir, index_path)
    try:
        from molgate_module3_engineA import build_engine_entry, load_master_index_file

        master = load_master_index_file(resolved)
        return build_engine_entry(master, target, session_dir)
    except (FileNotFoundError, KeyError, OSError):
        out: dict[str, Any] = {}
        for k in ("cocrystal_het", "het_id"):
            v = prep.get(k)
            if isinstance(v, str) and v.strip():
                out[k] = v.strip()
        ov = prep.get("override")
        if isinstance(ov, dict):
            out["override"] = ov
        return out


def get_primary_het_ids(entry: dict[str, Any], prep: dict[str, Any]) -> set[str]:
    """Master Index / prep 標記之共晶配體 ID（het_id 或 cocrystal_het）。"""
    s: set[str] = set()
    for src in (entry, prep):
        if not isinstance(src, dict):
            continue
        for key in ("het_id", "cocrystal_het"):
            v = src.get(key)
            if isinstance(v, str) and v.strip():
                s.add(v.strip().upper())
            elif isinstance(v, list):
                for x in v:
                    if isinstance(x, str) and x.strip():
                        s.add(x.strip().upper())
    return s


def get_keep_cofactors(entry: dict[str, Any]) -> set[str]:
    """除非列於此，否則套用 COFACTOR_BLACKLIST。可寫在條目或 override.keep_cofactors。"""
    s: set[str] = set()
    s |= get_primary_het_ids(entry, {})
    for key in ("keep_cofactors", "retain_cofactors", "cofactor_whitelist"):
        v = entry.get(key)
        if isinstance(v, list):
            s |= {x.strip().upper() for x in v if isinstance(x, str) and x.strip()}
    ov = entry.get("override")
    if isinstance(ov, dict):
        for key in ("keep_cofactors", "retain_cofactors", "cofactor_whitelist"):
            v = ov.get(key)
            if isinstance(v, list):
                s |= {x.strip().upper() for x in v if isinstance(x, str) and x.strip()}
    return s


def sort_stage21_candidates(cands: list[dict[str, Any]]) -> list[dict[str, Any]]:
    targets = [c for c in cands if c.get("ligand_role") == "target"]
    rest = [c for c in cands if c.get("ligand_role") != "target"]
    targets.sort(key=lambda c: (c["chain"], c["resseq"]))
    rest.sort(key=lambda c: (-float(c["mw_estimate"]), c["chain"], c["resseq"]))
    return targets + rest


def get_het_atoms(structure, resname: str, chain_id: str, resseq: int) -> list:
    atoms = []
    for model in structure:
        for chain in model:
            if chain.id != chain_id:
                continue
            for res in chain:
                if (res.get_resname().strip() == resname and
                        res.get_id()[1] == resseq):
                    atoms.extend(res.get_atoms())
    return atoms

def estimate_mw_from_pdb(resname: str) -> float:
    """
    從 PDB HETATM 估算 MW。
    優先用 RDKit CCD 查詢，失敗則用原子數粗估。
    """
    try:
        from rdkit.Chem import AllChem, Descriptors
        from rdkit.Chem.rdmolfiles import MolFromPDBBlock
        # 嘗試用 CCD SMILES 查詢（需要網路或本地 CCD）
        # 這裡用原子計數粗估作為 fallback
    except ImportError:
        pass

    # 粗估：平均有機分子每原子約 13 Da
    # 實際上我們用 PDB 結構裡的原子數來估
    return None  # 由呼叫端用原子數估算

# ── Stage 21：HETATM 掃描過濾 ─────────────────────────
def stage21_hetatm(
    session_dir: Path,
    prep: dict[str, Any],
    index_path: Path | None,
) -> dict[str, Any]:
    print("\n[Stage 21] HETATM 掃描過濾...")

    idx_entry = load_target_index_entry(session_dir, prep, index_path)
    primary_ids = get_primary_het_ids(idx_entry, prep)
    keep_cof = get_keep_cofactors(idx_entry)

    print(f"  → Master Index 目標 het_id / cocrystal_het：{sorted(primary_ids) or '（無）'}")
    if keep_cof:
        print(f"  → 輔因子保留（keep_cofactors）：{sorted(keep_cof)}")

    cleaned_path = load_cleaned_pdb(session_dir)
    parser = PDBParser(QUIET=True)
    structure = parser.get_structure("protein", str(cleaned_path))

    # 掃描所有 HETATM 殘基
    raw_het = []
    seen = set()

    for model in structure:
        for chain in model:
            for res in chain:
                hetflag, resseq, _ = res.get_id()
                resname = res.get_resname().strip().upper()

                # 只處理 HETATM（hetflag 非空白）
                if not hetflag.strip():
                    continue

                key = (resname, chain.id, resseq)
                if key in seen:
                    continue
                seen.add(key)

                atoms = list(res.get_atoms())
                atom_count = len(atoms)

                # 用原子數估算 MW（粗估，平均有機分子每非氫原子約 13 Da）
                # 更準確的做法是查 CCD，這裡用保守估算
                heavy_atoms = [a for a in atoms if a.element and a.element != 'H']
                mw_estimate = len(heavy_atoms) * 13.0 if heavy_atoms else atom_count * 10.0

                raw_het.append({
                    "resname": resname,
                    "chain": chain.id,
                    "resseq": resseq,
                    "atom_count": atom_count,
                    "heavy_atom_count": len(heavy_atoms),
                    "mw_estimate": round(mw_estimate, 1),
                })

    print(f"  ✓ 掃描到 HETATM 殘基：{len(raw_het)} 個")

    # 過濾：白名單（Index het）優先 → 糖 → 輔因子 → 緩衝/離子 → MW
    excluded_blacklist: list[str] = []
    excluded_glycan: list[str] = []
    excluded_cofactor: list[str] = []
    excluded_mw: list[str] = []
    candidates: list[dict[str, Any]] = []

    for het in raw_het:
        resname = het["resname"]
        mw = float(het["mw_estimate"])
        rsn = resname.upper()

        if rsn in primary_ids:
            row = dict(het)
            row["ligand_role"] = "target"
            candidates.append(row)
            continue

        if rsn in GLYCAN_BLACKLIST:
            excluded_glycan.append(rsn)
            continue

        if rsn in COFACTOR_BLACKLIST and rsn not in keep_cof:
            excluded_cofactor.append(rsn)
            continue

        if rsn in BLACKLIST:
            excluded_blacklist.append(rsn)
            continue

        if mw < MW_CUTOFF:
            excluded_mw.append(f"{resname}({mw:.0f}Da)")
            continue

        row = dict(het)
        row["ligand_role"] = "candidate"
        candidates.append(row)

    candidates = sort_stage21_candidates(candidates)

    print(f"  ✓ 目標配體（Index 白名單）：{sum(1 for c in candidates if c.get('ligand_role') == 'target')} 個")
    print(f"  ✓ 糖基化殘基排除：{len(excluded_glycan)} 個")
    print(f"  ✓ 輔因子排除：{len(excluded_cofactor)} 個")
    print(f"  ✓ 緩衝/離子/溶劑黑名單排除：{len(excluded_blacklist)} 個")
    print(f"  ✓ MW < {MW_CUTOFF} Da 排除：{len(excluded_mw)} 個")
    print(f"  ✓ 候選配體（目標 + 其餘）：{len(candidates)} 個")

    for c in candidates:
        tag = "【目標】" if c.get("ligand_role") == "target" else ""
        print(
            f"    {tag}{c['resname']} chain {c['chain']} res {c['resseq']} "
            f"(~{c['mw_estimate']:.0f} Da, {c['heavy_atom_count']} heavy atoms)"
        )

    if not candidates:
        status = "WARN"
        print(f"  ⚠ WARN: 無候選配體，可能是 apo 結構或 HET 已被移除")
    else:
        status = "PASS"

    print(f"  → Status: {status}")

    result = {
        "status": status,
        "raw_hetatm_count": len(raw_het),
        "excluded_glycan_count": len(excluded_glycan),
        "excluded_cofactor_count": len(excluded_cofactor),
        "excluded_blacklist_count": len(excluded_blacklist),
        "excluded_mw_count": len(excluded_mw),
        "candidate_count": len(candidates),
        "candidates": candidates,
        "mw_cutoff": MW_CUTOFF,
        "filter_rules": {
            "primary_het_ids": sorted(primary_ids),
            "keep_cofactors": sorted(keep_cof),
            "glycan_blacklist_size": len(GLYCAN_BLACKLIST),
            "cofactor_blacklist_size": len(COFACTOR_BLACKLIST),
        },
        "module3_prep": {
            "pdb_id": (prep.get("pdb_id") or "").strip().upper(),
            "target": prep.get("target"),
            "ready_for_engine_e": bool(prep.get("ready_for_engine_e")),
        },
    }
    return save_stage(session_dir, 21, "hetatm", result)

# ── 主流程 ────────────────────────────────────────────
def run_engine_a(
    session_dir: Path,
    *,
    target: str | None = None,
    index_path: Path | None = None,
    skip_ready_check: bool = False,
    skip_manifest_target_check: bool = False,
):
    print("=" * 55)
    print("  MolGate — Module 4 Engine A")
    print("  Stage 21: HETATM 掃描過濾")
    print("=" * 55)

    prep = validate_module3_prerequisites(
        session_dir,
        target,
        skip_ready_check=skip_ready_check,
        skip_manifest_target_check=skip_manifest_target_check,
    )
    pdb_id = (prep.get("pdb_id") or "").strip().upper()
    tgt = (prep.get("target") or target or "").strip().upper()
    print(f"\n  Module 3 銜接：PDB={pdb_id}" + (f" | Target={tgt}" if tgt else ""))

    try:
        result = stage21_hetatm(session_dir, prep, index_path)
    except FileNotFoundError as e:
        print(f"❌ {e}")
        sys.exit(1)
    except Exception as e:
        print(f"❌ Pipeline stopped: {e}")
        sys.exit(1)

    print("\n" + "=" * 55)
    if result["status"] == "WARN":
        print("  Engine A 完成（WARN — 無候選配體）")
        print("  → 若為 apo 結構，請在 Master Index 標記 pocket_rule: manual_coords")
    else:
        print(f"  Engine A 完成 ✓ — {result['candidate_count']} 個候選配體")
    print(f"  輸出：{session_dir / 'stage_21_hetatm.json'}")
    print("=" * 55)

    return result

# ── CLI ──────────────────────────────────────────────
if __name__ == "__main__":
    ap = argparse.ArgumentParser(description="MolGate Module 4 Engine A")
    ap.add_argument("--session-dir", required=True, help="Session 目錄路徑")
    ap.add_argument(
        "--target",
        default=None,
        help="標靶 ID（建議與 Module 3 一致；若提供則與 prep_decisions / manifest 交叉驗證）",
    )
    ap.add_argument(
        "--index",
        default=None,
        help="master_index.json 路徑（預設：session_meta / manifest 或 MolGate 目錄下）",
    )
    ap.add_argument(
        "--skip-ready-check",
        action="store_true",
        help="不檢查 ready_for_engine_e（僅除錯用）",
    )
    ap.add_argument(
        "--skip-manifest-target-check",
        action="store_true",
        help="不與 molgate_run_manifest.json 的 target 交叉驗證",
    )
    args = ap.parse_args()

    sd = Path(args.session_dir)
    if not sd.is_dir():
        print(f"❌ Session 目錄不存在：{sd}")
        sys.exit(1)

    idx = Path(args.index) if args.index else None
    run_engine_a(
        sd,
        target=args.target,
        index_path=idx,
        skip_ready_check=args.skip_ready_check,
        skip_manifest_target_check=args.skip_manifest_target_check,
    )