CupOFMoCA-PROTAC/src/protac/KRASVHL.py

"""
KRAS-VHL PROTAC Builder
Target:    KRAS (aminocyanothiophene-tetrahydroquinoline-isoxazole-pyrimidine
           warhead, ACBI3 family)
E3 Ligase: VHL (VH032 hydroxyproline/thiazolyl-benzyl binder)
Reference: ACBI3 (Astex/Boehringer Ingelheim); DegradMaster case_study_2

Warhead (ACBI3-derived KRAS binder, MW ≈ 478 Da):
    Aminocyanothiophene fused to a tetrahydroquinoline; the thiophene-fused
    ring is decorated with an isoxazole–pyrimidine. The pyrimidine bears a
    2-methyl-1,4-diazepane whose distal N-H is the linker attachment point
    (analogous to the gefitinib piperazine NH in EGFR-VHL).

E3 Ligand (VH032, MW ≈ 431 Da):
    Hydroxyproline + tert-leucine + 4-(4-methylthiazol-5-yl)benzyl cap.
    Free NH2 on tert-leucine alpha-amine is the reactive group.

NOTE on chemistry: ACBI3 itself uses a click-derived 1,2,3-triazole linker
joining the warhead alkyl chain to a VHL-side azide-bearing fragment. For
generation/scoring purposes we mirror EGFR-VHL chemistry (diacid linker +
2× amide bond formation), since the diazepane NH and VH032 NH2 are both
direct amide-coupling sites and this preserves a single linker chemistry
across all three VHL targets in the study.
"""

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import RDLogger
from typing import Optional

try:
    from .protac_builder import PROTAC_Builder
except ImportError:
    from protac_builder import PROTAC_Builder

RDLogger.DisableLog("rdApp.error")


# ---------------------------------------------------------------------------
# Component SMILES — derived from ACBI3 (DegradMaster case_study_2, c0 reference)
# ---------------------------------------------------------------------------

# ACBI3-derived KRAS warhead with diazepane-NH (linker stripped)
KRAS_WARHEAD_SMILES = (
    "C[C@H]1CNCCCN1c1nccc(-c2noc([C@@]3(C)CCCc4sc(N)c(C#N)c43)n2)n1"
)

# VH032: same VHL ligand used by EGFR-VHL and BRD9-VHL modules
VHL_E3_SMILES = "N[C@H](C(=O)N1C[C@H](O)C[C@H]1C(=O)NCc1ccc(-c2scnc2C)cc1)C(C)(C)C"

# Pre-built Mol objects
KRAS_WARHEAD_MOL = Chem.MolFromSmiles(KRAS_WARHEAD_SMILES)
VHL_E3_MOL       = Chem.MolFromSmiles(VHL_E3_SMILES)

# SMARTS patterns
_COOH_SMARTS           = "[CX3](=[OX1])[OX2H1]"
_NH2_SMARTS            = "[NX3;H2]"
_NH_SMARTS             = "[NX3;H1;!$(NC=O);!$(NS=O);!$(NC#N)]"
_AMIDE_SMARTS          = "[NX3H0,NX3H1][CX3](=[OX1])[#6]"
# Isoxazole–pyrimidine biaryl (signature of the ACBI3 KRAS warhead)
_ISOXAZOLE_PYRIMIDINE_SMARTS = "c1onc(n1)-c1ccnc(n1)"
_HYDROXYPROLINE_SMARTS = "[C@@H]1C[C@H](O)CN1"
_TERMINAL_CHO_SMARTS   = "[CX3H1](=[OX1])"
_TERMINAL_CH3_SMARTS   = "[CX4H3]"


