"""Redock a co-crystal ligand into its own receptor and measure the RMSD. The standard validation for a docking setup: take a solved complex, throw the ligand's pose away, dock it back, and see whether the top pose reproduces the crystal. Under 2.0 A is the accepted success criterion. Two things make this an honest test rather than a rehearsal. The ligand's bond orders come from the CCD block through `mol_from_block`, the same path the fingerprints use, so the molecule docked is the molecule the project believes it is. That path carries the 5 September bond order fix. The input conformer is generated fresh with ETKDG and the crystal coordinates are never given to the docker. Only the BOX comes from the crystal ligand, which is what a real screen would take from a known site. Feeding the deposited pose would let the optimizer start on the answer. """ import argparse, json, os, subprocess, sys, time import numpy as np HERE = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) sys.path.insert(0, os.path.join(HERE, "code")) WORK = os.path.join(HERE, "work") DOCK = os.path.join(HERE, "docking") ENV = os.path.join(DOCK, "env") import ligand_fingerprints as LF import ccd as CCD from rdkit import Chem, RDLogger from rdkit.Chem import AllChem, rdMolAlign RDLogger.DisableLog("rdApp.*") SKIP = {"HOH", "DOD"} def log(m): sys.stderr.write("[%s] %s\n" % (time.strftime("%H:%M:%S"), m)) sys.stderr.flush() def write_receptor_pdb(atoms, lig_key, path): """Protein only: no waters, no the ligand itself, no other heteroatoms.""" n = 0 with open(path, "w") as fh: for a in atoms: if a["resname"] in SKIP or a["el"] in ("H", "D"): continue if (a["chain"], a["resnum"], a["resname"]) == lig_key: continue n += 1 nm = a["name"] fh.write("%-6s%5d %s %3s %s%4s %8.3f%8.3f%8.3f 1.00 0.00" " %2s\n" % ("ATOM" if a["group"] == "ATOM" else "HETATM", n, (" %-3s" % nm) if len(nm) < 4 else nm, a["resname"][:3], (a["chain"] or "A")[:1], str(a["resnum"])[:4], a["xyz"][0], a["xyz"][1], a["xyz"][2], a["el"][:2].rjust(2))) fh.write("END\n") return n BLOCKS = os.path.join(WORK, "ccd_blocks_all.jsonl") def crystal_mol(comp_id, lig_atoms): """The ligand with CCD bond orders and its deposited coordinates.""" blocks = CCD.load(BLOCKS if os.path.exists(BLOCKS) else CCD.BLOCKS) if comp_id not in blocks: blocks = CCD.load() blk = blocks.get(comp_id) if not blk: raise RuntimeError("no CCD block for %s" % comp_id) mol, names = LF.mol_from_block(blk) if mol is None: raise RuntimeError("CCD block would not build") coords = {a["name"]: a["xyz"] for a in lig_atoms} if sum(1 for nm in names if nm in coords) < mol.GetNumAtoms() * 0.8: raise RuntimeError("crystal coordinates do not cover the component") conf = Chem.Conformer(mol.GetNumAtoms()) keep = [] for i, nm in enumerate(names): if nm in coords: conf.SetAtomPosition(i, coords[nm]); keep.append(i) else: conf.SetAtomPosition(i, (0.0, 0.0, 0.0)) m = Chem.Mol(mol); m.RemoveAllConformers(); m.AddConformer(conf, assignId=True) return m, names, keep def main(): ap = argparse.ArgumentParser() ap.add_argument("--pdb", required=True) ap.add_argument("--comp", required=True) ap.add_argument("--chain", default=None) ap.add_argument("--resnum", default=None) ap.add_argument("--exhaustiveness", type=int, default=16) ap.add_argument("--cpu", type=int, default=20) ap.add_argument("--pad", type=float, default=6.0) ap.add_argument("--seed", type=int, default=20260906) ap.add_argument("--outdir", default=os.path.join(DOCK, "runs")) a = ap.parse_args() os.makedirs(a.outdir, exist_ok=True) tag = "%s_%s" % (a.pdb.upper(), a.comp.upper()) d = os.path.join(a.outdir, tag) os.makedirs(d, exist_ok=True) atoms = LF.parse_atoms(LF.fetch_cif(a.pdb)) cands = {} for at in atoms: if at["resname"] == a.comp.upper() and at["el"] not in ("H", "D"): cands.setdefault((at["chain"], at["resnum"], at["resname"]), []).append(at) if a.chain: cands = {k: v for k, v in cands.items() if k[0] == a.chain} if a.resnum: cands = {k: v for k, v in cands.items() if k[1] == str(a.resnum)} if not cands: print("ligand copy not found"); return lig_key = sorted(cands, key=lambda k: -len(cands[k]))[0] lig_atoms = cands[lig_key] log("%s %s copy %s, %d heavy atoms" % (a.pdb, a.comp, lig_key, len(lig_atoms))) nrec = write_receptor_pdb(atoms, lig_key, os.path.join(d, "rec.pdb")) log("receptor atoms %d" % nrec) xtal, names, keep = crystal_mol(a.comp.upper(), lig_atoms) X = xtal.GetConformer().GetPositions()[keep] cen = X.mean(axis=0) ext = X.max(axis=0) - X.min(axis=0) box = np.maximum(ext + 2 * a.pad, 16.0) log("box centre %s size %s" % (np.round(cen, 2), np.round(box, 1))) # a FRESH conformer, so no crystal geometry reaches the docker fresh = Chem.Mol(xtal) fresh.RemoveAllConformers() fresh = Chem.AddHs(fresh) ps = AllChem.ETKDGv3(); ps.randomSeed = a.seed if AllChem.EmbedMolecule(fresh, ps) != 0: print("could not embed a fresh conformer"); return AllChem.MMFFOptimizeMolecule(fresh, maxIters=500) Chem.MolToMolFile(fresh, os.path.join(d, "lig_fresh.sdf")) py = os.path.join(ENV, "bin", "python") def run(cmd, name): r = subprocess.run(cmd, cwd=d, capture_output=True, text=True) open(os.path.join(d, name + ".log"), "w").write(r.stdout + r.stderr) return r r = run([py, os.path.join(ENV, "bin", "mk_prepare_ligand.py"), "-i", "lig_fresh.sdf", "-o", "lig.pdbqt"], "prep_lig") if not os.path.exists(os.path.join(d, "lig.pdbqt")): print("ligand prep failed:\n", r.stdout[-800:], r.stderr[-800:]); return r = run([py, os.path.join(ENV, "bin", "mk_prepare_receptor.py"), "--read_pdb", "rec.pdb", "-o", "rec", "-p", "--box_center", "%.3f" % cen[0], "%.3f" % cen[1], "%.3f" % cen[2], "--box_size", "%.1f" % box[0], "%.1f" % box[1], "%.1f" % box[2], "-x", "--forgive_extra_bonds"], "prep_rec") rec = os.path.join(d, "rec.pdbqt") if not os.path.exists(rec): print("receptor prep failed:\n", r.stdout[-1500:], r.stderr[-1500:]); return t0 = time.time() r = run([os.path.join(ENV, "bin", "vina"), "--receptor", "rec.pdbqt", "--ligand", "lig.pdbqt", "--center_x", "%.3f" % cen[0], "--center_y", "%.3f" % cen[1], "--center_z", "%.3f" % cen[2], "--size_x", "%.1f" % box[0], "--size_y", "%.1f" % box[1], "--size_z", "%.1f" % box[2], "--exhaustiveness", str(a.exhaustiveness), "--cpu", str(a.cpu), "--seed", str(a.seed), "--out", "docked.pdbqt"], "vina") dt = time.time() - t0 if not os.path.exists(os.path.join(d, "docked.pdbqt")): print("docking failed:\n", r.stdout[-1500:], r.stderr[-1500:]); return # meeko lives only in docking/env, deliberately, so the conversion runs # THERE and comes back as an SDF this interpreter can read conv = os.path.join(d, "convert.py") open(conv, "w").write( "from meeko import PDBQTMolecule, RDKitMolCreate\n" "from rdkit import Chem\n" "pm = PDBQTMolecule.from_file('docked.pdbqt', skip_typing=True)\n" "ms = RDKitMolCreate.from_pdbqt_mol(pm, only_cluster_leads=False)\n" "w = Chem.SDWriter('docked.sdf')\n" "m = ms[0]\n" "for c in m.GetConformers():\n" " w.write(m, confId=c.GetId())\n" "w.close()\n") r = run([py, "convert.py"], "convert") sdf = os.path.join(d, "docked.sdf") if not os.path.exists(sdf): print("pose conversion failed:\n", r.stdout[-800:], r.stderr[-800:]); return sup = [m for m in Chem.SDMolSupplier(sdf, removeHs=False) if m is not None] mol = Chem.Mol(sup[0]) mol.RemoveAllConformers() for m in sup: mol.AddConformer(m.GetConformer(), assignId=True) scores = [] for l in open(os.path.join(d, "docked.pdbqt")): if l.startswith("REMARK VINA RESULT"): scores.append(float(l.split()[3])) ref = Chem.RemoveHs(Chem.Mol(xtal)) out = [] for k in range(mol.GetNumConformers()): probe = Chem.RemoveHs(Chem.Mol(mol)) try: rms = rdMolAlign.CalcRMS(probe, ref, prbId=k, refId=0) except Exception: rms = float("nan") out.append(rms) print("\nRedocking %s ligand %s with AutoDock Vina 1.2.7" % (a.pdb.upper(), a.comp.upper())) print("receptor atoms %d, ligand heavy atoms %d, rotatable bonds %d" % (nrec, ref.GetNumHeavyAtoms(), Chem.rdMolDescriptors.CalcNumRotatableBonds(ref))) print("box %s A, exhaustiveness %d, %d cpu, %.1f s\n" % (np.round(box, 1).tolist(), a.exhaustiveness, a.cpu, dt)) print("%-5s %12s %10s" % ("pose", "vina score", "RMSD")) for k, rms in enumerate(out): print("%-5d %12.2f %10.2f" % (k + 1, scores[k] if k < len(scores) else float("nan"), rms)) best = float(np.nanmin(out)) if out else float("nan") print("\ntop pose RMSD %.2f A, best of %d poses %.2f A" % (out[0], len(out), best)) print("success criterion for redocking is under 2.00 A on the TOP pose") json.dump({"pdb": a.pdb.upper(), "comp": a.comp.upper(), "receptor_atoms": nrec, "seconds": dt, "scores": scores, "rmsd": out, "top_rmsd": out[0], "best_rmsd": best}, open(os.path.join(d, "result.json"), "w"), indent=1) if __name__ == "__main__": main()