Loading deepchem/feat/grid_featurizers.py +246 −2 Original line number Diff line number Diff line Loading @@ -5,6 +5,12 @@ import logging from deepchem.utils.rdkit_util import get_partial_charge from deepchem.feat import ComplexFeaturizer from deepchem.utils.voxel_utils import voxelize from deepchem.utils.rdkit_util import compute_salt_bridges from deepchem.utils.rdkit_util import compute_binding_pocket_cation_pi from deepchem.utils.voxel_utils import convert_atom_to_voxel from deepchem.utils.voxel_utils import convert_atom_pair_to_voxel from deepchem.utils.rdkit_util import compute_pairwise_distances from deepchem.utils.rdkit_util import compute_pi_stack logger = logging.getLogger(__name__) Loading @@ -28,7 +34,7 @@ class ChargeVoxelizer(ComplexFeaturizer): constitutent molecules, compute the partial (Gasteiger charge) on each molecule. For each atom, localize this partial charge in the voxel in which it originated to create a local charge array. Sum contributes to get an effective a local charge array. Sum contributions to get an effective charge at each voxel. Creates a tensor output of shape (box_width, box_width, Loading @@ -47,7 +53,6 @@ class ChargeVoxelizer(ComplexFeaturizer): voxel_width: float, optional (default 1.0) Size of a 3D voxel in a grid. """ self.cutoff = cutoff self.box_width = box_width self.voxel_width = voxel_width self.voxels_per_edge = int(self.box_width / self.voxel_width) Loading Loading @@ -79,3 +84,242 @@ class ChargeVoxelizer(ComplexFeaturizer): for xyz, mol in ((prot_xyz, prot_rdk), (lig_xyz, lig_rdk)) ]) ] class SaltBridgeVoxelizer(ComplexFeaturizer): """Localize salt bridges between atoms in macromolecular complexes. Given a macromolecular compelx made up of multiple constitutent molecules, compute salt bridges between atoms in the macromolecular complex. For each atom, localize this salt bridge in the voxel in which it originated to create a local salt bridge array. Note that if atoms in two different voxels interact in a salt-bridge, the interaction is double counted in both voxels. Creates a tensor output of shape (box_width, box_width, box_width, 1) for each macromolecular complex that computes the number of salt bridges at each voxel. """ def __init__(self, cutoff=5.0, box_width=16.0, voxel_width=1.0): """ Parameters ---------- cutoff: float, optional (default 5.0) The distance in angstroms within which atoms must be to be considered for a salt bridge between them. box_width: float, optional (default 16.0) Size of a box in which voxel features are calculated. Box is centered on a ligand centroid. voxel_width: float, optional (default 1.0) Size of a 3D voxel in a grid. """ self.cutoff = cutoff self.box_width = box_width self.voxel_width = voxel_width self.voxels_per_edge = int(self.box_width / self.voxel_width) def _featurize_complex(self, mol, protein): """ Compute featurization for a single mol/protein complex Parameters ---------- mol: object Representation of the molecule protein: object Representation of the protein """ (lig_xyz, lig_rdk), (prot_xyz, prot_rdk) = mol, protein distances = compute_pairwise_distances(prot_xyz, lig_xyz) return [ voxelize( convert_atom_pair_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, (prot_xyz, lig_xyz), feature_list=compute_salt_bridges( prot_rdk, lig_rdk, distances, cutoff=self.cutoff), nb_channel=1) ] class CationPiVoxelizer(ComplexFeaturizer): """Localize cation-Pi interactions between atoms in macromolecular complexes. Given a macromolecular compelx made up of multiple constitutent molecules, compute cation-Pi between atoms in the macromolecular complex. For each atom, localize this salt bridge in the voxel in which it originated to create a local cation-pi array. Creates a tensor output of shape (box_width, box_width, box_width, 1) for each macromolecular complex that computes the number of cation-pi interactions at each voxel. """ def __init__(self, distance_cutoff=6.5, angle_cutoff=30.0, box_width=16.0, voxel_width=1.0): """ Parameters ---------- distance_cutoff: float, optional (default 6.5) The distance in angstroms within which atoms must be to be considered for a cation-pi interaction between them. angle_cutoff: float, optional (default 30.0) Angle cutoff. Max allowed deviation from the ideal (0deg) angle between ring normal and vector pointing from ring center to cation (in degrees). box_width: float, optional (default 16.0) Size of a box in which voxel features are calculated. Box is centered on a ligand centroid. voxel_width: float, optional (default 1.0) Size of a 3D voxel in a grid. """ self.distance_cutoff = distance_cutoff self.angle_cutoff = angle_cutoff self.box_width = box_width self.voxel_width = voxel_width self.voxels_per_edge = int(self.box_width / self.voxel_width) def _featurize_complex(self, mol, protein): """ Compute featurization for a single mol/protein complex Parameters ---------- mol: object Representation of the molecule protein: object Representation of the protein """ (lig_xyz, lig_rdk), (prot_xyz, prot_rdk) = mol, protein distances = compute_pairwise_distances(prot_xyz, lig_xyz) return [ sum([ voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, xyz, feature_dict=cation_pi_dict, nb_channel=1) for xyz, cation_pi_dict in zip( (prot_xyz, lig_xyz), compute_binding_pocket_cation_pi( prot_rdk, lig_rdk, dist_cutoff=self.distance_cutoff, angle_cutoff=self.angle_cutoff, )) ]) ] class PiStackVoxelizer(ComplexFeaturizer): """Localize Pi stacking interactions between atoms in macromolecular complexes. Given a macromolecular compelx made up of multiple constitutent molecules, compute pi-stacking interactions between atoms in the macromolecular complex. For each atom, localize this salt bridge in the voxel in which it originated to create a local pi-stacking array. Creates a tensor output of shape (box_width, box_width, box_width, 2) for each macromolecular complex. Each voxel has 2 fields, with the first tracking the number of pi-pi parallel interactions, and the second tracking the number of pi-T interactions. """ def __init__(self, distance_cutoff=4.4, angle_cutoff=30.0, box_width=16.0, voxel_width=1.0): """ Parameters ---------- distance_cutoff: float, optional (default 4.4) The distance in angstroms within which atoms must be to be considered for a cation-pi interaction between them. angle_cutoff: float, optional (default 30.0) Angle cutoff. Max allowed deviation from the ideal (0 deg) angle between ring normal and vector pointing from ring center to other ring center (in degrees). box_width: float, optional (default 16.0) Size of a box in which voxel features are calculated. Box is centered on a ligand centroid. voxel_width: float, optional (default 1.0) Size of a 3D voxel in a grid. """ self.distance_cutoff = distance_cutoff self.angle_cutoff = angle_cutoff self.box_width = box_width self.voxel_width = voxel_width self.voxels_per_edge = int(self.box_width / self.voxel_width) def _featurize_complex(self, mol, protein): """ Compute featurization for a single mol/protein complex Parameters ---------- mol: object Representation of the molecule protein: object Representation of the protein """ (lig_xyz, lig_rdk), (prot_xyz, prot_rdk) = mol, protein distances = compute_pairwise_distances(prot_xyz, lig_xyz) protein_pi_t, protein_pi_parallel, ligand_pi_t, ligand_pi_parallel = ( compute_pi_stack( prot_rdk, lig_rdk, distances, dist_cutoff=self.distance_cutoff, angle_cutoff=self.angle_cutoff)) pi_parallel_tensor = voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, prot_xyz, feature_dict=protein_pi_parallel, nb_channel=1) pi_parallel_tensor += voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, lig_xyz, feature_dict=ligand_pi_parallel, nb_channel=1) pi_t_tensor = voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, prot_xyz, feature_dict=protein_pi_t, nb_channel=1) pi_t_tensor += voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, lig_xyz, feature_dict=ligand_pi_t, nb_channel=1) return [pi_parallel_tensor, pi_t_tensor] deepchem/feat/rdkit_grid_featurizer.py +34 −185 Original line number Diff line number Diff line Loading @@ -16,6 +16,11 @@ from deepchem.feat.contact_fingerprints import ContactCircularFingerprint from deepchem.feat.contact_fingerprints import ContactCircularVoxelizer from deepchem.feat.splif_fingerprint import SplifFingerprint from deepchem.feat.splif_fingerprint import SplifVoxelizer from deepchem.feat.grid_featurizers import compute_charge_dictionary from deepchem.feat.grid_featurizers import ChargeVoxelizer from deepchem.feat.grid_featurizers import SaltBridgeVoxelizer from deepchem.feat.grid_featurizers import CationPiVoxelizer from deepchem.feat.grid_featurizers import PiStackVoxelizer from deepchem.utils.rdkit_util import load_molecule from deepchem.utils.rdkit_util import compute_centroid from deepchem.utils.rdkit_util import subtract_centroid Loading @@ -42,67 +47,15 @@ logger = logging.getLogger(__name__) http://stackoverflow.com/questions/2827393/angles-between-two-n-dimensional-vectors-in-python """ sybyl_types = [ "C3", "C2", "C1", "Cac", "Car", "N3", "N3+", "Npl", "N2", "N1", "Ng+", "Nox", "Nar", "Ntr", "Nam", "Npl3", "N4", "O3", "O-", "O2", "O.co2", "O.spc", "O.t3p", "S3", "S3+", "S2", "So2", "Sox" "Sac" "SO", "P3", "P", "P3+", "F", "Cl", "Br", "I" ] # TODO(rbharath): Consider this comment on rdkit forums https://github.com/rdkit/rdkit/issues/1590 about sybyl featurization FLAT_FEATURES = ['ecfp_ligand', 'ecfp_hashed', 'splif_hashed', 'hbond_count'] VOXEL_FEATURES = [ 'ecfp', 'splif', 'sybyl', 'salt_bridge', 'charge', 'hbond', 'pi_stack', 'cation_pi' 'ecfp', 'splif', 'salt_bridge', 'charge', 'hbond', 'pi_stack', 'cation_pi' ] def hash_sybyl(sybyl, sybyl_types): return (sybyl_types.index(sybyl)) def compute_all_sybyl(mol, indices=None): """Computes Sybyl atom types for atoms in molecule.""" raise NotImplementedError("This function is not implemented yet") def featurize_binding_pocket_sybyl(protein_xyz, protein, ligand_xyz, ligand, pairwise_distances=None, cutoff=7.0): """Computes Sybyl dicts for ligand and binding pocket of the protein. TODO(rbharath): Consider this comment on rdkit forums https://github.com/rdkit/rdkit/issues/1590 Parameters ---------- protein_xyz: np.ndarray Of shape (N_protein_atoms, 3) protein: Rdkit Molecule Contains more metadata. ligand_xyz: np.ndarray Of shape (N_ligand_atoms, 3) ligand: Rdkit Molecule Contains more metadata pairwise_distances: np.ndarray Array of pairwise protein-ligand distances (Angstroms) cutoff: float Cutoff distance for contact consideration. """ features_dict = {} if pairwise_distances is None: pairwise_distances = compute_pairwise_distances(protein_xyz, ligand_xyz) contacts = np.nonzero((pairwise_distances < cutoff)) protein_atoms = set([int(c) for c in contacts[0].tolist()]) protein_sybyl_dict = compute_all_sybyl(protein, indices=protein_atoms) ligand_sybyl_dict = compute_all_sybyl(ligand) return (protein_sybyl_dict, ligand_sybyl_dict) def compute_hbonds_in_range(protein, protein_xyz, ligand, ligand_xyz, pairwise_distances, hbond_dist_bin, hbond_angle_cutoff): Loading Loading @@ -201,7 +154,7 @@ class RdkitGridFeaturizer(ComplexFeaturizer): hbond_angle_cutoffs: [5, 50, 90] splif_contact_bins: [(0, 2.0), (2.0, 3.0), (3.0, 4.5)] ecfp_cutoff: 4.5 sybyl_cutoff: 7.0 #sybyl_cutoff: 7.0 salt_bridges_cutoff: 5.0 pi_stack_dist_cutoff: 4.4 pi_stack_angle_cutoff: 30.0 Loading @@ -212,16 +165,11 @@ class RdkitGridFeaturizer(ComplexFeaturizer): # list of features that require sanitized molecules require_sanitized = ['pi_stack', 'cation_pi', 'ecfp_ligand'] # not implemented featurization types not_implemented = ['sybyl'] self.sanitize = sanitize self.flatten = flatten self.ecfp_degree = ecfp_degree self.ecfp_power = ecfp_power self.splif_power = splif_power self.nb_rotations = nb_rotations # default values Loading @@ -230,7 +178,6 @@ class RdkitGridFeaturizer(ComplexFeaturizer): 'hbond_angle_cutoffs': [5, 50, 90], 'splif_contact_bins': [(0, 2.0), (2.0, 3.0), (3.0, 4.5)], 'ecfp_cutoff': 4.5, 'sybyl_cutoff': 7.0, 'salt_bridges_cutoff': 5.0, 'pi_stack_dist_cutoff': 4.4, 'pi_stack_angle_cutoff': 30.0, Loading Loading @@ -259,7 +206,7 @@ class RdkitGridFeaturizer(ComplexFeaturizer): ignored_features = [] if self.sanitize is False: ignored_features += require_sanitized ignored_features += not_implemented #ignored_features += not_implemented # Intantiate Featurizers self.ecfp_featurizer = CircularFingerprint( Loading @@ -284,6 +231,22 @@ class RdkitGridFeaturizer(ComplexFeaturizer): size=2**self.splif_power, box_width=self.box_width, voxel_width=self.voxel_width) self.charge_voxelizer = ChargeVoxelizer( box_width=self.box_width, voxel_width=self.voxel_width) self.salt_bridge_voxelizer = SaltBridgeVoxelizer( cutoff=self.cutoffs['salt_bridges_cutoff'], box_width=self.box_width, voxel_width=self.voxel_width) self.cation_pi_voxelizer = CationPiVoxelizer( distance_cutoff=self.cutoffs['cation_pi_dist_cutoff'], angle_cutoff=self.cutoffs['cation_pi_angle_cutoff'], box_width=self.box_width, voxel_width=self.voxel_width) self.pi_stack_voxelizer = PiStackVoxelizer( distance_cutoff=self.cutoffs['pi_stack_dist_cutoff'], angle_cutoff=self.cutoffs['pi_stack_angle_cutoff'], box_width=self.box_width, voxel_width=self.voxel_width) # parse provided feature types for feature_type in feature_types: Loading @@ -291,10 +254,6 @@ class RdkitGridFeaturizer(ComplexFeaturizer): logger.warning('sanitize is set to False, %s feature will be ignored' % feature_type) continue if feature_type in not_implemented: logger.warning('%s feature is not implemented yet and will be ignored' % feature_type) continue if feature_type in FLAT_FEATURES: self.feature_types.append((True, feature_type)) Loading Loading @@ -354,57 +313,12 @@ class RdkitGridFeaturizer(ComplexFeaturizer): if feature_name == 'splif': return self.splif_voxelizer._featurize_complex((lig_xyz, lig_rdk), (prot_xyz, prot_rdk)) if feature_name == 'sybyl': return [ voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, lambda x: hash_sybyl(x, sybyl_types=sybyl_types), xyz, feature_dict=sybyl_dict, nb_channel=len(sybyl_types)) for xyz, sybyl_dict in zip((prot_xyz, lig_xyz), featurize_binding_pocket_sybyl( prot_xyz, prot_rdk, lig_xyz, lig_rdk, distances, cutoff=self.cutoffs['sybyl_cutoff'])) ] if feature_name == 'salt_bridge': return [ voxelize( convert_atom_pair_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, (prot_xyz, lig_xyz), feature_list=compute_salt_bridges( prot_rdk, lig_rdk, distances, cutoff=self.cutoffs['salt_bridges_cutoff']), nb_channel=1) ] return self.salt_bridge_voxelizer._featurize_complex((lig_xyz, lig_rdk), (prot_xyz, prot_rdk)) if feature_name == 'charge': return [ sum([ voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, xyz, feature_dict=compute_charge_dictionary(mol), nb_channel=1, dtype="np.float16") for xyz, mol in ((prot_xyz, prot_rdk), (lig_xyz, lig_rdk)) ]) ] return self.charge_voxelizer._featurize_complex((lig_xyz, lig_rdk), (prot_xyz, prot_rdk)) if feature_name == 'hbond': return [ voxelize( Loading @@ -414,34 +328,16 @@ class RdkitGridFeaturizer(ComplexFeaturizer): self.voxel_width, None, (prot_xyz, lig_xyz), feature_list=hbond_list, nb_channel=2**0) for hbond_list in compute_hydrogen_bonds( nb_channel=1) for hbond_list in compute_hydrogen_bonds( prot_xyz, prot_rdk, lig_xyz, lig_rdk, distances, self.cutoffs[ 'hbond_dist_bins'], self.cutoffs['hbond_angle_cutoffs']) ] if feature_name == 'pi_stack': return self.voxelize_pi_stack(prot_xyz, prot_rdk, lig_xyz, lig_rdk, distances) return self.pi_stack_voxelizer._featurize_complex((lig_xyz, lig_rdk), (prot_xyz, prot_rdk)) if feature_name == 'cation_pi': return [ sum([ voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, xyz, feature_dict=cation_pi_dict, nb_channel=1) for xyz, cation_pi_dict in zip( (prot_xyz, lig_xyz), compute_binding_pocket_cation_pi( prot_rdk, lig_rdk, dist_cutoff=self.cutoffs['cation_pi_dist_cutoff'], angle_cutoff=self.cutoffs['cation_pi_angle_cutoff'], )) ]) ] return self.cation_pi_voxelizer._featurize_complex((lig_xyz, lig_rdk), (prot_xyz, prot_rdk)) raise ValueError('Unknown feature type "%s"' % feature_name) def _featurize(self, mol_pdb_file, protein_pdb_file): Loading Loading @@ -524,50 +420,3 @@ class RdkitGridFeaturizer(ComplexFeaturizer): features = np.squeeze(np.array(list(features_dict.values()))) return features def voxelize_pi_stack(self, prot_xyz, prot_rdk, lig_xyz, lig_rdk, distances): protein_pi_t, protein_pi_parallel, ligand_pi_t, ligand_pi_parallel = ( compute_pi_stack( prot_rdk, lig_rdk, distances, dist_cutoff=self.cutoffs['pi_stack_dist_cutoff'], angle_cutoff=self.cutoffs['pi_stack_angle_cutoff'])) pi_parallel_tensor = voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, prot_xyz, feature_dict=protein_pi_parallel, nb_channel=1) pi_parallel_tensor += voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, lig_xyz, feature_dict=ligand_pi_parallel, nb_channel=1) pi_t_tensor = voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, prot_xyz, feature_dict=protein_pi_t, nb_channel=1) pi_t_tensor += voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, lig_xyz, feature_dict=ligand_pi_t, nb_channel=1) return [pi_parallel_tensor, pi_t_tensor] deepchem/feat/tests/test_rdkit_grid_features.py +10 −8 Original line number Diff line number Diff line Loading @@ -11,6 +11,8 @@ np.random.seed(123) from deepchem.utils import rdkit_util from deepchem.utils import hash_utils from deepchem.utils import voxel_utils from deepchem.feat.contact_fingerprints import featurize_binding_pocket_ecfp from deepchem.feat.splif_fingerprint import featurize_splif from deepchem.feat import rdkit_grid_featurizer as rgf from deepchem.feat.splif_fingerprint import compute_splif_features_in_range Loading Loading @@ -124,13 +126,13 @@ class TestFeaturizationFunctions(unittest.TestCase): distance = rdkit_util.compute_pairwise_distances(prot_xyz, lig_xyz) # check if results are the same if we provide precomputed distances prot_dict, lig_dict = rgf.featurize_binding_pocket_ecfp( prot_dict, lig_dict = featurize_binding_pocket_ecfp( prot_xyz, prot_rdk, lig_xyz, lig_rdk, ) prot_dict_dist, lig_dict_dist = rgf.featurize_binding_pocket_ecfp( prot_dict_dist, lig_dict_dist = featurize_binding_pocket_ecfp( prot_xyz, prot_rdk, lig_xyz, lig_rdk, pairwise_distances=distance) # ...but first check if we actually got two dicts self.assertIsInstance(prot_dict, dict) Loading @@ -140,14 +142,14 @@ class TestFeaturizationFunctions(unittest.TestCase): self.assertEqual(lig_dict, lig_dict_dist) # check if we get less features with smaller distance cutoff prot_dict_d2, lig_dict_d2 = rgf.featurize_binding_pocket_ecfp( prot_dict_d2, lig_dict_d2 = featurize_binding_pocket_ecfp( prot_xyz, prot_rdk, lig_xyz, lig_rdk, cutoff=2.0, ) prot_dict_d6, lig_dict_d6 = rgf.featurize_binding_pocket_ecfp( prot_dict_d6, lig_dict_d6 = featurize_binding_pocket_ecfp( prot_xyz, prot_rdk, lig_xyz, Loading @@ -161,7 +163,7 @@ class TestFeaturizationFunctions(unittest.TestCase): self.assertGreaterEqual(len(lig_dict_d6), len(lig_dict)) # check if using different ecfp_degree changes anything prot_dict_e3, lig_dict_e3 = rgf.featurize_binding_pocket_ecfp( prot_dict_e3, lig_dict_e3 = featurize_binding_pocket_ecfp( prot_xyz, prot_rdk, lig_xyz, Loading Loading @@ -209,7 +211,7 @@ class TestFeaturizationFunctions(unittest.TestCase): bins = [(1, 2), (2, 3)] dicts = rgf.featurize_splif( dicts = featurize_splif( prot_xyz, prot_rdk, lig_xyz, Loading Loading @@ -329,7 +331,7 @@ class TestRdkitGridFeaturizer(unittest.TestCase): 'hbond_angle_cutoffs': [5, 90], 'splif_contact_bins': [(0, 3.5), (3.5, 6.0)], 'ecfp_cutoff': 5.0, 'sybyl_cutoff': 3.0, #'sybyl_cutoff': 3.0, 'salt_bridges_cutoff': 4.0, 'pi_stack_dist_cutoff': 5.0, 'pi_stack_angle_cutoff': 15.0, Loading Loading @@ -357,7 +359,7 @@ class TestRdkitGridFeaturizer(unittest.TestCase): prot_xyz = rgf.subtract_centroid(prot_xyz, centroid) lig_xyz = rgf.subtract_centroid(lig_xyz, centroid) prot_ecfp_dict, lig_ecfp_dict = rgf.featurize_binding_pocket_ecfp( prot_ecfp_dict, lig_ecfp_dict = featurize_binding_pocket_ecfp( prot_xyz, prot_rdk, lig_xyz, lig_rdk) box_w = 20 Loading deepchem/utils/rdkit_util.py +119 −69 File changed.Preview size limit exceeded, changes collapsed. Show changes deepchem/feat/contact_fingerprints.py +1 −1 File changed.Contains only whitespace changes. Show changes Loading
deepchem/feat/grid_featurizers.py +246 −2 Original line number Diff line number Diff line Loading @@ -5,6 +5,12 @@ import logging from deepchem.utils.rdkit_util import get_partial_charge from deepchem.feat import ComplexFeaturizer from deepchem.utils.voxel_utils import voxelize from deepchem.utils.rdkit_util import compute_salt_bridges from deepchem.utils.rdkit_util import compute_binding_pocket_cation_pi from deepchem.utils.voxel_utils import convert_atom_to_voxel from deepchem.utils.voxel_utils import convert_atom_pair_to_voxel from deepchem.utils.rdkit_util import compute_pairwise_distances from deepchem.utils.rdkit_util import compute_pi_stack logger = logging.getLogger(__name__) Loading @@ -28,7 +34,7 @@ class ChargeVoxelizer(ComplexFeaturizer): constitutent molecules, compute the partial (Gasteiger charge) on each molecule. For each atom, localize this partial charge in the voxel in which it originated to create a local charge array. Sum contributes to get an effective a local charge array. Sum contributions to get an effective charge at each voxel. Creates a tensor output of shape (box_width, box_width, Loading @@ -47,7 +53,6 @@ class ChargeVoxelizer(ComplexFeaturizer): voxel_width: float, optional (default 1.0) Size of a 3D voxel in a grid. """ self.cutoff = cutoff self.box_width = box_width self.voxel_width = voxel_width self.voxels_per_edge = int(self.box_width / self.voxel_width) Loading Loading @@ -79,3 +84,242 @@ class ChargeVoxelizer(ComplexFeaturizer): for xyz, mol in ((prot_xyz, prot_rdk), (lig_xyz, lig_rdk)) ]) ] class SaltBridgeVoxelizer(ComplexFeaturizer): """Localize salt bridges between atoms in macromolecular complexes. Given a macromolecular compelx made up of multiple constitutent molecules, compute salt bridges between atoms in the macromolecular complex. For each atom, localize this salt bridge in the voxel in which it originated to create a local salt bridge array. Note that if atoms in two different voxels interact in a salt-bridge, the interaction is double counted in both voxels. Creates a tensor output of shape (box_width, box_width, box_width, 1) for each macromolecular complex that computes the number of salt bridges at each voxel. """ def __init__(self, cutoff=5.0, box_width=16.0, voxel_width=1.0): """ Parameters ---------- cutoff: float, optional (default 5.0) The distance in angstroms within which atoms must be to be considered for a salt bridge between them. box_width: float, optional (default 16.0) Size of a box in which voxel features are calculated. Box is centered on a ligand centroid. voxel_width: float, optional (default 1.0) Size of a 3D voxel in a grid. """ self.cutoff = cutoff self.box_width = box_width self.voxel_width = voxel_width self.voxels_per_edge = int(self.box_width / self.voxel_width) def _featurize_complex(self, mol, protein): """ Compute featurization for a single mol/protein complex Parameters ---------- mol: object Representation of the molecule protein: object Representation of the protein """ (lig_xyz, lig_rdk), (prot_xyz, prot_rdk) = mol, protein distances = compute_pairwise_distances(prot_xyz, lig_xyz) return [ voxelize( convert_atom_pair_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, (prot_xyz, lig_xyz), feature_list=compute_salt_bridges( prot_rdk, lig_rdk, distances, cutoff=self.cutoff), nb_channel=1) ] class CationPiVoxelizer(ComplexFeaturizer): """Localize cation-Pi interactions between atoms in macromolecular complexes. Given a macromolecular compelx made up of multiple constitutent molecules, compute cation-Pi between atoms in the macromolecular complex. For each atom, localize this salt bridge in the voxel in which it originated to create a local cation-pi array. Creates a tensor output of shape (box_width, box_width, box_width, 1) for each macromolecular complex that computes the number of cation-pi interactions at each voxel. """ def __init__(self, distance_cutoff=6.5, angle_cutoff=30.0, box_width=16.0, voxel_width=1.0): """ Parameters ---------- distance_cutoff: float, optional (default 6.5) The distance in angstroms within which atoms must be to be considered for a cation-pi interaction between them. angle_cutoff: float, optional (default 30.0) Angle cutoff. Max allowed deviation from the ideal (0deg) angle between ring normal and vector pointing from ring center to cation (in degrees). box_width: float, optional (default 16.0) Size of a box in which voxel features are calculated. Box is centered on a ligand centroid. voxel_width: float, optional (default 1.0) Size of a 3D voxel in a grid. """ self.distance_cutoff = distance_cutoff self.angle_cutoff = angle_cutoff self.box_width = box_width self.voxel_width = voxel_width self.voxels_per_edge = int(self.box_width / self.voxel_width) def _featurize_complex(self, mol, protein): """ Compute featurization for a single mol/protein complex Parameters ---------- mol: object Representation of the molecule protein: object Representation of the protein """ (lig_xyz, lig_rdk), (prot_xyz, prot_rdk) = mol, protein distances = compute_pairwise_distances(prot_xyz, lig_xyz) return [ sum([ voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, xyz, feature_dict=cation_pi_dict, nb_channel=1) for xyz, cation_pi_dict in zip( (prot_xyz, lig_xyz), compute_binding_pocket_cation_pi( prot_rdk, lig_rdk, dist_cutoff=self.distance_cutoff, angle_cutoff=self.angle_cutoff, )) ]) ] class PiStackVoxelizer(ComplexFeaturizer): """Localize Pi stacking interactions between atoms in macromolecular complexes. Given a macromolecular compelx made up of multiple constitutent molecules, compute pi-stacking interactions between atoms in the macromolecular complex. For each atom, localize this salt bridge in the voxel in which it originated to create a local pi-stacking array. Creates a tensor output of shape (box_width, box_width, box_width, 2) for each macromolecular complex. Each voxel has 2 fields, with the first tracking the number of pi-pi parallel interactions, and the second tracking the number of pi-T interactions. """ def __init__(self, distance_cutoff=4.4, angle_cutoff=30.0, box_width=16.0, voxel_width=1.0): """ Parameters ---------- distance_cutoff: float, optional (default 4.4) The distance in angstroms within which atoms must be to be considered for a cation-pi interaction between them. angle_cutoff: float, optional (default 30.0) Angle cutoff. Max allowed deviation from the ideal (0 deg) angle between ring normal and vector pointing from ring center to other ring center (in degrees). box_width: float, optional (default 16.0) Size of a box in which voxel features are calculated. Box is centered on a ligand centroid. voxel_width: float, optional (default 1.0) Size of a 3D voxel in a grid. """ self.distance_cutoff = distance_cutoff self.angle_cutoff = angle_cutoff self.box_width = box_width self.voxel_width = voxel_width self.voxels_per_edge = int(self.box_width / self.voxel_width) def _featurize_complex(self, mol, protein): """ Compute featurization for a single mol/protein complex Parameters ---------- mol: object Representation of the molecule protein: object Representation of the protein """ (lig_xyz, lig_rdk), (prot_xyz, prot_rdk) = mol, protein distances = compute_pairwise_distances(prot_xyz, lig_xyz) protein_pi_t, protein_pi_parallel, ligand_pi_t, ligand_pi_parallel = ( compute_pi_stack( prot_rdk, lig_rdk, distances, dist_cutoff=self.distance_cutoff, angle_cutoff=self.angle_cutoff)) pi_parallel_tensor = voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, prot_xyz, feature_dict=protein_pi_parallel, nb_channel=1) pi_parallel_tensor += voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, lig_xyz, feature_dict=ligand_pi_parallel, nb_channel=1) pi_t_tensor = voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, prot_xyz, feature_dict=protein_pi_t, nb_channel=1) pi_t_tensor += voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, lig_xyz, feature_dict=ligand_pi_t, nb_channel=1) return [pi_parallel_tensor, pi_t_tensor]
deepchem/feat/rdkit_grid_featurizer.py +34 −185 Original line number Diff line number Diff line Loading @@ -16,6 +16,11 @@ from deepchem.feat.contact_fingerprints import ContactCircularFingerprint from deepchem.feat.contact_fingerprints import ContactCircularVoxelizer from deepchem.feat.splif_fingerprint import SplifFingerprint from deepchem.feat.splif_fingerprint import SplifVoxelizer from deepchem.feat.grid_featurizers import compute_charge_dictionary from deepchem.feat.grid_featurizers import ChargeVoxelizer from deepchem.feat.grid_featurizers import SaltBridgeVoxelizer from deepchem.feat.grid_featurizers import CationPiVoxelizer from deepchem.feat.grid_featurizers import PiStackVoxelizer from deepchem.utils.rdkit_util import load_molecule from deepchem.utils.rdkit_util import compute_centroid from deepchem.utils.rdkit_util import subtract_centroid Loading @@ -42,67 +47,15 @@ logger = logging.getLogger(__name__) http://stackoverflow.com/questions/2827393/angles-between-two-n-dimensional-vectors-in-python """ sybyl_types = [ "C3", "C2", "C1", "Cac", "Car", "N3", "N3+", "Npl", "N2", "N1", "Ng+", "Nox", "Nar", "Ntr", "Nam", "Npl3", "N4", "O3", "O-", "O2", "O.co2", "O.spc", "O.t3p", "S3", "S3+", "S2", "So2", "Sox" "Sac" "SO", "P3", "P", "P3+", "F", "Cl", "Br", "I" ] # TODO(rbharath): Consider this comment on rdkit forums https://github.com/rdkit/rdkit/issues/1590 about sybyl featurization FLAT_FEATURES = ['ecfp_ligand', 'ecfp_hashed', 'splif_hashed', 'hbond_count'] VOXEL_FEATURES = [ 'ecfp', 'splif', 'sybyl', 'salt_bridge', 'charge', 'hbond', 'pi_stack', 'cation_pi' 'ecfp', 'splif', 'salt_bridge', 'charge', 'hbond', 'pi_stack', 'cation_pi' ] def hash_sybyl(sybyl, sybyl_types): return (sybyl_types.index(sybyl)) def compute_all_sybyl(mol, indices=None): """Computes Sybyl atom types for atoms in molecule.""" raise NotImplementedError("This function is not implemented yet") def featurize_binding_pocket_sybyl(protein_xyz, protein, ligand_xyz, ligand, pairwise_distances=None, cutoff=7.0): """Computes Sybyl dicts for ligand and binding pocket of the protein. TODO(rbharath): Consider this comment on rdkit forums https://github.com/rdkit/rdkit/issues/1590 Parameters ---------- protein_xyz: np.ndarray Of shape (N_protein_atoms, 3) protein: Rdkit Molecule Contains more metadata. ligand_xyz: np.ndarray Of shape (N_ligand_atoms, 3) ligand: Rdkit Molecule Contains more metadata pairwise_distances: np.ndarray Array of pairwise protein-ligand distances (Angstroms) cutoff: float Cutoff distance for contact consideration. """ features_dict = {} if pairwise_distances is None: pairwise_distances = compute_pairwise_distances(protein_xyz, ligand_xyz) contacts = np.nonzero((pairwise_distances < cutoff)) protein_atoms = set([int(c) for c in contacts[0].tolist()]) protein_sybyl_dict = compute_all_sybyl(protein, indices=protein_atoms) ligand_sybyl_dict = compute_all_sybyl(ligand) return (protein_sybyl_dict, ligand_sybyl_dict) def compute_hbonds_in_range(protein, protein_xyz, ligand, ligand_xyz, pairwise_distances, hbond_dist_bin, hbond_angle_cutoff): Loading Loading @@ -201,7 +154,7 @@ class RdkitGridFeaturizer(ComplexFeaturizer): hbond_angle_cutoffs: [5, 50, 90] splif_contact_bins: [(0, 2.0), (2.0, 3.0), (3.0, 4.5)] ecfp_cutoff: 4.5 sybyl_cutoff: 7.0 #sybyl_cutoff: 7.0 salt_bridges_cutoff: 5.0 pi_stack_dist_cutoff: 4.4 pi_stack_angle_cutoff: 30.0 Loading @@ -212,16 +165,11 @@ class RdkitGridFeaturizer(ComplexFeaturizer): # list of features that require sanitized molecules require_sanitized = ['pi_stack', 'cation_pi', 'ecfp_ligand'] # not implemented featurization types not_implemented = ['sybyl'] self.sanitize = sanitize self.flatten = flatten self.ecfp_degree = ecfp_degree self.ecfp_power = ecfp_power self.splif_power = splif_power self.nb_rotations = nb_rotations # default values Loading @@ -230,7 +178,6 @@ class RdkitGridFeaturizer(ComplexFeaturizer): 'hbond_angle_cutoffs': [5, 50, 90], 'splif_contact_bins': [(0, 2.0), (2.0, 3.0), (3.0, 4.5)], 'ecfp_cutoff': 4.5, 'sybyl_cutoff': 7.0, 'salt_bridges_cutoff': 5.0, 'pi_stack_dist_cutoff': 4.4, 'pi_stack_angle_cutoff': 30.0, Loading Loading @@ -259,7 +206,7 @@ class RdkitGridFeaturizer(ComplexFeaturizer): ignored_features = [] if self.sanitize is False: ignored_features += require_sanitized ignored_features += not_implemented #ignored_features += not_implemented # Intantiate Featurizers self.ecfp_featurizer = CircularFingerprint( Loading @@ -284,6 +231,22 @@ class RdkitGridFeaturizer(ComplexFeaturizer): size=2**self.splif_power, box_width=self.box_width, voxel_width=self.voxel_width) self.charge_voxelizer = ChargeVoxelizer( box_width=self.box_width, voxel_width=self.voxel_width) self.salt_bridge_voxelizer = SaltBridgeVoxelizer( cutoff=self.cutoffs['salt_bridges_cutoff'], box_width=self.box_width, voxel_width=self.voxel_width) self.cation_pi_voxelizer = CationPiVoxelizer( distance_cutoff=self.cutoffs['cation_pi_dist_cutoff'], angle_cutoff=self.cutoffs['cation_pi_angle_cutoff'], box_width=self.box_width, voxel_width=self.voxel_width) self.pi_stack_voxelizer = PiStackVoxelizer( distance_cutoff=self.cutoffs['pi_stack_dist_cutoff'], angle_cutoff=self.cutoffs['pi_stack_angle_cutoff'], box_width=self.box_width, voxel_width=self.voxel_width) # parse provided feature types for feature_type in feature_types: Loading @@ -291,10 +254,6 @@ class RdkitGridFeaturizer(ComplexFeaturizer): logger.warning('sanitize is set to False, %s feature will be ignored' % feature_type) continue if feature_type in not_implemented: logger.warning('%s feature is not implemented yet and will be ignored' % feature_type) continue if feature_type in FLAT_FEATURES: self.feature_types.append((True, feature_type)) Loading Loading @@ -354,57 +313,12 @@ class RdkitGridFeaturizer(ComplexFeaturizer): if feature_name == 'splif': return self.splif_voxelizer._featurize_complex((lig_xyz, lig_rdk), (prot_xyz, prot_rdk)) if feature_name == 'sybyl': return [ voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, lambda x: hash_sybyl(x, sybyl_types=sybyl_types), xyz, feature_dict=sybyl_dict, nb_channel=len(sybyl_types)) for xyz, sybyl_dict in zip((prot_xyz, lig_xyz), featurize_binding_pocket_sybyl( prot_xyz, prot_rdk, lig_xyz, lig_rdk, distances, cutoff=self.cutoffs['sybyl_cutoff'])) ] if feature_name == 'salt_bridge': return [ voxelize( convert_atom_pair_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, (prot_xyz, lig_xyz), feature_list=compute_salt_bridges( prot_rdk, lig_rdk, distances, cutoff=self.cutoffs['salt_bridges_cutoff']), nb_channel=1) ] return self.salt_bridge_voxelizer._featurize_complex((lig_xyz, lig_rdk), (prot_xyz, prot_rdk)) if feature_name == 'charge': return [ sum([ voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, xyz, feature_dict=compute_charge_dictionary(mol), nb_channel=1, dtype="np.float16") for xyz, mol in ((prot_xyz, prot_rdk), (lig_xyz, lig_rdk)) ]) ] return self.charge_voxelizer._featurize_complex((lig_xyz, lig_rdk), (prot_xyz, prot_rdk)) if feature_name == 'hbond': return [ voxelize( Loading @@ -414,34 +328,16 @@ class RdkitGridFeaturizer(ComplexFeaturizer): self.voxel_width, None, (prot_xyz, lig_xyz), feature_list=hbond_list, nb_channel=2**0) for hbond_list in compute_hydrogen_bonds( nb_channel=1) for hbond_list in compute_hydrogen_bonds( prot_xyz, prot_rdk, lig_xyz, lig_rdk, distances, self.cutoffs[ 'hbond_dist_bins'], self.cutoffs['hbond_angle_cutoffs']) ] if feature_name == 'pi_stack': return self.voxelize_pi_stack(prot_xyz, prot_rdk, lig_xyz, lig_rdk, distances) return self.pi_stack_voxelizer._featurize_complex((lig_xyz, lig_rdk), (prot_xyz, prot_rdk)) if feature_name == 'cation_pi': return [ sum([ voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, xyz, feature_dict=cation_pi_dict, nb_channel=1) for xyz, cation_pi_dict in zip( (prot_xyz, lig_xyz), compute_binding_pocket_cation_pi( prot_rdk, lig_rdk, dist_cutoff=self.cutoffs['cation_pi_dist_cutoff'], angle_cutoff=self.cutoffs['cation_pi_angle_cutoff'], )) ]) ] return self.cation_pi_voxelizer._featurize_complex((lig_xyz, lig_rdk), (prot_xyz, prot_rdk)) raise ValueError('Unknown feature type "%s"' % feature_name) def _featurize(self, mol_pdb_file, protein_pdb_file): Loading Loading @@ -524,50 +420,3 @@ class RdkitGridFeaturizer(ComplexFeaturizer): features = np.squeeze(np.array(list(features_dict.values()))) return features def voxelize_pi_stack(self, prot_xyz, prot_rdk, lig_xyz, lig_rdk, distances): protein_pi_t, protein_pi_parallel, ligand_pi_t, ligand_pi_parallel = ( compute_pi_stack( prot_rdk, lig_rdk, distances, dist_cutoff=self.cutoffs['pi_stack_dist_cutoff'], angle_cutoff=self.cutoffs['pi_stack_angle_cutoff'])) pi_parallel_tensor = voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, prot_xyz, feature_dict=protein_pi_parallel, nb_channel=1) pi_parallel_tensor += voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, lig_xyz, feature_dict=ligand_pi_parallel, nb_channel=1) pi_t_tensor = voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, prot_xyz, feature_dict=protein_pi_t, nb_channel=1) pi_t_tensor += voxelize( convert_atom_to_voxel, self.voxels_per_edge, self.box_width, self.voxel_width, None, lig_xyz, feature_dict=ligand_pi_t, nb_channel=1) return [pi_parallel_tensor, pi_t_tensor]
deepchem/feat/tests/test_rdkit_grid_features.py +10 −8 Original line number Diff line number Diff line Loading @@ -11,6 +11,8 @@ np.random.seed(123) from deepchem.utils import rdkit_util from deepchem.utils import hash_utils from deepchem.utils import voxel_utils from deepchem.feat.contact_fingerprints import featurize_binding_pocket_ecfp from deepchem.feat.splif_fingerprint import featurize_splif from deepchem.feat import rdkit_grid_featurizer as rgf from deepchem.feat.splif_fingerprint import compute_splif_features_in_range Loading Loading @@ -124,13 +126,13 @@ class TestFeaturizationFunctions(unittest.TestCase): distance = rdkit_util.compute_pairwise_distances(prot_xyz, lig_xyz) # check if results are the same if we provide precomputed distances prot_dict, lig_dict = rgf.featurize_binding_pocket_ecfp( prot_dict, lig_dict = featurize_binding_pocket_ecfp( prot_xyz, prot_rdk, lig_xyz, lig_rdk, ) prot_dict_dist, lig_dict_dist = rgf.featurize_binding_pocket_ecfp( prot_dict_dist, lig_dict_dist = featurize_binding_pocket_ecfp( prot_xyz, prot_rdk, lig_xyz, lig_rdk, pairwise_distances=distance) # ...but first check if we actually got two dicts self.assertIsInstance(prot_dict, dict) Loading @@ -140,14 +142,14 @@ class TestFeaturizationFunctions(unittest.TestCase): self.assertEqual(lig_dict, lig_dict_dist) # check if we get less features with smaller distance cutoff prot_dict_d2, lig_dict_d2 = rgf.featurize_binding_pocket_ecfp( prot_dict_d2, lig_dict_d2 = featurize_binding_pocket_ecfp( prot_xyz, prot_rdk, lig_xyz, lig_rdk, cutoff=2.0, ) prot_dict_d6, lig_dict_d6 = rgf.featurize_binding_pocket_ecfp( prot_dict_d6, lig_dict_d6 = featurize_binding_pocket_ecfp( prot_xyz, prot_rdk, lig_xyz, Loading @@ -161,7 +163,7 @@ class TestFeaturizationFunctions(unittest.TestCase): self.assertGreaterEqual(len(lig_dict_d6), len(lig_dict)) # check if using different ecfp_degree changes anything prot_dict_e3, lig_dict_e3 = rgf.featurize_binding_pocket_ecfp( prot_dict_e3, lig_dict_e3 = featurize_binding_pocket_ecfp( prot_xyz, prot_rdk, lig_xyz, Loading Loading @@ -209,7 +211,7 @@ class TestFeaturizationFunctions(unittest.TestCase): bins = [(1, 2), (2, 3)] dicts = rgf.featurize_splif( dicts = featurize_splif( prot_xyz, prot_rdk, lig_xyz, Loading Loading @@ -329,7 +331,7 @@ class TestRdkitGridFeaturizer(unittest.TestCase): 'hbond_angle_cutoffs': [5, 90], 'splif_contact_bins': [(0, 3.5), (3.5, 6.0)], 'ecfp_cutoff': 5.0, 'sybyl_cutoff': 3.0, #'sybyl_cutoff': 3.0, 'salt_bridges_cutoff': 4.0, 'pi_stack_dist_cutoff': 5.0, 'pi_stack_angle_cutoff': 15.0, Loading Loading @@ -357,7 +359,7 @@ class TestRdkitGridFeaturizer(unittest.TestCase): prot_xyz = rgf.subtract_centroid(prot_xyz, centroid) lig_xyz = rgf.subtract_centroid(lig_xyz, centroid) prot_ecfp_dict, lig_ecfp_dict = rgf.featurize_binding_pocket_ecfp( prot_ecfp_dict, lig_ecfp_dict = featurize_binding_pocket_ecfp( prot_xyz, prot_rdk, lig_xyz, lig_rdk) box_w = 20 Loading
deepchem/utils/rdkit_util.py +119 −69 File changed.Preview size limit exceeded, changes collapsed. Show changes
deepchem/feat/contact_fingerprints.py +1 −1 File changed.Contains only whitespace changes. Show changes
