Failed to generate sdf from SMILES

[ghliu521]

Dear Forlilab team,
Thanks for the great work on this package, but I encountered with some issues while using this tool. Some SMILES can't be processes:

"JS5": "C1CC(CNC1)CNCCO[C@@H]2[C@@H]([C@H](O[C@H]2O[C@@H]3[C@H]([C@@H](C[C@@H]([C@H]3O[C@@H]4[C@@H]([C@H]([C@@H]([C@H](O4)CO)O)O)N)N)N)O)CO)O[C@@H]5[C@@H]([C@H]([C@@H]([C@@H](O5)CN)O)O)N", 

and the error infomation is:

Generating ligands_1/JS5_scrubbed.sdf...
/home/ghliu/miniconda3/envs/AutoDock/bin/scrub.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
import('pkg_resources').require('molscrub==0.1.1')
{'INITIAL_COORDS': 1200}
Scrub completed.
Summary of what happened:
Input molecules supplied: 1
mols processed: 0, skipped by rdkit: 0, failed: 1

Other failed SIMILES with similar errors are listed below :

"JS6": "c1ccc(cc1)CCNCCO[C@@H]2[C@@H]([C@H](O[C@H]2O[C@@H]3[C@H]([C@@H](C[C@@H]([C@H]3O[C@@H]4[C@@H]([C@H]([C@@H]([C@H](O4)CO)O)O)N)N)N)O)CO)O[C@@H]5[C@@H]([C@H]([C@@H]([C@@H](O5)CN)O)O)N",
 "G0B": "C1[C@H]([C@@H]([C@H]([C@H]([C@H]1NC(=[NH2+])N)O[C@@H]2[C@H]([C@@H]([C@H]([C@@H](O2)CNC(=[NH2+])N)O)O)NC(=[NH2+])N)O[C@H]3[C@H]([C@@H]([C@H](O3)CO)O[C@@H]4[C@H]([C@@H]([C@H]([C@@H](O4)CNC(=[NH2+])N)O)O)NC(=[NH2+])N)O)O)NC(=[NH2+])N",
"74G": "C[n+]1ccc(c2c1cccc2)C[C@H]3N(c4ccccc4S3)CC(=O)NCCOCCOCCOCCNC(=O)CCCC[C@H]5[C@@H]6[C@H](CS5)NC(=O)N6",
"OIY": "C1[C@H]([C@@H]2[C@@H](C(=O)N1)N=C(N2)N[C@H]3[C@@H]([C@@H]([C@H]([C@H](O3)CO)OC(=O)N)O)NC(=O)C[C@H](CCCNC(=O)C[C@H](CCCNC(=O)C[C@H](CCCN)N)N)N)O",
 "VLR": "CN1C=CC(=Cc2[n+](c3ccccc3s2)CC(=O)NCCOCCOCCOCCNC(=O)CCCC[C@@H]4[C@@H]5[C@@H](CS4)NC(=O)N5)c6c1ccc(c6)c7ccc(cc7)OC",
"LHA": "c1ccc(cc1)CCNCCO[C@@H]2[C@@H]([C@H](O[C@H]2O[C@@H]3[C@H]([C@@H](C[C@@H]([C@H]3O[C@@H]4[C@@H]([C@H]([C@@H]([C@H](O4)CO)O)O)N)N)NC(=O)[C@H](CCN)O)O)CO)O[C@@H]5[C@@H]([C@H]([C@@H]([C@@H](O5)CN)O)O)N",

It will be greatly appreciated if you could have a look.

[ghliu521]

For
"A1IY6": "c1nc2c(nc(nc2n1[C@H]3C@@HO)N)N", even though the sdf file seemed to be successfully generated with output:

Generating ligands_1/A1IY6_scrubbed.sdf...
/home/ghliu/miniconda3/envs/AutoDock/bin/scrub.py:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
import('pkg_resources').require('molscrub==0.1.1')
Scrub completed.
Summary of what happened:
Input molecules supplied: 1
mols processed: 1, skipped by rdkit: 0, failed: 0
nr isomers (tautomers and acid/base conjugates): 1 (avg. 1.000 per mol)
nr conformers: 1 (avg. 1.000 per isomer, 1.000 per mol)

But when I used the generated sdf file for mk_prepare_ligand.py in Meeko, an error occured as below:

Input molecules processed: 0, skipped: 0
PDBQT files written: 0
PDBQT files not written due to error: 1
Input molecules with errors: 1
No duplicate molecule molecule names were found
No PDBQT files were written due to errors!

[diogomart]

They are big molecules, the ETKDG step fails. There is no solution at the moment besides trying multiple times, ETKDG might randomly succeed by chance.

Can you post the SDF of the molecule that fails with meeko?

[ghliu521]

big molecules