class KRASVHL_PROTAC_Builder(PROTAC_Builder):
    """
    PROTAC builder: ACBI3-derived KRAS warhead + VH032 VHL E3.

    Linker requirement: two COOH groups (diacid), one per end.
    Warhead NH (diazepane) + linker COOH → tertiary amide.
    E3 NH2 (tert-leucine) + linker COOH → secondary amide.
    """

    def __init__(self, warhead: Optional[Chem.Mol] = None,
                 e3: Optional[Chem.Mol] = None):
        super().__init__(
            warhead=warhead if warhead is not None else KRAS_WARHEAD_MOL,
            e3=e3 if e3 is not None else VHL_E3_MOL,
            warhead_name="ACBI3-diazepane (KRAS)",
            e3_name="VHL binder NH2 (VH032)",
        )
        self._cooh_pat       = Chem.MolFromSmarts(_COOH_SMARTS)
        self._nh2_pat        = Chem.MolFromSmarts(_NH2_SMARTS)
        self._nh_pat         = Chem.MolFromSmarts(_NH_SMARTS)
        self._amide_pat      = Chem.MolFromSmarts(_AMIDE_SMARTS)
        self._iso_pyrim      = Chem.MolFromSmarts(_ISOXAZOLE_PYRIMIDINE_SMARTS)
        self._hyp            = Chem.MolFromSmarts(_HYDROXYPROLINE_SMARTS)
        self._terminal_cho   = Chem.MolFromSmarts(_TERMINAL_CHO_SMARTS)
        self._terminal_ch3   = Chem.MolFromSmarts(_TERMINAL_CH3_SMARTS)

        # --- Amide bond reactions (diacid linker path) ---
        self._amide_rxn_nh2 = AllChem.ReactionFromSmarts(
            "[NX3H2:1].[CX3:2](=[OX1:3])[OX2H1]>>[NX3H1:1][CX3:2](=[OX1:3])"
        )
        self._amide_rxn_nh = AllChem.ReactionFromSmarts(
            "[NX3H1;!$(NC=O);!$(NS=O):1].[CX3:2](=[OX1:3])[OX2H1]>>[NX3H0:1][CX3:2](=[OX1:3])"
        )

        # --- C-N alkylation reactions (alkyl/PEG linker path) ---
        self._cn_rxn_nh2 = AllChem.ReactionFromSmarts(
            "[NX3H2:1].[CX4H3:2]>>[NX3H1:1][CX4H2:2]"
        )
        self._cn_rxn_nh = AllChem.ReactionFromSmarts(
            "[NX3H1;!$(NC=O);!$(NS=O):1].[CX4H3:2]>>[NX3H0:1][CX4H2:2]"
        )

    # ------------------------------------------------------------------
    # Verification
    # ------------------------------------------------------------------

    def _verify_warhead_reactive_group(self):
        pat = Chem.MolFromSmarts(_NH_SMARTS)
        if self.warhead_mol.HasSubstructMatch(pat):
            print(f"✓ {self.warhead_name}: NH reactive group confirmed (diazepane)")
        else:
            print(f"⚠ {self.warhead_name}: no reactive NH detected — check warhead")

    def _verify_e3_reactive_group(self):
        pat = Chem.MolFromSmarts(_NH2_SMARTS)
        if self.e3_mol.HasSubstructMatch(pat):
            print(f"✓ {self.e3_name}: NH2 reactive group confirmed (tert-leucine)")
        else:
            print(f"⚠ {self.e3_name}: no reactive NH2 detected — check E3")

    # ------------------------------------------------------------------
    # Linker attachment site detection
    # ------------------------------------------------------------------

    def _reactive_cooh_atoms(self, linker_mol: Chem.Mol) -> list:
        results = []
        seen = set()

        for m in linker_mol.GetSubstructMatches(self._cooh_pat):
            c_idx = m[0]
            if c_idx not in seen:
                seen.add(c_idx)
                results.append({"idx": c_idx, "restored": False})

        for m in linker_mol.GetSubstructMatches(self._terminal_cho):
            c_idx = m[0]
            if c_idx in seen:
                continue

            rw = Chem.RWMol(linker_mol)
            c_atom = rw.GetAtomWithIdx(c_idx)
            c_atom.SetNumExplicitHs(0)
            c_atom.SetNoImplicit(True)
            o_idx = rw.AddAtom(Chem.Atom("O"))
            rw.GetAtomWithIdx(o_idx).SetNumExplicitHs(1)
            rw.AddBond(c_idx, o_idx, Chem.BondType.SINGLE)
            rw.GetAtomWithIdx(c_idx).SetAtomMapNum(99)
            try:
                Chem.SanitizeMol(rw)
            except Exception:
                continue

            step1_rxn = AllChem.ReactionFromSmarts(
                "[NX3H2:1].[CX3:99](=[OX1:3])[OX2H1]>>[NX3H1:1][CX3:99](=[OX1:3])"
            )
            if step1_rxn.RunReactants((self.e3_mol, rw.GetMol())):
                seen.add(c_idx)
                results.append({"idx": c_idx, "restored": True})

        return results

    def _reactive_terminal_c_atoms(self, linker_mol: Chem.Mol) -> list:
        results = []
        seen = set()

        for m in linker_mol.GetSubstructMatches(self._terminal_ch3):
            c_idx = m[0]
            if c_idx in seen:
                continue

            rw = Chem.RWMol(linker_mol)
            rw.GetAtomWithIdx(c_idx).SetAtomMapNum(98)
            try:
                Chem.SanitizeMol(rw)
            except Exception:
                continue

            trial_rxn = AllChem.ReactionFromSmarts(
                "[NX3H2:1].[CX4H3:98]>>[NX3H1:1][CX4H2:98]"
            )
            if trial_rxn.RunReactants((self.e3_mol, rw.GetMol())):
                seen.add(c_idx)
                results.append({"idx": c_idx, "type": "alkyl"})

        return results

    def find_linker_attachment_sites(self, linker_mol: Chem.Mol,
                                      amide_only: bool = False) -> list:
        amide_cands  = self._reactive_cooh_atoms(linker_mol)
        alkyl_cands  = [] if amide_only else self._reactive_terminal_c_atoms(linker_mol)

        sites = []
        seen  = set()

        def _add(w, e, wtype, etype):
            key = (w["idx"], e["idx"])
            if key not in seen and w["idx"] != e["idx"]:
                seen.add(key)
                sites.append({
                    "warhead_site":     w["idx"],
                    "warhead_type":     wtype,
                    "warhead_restored": w.get("restored", False),
                    "e3_site":          e["idx"],
                    "e3_type":          etype,
                    "e3_restored":      e.get("restored", False),
                })

        direct   = [c for c in amide_cands if not c["restored"]]
        restored = [c for c in amide_cands if     c["restored"]]
        for w in direct:
            for e in direct:    _add(w, e, "cooh", "cooh")
        for w in direct:
            for e in restored:  _add(w, e, "cooh", "cooh")
            for e in direct:    _add(e, w, "cooh", "cooh")
        for w in restored:
            for e in restored:  _add(w, e, "cooh", "cooh")

        if not amide_only:
            for w in alkyl_cands:
                for e in alkyl_cands:
                    _add(w, e, "alkyl", "alkyl")

            for w in amide_cands:
                for e in alkyl_cands:
                    _add(w, e, "cooh", "alkyl")
                    _add(e, w, "alkyl", "cooh")

        return sites

    # ------------------------------------------------------------------
    # Assembly
    # ------------------------------------------------------------------

    def _restore_cho_to_cooh(self, rw: Chem.RWMol, c_idx: int):
        c_atom = rw.GetAtomWithIdx(c_idx)
        c_atom.SetNumExplicitHs(0)
        c_atom.SetNoImplicit(True)
        o_idx = rw.AddAtom(Chem.Atom("O"))
        rw.GetAtomWithIdx(o_idx).SetNumExplicitHs(1)
        rw.AddBond(c_idx, o_idx, Chem.BondType.SINGLE)

    def _assemble_protac(self, linker_mol: Chem.Mol, warhead_site_idx: int,
                         e3_site_idx: int,
                         warhead_type: str = "cooh", e3_type: str = "cooh",
                         warhead_restored: bool = False,
                         e3_restored: bool = False) -> Optional[str]:
        rw = Chem.RWMol(linker_mol)

        if e3_type == "cooh" and e3_restored:
            self._restore_cho_to_cooh(rw, e3_site_idx)
        if warhead_type == "cooh" and warhead_restored:
            self._restore_cho_to_cooh(rw, warhead_site_idx)
        try:
            Chem.SanitizeMol(rw)
        except Exception:
            return None

        rw.GetAtomWithIdx(e3_site_idx).SetAtomMapNum(99)
        tagged_linker = rw.GetMol()

        if e3_type == "cooh":
            step1_rxn = AllChem.ReactionFromSmarts(
                "[NX3H2:1].[CX3:99](=[OX1:3])[OX2H1]>>[NX3H1:1][CX3:99](=[OX1:3])"
            )
        else:
            step1_rxn = AllChem.ReactionFromSmarts(
                "[NX3H2:1].[CX4H3:99]>>[NX3H1:1][CX4H2:99]"
            )

        products1 = step1_rxn.RunReactants((self.e3_mol, tagged_linker))
        if not products1:
            return None

        intermediate = None
        for prod_tuple in products1:
            try:
                mol = prod_tuple[0]
                Chem.SanitizeMol(mol)
                intermediate = mol
                break
            except Exception:
                continue
        if intermediate is None:
            return None

        rw2 = Chem.RWMol(intermediate)
        for atom in rw2.GetAtoms():
            if atom.GetAtomMapNum() == 99:
                atom.SetAtomMapNum(0)
        intermediate = rw2.GetMol()

        if warhead_type == "cooh":
            step2_rxn = self._amide_rxn_nh
        else:
            step2_rxn = self._cn_rxn_nh

        products2 = step2_rxn.RunReactants((self.warhead_mol, intermediate))
        if not products2:
            return None

        for prod_tuple in products2:
            try:
                protac_mol = prod_tuple[0]
                Chem.SanitizeMol(protac_mol)
                return Chem.MolToSmiles(protac_mol)
            except Exception:
                continue

        return None

    # ------------------------------------------------------------------
    # Validation
    # ------------------------------------------------------------------

    def check_attachment_site_compatibility(self, protac_mol: Chem.Mol) -> bool:
        has_amide  = len(protac_mol.GetSubstructMatches(self._amide_pat)) >= 2
        no_cooh    = len(protac_mol.GetSubstructMatches(self._cooh_pat)) == 0
        has_cores  = self.has_warhead_core(protac_mol) and self.has_e3_core(protac_mol)
        return (has_amide and no_cooh) or has_cores

    def has_warhead_core(self, protac_mol: Chem.Mol) -> bool:
        return protac_mol.HasSubstructMatch(self._iso_pyrim)

    def has_e3_core(self, protac_mol: Chem.Mol) -> bool:
        return protac_mol.HasSubstructMatch(self._hyp)

    def has_both_amide_bonds(self, protac_mol: Chem.Mol) -> bool:
        return len(protac_mol.GetSubstructMatches(self._amide_pat)) >= 2

    def no_unreacted_cooh(self, protac_mol: Chem.Mol) -> bool:
        return len(protac_mol.GetSubstructMatches(self._cooh_pat)) == 0

    def _basic_validation(self, protac_mol: Chem.Mol) -> dict:
        checks = {
            "has_warhead_core":      self.has_warhead_core(protac_mol),
            "has_e3_core":           self.has_e3_core(protac_mol),
            "has_both_amide_bonds":  self.has_both_amide_bonds(protac_mol),
            "no_unreacted_cooh":     self.no_unreacted_cooh(protac_mol),
            "attachment_site_valid": self.check_attachment_site_compatibility(protac_mol),
        }
        checks["_passed"] = all([
            checks["has_warhead_core"],
            checks["has_e3_core"],
            checks["has_both_amide_bonds"],
            checks["no_unreacted_cooh"],
        ])
        return checks


