"""
Create polymer descriptor and molecule features using rdkit.
Copyright Schrodinger, LLC. All rights reserved.
"""
import numpy
import pandas
from rdkit import Chem
from rdkit.Chem import rdMolDescriptors as rdMol
from rdkit.Chem import rdFingerprintGenerator
from rdkit.Chem.Lipinski import RotatableBondSmarts
from rdkit.Chem.rdMolDescriptors import CalcMolFormula
from rdkit.ML.Descriptors import MoleculeDescriptors
from matminer.featurizers.composition import ElementProperty
from matminer.featurizers.conversions import StrToComposition
SMILES = 'smiles'
MOL = 'mol'
# Descriptors type
POLYMER_DESCRIPTOR = 'polymer_descriptor'
RDKIT_DESCRIPTOR = 'rdkit_descriptor'
MORGAN = 'Morgan'
MATMINER = 'Matminer'
[docs]class FeaturizeError(RuntimeError):
pass
[docs]class Featurize:
"""
Create features from the structures
"""
[docs] def __init__(self,
dataframe,
smile_col=SMILES,
feature_types=(RDKIT_DESCRIPTOR, MORGAN),
return_float64=False,
logger=None):
"""
Initialize FeaturizeStructure class
:type dataframe: pandas.DataFrame
:param dataframe: pandas dataframe containing structures
:type smile_col: str
:param smile_col: Name of smiles column.
:type feature_types: list
:param feature_types: list of features to be calculated for given
structures
:type return_float64: bool
:param return_float64:
:type logger: `logging.Logger`
:param logger: The logger to log with
"""
self.dataframe = dataframe
self.feature_types = feature_types
self.smile_col = smile_col
self.logger = logger
self.feature_cols = []
self.return_float64 = return_float64
[docs] def getFeatures(self):
"""
Create features from SMILES pattern of dataset
:rtype: pandas.DataFrame
:return: Dataframe containing structure, SMILES and full feature vector
"""
# Generate RDKIT mol
self.createMols()
# Generate feature vector
self.feature_cols = []
if RDKIT_DESCRIPTOR in self.feature_types:
self.createRdkitDescriptors()
if MATMINER in self.feature_types:
self.createMatMinerDescriptors()
if MORGAN in self.feature_types:
self.createMorganDescriptors()
if POLYMER_DESCRIPTOR in self.feature_types:
self.createPolymerDescriptors()
self.clearDataFrame()
return self.dataframe, self.feature_cols
[docs] def clearDataFrame(self):
"""
Clear data
"""
self.dataframe[self.feature_cols] = self.dataframe[
self.feature_cols].fillna(0)
self.dataframe[self.feature_cols] = self.dataframe[
self.feature_cols].replace([numpy.inf, -numpy.inf], 0)
if self.return_float64:
self.dataframe[self.dataframe.select_dtypes(
numpy.float64).columns] = self.dataframe.select_dtypes(
numpy.float64).astype(numpy.float32)
[docs] def createMols(self):
"""
Create rdkit molecules from SMILES
"""
nmol = len(self.dataframe.index)
if self.logger:
self.logger.info(f"Creating descriptor for {nmol} molecules")
# Importing pandas tool cause issue. See SHARED-8652
from rdkit.Chem import PandasTools
PandasTools.AddMoleculeColumnToFrame(self.dataframe,
smilesCol=self.smile_col,
molCol=MOL)
PandasTools.UninstallPandasTools()
self.dataframe.dropna(subset=[MOL], axis='rows', inplace=True)
self.dataframe.reset_index(drop=True, inplace=True)
# Check if all molecules are properly created
if len(self.dataframe.index) != nmol and self.logger:
self.logger.info("Could not generate rdkit molecules from SMILES "
"for some inputs. Dropping "
f"{nmol - len(self.dataframe.index)} molecules")
if self.dataframe.empty:
raise FeaturizeError("No molecule is generated. Exiting now.")
[docs] def createRdkitDescriptors(self):
"""
Create RDKit descriptors for the molecules
"""
descriptors_names = [x[0] for x in Chem.Descriptors._descList]
calculator = MoleculeDescriptors.MolecularDescriptorCalculator(
descriptors_names)
descript_val = numpy.stack(
[calculator.CalcDescriptors(mol) for mol in self.dataframe[MOL]])
# set Nan and Inf to zero
descript_val[numpy.isnan(descript_val)] = 0
descript_val[numpy.isinf(descript_val)] = 0
rdkit_descr_name = [
'rdkit_descr_' + descript for descript in descriptors_names
]
descr_df = pandas.DataFrame(descript_val, columns=rdkit_descr_name)
self.dataframe = pandas.concat((self.dataframe, descr_df), axis=1)
self.feature_cols.extend(rdkit_descr_name)
# Contiguous bond group list
self.dataframe['max_contiguous'] = self.dataframe[MOL].apply(
self.findBondGroups)
self.feature_cols.append('max_contiguous')
[docs] def findBondGroups(self, mol):
"""
Find largest size of contiguous rotatable bonds present in the system
:param mol: rdkit molecule
:type mol: `rdkit.Chem.mol`
:rtype: int
:return: length of largest contiguous rotatable bonds present
"""
# get atom pair with rotatable bond
rot_atom_pairs = mol.GetSubstructMatches(RotatableBondSmarts)
# Extract bond id from atom pairs
rot_bond_set = set(
[mol.GetBondBetweenAtoms(*ap).GetIdx() for ap in rot_atom_pairs])
# Note: Partition rot_bond_set based on connectivity. All connected
# rotational bonds are placed in one set. The scheme pops bonds from
# the pre-existing set of rotatable bond. Checks if, atoms in bonds are
# also connected to any other atom by a rotatable bond which exists
# in pre-existing dataset but not in
# result set yet.
# TODO: Use networkx.connected_component to simplify this loop.
rot_bond_groups = []
while rot_bond_set:
idx = rot_bond_set.pop()
connected_bond_set = {idx}
stack = [idx]
while stack:
idx = stack.pop()
bond = mol.GetBondWithIdx(idx)
for atom in (bond.GetBeginAtom(), bond.GetEndAtom()):
connected_bond = [
b.GetIdx()
for b in atom.GetBonds()
if b.GetIdx() in rot_bond_set and
b.GetIdx() not in connected_bond_set
]
connected_bond_set.update(connected_bond)
stack.extend(connected_bond)
# Remove all the bonds from pre-existing dataset as these are
# already found as connected.
rot_bond_set.difference_update(connected_bond_set)
rot_bond_groups.append(len(connected_bond_set))
return max(rot_bond_groups) if rot_bond_groups else 0
[docs] def createMatMinerDescriptors(self):
"""
Create Mat-Miner descriptor from the rdkit molecule
"""
# Clean up molecular formula
list_formula = [CalcMolFormula(mol) for mol in self.dataframe[MOL]]
list_formula = [form.replace('+2', '') for form in list_formula]
list_formula = [form.replace('+', '') for form in list_formula]
df_mat = pandas.DataFrame()
df_mat['formula'] = list_formula
featurizer = StrToComposition()
featurizer.set_n_jobs(1)
df_mat = featurizer.featurize_dataframe(df_mat,
"formula",
pbar=False,
ignore_errors=True)
ep_feat = ElementProperty.from_preset(preset_name="magpie")
ep_feat.set_n_jobs(1)
df_mat = ep_feat.featurize_dataframe(df_mat,
col_id="composition",
pbar=False,
ignore_errors=True)
magpie_features = [col for col in df_mat.columns if 'Magpie' in col]
self.dataframe[magpie_features] = df_mat[magpie_features]
self.feature_cols.extend(magpie_features)
[docs] def createMorganDescriptors(self):
"""
Create Morgan descriptor from the rdkit molecule
"""
morgan_fp_gen = rdFingerprintGenerator.GetMorganGenerator(
includeChirality=False,
radius=2,
fpSize=1000,
countSimulation=False)
self.dataframe['morgan_fp'] = self.dataframe[MOL].apply(
self.getMorganCount, morgan_fp_generator=morgan_fp_gen)
x_df_morgan = pandas.DataFrame(list(self.dataframe['morgan_fp'].values))
x_df_morgan.fillna(0, inplace=True)
self.dataframe = pandas.concat([self.dataframe, x_df_morgan], axis=1)
self.feature_cols.extend(x_df_morgan.columns)
[docs] def getMorganCount(self, mol, morgan_fp_generator):
"""
Generate morgan fingerprint for one molecule
:param mol: rdkit molecule
:type mol: `rdkit.Chem.mol`
:type morgan_fp_generator: class object
:param morgan_fp_generator: Morgan fingerprint generator object
:return: Dictionary with key as fingerprint name and value as
fingerprint
:rtype: Dict
"""
finger_print = numpy.array(morgan_fp_generator.GetFingerprint(mol))
fp_dict = {
'Morgan_' + str(i): finger_print[i]
for i in range(len(finger_print))
}
return fp_dict
[docs] def createPolymerDescriptors(self):
"""
Create polymer descriptors in the pandas dataframe.
"""
def getPolymerDescriptors(smiles):
"""
Generate polymer descriptor based on SMILES
:rtype: dict
:return: dictionary of feature name and value for input smiles
"""
descript = PolymerDescriptor(smiles)
try:
return descript.getFeatures()
except PolymerDescriptorError as err:
raise FeaturizeError(str(err))
polymer_descr = sorted(PolymerDescriptor.DESCRIPTORS)
dataframe = self.dataframe[self.smile_col].apply(getPolymerDescriptors)
x_df_polymer = pandas.DataFrame(list(dataframe.values))
self.dataframe = pandas.concat([self.dataframe, x_df_polymer], axis=1)
self.dataframe.dropna(subset=polymer_descr, axis='rows', inplace=True)
self.dataframe.reset_index(drop=True, inplace=True)
self.feature_cols.extend(polymer_descr)
[docs]class PolymerDescriptorError(RuntimeError):
pass
[docs]class PolymerDescriptor:
"""
Class to create the polymer descriptor
"""
FLOAT_BASE = 'r_matsci'
DESCRIPTORS = {
f'{FLOAT_BASE}_Ring_Atoms_Fraction': 'ringAtomsFraction',
f'{FLOAT_BASE}_Double_Ring_Atoms_Fraction': 'doubleRingAtomsFraction',
f'{FLOAT_BASE}_Triple_Ring_Atoms_Fraction': 'tripleRingAtomsFraction',
f'{FLOAT_BASE}_Backbone_Atoms_Fraction': 'backboneAtomsFraction',
f'{FLOAT_BASE}_Rotatable_Bonds_Fraction': 'rotatableBondsFraction',
f'{FLOAT_BASE}_Sp3_Atoms_Fraction': 'sp3AtomsFraction',
}
[docs] def __init__(self, smiles):
"""
Initialize polymer descriptor class
:param `rdkit.Chem.rdchem.Mol` mol: rdkit mol object of polymer
"""
self.monomer_smiles = smiles
self.dimer = None
self.n_dimer_atoms = None
self.dimer_rings = None
self.features = {}
[docs] def getFeatures(self):
"""
Get descriptor for given polymer
:return dict: Dictionary with key as name of the descriptor and value
as descriptor
"""
self.run()
return self.features
[docs] def run(self):
"""
Generate polymer descriptors for the structure
"""
try:
self.dimer = create_oligomer(self.monomer_smiles, 2)
except PolymerDescriptorError as err:
raise PolymerDescriptorError(str(err))
self.n_dimer_atoms = self.dimer.GetNumAtoms() - 2 # Exclude heavy atoms
ring_info = self.dimer.GetRingInfo()
self.dimer_rings = ring_info.AtomRings()
self.getDoubleAndTripleRingAtoms()
for prop_name, func_name, in self.DESCRIPTORS.items():
self.features[prop_name] = getattr(self, func_name)()
[docs] def getDoubleAndTripleRingAtoms(self):
"""
Calculate the number of atoms in double and triple fused rings
"""
def merge_sets(sets):
"""
Merge sets that have common elements
"""
for idx_1, set_1 in enumerate(sets):
for set_2 in sets[idx_1 + 1:]:
if not set_1.isdisjoint(set_2):
set_1.update(set_2)
sets.remove(set_2)
merge_sets(sets)
return
# Get neighbors and merge them together
ring_sets = [set(x) for x in self.dimer_rings]
neighbor_groups = []
for idx_1, ring_1 in enumerate(ring_sets):
neighbors = {idx_1}
start = idx_1 + 1
for idx_2, ring_2 in enumerate(ring_sets[start:], start):
if not ring_1.isdisjoint(ring_2):
neighbors.add(idx_2)
neighbor_groups.append(neighbors)
merge_sets(neighbor_groups)
# Get atoms for 2- and 3-ring neighbors
double_rings_atoms = []
triple_rings_atoms = []
for neighbors in neighbor_groups:
size = len(neighbors)
if size in (2, 3):
atoms = []
for ring_idx in neighbors:
atoms.extend(ring_sets[ring_idx])
if size == 2:
double_rings_atoms += atoms
elif size == 3:
triple_rings_atoms += atoms
self.double_rings_atoms = set(double_rings_atoms)
self.triple_rings_atoms = set(triple_rings_atoms)
[docs] def ringAtomsFraction(self):
"""
Get the fraction of atoms in rings in the polymer
:rtype: float
:return: The fraction of atoms in rings
"""
ring_atoms = {atom for ring in self.dimer_rings for atom in ring}
return len(ring_atoms) / self.n_dimer_atoms
[docs] def doubleRingAtomsFraction(self):
"""
Get the fraction of atoms in double fused ring systems in the polymer
:rtype: float
:return: The fraction of atoms in double fused ring systems
"""
return len(self.double_rings_atoms) / self.n_dimer_atoms
[docs] def tripleRingAtomsFraction(self):
"""
Get the fraction of atoms in triple fused ring systems in the polymer
:rtype: float
:return: The fraction of atoms in triple fused ring systems
"""
return len(self.triple_rings_atoms) / self.n_dimer_atoms
[docs] def backboneAtomsFraction(self):
"""
Get the fraction of backbone atoms in the polymer
:rtype: float
:return: The fraction of backbone atoms
"""
pattern = Chem.MolFromSmarts('[At]')
at_indexes = self.dimer.GetSubstructMatches(pattern)
head = at_indexes[0][0]
tail = at_indexes[-1][0]
original_backbone = Chem.rdmolops.GetShortestPath(
self.dimer, head, tail)
new_backbone = set(original_backbone)
for ring_atoms in self.dimer_rings:
if not set(ring_atoms).isdisjoint(original_backbone):
new_backbone.update(ring_atoms)
backbone_norm = (len(new_backbone) - 2) / self.n_dimer_atoms
return backbone_norm
[docs] def rotatableBondsFraction(self):
"""
Get the fraction of rotatable bonds in the polymer
:rtype: float
:return: The fraction of rotatable bonds
"""
monomer = Chem.MolFromSmiles(self.monomer_smiles)
rep_unit_len = self.dimer.GetNumAtoms() - monomer.GetNumAtoms()
dimer_rot_bonds = rdMol.CalcNumRotatableBonds(self.dimer)
monomer_rot_bonds = rdMol.CalcNumRotatableBonds(monomer)
rot_bonds_frac = (dimer_rot_bonds - monomer_rot_bonds) / rep_unit_len
return rot_bonds_frac
[docs] def sp3AtomsFraction(self):
"""
Get the fraction of sp3 atoms in the polymer
:rtype: float
:return: The fraction of sp3 atoms bonds
"""
return rdMol.CalcFractionCSP3(self.dimer)
[docs]def create_oligomer(smiles, monomers):
"""
Create an oligomer given the monomer SMILES and the number of monomer
repetitions. The head and tail of the monomer should be denoted by the atom
[At].
:param str smiles: The SMILES string of the monomer
:param int monomers: The number of monomers to repeat
:rtype: `Chem.rdchem.Mol`
:return: The oligomer as `Chem.rdchem.Mol`
"""
smiles_chain = '.'.join([smiles] * monomers)
mol = Chem.MolFromSmiles(smiles_chain)
if mol is None:
raise PolymerDescriptorError(f"Error in parsing SMILES: {smiles}. "
"Please check the input structure.")
pattern = Chem.MolFromSmarts('[At]')
at_indexes = mol.GetSubstructMatches(pattern)
if len(at_indexes) != 2 * monomers:
raise PolymerDescriptorError('Unexpected number of [At] indexes in '
f'oligomer:\n Expected {2 * monomers}, '
f'got {len(at_indexes)}.')
all_neighbor_indexes = [
mol.GetAtomWithIdx(at_index[0]).GetNeighbors()[0].GetIdx()
for at_index in at_indexes
]
# Add bonds
editable_mol = Chem.EditableMol(mol)
for idx in range(1, 2 * monomers - 1, 2):
editable_mol.AddBond(all_neighbor_indexes[idx],
all_neighbor_indexes[idx + 1],
order=Chem.rdchem.BondType.SINGLE)
# Remove all heavy atoms except the first and last ones
for idx in range(2, 2 * monomers):
editable_mol.RemoveAtom(at_indexes[-idx][0])
new_mol = editable_mol.GetMol()
# Alternative to sanitizing is converting to smiles and back to mol
Chem.SanitizeMol(new_mol)
return new_mol