Thank you for your reply. Here's the A1IY6 SDF file:

     RDKit          3D

 50 52  0  0  0  0  0  0  0  0999 V2000
    5.2291    0.4158   -0.6228 C   0  0  1  0  0  0  0  0  0  0  0  0
    6.3946   -0.0645    0.0753 N   0  0  0  0  0  0  0  0  0  0  0  0
    4.1085    0.4735    0.3403 C   0  0  0  0  0  0  0  0  0  0  0  0
    5.5354    1.8301   -1.0372 C   0  0  0  0  0  0  0  0  0  0  0  0
    3.6588   -0.8480    0.9747 C   0  0  0  0  0  0  0  0  0  0  0  0
    5.4438    2.0797   -2.2562 O   0  0  0  0  0  0  0  0  0  0  0  0
    5.9184    2.8461   -0.0836 N   0  0  0  0  0  0  0  0  0  0  0  0
    2.3187   -0.4376    2.2131 Se  0  0  0  0  0  2  0  0  0  0  0  0
    0.6311   -1.2290    2.2740 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.0181   -1.6508    0.9845 C   0  0  2  0  0  0  0  0  0  0  0  0
   -0.2097   -0.5028    0.0067 C   0  0  2  0  0  0  0  0  0  0  0  0
   -1.2744   -2.1857    1.1789 O   0  0  0  0  0  0  0  0  0  0  0  0
   -1.3846   -0.9304   -0.7912 C   0  0  2  0  0  0  0  0  0  0  0  0
    0.9095   -0.3507   -0.8116 O   0  0  0  0  0  0  0  0  0  0  0  0
   -2.1087   -1.7900    0.2272 C   0  0  2  0  0  0  0  0  0  0  0  0
   -1.0120   -1.7200   -1.8861 O   0  0  0  0  0  0  0  0  0  0  0  0
   -3.4761   -1.4944    0.3161 N   0  0  0  0  0  0  0  0  0  0  0  0
   -3.9984   -0.2411    0.2345 C   0  0  0  0  0  0  0  0  0  0  0  0
   -4.5220   -2.3741    0.3416 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.5130    1.0097    0.1777 N   0  0  0  0  0  0  0  0  0  0  0  0
   -5.3912   -0.3727    0.2000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.6863   -1.6536    0.2647 N   0  0  0  0  0  0  0  0  0  0  0  0
   -4.2681    2.0855    0.0909 C   0  0  0  0  0  0  0  0  0  0  0  0
   -6.2132    0.7425    0.1056 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.6070    1.9468    0.0545 N   0  0  0  0  0  0  0  0  0  0  0  0
   -3.8415    3.4487    0.0007 N   0  0  0  0  0  0  0  0  0  0  0  0
   -7.5956    0.5896    0.0318 N   0  0  0  0  0  0  0  0  0  0  0  0
    5.1213   -0.2124   -1.5480 H   0  0  0  0  0  0  0  0  0  0  0  0
    7.1631    0.6646    0.1725 H   0  0  0  0  0  0  0  0  0  0  0  0
    6.8802   -0.8712   -0.3737 H   0  0  0  0  0  0  0  0  0  0  0  0
    6.2267   -0.2982    1.0871 H   0  0  0  0  0  0  0  0  0  0  0  0
    3.2294    0.8738   -0.2545 H   0  0  0  0  0  0  0  0  0  0  0  0
    4.3817    1.1645    1.1609 H   0  0  0  0  0  0  0  0  0  0  0  0
    4.5647   -1.4057    1.3735 H   0  0  0  0  0  0  0  0  0  0  0  0
    3.3265   -1.5144    0.1479 H   0  0  0  0  0  0  0  0  0  0  0  0
    5.3026    3.6843    0.0225 H   0  0  0  0  0  0  0  0  0  0  0  0
    6.7337    2.7612    0.5227 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.0946   -0.5473    2.8344 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.6145   -2.1605    2.9489 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.5833   -2.4198    0.4209 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.4762    0.4227    0.5156 H   0  0  0  0  0  0  0  0  0  0  0  0
   -2.0423   -0.1227   -1.1317 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.6234   -0.1312   -1.6959 H   0  0  0  0  0  0  0  0  0  0  0  0
   -2.1956   -2.8560   -0.3646 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.8381   -2.1322   -2.2741 H   0  0  0  0  0  0  0  0  0  0  0  0
   -4.4161   -3.4518    0.4260 H   0  0  0  0  0  0  0  0  0  0  0  0
   -3.0446    3.7413    0.5740 H   0  0  0  0  0  0  0  0  0  0  0  0
   -4.2971    4.1429   -0.6261 H   0  0  0  0  0  0  0  0  0  0  0  0
   -8.1792    1.0554   -0.7456 H   0  0  0  0  0  0  0  0  0  0  0  0
   -8.1952   -0.0100    0.6688 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0
  3  1  1  0
  1  4  1  0
  5  3  1  0
  4  6  2  0
  4  7  1  0
  8  5  1  0
  9  8  1  0
 10  9  1  0
 11 10  1  0
 10 12  1  0
 13 11  1  0
 11 14  1  0
 12 15  1  0
 15 13  1  0
 13 16  1  0
 17 15  1  0
 18 17  1  0
 17 19  1  0
 20 18  1  0
 18 21  2  0
 19 22  2  0
 23 20  2  0
 22 21  1  0
 21 24  1  0
 25 23  1  0
 23 26  1  0
 24 25  2  0
 24 27  1  0
  1 28  1  6
  2 29  1  0
  2 30  1  0
  2 31  1  0
  3 32  1  0
  3 33  1  0
  5 34  1  0
  5 35  1  0
  7 36  1  0
  7 37  1  0
  9 38  1  0
  9 39  1  0
 10 40  1  6
 11 41  1  1
 13 42  1  6
 14 43  1  0
 15 44  1  6
 16 45  1  0
 19 46  1  0
 26 47  1  0
 26 48  1  0
 27 49  1  0
 27 50  1  0