# ---------------------------------------------------------------------------
# Standard interface aliases (used by defog_guided.py / batch_report.py)
# ---------------------------------------------------------------------------

WARHEAD_MOL        = KRAS_WARHEAD_MOL
E3_MOL             = VHL_E3_MOL
WARHEAD_SMILES     = KRAS_WARHEAD_SMILES
E3_SMILES          = VHL_E3_SMILES
Builder            = KRASVHL_PROTAC_Builder
WH_REACTIVE_SMARTS = _NH_SMARTS
E3_REACTIVE_SMARTS = _NH2_SMARTS
WARHEAD_PAT        = Chem.MolFromSmarts(_ISOXAZOLE_PYRIMIDINE_SMARTS)  # KRAS warhead signature
E3_PAT             = Chem.MolFromSmarts("N1CC(O)CC1")                   # hydroxyproline (VHL)


# ---------------------------------------------------------------------------
# Self-test
# ---------------------------------------------------------------------------

if __name__ == "__main__":
    from rdkit import RDLogger
    RDLogger.DisableLog("rdApp.error")

    print("=" * 70)
    print("KRAS-VHL PROTAC Builder — Self-Test")
    print(f"Warhead : KRAS (ACBI3-diazepane, NH reactive)")
    print(f"E3      : VHL  (VH032, NH2 reactive)")
    print(f"Chemistry: diacid linker + 2× amide bond formation")
    print("=" * 70)

    builder = KRASVHL_PROTAC_Builder()

    test_linkers = [
        ("OC(=O)CCCCCCCCC(=O)O",   "Sebacic acid (C10 diacid)"),
        ("OC(=O)CCOCCOCCOC(=O)O",  "PEG-diacid (valid)"),
        ("NCCOCCOCCN",              "Diamine — no COOH, should fail"),
        ("OC(=O)CCCC",             "Single COOH — should fail"),
    ]

    print()
    for smiles, desc in test_linkers:
        linker_mol = Chem.MolFromSmiles(smiles)
        print(f"Linker : {desc}")
        if linker_mol is None:
            print("  [SKIP] unparseable SMILES")
            print("-" * 70)
            continue

        cooh_count = len(linker_mol.GetSubstructMatches(Chem.MolFromSmarts(_COOH_SMARTS)))
        print(f"COOH sites: {cooh_count}")

        try:
            protac_smi = builder.assemble(linker_mol, return_all=False)
            protac_mol = Chem.MolFromSmiles(protac_smi)
            passed, checks = builder.validate_protac(protac_mol)
            tag = "VALID ✓" if passed else "ASSEMBLED but INVALID ✗"
            print(f"Assembly: {tag}")
            for k, v in checks.items():
                print(f"  {'✓' if v else '✗'} {k}")
            props = builder.get_protac_properties(protac_mol)
            print(f"  MW: {props['mw']} Da  LogP: {props['logp']}  "
                  f"HeavyAtoms: {props['heavy_atoms']}")
        except Exception as e:
            print(f"Assembly: FAILED — {e}")

        print("-" * 70)