Voxel descriptors - problem with assigning atom types to SmallMol

[larry0x]

Hi devs,

Before you get the documentations ready, I tried to figure out how to generate voxel descriptors for my small molecules myself. I used this molecule from the PDB as a test. Below is the code:

from rdkit import Chem

from moleculekit.smallmol.smallmol import SmallMol
from moleculekit.tools.voxeldescriptors import getVoxelDescriptors

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# read molecule
suppl = Chem.SDMolSupplier("../data/GDP_model.sdf")
mol = SmallMol(suppl[0])

# generate voxels
features, _, _ = getVoxelDescriptors(mol,
                                     boxsize=(16, 16, 16),
                                     voxelsize=0.5,
                                     center=mol.getCenter())
features = features.reshape(32, 32, 32, 8)

# visualize voxels
channel_names = ["hydrophobic", "aromatic", "hbond_acceptor", "hbond_donor",
                 "positive_ionizable", "negative_ionizable", "metal", "occupancies"]
threshold = 0.9
features = (features >= threshold).astype(float)

fig = plt.figure(figsize=(16, 8))

for i in range(8):
    ax = fig.add_subplot(241+i, projection="3d")
    ax.voxels(features[:, :, :, i], facecolors="red", edgecolor="k")
    ax.set_title(channel_names[i])

fig.tight_layout()
plt.show()

Most of the channels generated by this code look fine, but hydrophobic and negative_ionizable channels are completely empty. This does not make sense chemically, since at least some carbon atoms on the guanine and ribose rings should be considered hydrophobic, and the hydroxyls on the phosphates are definitively negatively ionizable.

I noticed when a SmallMol object is passed into getVoxelDescriptors, the function calls _getPropertiesRDkit, which in turn calls some RDKit functions to assign properties to atoms. I tried to print the channels array generated by this function, and it looked like all atoms are assigned zeros for hydrophobic and negative_ionizable channels.

Did I do this wrong, or is this an issue with RDKit? Thank you

[stefdoerr]

Hi Larry, I'll look into it the coming week. Thanks for the report :)

[stefdoerr]

From the code side it looks fine. If you have VMD use this code to visualize the features in VMD, and also use SmallMolLib if you want.

from moleculekit.smallmol.smallmol import SmallMol
from moleculekit.smallmol.smallmollib import SmallMolLib
from moleculekit.tools.voxeldescriptors import getVoxelDescriptors, viewVoxelFeatures
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# read molecule
suppl = SmallMolLib("GDP_ideal.sdf")
mol = suppl[0]

# generate voxels
features, centers, N = getVoxelDescriptors(mol,
                                     boxsize=(18, 16, 16),
                                     voxelsize=0.5,
                                     center=mol.getCenter())
viewVoxelFeatures(features, centers, N, voxelsize=0.5)
mol.view(style='Licorice')

As for the hydrophobicity and negative ionizable, @cuzzo87 will have to chime in as I'm no chemist. I assume it's an issue of rdkit though as you said

[stefdoerr]

According to @cuzzo87 it's an rdkit issue since it can only return one descriptor/feature per atom. This means that an atom can be either hydrophobic or aromatic but not both. This is not something we can correct easily. But from an ML point of view it might not matter in the end as long as it's consistently putting atoms into one category.

[larry0x]

Hi @stefdoerr @cuzzo87 thanks for looking into into this issue. I agree this should not be a problem for ML.

[stefdoerr]

I will close this as we are not planning to modify rdkit. Thanks for the input!