M  CHG  1   2   1
M  END
>  <ScrubInfo>  (1) 
{"isomerGroup": 0, "isomerId": 0, "confId": 0, "nr_conformers:": 1, "nr_isomers:": 1}

$$$$

[diogomart]

Thanks! AutoDock has no parameters for Selenium. I recommend typing Se atoms with type "S" (Sulfur, not H-bond acceptor).

Docs: https://meeko.readthedocs.io/en/release-doc/lig_prep_advanced.html#modifying-atom-types

Example:

mk_prepare_ligand.py -i x.sdf -aa '[{"smarts": "[Se]", "atype": "S"}]'

You'll get the Selenium back after running mk_export.py on docked poses.

[ghliu521]

Thanks. But I am a little confused. Why will I get the Selenium back after running mk_export.py on docked poses. Should I by myself do something to achieve this like by setting some paramers while running mk_export.py ? Except that, I followed your suggestion typing Se atoms with type "S" :
! python3 /home/ghliu/miniconda3/envs/AutoDock/bin/mk_prepare_ligand.py -i "{ligandSDF}" --multimol_outdir "{lig_folder}" --multimol_prefix "{lig}" -aa '[{"smarts": "[Se]", "atype": "S"}]'
but I got another error as below:
*ERROR* Format [] not in supported formats [sdf/mol2/mol]

Why it can't recognise the file format properly? Look forward to your reply.

[diogomart]

The filename of the ligand must end in .sdf, .mol or `.mol2.

mk_export.py reads a SMILES string from the PDBQT remarks, and that string defines the molecule, with Selenium. mk_export.py ignores the atom types.

[ghliu521]

The filename of the ligand must end in .sdf, .mol or `.mol2.

mk_export.py reads a SMILES string from the PDBQT remarks, and that string defines the molecule, with Selenium. mk_export.py ignores the atom types.

Yes my filename is ending with .sdf, here's my code:

ligandSDF = f"{ligandName}_scrubbed.sdf"
print(f"Generating {ligandSDF}...")
if not os.path.exists(ligandSDF):
    ! python3 /home/ghliu/miniconda3/envs/AutoDock/bin/scrub.py "{ligand}" -o {ligandSDF} --ph {pH} {args}

# Prepare ligand PDBQT
if not os.path.exists(ligandPDBQT):
    # ligandSDF = os.path.join('ligands', f"{lig.lower()}_ideal.sdf")
    ! python3 /home/ghliu/miniconda3/envs/AutoDock/bin/mk_prepare_ligand.py -i "{ligandSDF}" --multimol_outdir "{lig_folder}" --multimol_prefix "{lig}" -aa '[{"smarts": "[Se]", "atype": "S"}]

`

[diogomart]

I'm not familiar with that python interpreter. Maybe an f is missing before `"{ligandSDF}" or something similar. The error was raised because the ligand filename extension was not ".sdf".

Here's a suggestion: you could pass rdkit molecule within python without using the filesystem:

from rdkit import Chem
from meeko import MoleculePreparation
from meeko import PDBQTMolecule
from meeko import RDKitMolCreate
from scrubber import Scrub


# input
ligand_smiles = "c1c(Cl)ccn1CCO"

# embed ligand
scrub = Scrub()
rdmol = Chem.MolFromSmiles(ligand_smiles)
isomers = scrub(rdmol)

# consider only the first protonation state
# each molecule may have more than one conformer
# meeko will consider only the first conformer, need to change meeko
mol = isomers[0]

# prepare ligand pdbqt
meeko_prep = MoleculePreparation()
meeko_prep.prepare(mol)
lig_pdbqt = meeko_prep.write_pdbqt_string()

print(lig_pdbqt)

You can even do the docking from vina directly using its Python bindings

[ghliu521]

Thank you anyway. Maybe I should put this error under the Meeko issue area. So anyone who have encountered the same error can provide an solution.