Fix scrambled-molecule docking bug + populate Layer 1.7 blind validation#11
Open
Rockman6 wants to merge 5 commits into
Open
Fix scrambled-molecule docking bug + populate Layer 1.7 blind validation#11Rockman6 wants to merge 5 commits into
Rockman6 wants to merge 5 commits into
Conversation
…lecule bug The pose→molecule reconstruction assumed Meeko "preserves SMILES atom order". It does NOT — Meeko reorders atoms during ligand prep (biotin SMILES order OCOCCCCC… becomes PDBQT order CCCOOCCC…). Both consumers that need a chemically-correct 3D mol — _build_pose_mol (PoseBusters) and pose_rmsd_symmetry_aware (crystal RMSD) — mapped the pose's PDBQT-order coordinates onto a MolFromSmiles template BY INDEX, so they built a topologically scrambled molecule with physically impossible 3–5 Å "bonds" (12 of 17 heavy-atom bonds wrong on biotin). Impact: this silently corrupted BOTH pose-recovery RMSD and PoseBusters validity — the two headline numbers of the docking layer. On 1stp/biotin the scramble reported top-1 RMSD 2.06 Å (borderline fail) and 0% PoseBusters; with correct reconstruction the SAME pose is 0.69 Å (excellent) and passes PoseBusters. Good poses were being scored as bad. Fix: reconstruct each pose with Meeko's own reverse conversion (PDBQTMolecule + RDKitMolCreate.from_pdbqt_mol), which carries the true atom mapping + hydrogens, and store it on DockingPose.rdkit_molblock. _build_pose_mol and pose_rmsd now use that; the legacy template path is kept only as a fallback for poses lacking a molblock. This also makes the separate OpenMM pose-refine step unnecessary for validity. Regression test docks biotin (a ligand Meeko reorders) and asserts every reconstructed heavy-atom bond is physically plausible. Full tests/dock suite: 41 passed.
…st results Populates the previously-empty Layer 1.7 (blind validation) with a real, reproducible benchmark — the anti-circular-validation gate that was missing. scripts/run_blind_dock_bench.py fetches structures AND authoritative ligand SMILES from RCSB at run time (never from memory or our own inputs), auto-picks the drug ligand (skipping waters/ions/ buffers), re-docks blind, and scores top-1 RMSD (criterion #1, gate >= 75%) + PoseBusters physical validity (criterion #3, gate >= 95%). Bad entries are reported and excluded from the denominators, never silently counted. First run on a 15-cocrystal diverse set (benchmarks/pdbbind/blind_set.yaml), after the Meeko pose-reconstruction fix: - PoseBusters validity: 100% (10/10 scored) -> PASSES #3. - Pose recovery (top-1 <= 2.0 A): 70% (7/10) -> just under #1's 75% gate. 5 entries failed the crystal-RMSD computation (docking + PB succeeded; only the RMSD-vs-crystal atom match failed) and are excluded — a known robustness gap to fix next. 3 more are borderline top-1 but good top-3 (1m17, 1fkg, 1hpx). This is the first honest, external measurement of criteria #1/#3. Structures cached under benchmarks/pdbbind/structures/ (gitignored, re-fetched on demand).
The crystal-RMSD path dropped 5/15 blind-set entries to None (docking +
PoseBusters succeeded; only the RMSD-vs-crystal atom match failed):
- multiple ligand copies in the crystal (imatinib in 2hyy appears 4x,
raloxifene in 1err 2x) -> extract_hetatm_ligand grabbed all of them,
so the atom count was a multiple of one ligand;
- crystal missing a disordered atom, or the deposited SMILES carrying
an extra (3ptb 10 vs 9, 1dwd 38 vs 37) -> exact-count match refused;
- a subtle topology mismatch (1hvr) even at equal counts.
Fixes in pose_rmsd_symmetry_aware:
- extract_hetatm_ligand(first_copy_only=True): one ligand copy (first
chain/resSeq, blank/'A' altloc) for RMSD references.
- a robust matcher that prefers the Meeko-reconstructed pose mol and
falls back to an MCS common-core in-place RMSD when atom counts
differ, so near-misses score instead of vanishing. (Ring perception
is initialised before FMCS to avoid the sanitize=False RingInfo
error.)
- the legacy exact-count template path is kept only as a last resort.
Result: all 15/15 blind cocrystals now score. Harness also reports top-3
recovery (sampling power, ranking-independent) alongside top-1.
Honest Layer-1.7 numbers on the 15-cocrystal diverse blind set:
- PoseBusters validity: 100% -> PASSES criterion #3 (>= 95%).
- Pose recovery top-1 <= 2 A: 67% -> under criterion #1 (>= 75%).
- Pose recovery top-3 <= 2 A: 87% -> Vina SAMPLES the right pose
13/15; its SCORING just doesn't rank it #1. Criterion #1's shortfall
is a ranking limitation, the same root as criterion #2 (kinase
ranking) — which is exactly the gap FEP is meant to close.
tests/dock: 43 passed.
…tals) run_one derived the box centre from extract_hetatm_ligand over ALL HETATM copies of the ligand. For crystals with the ligand in several chains (imatinib in 2hyy = 4 copies spread across the structure) the centroid landed BETWEEN the binding sites, so Vina docked into empty space and the "pose recovery" for those entries was an artefact of the wrong box, not of the docking. Use first_copy_only for the centre. Effect on the 15-cocrystal blind set: 2hyy 2.77 -> 0.80 A (fail -> pass), 1err 0.83 -> 0.60 A. Updated honest numbers: - PoseBusters validity: 100% -> PASSES criterion #3. - Pose recovery top-1 <= 2A: 73% (11/15) -> just under criterion #1. - Pose recovery top-3 <= 2A: 87% (13/15) -> sampling power. Remaining top-1 misses: 1m17 + 1hpx are ranking failures (correct pose at rank 2-3), 1fkg + 1dwd are genuine sampling near-misses (top-3 best 2.26 / 2.34 A). Confirmed separately that UFF-strain re-ranking does NOT recover the ranking failures (the correct pose is not lower-strain), so top-1 is a Vina-scoring limit, not fixable by a strain gate — the same root as criterion #2 (kinase ranking), which is what FEP addresses.
…ed mols) README, BENCHMARKS, and campaign1_status quoted "2/3 = 67% pose recovery" and per-system RMSDs (1STP 2.02 Å, 1M17 4.23 Å) that were computed on molecules scrambled by the Meeko atom-order bug. Replace with the honest blind-set numbers on correct molecules: PoseBusters 100% (criterion #3 PASS), pose recovery top-3 87% / top-1 73%. Note the fix and the top-1/top-3 (scoring-vs-sampling) story throughout.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
The bug (load-bearing)
The docking layer assumed Meeko "preserves SMILES atom order." It does not — Meeko reorders atoms during ligand prep (biotin: SMILES
OCOCCCCC…→ PDBQTCCCOOCCC…). Both consumers that need a correct 3D molecule —_build_pose_mol(PoseBusters) andpose_rmsd_symmetry_aware(crystal RMSD) — mapped pose coordinates onto aMolFromSmilestemplate by index, producing a topologically scrambled molecule with physically impossible 3–5 Å "bonds" (12/17 wrong on biotin).Impact: this silently corrupted both headline numbers of Layer 1.3a — pose recovery and PoseBusters validity. Biotin/1STP was reported at 2.06 Å (borderline fail) when the real redock is 0.69 Å (excellent). Good poses were scored as bad.
The fixes (5 commits)
2d222e4— reconstruct poses via Meeko's reverse conversion (PDBQTMolecule+RDKitMolCreate.from_pdbqt_mol), storing a correct mol block on each pose. 41/41 dock tests pass; regression test guards against the scramble returning.1e9e845— new blind-validation harness (scripts/run_blind_dock_bench.py): fetches structures and authoritative ligand SMILES from RCSB at run time (not from memory — so it's not circular on our inputs), auto-picks the drug ligand, re-docks blind, scores RMSD + PoseBusters. Populates the previously-empty Layer 1.7.7e52547— robust crystal-RMSD: one-copy extraction (multi-copy crystals) + MCS common-core matching (off-by-one atom counts) + ring-info init. All 15/15 blind entries now score.0c8f265— centre the docking box on one ligand copy (multi-copy crystals were centred between binding sites → docked into empty space). 2hyy 2.77 → 0.80 Å.docs— correct the README/BENCHMARKS/status numbers that were computed on scrambled molecules.Honest Layer-1.7 results (15 diverse cocrystals)
The top-1/top-3 gap is the finding: Vina's sampling finds the correct pose 87% of the time; its scoring ranks it #1 only 73%. A UFF-strain re-rank does not recover the ranking failures. This is the same limitation as criterion #2 (kinase ranking) — Vina scoring — which is exactly what FEP replaces.
tests/dock: 43 passed. Separate from the FEP audit branch (PR #10).