|
| 1 | +# pylint: disable=E0401,E1101,I1101 |
| 2 | +"""Filter molecules with kekulization errors. |
| 3 | +Handle molecules with sanitization errors by assign formal charge based on valence.""" |
| 4 | +# https://depth-first.com/articles/2020/02/10/a-comprehensive-treatment-of-aromaticity-in-the-smiles-language/ |
| 5 | +import argparse |
| 6 | +from pathlib import Path |
| 7 | +from typing import Dict, Tuple |
| 8 | + |
| 9 | +import rdkit |
| 10 | +from rdkit import Chem |
| 11 | +from rdkit.Chem import AllChem |
| 12 | + |
| 13 | +parser = argparse.ArgumentParser() |
| 14 | +parser.add_argument('--input_small_mol_ligand', type=str, help='Ligand file') |
| 15 | +args = parser.parse_args() |
| 16 | +input_small_mol_ligand = Path(args.input_small_mol_ligand).name |
| 17 | + |
| 18 | + |
| 19 | +def adjust_formal_charges(molecule: Chem.SDMolSupplier) -> Chem.SDMolSupplier: |
| 20 | + """Sometimes input structures do not have correct formal charges corresponding |
| 21 | + to bond order topology. So choose to trust bond orders assigned and generate formal |
| 22 | + charges based on that. |
| 23 | + Explicit valence determined what the formal charge should be from dictionary of valence |
| 24 | + to formal charge for that atomic number. Special case if atom == carbon or nitrogen |
| 25 | + and if neighbors contain nitrogen, oyxgen or sulfur (polarizable atoms) then if carbon |
| 26 | + and explicit valence only 3, give formal charge of +1 (more stable then -1 case). |
| 27 | +
|
| 28 | + Args: |
| 29 | + molecule (Chem.SDMolSupplier): The rdkit molecule object |
| 30 | +
|
| 31 | + Returns: |
| 32 | + Chem.SDMolSupplier: Molecule object with adjusted formal charges |
| 33 | + """ |
| 34 | + # 1=H, 5=B, 6=C, 7=N, 8=O, 15=P, 16=S, 17=Cl, 9=F, 35=Br, 53=I |
| 35 | + atomicnumtoformalchg: Dict[int, Dict[int, int]] = {1: {2: 1}, 5: {4: 1}, 6: {3: -1}, 7: {2: -1, 4: 1}, |
| 36 | + 8: {1: -1, 3: 1}, 15: {4: 1}, 16: {1: -1, 3: 1, 5: -1}, |
| 37 | + 17: {0: -1, 4: 3}, 9: {0: -1}, 35: {0: -1}, 53: {0: -1}} |
| 38 | + for atom in molecule.GetAtoms(): |
| 39 | + atomnum = atom.GetAtomicNum() |
| 40 | + val = atom.GetExplicitValence() |
| 41 | + if atomicnumtoformalchg.get(atomnum) is None: |
| 42 | + continue |
| 43 | + valtochg = atomicnumtoformalchg[atomnum] |
| 44 | + chg = valtochg.get(val, 0) |
| 45 | + # special case of polar neighbors surrounding carbon or nitrogen |
| 46 | + # https://docs.eyesopen.com/toolkits/cpp/oechemtk/valence.html#openeye-charge-model |
| 47 | + # |
| 48 | + # See Section 6: Factors That Stabilize Carbocations – Adjacent Lone Pairs |
| 49 | + # https://www.masterorganicchemistry.com/2011/03/11/3-factors-that-stabilize-carbocations/#six |
| 50 | + polneighb = False |
| 51 | + if atomnum in (6, 7): |
| 52 | + for natom in atom.GetNeighbors(): |
| 53 | + natomicnum = natom.GetAtomicNum() |
| 54 | + if natomicnum in (7, 8, 16): |
| 55 | + polneighb = True |
| 56 | + if polneighb and val == 3 and atomnum == 6: |
| 57 | + chg = 1 |
| 58 | + |
| 59 | + atom.SetFormalCharge(chg) |
| 60 | + return molecule |
| 61 | + |
| 62 | + |
| 63 | +def generate_conformer(molecule: Chem.SDMolSupplier) -> None: |
| 64 | + """ Generate conformer for molecule. Sometimes rdkit embedding can fail, |
| 65 | + so use random coordinates if failed at first. |
| 66 | +
|
| 67 | + Args: |
| 68 | + molecule (Chem.SDMolSupplier): _description_ |
| 69 | + """ |
| 70 | + ps = AllChem.ETKDGv2() |
| 71 | + confid = AllChem.EmbedMolecule(molecule, ps) |
| 72 | + if confid == -1: |
| 73 | + ps.useRandomCoords = True |
| 74 | + AllChem.EmbedMolecule(molecule, ps) |
| 75 | + # only want to catch error Bad Conformation Id, dont need to spend time optimizing |
| 76 | + AllChem.MMFFOptimizeMolecule(molecule, confId=0, maxIters=1) |
| 77 | + |
| 78 | + |
| 79 | +def is_valid_ligand(molecule: Chem.SDMolSupplier) -> Tuple[bool, bool, Chem.SDMolSupplier]: |
| 80 | + """Check for sanitization errors, attempt to fix formal charges/valence consistency errors. |
| 81 | + DiffDock uses rdkit to generate a seed conformation that will sometimes crash, so generating |
| 82 | + conformations here to catch that error and prevent DiffDock from running that ligand. |
| 83 | +
|
| 84 | + Args: |
| 85 | + mol (Chem.SDMolSupplier): The rdkit small molecule object |
| 86 | +
|
| 87 | + Returns: |
| 88 | + bool: if ligand is valid |
| 89 | + bool: if symlink should be used |
| 90 | + Chem.SDMolSupplier: molecule object |
| 91 | + """ |
| 92 | + valid_lig = True |
| 93 | + use_symlink = True |
| 94 | + try: |
| 95 | + Chem.SanitizeMol(molecule) |
| 96 | + generate_conformer(molecule) |
| 97 | + except rdkit.Chem.rdchem.KekulizeException as e: |
| 98 | + valid_lig = False |
| 99 | + # Not handling kekulization error now so just remove file to prevent DiffDock execution |
| 100 | + except rdkit.Chem.rdchem.MolSanitizeException as e: |
| 101 | + # can also be explicit valence error (i.e.) formal charge not consistent with bond topology |
| 102 | + # choose to trust bond topology around atom and add formal charge based on that |
| 103 | + molecule = adjust_formal_charges(molecule) |
| 104 | + use_symlink = False |
| 105 | + except ValueError: |
| 106 | + # assuming this is Bad Conformer Id error from generate_conformer |
| 107 | + valid_lig = False |
| 108 | + except Exception: # pylint: disable=broad-exception-caught |
| 109 | + # catch *all* exceptions rdkit can throw |
| 110 | + valid_lig = False |
| 111 | + |
| 112 | + return valid_lig, use_symlink, molecule |
| 113 | + |
| 114 | + |
| 115 | +mol: Chem.SDMolSupplier = Chem.SDMolSupplier(args.input_small_mol_ligand, sanitize=False, removeHs=False)[0] |
| 116 | + |
| 117 | +valid_ligand, symlink, rdkit_mol = is_valid_ligand(mol) |
| 118 | + |
| 119 | +if not symlink and valid_ligand: |
| 120 | + with Chem.SDWriter(input_small_mol_ligand) as w: |
| 121 | + w.write(rdkit_mol) |
| 122 | + |
| 123 | +with open('valid.txt', 'w', encoding="utf-8") as f: |
| 124 | + f.write(str(valid_ligand)) |
0 commit comments