Loading deepchem/feat/atomic_coordinates.py +90 −26 Original line number Diff line number Diff line Loading @@ -8,6 +8,7 @@ from deepchem.feat import MolecularFeaturizer from deepchem.feat import ComplexFeaturizer from deepchem.utils import rdkit_util, pad_array from deepchem.utils.rdkit_util import MoleculeLoadException from deepchem.utils.rdkit_util import reduce_molecular_complex_to_contacts logger = logging.getLogger(__name__) Loading Loading @@ -58,7 +59,19 @@ class AtomicCoordinates(MolecularFeaturizer): def _compute_neighbor_list(coords, neighbor_cutoff, max_num_neighbors, periodic_box_size): """Computes a neighbor list from atom coordinates.""" """Computes a neighbor list from atom coordinates. Parameters ---------- coords: Numpy array Of shape (N, 3) with all the atoms in this systems neighbor_cutoff: float The neighbor cutoff in angstroms max_num_neighbors: int The maximum number of neighbors per atom periodic_box_size: tuple With (x, y, z) box sizes per dimension """ N = coords.shape[0] import mdtraj traj = mdtraj.Trajectory(coords.reshape((1, N, 3)), None) Loading Loading @@ -89,6 +102,11 @@ def _compute_neighbor_list(coords, neighbor_cutoff, max_num_neighbors, def get_coords(mol): """ Gets coordinates in Angstrom for RDKit mol. Parameters ---------- mol: rdkit mol Molecule to get coordinates for """ N = mol.GetNumAtoms() coords = np.zeros((N, 3)) Loading Loading @@ -182,7 +200,7 @@ class NeighborListComplexAtomicCoordinates(ComplexFeaturizer): self.dtype = object self.coordinates_featurizer = AtomicCoordinates() def _featurize(self, mol_pdb_file, protein_pdb_file): def _featurize(self, molecular_complex): """ Compute neighbor list for complex. Loading @@ -191,9 +209,13 @@ class NeighborListComplexAtomicCoordinates(ComplexFeaturizer): molecular_complex: Object Some representation of a molecular complex. """ mol_coords, ob_mol = rdkit_util.load_molecule(mol_pdb_file) protein_coords, protein_mol = rdkit_util.load_molecule(protein_pdb_file) system_coords = rdkit_util.merge_molecules_xyz([mol_coords, protein_coords]) fragments = rdkit_util.load_complex(molecular_complex, add_hydrogens=False) #mol_coords, ob_mol = rdkit_util.load_molecule(mol_pdb_file) #protein_coords, protein_mol = rdkit_util.load_molecule(protein_pdb_file) coords = [frag[0] for frag in fragments] mols = [frag[1] for frag in fragments] #system_coords = rdkit_util.merge_molecules_xyz(mol_coords, protein_coords) system_coords = rdkit_util.merge_molecules_xyz(coords) system_neighbor_list = _compute_neighbor_list( system_coords, self.neighbor_cutoff, self.max_num_neighbors, None) Loading Loading @@ -223,31 +245,37 @@ class AtomicConvFeaturizer(ComplexFeaturizer): """ def __init__(self, frag_num_atoms, complex_num_atoms, frag_max_atoms, max_num_neighbors, neighbor_cutoff, reduce_to_contacts=True, cutoff=4.5, strip_hydrogens=True): """Initialize an AtomicConvFeaturizer object. Parameters ---------- frag_num_atoms: list[int] List of the number of atoms in each fragment. frag_max_atoms: int or list[int] List of the max number of atoms in each fragment. If int, assumes you're setting the same bound for all fragments. max_num_neighbors: int The maximum number of neighbors allowed neighbor_cutoff: float The distance in angstroms after which neighbors are cutoff. reduce_to_contacts: bool, optional If True, reduce the atoms in the complex to those near a contact region. cutoff: float The cutoff distance in angstroms. Only used if `reduce_to_contacts` is `True`. strip_hydrogens: bool, optional If true, remove hydrogens before featurizing. If True, remove hydrogens before featurizing. """ # TODO(rbharath): extend to more fragments if len(frag_num_atoms) != 2: raise ValueError("Currently only supports two fragments") self.frag_num_atoms = frag_num_atoms self.complex_num_atoms = sum(frag_num_atoms) self.frag_max_atoms = frag_max_atoms self.max_num_neighbors = max_num_neighbors self.neighbor_cutoff = neighbor_cutoff self.reduce_to_contacts = reduce_to_contacts self.cutoff = cutoff self.strip_hydrogens = strip_hydrogens def _featurize(self, molecular_complex): Loading @@ -258,6 +286,13 @@ class AtomicConvFeaturizer(ComplexFeaturizer): molecular_complex: Object Some representation of a molecular complex. """ # If our upper bound is an int, expand it to a list for each # fragment if isinstance(self.frag_max_atoms, int): frag_max_atoms = [self.frag_max_atoms] * len(molecular_complex) else: frag_max_atoms = self.frag_max_atoms complex_max_atoms = sum(frag_max_atoms) frag_coords = [] frag_mols = [] try: Loading @@ -267,25 +302,28 @@ class AtomicConvFeaturizer(ComplexFeaturizer): except MoleculeLoadException: logging.warning("This molecule cannot be loaded by Rdkit. Returning None") return None mols = [frag[1] for frag in fragments] system_mol = rdkit_util.merge_molecules(mols) system_coords = rdkit_util.get_xyz_from_mol(system_mol) if self.strip_hydrogens: fragments = [ rdkit_util.strip_hydrogens(frag[0], frag[1]) for frag in fragments ] system_coords, system_mol = rdkit_util.strip_hydrogens( system_coords, system_mol) if self.reduce_to_contacts: fragments = reduce_molecular_complex_to_contacts(fragments, self.cutoff) coords = [frag[0] for frag in fragments] mols = [frag[1] for frag in fragments] #system_mol = rdkit_util.merge_molecules(mols) system_mol = rdkit_util.merge_molecular_fragments(mols) #system_coords = rdkit_util.get_xyz_from_mol(system_mol) system_coords = rdkit_util.merge_molecules_xyz(coords) try: frag_inputs = [ self.featurize_mol(frag[0], frag[1], frag_num_atoms) for (frag, frag_num_atoms) in zip(fragments, self.frag_num_atoms) frag_outputs = [ self.featurize_mol(frag[0], frag[1], frag_max) for (frag, frag_max) in zip(fragments, frag_max_atoms) ] system_outputs = self.featurize_mol(system_coords, system_mol, self.complex_num_atoms) complex_max_atoms) except ValueError as e: logging.warning( "max_atoms was set too low. Some complexes too large and skipped") Loading @@ -294,6 +332,21 @@ class AtomicConvFeaturizer(ComplexFeaturizer): return frag_outputs, system_outputs def get_Z_matrix(self, mol, max_atoms): """Helper function to make the matrix of atomic-numbers Parameters ---------- mol: rdkit mol The molecule to featurize. max_atoms: int The max number of atoms allowed. Returns ------- Numpy array of shape `(max_atoms,)`. The first `len(mol.GetAtoms())` entries will contain the atomic numbers with the remaining entries zero padded. """ if len(mol.GetAtoms()) > max_atoms: raise ValueError("A molecule is larger than permitted by max_atoms. " "Increase max_atoms and try again.") Loading @@ -301,7 +354,18 @@ class AtomicConvFeaturizer(ComplexFeaturizer): np.array([atom.GetAtomicNum() for atom in mol.GetAtoms()]), max_atoms) def featurize_mol(self, coords, mol, max_num_atoms): logging.info("Featurizing molecule of size: %d", len(mol.GetAtoms())) """Helper function to featurize each molecule in complex. Parameters ---------- coords: Numpy array Shape `(N, 3)` for this molecule mol: rdkit mol Rdkit mol corresponding to `coords` max_num_atoms: int Max number of atoms for this molecules. """ logger.info("Featurizing molecule of size: %d", len(mol.GetAtoms())) neighbor_list = _compute_neighbor_list(coords, self.neighbor_cutoff, self.max_num_neighbors, None) z = self.get_Z_matrix(mol, max_num_atoms) Loading deepchem/feat/contact_fingerprints.py +0 −1 Original line number Diff line number Diff line Loading @@ -40,7 +40,6 @@ def featurize_contacts_ecfp(frag1, pairwise_distances = compute_pairwise_distances(frag1[0], frag2[0]) # contacts is of form (x_coords, y_coords), a tuple of 2 lists contacts = np.nonzero((pairwise_distances < cutoff)) #protein_atoms = set([int(c) for c in contacts[0].tolist()]) # contacts[0] is the x_coords, that is the frag1 atoms that have # nonzero contact. frag1_atoms = set([int(c) for c in contacts[0].tolist()]) Loading deepchem/utils/rdkit_util.py +311 −0 Original line number Diff line number Diff line Loading @@ -433,6 +433,317 @@ def merge_molecules(molecules): combined = rdmolops.CombineMols(combined, nextmol) return combined def merge_molecular_fragments(molecules): """Helper method to merge two molecular fragments. Parameters ---------- molecules: list List of `MolecularFragment` objects. Returns ------- merged: `MolecularFragment` """ if len(molecules) == 0: return None if len(molecules) == 1: return molecules[0] else: all_atoms = [] for mol_frag in molecules: all_atoms += mol_frag.GetAtoms() return MolecularFragment(all_atoms) def strip_hydrogens(coords, mol): """Strip the hydrogens from input molecule Parameters ---------- coords: Numpy ndarray Must be of shape (N, 3) and correspond to coordinates of mol. mol: Rdkit mol or `MolecularFragment` The molecule to strip Returns ------- A tuple of (coords, mol_frag) where coords is a Numpy array of coordinates with hydrogen coordinates. mol_frag is a `MolecularFragment`. """ mol_atoms = mol.GetAtoms() atomic_numbers = [atom.GetAtomicNum() for atom in mol_atoms] atom_indices_to_keep = [ind for (ind, atomic_number) in enumerate(atomic_numbers) if (atomic_number != 1)] return get_mol_subset(coords, mol, atom_indices_to_keep) class MolecularFragment(object): """A class that represents a fragment of a molecule. It's often convenient to represent a fragment of a molecule. For example, if two molecules form a molecular complex, it may be useful to create two fragments which represent the subsets of each molecule that's close to the other molecule (in the contact region). """ def __init__(self, atoms): """Initialize this object. Parameters ---------- atoms: list Each entry in this list should be an rdkit Atom object. """ self.atoms = [AtomShim(x) for x in atoms] #self.atoms = atoms def GetAtoms(self): """Returns the list of atoms Returns ------- list of atoms in this fragment. """ return self.atoms class AtomShim(object): """This is a shim object wrapping an atom. We use this class instead of raw RDKit atoms since manipulating a large number of rdkit Atoms seems to result in segfaults. Wrapping the basic information in an AtomShim seems to avoid issues. """ def __init__(self, atomic_num): """Initialize this object Parameters ---------- atomic_num: int Atomic number for this atom. """ self.atomic_num = atomic_num def GetAtomicNum(self): return self.atomic_num def get_mol_subset(coords, mol, atom_indices_to_keep): """Strip a subset of the atoms in this molecule Parameters ---------- coords: Numpy ndarray Must be of shape (N, 3) and correspond to coordinates of mol. mol: Rdkit mol The molecule to strip atom_indices_to_keep: list List of the indices of the atoms to keep. Each index is a unique number between `[0, N)`. Returns ------- A tuple of (coords, mol_frag) where coords is a Numpy array of coordinates with hydrogen coordinates. mol_frag is a `MolecularFragment`. """ indexes_to_keep = [] atoms_to_keep = [] atoms = list(mol.GetAtoms()) for index in atom_indices_to_keep: indexes_to_keep.append(index) #atomic_numbers.append(atom.GetAtomicNum()) atoms_to_keep.append(atoms[index]) #mol = MolecularFragment(atomic_numbers) mol_frag = MolecularFragment(atoms_to_keep) coords = coords[indexes_to_keep] return coords, mol_frag def reduce_molecular_complex_to_contacts(fragments, cutoff=4.5): """Reduce a molecular complex to only those atoms near a contact. Molecular complexes can get very large. This can make it unwieldy to compute functions on them. To improve memory usage, it can be very useful to trim out atoms that aren't close to contact regions. This function takes in a molecular complex and returns a new molecular complex representation that contains only contact atoms. The contact atoms are computed by calling `get_contact_atom_indices` under the hood. Parameters ---------- fragments: List As returned by `rdkit_util.load_complex`, a list of tuples of `(coords, mol)` where `coords` is a `(N_atoms, 3)` array and `mol` is the rdkit molecule object. cutoff: float The cutoff distance in angstroms. Returns ------- A list of length `len(molecular_complex)`. Each entry in this list is a tuple of `(coords, MolecularShim)`. The coords is stripped down to `(N_contact_atoms, 3)` where `N_contact_atoms` is the number of contact atoms for this complex. `MolecularShim` is used since it's tricky to make a RDKit sub-molecule. """ atoms_to_keep = get_contact_atom_indices(fragments, cutoff) reduced_complex = [] for frag, keep in zip(fragments, atoms_to_keep): contact_frag = get_mol_subset(frag[0], frag[1], keep) reduced_complex.append(contact_frag) return reduced_complex def get_contact_atom_indices(fragments, cutoff=4.5): """Compute that atoms close to contact region. Molecular complexes can get very large. This can make it unwieldy to compute functions on them. To improve memory usage, it can be very useful to trim out atoms that aren't close to contact regions. This function computes pairwise distances between all pairs of molecules in the molecular complex. If an atom is within cutoff distance of any atom on another molecule in the complex, it is regarded as a contact atom. Otherwise it is trimmed. Parameters ---------- fragments: List As returned by `rdkit_util.load_complex`, a list of tuples of `(coords, mol)` where `coords` is a `(N_atoms, 3)` array and `mol` is the rdkit molecule object. cutoff: float The cutoff distance in angstroms. Returns ------- A list of length `len(molecular_complex)`. Each entry in this list is a list of atom indices from that molecule which should be kept, in sorted order. """ # indices to atoms to keep keep_inds = [set([]) for _ in fragments] for (ind1, ind2) in itertools.combinations(range(len(fragments)), 2): frag1, frag2 = fragments[ind1], fragments[ind2] pairwise_distances = compute_pairwise_distances(frag1[0], frag2[0]) # contacts is of form (x_coords, y_coords), a tuple of 2 lists contacts = np.nonzero((pairwise_distances < cutoff)) # contacts[0] is the x_coords, that is the frag1 atoms that have # nonzero contact. frag1_atoms = set([int(c) for c in contacts[0].tolist()]) # contacts[1] is the y_coords, the frag2 atoms with nonzero contacts frag2_atoms = set([int(c) for c in contacts[1].tolist()]) keep_inds[ind1] = keep_inds[ind1].union(frag1_atoms) keep_inds[ind2] = keep_inds[ind2].union(frag2_atoms) keep_inds = [sorted(list(keep)) for keep in keep_inds] return keep_inds # Now extract atoms #atoms_to_keep = [] #for i, frag_keep_inds in enumerate(keep_inds): # frag = fragments[i] # mol = frag[1] # atoms = mol.GetAtoms() # frag_keep = [atoms[keep_ind] for keep_ind in frag_keep_inds] # atoms_to_keep.append(frag_keep) #return atoms_to_keep def compute_centroid(coordinates): """Compute the x,y,z centroid of provided coordinates coordinates: np.ndarray Shape (N, 3), where N is number atoms. """ centroid = np.mean(coordinates, axis=0) return (centroid) def subtract_centroid(xyz, centroid): """Subtracts centroid from each coordinate. Subtracts the centroid, a numpy array of dim 3, from all coordinates of all atoms in the molecule """ xyz -= np.transpose(centroid) return (xyz) def compute_ring_center(mol, ring_indices): """Computes 3D coordinates of a center of a given ring. Parameters: ----------- mol: rdkit.rdchem.Mol Molecule containing a ring ring_indices: array-like Indices of atoms forming a ring Returns: -------- ring_centroid: np.ndarray Position of a ring center """ conformer = mol.GetConformer() ring_xyz = np.zeros((len(ring_indices), 3)) for i, atom_idx in enumerate(ring_indices): atom_position = conformer.GetAtomPosition(atom_idx) ring_xyz[i] = np.array(atom_position) ring_centroid = compute_centroid(ring_xyz) return ring_centroid def compute_ring_normal(mol, ring_indices): """Computes normal to a plane determined by a given ring. Parameters: ----------- mol: rdkit.rdchem.Mol Molecule containing a ring ring_indices: array-like Indices of atoms forming a ring Returns: -------- normal: np.ndarray Normal vector """ conformer = mol.GetConformer() points = np.zeros((3, 3)) for i, atom_idx in enumerate(ring_indices[:3]): atom_position = conformer.GetAtomPosition(atom_idx) points[i] = np.array(atom_position) v1 = points[1] - points[0] v2 = points[2] - points[0] normal = np.cross(v1, v2) return normal def rotate_molecules(mol_coordinates_list): """Rotates provided molecular coordinates. Pseudocode: 1. Generate random rotation matrix. This matrix applies a random transformation to any 3-vector such that, were the random transformation repeatedly applied, it would randomly sample along the surface of a sphere with radius equal to the norm of the given 3-vector cf. generate_random_rotation_matrix() for details 2. Apply R to all atomic coordinates. 3. Return rotated molecule Parameters ---------- mol_coordinates_list: list Elements of list must be (N_atoms, 3) shaped arrays """ from rdkit.Chem import rdmolops if len(molecules) == 0: return None elif len(molecules) == 1: return molecules[0] else: combined = molecules[0] for nextmol in molecules[1:]: combined = rdmolops.CombineMols(combined, nextmol) return combined def is_hydrogen_bond(protein_xyz, protein, Loading examples/pdbbind/featurize_atomic_conv.py +1 −3 Original line number Diff line number Diff line Loading @@ -6,13 +6,11 @@ from deepchem.molnet.load_function.pdbbind_datasets import get_pdbbind_molecular complex_files = get_pdbbind_molecular_complex_files(subset="core", version="v2015", interactions="protein-ligand", load_binding_pocket=False) core_subset = complex_files[:2] complex_num_atoms = 24070 # in total max_num_neighbors = 4 # Cutoff in angstroms neighbor_cutoff = 4 featurizer = dc.feat.AtomicConvFeaturizer( frag_num_atoms=[70 ,24000], complex_num_atoms=complex_num_atoms, frag_max_atoms=70, max_num_neighbors=max_num_neighbors, neighbor_cutoff=neighbor_cutoff) features, failures = featurizer.featurize_complexes( Loading examples/pdbbind/featurize_atomic_coordinates.py +0 −6 Original line number Diff line number Diff line Loading @@ -9,10 +9,4 @@ core_subset = ligand_files[:2] featurizer = dc.feat.AtomicCoordinates() features = featurizer.featurize(core_subset) print("features.shape") print(features.shape) print("features[0].shape") print(features[0].shape) print("features[1].shape") print(features[1].shape) assert features.shape == (2,) Loading
deepchem/feat/atomic_coordinates.py +90 −26 Original line number Diff line number Diff line Loading @@ -8,6 +8,7 @@ from deepchem.feat import MolecularFeaturizer from deepchem.feat import ComplexFeaturizer from deepchem.utils import rdkit_util, pad_array from deepchem.utils.rdkit_util import MoleculeLoadException from deepchem.utils.rdkit_util import reduce_molecular_complex_to_contacts logger = logging.getLogger(__name__) Loading Loading @@ -58,7 +59,19 @@ class AtomicCoordinates(MolecularFeaturizer): def _compute_neighbor_list(coords, neighbor_cutoff, max_num_neighbors, periodic_box_size): """Computes a neighbor list from atom coordinates.""" """Computes a neighbor list from atom coordinates. Parameters ---------- coords: Numpy array Of shape (N, 3) with all the atoms in this systems neighbor_cutoff: float The neighbor cutoff in angstroms max_num_neighbors: int The maximum number of neighbors per atom periodic_box_size: tuple With (x, y, z) box sizes per dimension """ N = coords.shape[0] import mdtraj traj = mdtraj.Trajectory(coords.reshape((1, N, 3)), None) Loading Loading @@ -89,6 +102,11 @@ def _compute_neighbor_list(coords, neighbor_cutoff, max_num_neighbors, def get_coords(mol): """ Gets coordinates in Angstrom for RDKit mol. Parameters ---------- mol: rdkit mol Molecule to get coordinates for """ N = mol.GetNumAtoms() coords = np.zeros((N, 3)) Loading Loading @@ -182,7 +200,7 @@ class NeighborListComplexAtomicCoordinates(ComplexFeaturizer): self.dtype = object self.coordinates_featurizer = AtomicCoordinates() def _featurize(self, mol_pdb_file, protein_pdb_file): def _featurize(self, molecular_complex): """ Compute neighbor list for complex. Loading @@ -191,9 +209,13 @@ class NeighborListComplexAtomicCoordinates(ComplexFeaturizer): molecular_complex: Object Some representation of a molecular complex. """ mol_coords, ob_mol = rdkit_util.load_molecule(mol_pdb_file) protein_coords, protein_mol = rdkit_util.load_molecule(protein_pdb_file) system_coords = rdkit_util.merge_molecules_xyz([mol_coords, protein_coords]) fragments = rdkit_util.load_complex(molecular_complex, add_hydrogens=False) #mol_coords, ob_mol = rdkit_util.load_molecule(mol_pdb_file) #protein_coords, protein_mol = rdkit_util.load_molecule(protein_pdb_file) coords = [frag[0] for frag in fragments] mols = [frag[1] for frag in fragments] #system_coords = rdkit_util.merge_molecules_xyz(mol_coords, protein_coords) system_coords = rdkit_util.merge_molecules_xyz(coords) system_neighbor_list = _compute_neighbor_list( system_coords, self.neighbor_cutoff, self.max_num_neighbors, None) Loading Loading @@ -223,31 +245,37 @@ class AtomicConvFeaturizer(ComplexFeaturizer): """ def __init__(self, frag_num_atoms, complex_num_atoms, frag_max_atoms, max_num_neighbors, neighbor_cutoff, reduce_to_contacts=True, cutoff=4.5, strip_hydrogens=True): """Initialize an AtomicConvFeaturizer object. Parameters ---------- frag_num_atoms: list[int] List of the number of atoms in each fragment. frag_max_atoms: int or list[int] List of the max number of atoms in each fragment. If int, assumes you're setting the same bound for all fragments. max_num_neighbors: int The maximum number of neighbors allowed neighbor_cutoff: float The distance in angstroms after which neighbors are cutoff. reduce_to_contacts: bool, optional If True, reduce the atoms in the complex to those near a contact region. cutoff: float The cutoff distance in angstroms. Only used if `reduce_to_contacts` is `True`. strip_hydrogens: bool, optional If true, remove hydrogens before featurizing. If True, remove hydrogens before featurizing. """ # TODO(rbharath): extend to more fragments if len(frag_num_atoms) != 2: raise ValueError("Currently only supports two fragments") self.frag_num_atoms = frag_num_atoms self.complex_num_atoms = sum(frag_num_atoms) self.frag_max_atoms = frag_max_atoms self.max_num_neighbors = max_num_neighbors self.neighbor_cutoff = neighbor_cutoff self.reduce_to_contacts = reduce_to_contacts self.cutoff = cutoff self.strip_hydrogens = strip_hydrogens def _featurize(self, molecular_complex): Loading @@ -258,6 +286,13 @@ class AtomicConvFeaturizer(ComplexFeaturizer): molecular_complex: Object Some representation of a molecular complex. """ # If our upper bound is an int, expand it to a list for each # fragment if isinstance(self.frag_max_atoms, int): frag_max_atoms = [self.frag_max_atoms] * len(molecular_complex) else: frag_max_atoms = self.frag_max_atoms complex_max_atoms = sum(frag_max_atoms) frag_coords = [] frag_mols = [] try: Loading @@ -267,25 +302,28 @@ class AtomicConvFeaturizer(ComplexFeaturizer): except MoleculeLoadException: logging.warning("This molecule cannot be loaded by Rdkit. Returning None") return None mols = [frag[1] for frag in fragments] system_mol = rdkit_util.merge_molecules(mols) system_coords = rdkit_util.get_xyz_from_mol(system_mol) if self.strip_hydrogens: fragments = [ rdkit_util.strip_hydrogens(frag[0], frag[1]) for frag in fragments ] system_coords, system_mol = rdkit_util.strip_hydrogens( system_coords, system_mol) if self.reduce_to_contacts: fragments = reduce_molecular_complex_to_contacts(fragments, self.cutoff) coords = [frag[0] for frag in fragments] mols = [frag[1] for frag in fragments] #system_mol = rdkit_util.merge_molecules(mols) system_mol = rdkit_util.merge_molecular_fragments(mols) #system_coords = rdkit_util.get_xyz_from_mol(system_mol) system_coords = rdkit_util.merge_molecules_xyz(coords) try: frag_inputs = [ self.featurize_mol(frag[0], frag[1], frag_num_atoms) for (frag, frag_num_atoms) in zip(fragments, self.frag_num_atoms) frag_outputs = [ self.featurize_mol(frag[0], frag[1], frag_max) for (frag, frag_max) in zip(fragments, frag_max_atoms) ] system_outputs = self.featurize_mol(system_coords, system_mol, self.complex_num_atoms) complex_max_atoms) except ValueError as e: logging.warning( "max_atoms was set too low. Some complexes too large and skipped") Loading @@ -294,6 +332,21 @@ class AtomicConvFeaturizer(ComplexFeaturizer): return frag_outputs, system_outputs def get_Z_matrix(self, mol, max_atoms): """Helper function to make the matrix of atomic-numbers Parameters ---------- mol: rdkit mol The molecule to featurize. max_atoms: int The max number of atoms allowed. Returns ------- Numpy array of shape `(max_atoms,)`. The first `len(mol.GetAtoms())` entries will contain the atomic numbers with the remaining entries zero padded. """ if len(mol.GetAtoms()) > max_atoms: raise ValueError("A molecule is larger than permitted by max_atoms. " "Increase max_atoms and try again.") Loading @@ -301,7 +354,18 @@ class AtomicConvFeaturizer(ComplexFeaturizer): np.array([atom.GetAtomicNum() for atom in mol.GetAtoms()]), max_atoms) def featurize_mol(self, coords, mol, max_num_atoms): logging.info("Featurizing molecule of size: %d", len(mol.GetAtoms())) """Helper function to featurize each molecule in complex. Parameters ---------- coords: Numpy array Shape `(N, 3)` for this molecule mol: rdkit mol Rdkit mol corresponding to `coords` max_num_atoms: int Max number of atoms for this molecules. """ logger.info("Featurizing molecule of size: %d", len(mol.GetAtoms())) neighbor_list = _compute_neighbor_list(coords, self.neighbor_cutoff, self.max_num_neighbors, None) z = self.get_Z_matrix(mol, max_num_atoms) Loading
deepchem/feat/contact_fingerprints.py +0 −1 Original line number Diff line number Diff line Loading @@ -40,7 +40,6 @@ def featurize_contacts_ecfp(frag1, pairwise_distances = compute_pairwise_distances(frag1[0], frag2[0]) # contacts is of form (x_coords, y_coords), a tuple of 2 lists contacts = np.nonzero((pairwise_distances < cutoff)) #protein_atoms = set([int(c) for c in contacts[0].tolist()]) # contacts[0] is the x_coords, that is the frag1 atoms that have # nonzero contact. frag1_atoms = set([int(c) for c in contacts[0].tolist()]) Loading
deepchem/utils/rdkit_util.py +311 −0 Original line number Diff line number Diff line Loading @@ -433,6 +433,317 @@ def merge_molecules(molecules): combined = rdmolops.CombineMols(combined, nextmol) return combined def merge_molecular_fragments(molecules): """Helper method to merge two molecular fragments. Parameters ---------- molecules: list List of `MolecularFragment` objects. Returns ------- merged: `MolecularFragment` """ if len(molecules) == 0: return None if len(molecules) == 1: return molecules[0] else: all_atoms = [] for mol_frag in molecules: all_atoms += mol_frag.GetAtoms() return MolecularFragment(all_atoms) def strip_hydrogens(coords, mol): """Strip the hydrogens from input molecule Parameters ---------- coords: Numpy ndarray Must be of shape (N, 3) and correspond to coordinates of mol. mol: Rdkit mol or `MolecularFragment` The molecule to strip Returns ------- A tuple of (coords, mol_frag) where coords is a Numpy array of coordinates with hydrogen coordinates. mol_frag is a `MolecularFragment`. """ mol_atoms = mol.GetAtoms() atomic_numbers = [atom.GetAtomicNum() for atom in mol_atoms] atom_indices_to_keep = [ind for (ind, atomic_number) in enumerate(atomic_numbers) if (atomic_number != 1)] return get_mol_subset(coords, mol, atom_indices_to_keep) class MolecularFragment(object): """A class that represents a fragment of a molecule. It's often convenient to represent a fragment of a molecule. For example, if two molecules form a molecular complex, it may be useful to create two fragments which represent the subsets of each molecule that's close to the other molecule (in the contact region). """ def __init__(self, atoms): """Initialize this object. Parameters ---------- atoms: list Each entry in this list should be an rdkit Atom object. """ self.atoms = [AtomShim(x) for x in atoms] #self.atoms = atoms def GetAtoms(self): """Returns the list of atoms Returns ------- list of atoms in this fragment. """ return self.atoms class AtomShim(object): """This is a shim object wrapping an atom. We use this class instead of raw RDKit atoms since manipulating a large number of rdkit Atoms seems to result in segfaults. Wrapping the basic information in an AtomShim seems to avoid issues. """ def __init__(self, atomic_num): """Initialize this object Parameters ---------- atomic_num: int Atomic number for this atom. """ self.atomic_num = atomic_num def GetAtomicNum(self): return self.atomic_num def get_mol_subset(coords, mol, atom_indices_to_keep): """Strip a subset of the atoms in this molecule Parameters ---------- coords: Numpy ndarray Must be of shape (N, 3) and correspond to coordinates of mol. mol: Rdkit mol The molecule to strip atom_indices_to_keep: list List of the indices of the atoms to keep. Each index is a unique number between `[0, N)`. Returns ------- A tuple of (coords, mol_frag) where coords is a Numpy array of coordinates with hydrogen coordinates. mol_frag is a `MolecularFragment`. """ indexes_to_keep = [] atoms_to_keep = [] atoms = list(mol.GetAtoms()) for index in atom_indices_to_keep: indexes_to_keep.append(index) #atomic_numbers.append(atom.GetAtomicNum()) atoms_to_keep.append(atoms[index]) #mol = MolecularFragment(atomic_numbers) mol_frag = MolecularFragment(atoms_to_keep) coords = coords[indexes_to_keep] return coords, mol_frag def reduce_molecular_complex_to_contacts(fragments, cutoff=4.5): """Reduce a molecular complex to only those atoms near a contact. Molecular complexes can get very large. This can make it unwieldy to compute functions on them. To improve memory usage, it can be very useful to trim out atoms that aren't close to contact regions. This function takes in a molecular complex and returns a new molecular complex representation that contains only contact atoms. The contact atoms are computed by calling `get_contact_atom_indices` under the hood. Parameters ---------- fragments: List As returned by `rdkit_util.load_complex`, a list of tuples of `(coords, mol)` where `coords` is a `(N_atoms, 3)` array and `mol` is the rdkit molecule object. cutoff: float The cutoff distance in angstroms. Returns ------- A list of length `len(molecular_complex)`. Each entry in this list is a tuple of `(coords, MolecularShim)`. The coords is stripped down to `(N_contact_atoms, 3)` where `N_contact_atoms` is the number of contact atoms for this complex. `MolecularShim` is used since it's tricky to make a RDKit sub-molecule. """ atoms_to_keep = get_contact_atom_indices(fragments, cutoff) reduced_complex = [] for frag, keep in zip(fragments, atoms_to_keep): contact_frag = get_mol_subset(frag[0], frag[1], keep) reduced_complex.append(contact_frag) return reduced_complex def get_contact_atom_indices(fragments, cutoff=4.5): """Compute that atoms close to contact region. Molecular complexes can get very large. This can make it unwieldy to compute functions on them. To improve memory usage, it can be very useful to trim out atoms that aren't close to contact regions. This function computes pairwise distances between all pairs of molecules in the molecular complex. If an atom is within cutoff distance of any atom on another molecule in the complex, it is regarded as a contact atom. Otherwise it is trimmed. Parameters ---------- fragments: List As returned by `rdkit_util.load_complex`, a list of tuples of `(coords, mol)` where `coords` is a `(N_atoms, 3)` array and `mol` is the rdkit molecule object. cutoff: float The cutoff distance in angstroms. Returns ------- A list of length `len(molecular_complex)`. Each entry in this list is a list of atom indices from that molecule which should be kept, in sorted order. """ # indices to atoms to keep keep_inds = [set([]) for _ in fragments] for (ind1, ind2) in itertools.combinations(range(len(fragments)), 2): frag1, frag2 = fragments[ind1], fragments[ind2] pairwise_distances = compute_pairwise_distances(frag1[0], frag2[0]) # contacts is of form (x_coords, y_coords), a tuple of 2 lists contacts = np.nonzero((pairwise_distances < cutoff)) # contacts[0] is the x_coords, that is the frag1 atoms that have # nonzero contact. frag1_atoms = set([int(c) for c in contacts[0].tolist()]) # contacts[1] is the y_coords, the frag2 atoms with nonzero contacts frag2_atoms = set([int(c) for c in contacts[1].tolist()]) keep_inds[ind1] = keep_inds[ind1].union(frag1_atoms) keep_inds[ind2] = keep_inds[ind2].union(frag2_atoms) keep_inds = [sorted(list(keep)) for keep in keep_inds] return keep_inds # Now extract atoms #atoms_to_keep = [] #for i, frag_keep_inds in enumerate(keep_inds): # frag = fragments[i] # mol = frag[1] # atoms = mol.GetAtoms() # frag_keep = [atoms[keep_ind] for keep_ind in frag_keep_inds] # atoms_to_keep.append(frag_keep) #return atoms_to_keep def compute_centroid(coordinates): """Compute the x,y,z centroid of provided coordinates coordinates: np.ndarray Shape (N, 3), where N is number atoms. """ centroid = np.mean(coordinates, axis=0) return (centroid) def subtract_centroid(xyz, centroid): """Subtracts centroid from each coordinate. Subtracts the centroid, a numpy array of dim 3, from all coordinates of all atoms in the molecule """ xyz -= np.transpose(centroid) return (xyz) def compute_ring_center(mol, ring_indices): """Computes 3D coordinates of a center of a given ring. Parameters: ----------- mol: rdkit.rdchem.Mol Molecule containing a ring ring_indices: array-like Indices of atoms forming a ring Returns: -------- ring_centroid: np.ndarray Position of a ring center """ conformer = mol.GetConformer() ring_xyz = np.zeros((len(ring_indices), 3)) for i, atom_idx in enumerate(ring_indices): atom_position = conformer.GetAtomPosition(atom_idx) ring_xyz[i] = np.array(atom_position) ring_centroid = compute_centroid(ring_xyz) return ring_centroid def compute_ring_normal(mol, ring_indices): """Computes normal to a plane determined by a given ring. Parameters: ----------- mol: rdkit.rdchem.Mol Molecule containing a ring ring_indices: array-like Indices of atoms forming a ring Returns: -------- normal: np.ndarray Normal vector """ conformer = mol.GetConformer() points = np.zeros((3, 3)) for i, atom_idx in enumerate(ring_indices[:3]): atom_position = conformer.GetAtomPosition(atom_idx) points[i] = np.array(atom_position) v1 = points[1] - points[0] v2 = points[2] - points[0] normal = np.cross(v1, v2) return normal def rotate_molecules(mol_coordinates_list): """Rotates provided molecular coordinates. Pseudocode: 1. Generate random rotation matrix. This matrix applies a random transformation to any 3-vector such that, were the random transformation repeatedly applied, it would randomly sample along the surface of a sphere with radius equal to the norm of the given 3-vector cf. generate_random_rotation_matrix() for details 2. Apply R to all atomic coordinates. 3. Return rotated molecule Parameters ---------- mol_coordinates_list: list Elements of list must be (N_atoms, 3) shaped arrays """ from rdkit.Chem import rdmolops if len(molecules) == 0: return None elif len(molecules) == 1: return molecules[0] else: combined = molecules[0] for nextmol in molecules[1:]: combined = rdmolops.CombineMols(combined, nextmol) return combined def is_hydrogen_bond(protein_xyz, protein, Loading
examples/pdbbind/featurize_atomic_conv.py +1 −3 Original line number Diff line number Diff line Loading @@ -6,13 +6,11 @@ from deepchem.molnet.load_function.pdbbind_datasets import get_pdbbind_molecular complex_files = get_pdbbind_molecular_complex_files(subset="core", version="v2015", interactions="protein-ligand", load_binding_pocket=False) core_subset = complex_files[:2] complex_num_atoms = 24070 # in total max_num_neighbors = 4 # Cutoff in angstroms neighbor_cutoff = 4 featurizer = dc.feat.AtomicConvFeaturizer( frag_num_atoms=[70 ,24000], complex_num_atoms=complex_num_atoms, frag_max_atoms=70, max_num_neighbors=max_num_neighbors, neighbor_cutoff=neighbor_cutoff) features, failures = featurizer.featurize_complexes( Loading
examples/pdbbind/featurize_atomic_coordinates.py +0 −6 Original line number Diff line number Diff line Loading @@ -9,10 +9,4 @@ core_subset = ligand_files[:2] featurizer = dc.feat.AtomicCoordinates() features = featurizer.featurize(core_subset) print("features.shape") print(features.shape) print("features[0].shape") print(features[0].shape) print("features[1].shape") print(features[1].shape) assert features.shape == (2,)
