Problem when assigning bonds from a PDB structure

[PatWalters]

I'm not sure what I'm doing wrong here.

from rdkit import Chem
from steamroll.steamroll import to_rdkit

pdb_block = """HETATM    2  C9  BZU A 555      -3.394   5.205  13.019  1.00  0.00           C
HETATM    3  C10 BZU A 555      -3.055   3.867  13.710  1.00  0.00           C
HETATM    4  C11 BZU A 555      -3.349   7.408  13.935  1.00  0.00           C
HETATM    5  C12 BZU A 555      -2.076   7.426  14.764  1.00  0.00           C
HETATM    6  C14 BZU A 555      -1.915   9.769  15.059  1.00  0.00           C
HETATM    7  C15 BZU A 555      -0.878  10.903  15.178  1.00  0.00           C
HETATM    8  O1A BZU A 555      -5.425   0.994  18.287  1.00  0.00           O
HETATM    9  O2A BZU A 555      -7.398   0.817  16.995  1.00  0.00           O
HETATM   10  N21 BZU A 555      -5.486  -0.545  16.461  1.00  0.00           N
HETATM   11  S1  BZU A 555      -5.966   0.750  16.997  1.00  0.00           S
HETATM   12  C4  BZU A 555      -4.321   2.152  15.292  1.00  0.00           C
HETATM   13  C5  BZU A 555      -4.165   3.345  14.642  1.00  0.00           C
HETATM   14  C6  BZU A 555      -5.216   4.200  14.931  1.00  0.00           C
HETATM   15  S2  BZU A 555      -6.379   3.501  15.979  1.00  0.00           S
HETATM   16  S7  BZU A 555      -5.295   5.827  14.368  1.00  0.00           S
HETATM   17  O3B BZU A 555      -6.160   5.873  13.244  1.00  0.00           O
HETATM   18  O4B BZU A 555      -5.858   6.680  15.373  1.00  0.00           O
HETATM   19  N8  BZU A 555      -3.848   6.065  14.090  1.00  0.00           N
HETATM   20  N16 BZU A 555      -2.776   2.895  12.649  1.00  0.00           N
HETATM   21  O13 BZU A 555      -1.307   8.596  14.533  1.00  0.00           O
HETATM   22  C17 BZU A 555      -1.564   3.011  11.841  1.00  0.00           C
HETATM   23  C18 BZU A 555      -1.463   1.927  10.754  1.00  0.00           C
"""

mol = Chem.MolFromPDBBlock(pdb_block)
crds = mol.GetConformer(0).GetPositions()
atomic_nums = [x.GetAtomicNum() for x in mol.GetAtoms()]
rd_mol = to_rdkit(atomic_nums,crds,charge=0)
print(Chem.MolToSmiles(rd_mol))

This fails to generate bonds or bond orders. Adding the SMILES to the input doesn't help.

[vandan-revanur]

Just for some more clarity the problem seems to be that the rdkit mol that the to_rdkit function from steamroll doesn't produce a proper/true rdkit mol.

It can be seen from some print statements I added in the snippet from @PatWalters

mol = Chem.MolFromPDBBlock(pdb_block)
print(mol)
>> <rdkit.Chem.rdchem.Mol object at 0x7b8bb41e6c00>

print(Chem.MolToSmiles(mol))

>> CCN[C@H]1CN(CCOCC)S(O)(O)[C@H](S)C1C.N[SH](O)O

crds = mol.GetConformer(0).GetPositions()
atomic_nums = [x.GetAtomicNum() for x in mol.GetAtoms()]

print(crds)
>> coords [[-3.394  5.205 13.019]
 [-3.055  3.867 13.71 ]
 [-3.349  7.408 13.935]
 [-2.076  7.426 14.764]
 [-1.915  9.769 15.059]
 [-0.878 10.903 15.178]
 [-5.425  0.994 18.287]
 [-7.398  0.817 16.995]
 [-5.486 -0.545 16.461]
 [-5.966  0.75  16.997]
 [-4.321  2.152 15.292]
 [-4.165  3.345 14.642]
 [-5.216  4.2   14.931]
 [-6.379  3.501 15.979]
 [-5.295  5.827 14.368]
 [-6.16   5.873 13.244]
 [-5.858  6.68  15.373]
 [-3.848  6.065 14.09 ]
 [-2.776  2.895 12.649]
 [-1.307  8.596 14.533]
 [-1.564  3.011 11.841]
 [-1.463  1.927 10.754]]

print(atomic_nums)
>> [6, 6, 6, 6, 6, 6, 8, 8, 7, 16, 6, 6, 6, 16, 16, 8, 8, 7, 7, 8, 6, 6]

# Here we obtain rdkit mol from steamroll

rd_mol = to_rdkit(atomic_nums,crds,charge=0)
print(rd_mol)

>> <rdkit.Chem.rdchem.Mol object at 0x7b8bb674df50>

print(Chem.MolToSmiles(rd_mol))

>> C.C.C.C.C.C.C.C.C.C.C.N.N.N.O.O.O.O.O.S.S.S

Notice that the Smiles of the rdkit mol obtained from steamroll is different from the rdkit mol directly constructed by rdkit.

[corinwagen]

Thanks for the bug report, Pat - we've used this much less in cases where the molecule is missing Hs (relevant for PDBs, obviously) and there's a relatively simple error which is causing issues here. I'll let you know when this is fixed.