Домашнее задание №3


Загрузим необходимые пакеты.

In [1]:
from rdkit import Chem
import pubchempy as pcp
from rdkit.Chem import Draw, rdMolDescriptors, Descriptors, Crippen
import itertools

Загрузим из PubChem SMILES-идентификаторы всех соединений, содержащих подструктуры изоморфные мезомерам азидной группы.

In [2]:
smiles_list = ['N=[N+]=[N-]', '[N-][N+]#N', '[N-]N=[N+]']
STEP_SIZE = 10000
output = []

for smiles in smiles_list:
    for i in itertools.count():
        try:
            result = pcp.get_properties(properties="CanonicalSMILES", 
                                    identifier=smiles, 
                                    namespace="smiles",
                                    searchtype="substructure", 
                                    RingsNotEmbedded=True, 
                                    MatchCharges=True, 
                                    listkey_count=STEP_SIZE, 
                                    listkey_start=i * STEP_SIZE)
            output += result
        except:
            break
In [3]:
output = list(set(map(lambda x: x['CanonicalSMILES'], output)))
print(len(output))

119612

Загрузим молекулу ибупрофена, заменим в ней изопропил на ацетилен и пометим в измененном ибупрофене сайт будущей реакции, то есть те два атома, к которым в ходе реакции CuAAC будут присоединяться новые связи.

In [4]:
def label_reaction_site(molecule):
    patt = Chem.MolFromSmarts("C#C")
    matching = molecule.GetSubstructMatch(patt)
    if molecule.GetAtomWithIdx(matching[0]).GetDegree() == 1:
        molecule.GetAtomWithIdx(matching[0]).SetAtomMapNum(1)
        molecule.GetAtomWithIdx(matching[1]).SetAtomMapNum(2)
    else:
        molecule.GetAtomWithIdx(matching[1]).SetAtomMapNum(1)
        molecule.GetAtomWithIdx(matching[0]).SetAtomMapNum(2)
        
ibuprophen = Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(C)C(=O)O')
altered_ibuprophen = Chem.MolFromSmiles('C#CCC1=CC=C(C=C1)C(C)C(=O)O')
label_reaction_site(altered_ibuprophen)
Draw.MolsToGridImage([ibuprophen,altered_ibuprophen], molsPerRow=2, subImgSize=(400,400), legends = ['Ibuprophen','Altered Ibuprophen'])
Out[4]:

Для каждого скачанного из PubChem соединения создается RDKit-овский Mol объект. Теперь мы можем работать с соединениями, как с графами. Для каждого соединения найдем все содержащиеся в нем азидные группы и для каждой такой группы осуществим реакцию CuAAC между данным соединением и измененным ибупрофеном по данной азидной группе.

In [5]:
def is_good_triple(molecule, triple):
    charges = list(map(lambda x: x.GetFormalCharge(), triple))
    first_bond_type = molecule.GetBondBetweenAtoms(triple[0].GetIdx(), triple[1].GetIdx()).GetBondType()
    second_bond_type = molecule.GetBondBetweenAtoms(triple[1].GetIdx(), triple[2].GetIdx()).GetBondType()
    if charges == [-1, 1, 0]:
        return first_bond_type == Chem.rdchem.BondType.DOUBLE and second_bond_type == Chem.rdchem.BondType.DOUBLE
    elif charges == [0, 1, -1]:
        return first_bond_type == Chem.rdchem.BondType.TRIPLE and second_bond_type == Chem.rdchem.BondType.SINGLE
    elif charges == [1, 0, -1]:
        return first_bond_type == Chem.rdchem.BondType.DOUBLE and second_bond_type == Chem.rdchem.BondType.SINGLE
    
    return False

def label_atoms(atoms, labels):
    for label, atom in zip(labels, atoms):
        atom.SetAtomMapNum(label)

def create_product(alt_ibuprophen, mol):
    product = Chem.CombineMols(alt_ibuprophen, mol)
    reaction_site = [0] * 5
    for atom in product.GetAtoms():
        if atom.GetAtomMapNum() > 0:
            reaction_site[atom.GetAtomMapNum() - 1] = atom.GetIdx()
    
    for i in range(2, 5):
        product.GetAtomWithIdx(reaction_site[i]).SetFormalCharge(0)
    edproduct = Chem.EditableMol(product)
    edproduct.AddBond(reaction_site[0], reaction_site[4], order=Chem.rdchem.BondType.SINGLE)
    edproduct.AddBond(reaction_site[1], reaction_site[2], order=Chem.rdchem.BondType.SINGLE)
    edproduct.RemoveBond(reaction_site[0], reaction_site[1])
    edproduct.AddBond(reaction_site[0], reaction_site[1], order=Chem.rdchem.BondType.DOUBLE)
    edproduct.RemoveBond(reaction_site[2], reaction_site[3])
    edproduct.AddBond(reaction_site[2], reaction_site[3], order=Chem.rdchem.BondType.DOUBLE)
    edproduct.RemoveBond(reaction_site[3], reaction_site[4])
    edproduct.AddBond(reaction_site[3], reaction_site[4], order=Chem.rdchem.BondType.SINGLE)
    
    product = edproduct.GetMol()
    for atom in product.GetAtoms():
        atom.SetAtomMapNum(0)
    
    return product
        
def main():
    for smiles in output:
        mol = Chem.MolFromSmiles(smiles)
        if mol is None:
            continue
        for atom in mol.GetAtoms():
            triple = []
            if atom.GetDegree() == 1 and atom.GetSymbol() == "N":
                neighbour_1 = atom.GetNeighbors()[0]
                if neighbour_1.GetDegree() == 2 and neighbour_1.GetSymbol() == "N":
                    neighbours = neighbour_1.GetNeighbors()
                    neighbour_2 = (neighbours[0] if neighbours[0] != neighbour_1 else neighbours[1])
                    if neighbour_2.GetDegree() == 2 and neighbour_2.GetSymbol() == "N":
                        triple.append(atom)
                        triple.append(neighbour_1)
                        triple.append(neighbour_2)

            if len(triple) == 3 and is_good_triple(mol, triple):
                label_atoms(triple, [3, 4, 5])
                products.append(create_product(altered_ibuprophen, mol))
                label_atoms(triple, [0, 0, 0])

products = []
main()

В силу того, что молекулярные графы некоторых соединений, скачанных из pubchem, обладают нетривиальной (то есть отличной от группы, содержащий лишь тождественный автоморфизм) группой автоморфизмов, то при присоединении измененного ибупрофена к разным азидным группам таких соединений могут получаться одинаковые продукты. Дубликаты следует отсеить.

In [6]:
import networkx as nx
import networkx.algorithms.isomorphism as iso

class WrappedMol:
    def __init__(self, molecule):
        self.rdkit_molecule = molecule
        self.networkx_graph = self.__rdkit_to_networkx(molecule)
        self.__edge_match = iso.numerical_edge_match('type', -1)
        self.__node_match = iso.categorical_node_match('label', '')
        
    def __eq__(self, other):
        return nx.is_isomorphic(self.networkx_graph, other.networkx_graph, 
                                node_match=self.__node_match, edge_match=self.__edge_match)
    
    def __hash__(self):
        return hash(Chem.MolToSmiles(self.rdkit_molecule))
    
    def __rdkit_to_networkx(self, rdkit_molecule):
        networkx_graph = nx.Graph()
        for index, atom in enumerate(rdkit_molecule.GetAtoms()):
            atom.SetAtomMapNum(index + 1)
        for atom in rdkit_molecule.GetAtoms():
            networkx_graph.add_node(atom.GetAtomMapNum(), label=atom.GetSymbol() + str(atom.GetFormalCharge()))
        for bond in rdkit_molecule.GetBonds():
            networkx_graph.add_edge(bond.GetBeginAtom().GetAtomMapNum(), bond.GetEndAtom().GetAtomMapNum(),
                                    type=bond.GetBondTypeAsDouble())
        for atom in rdkit_molecule.GetAtoms():
            atom.SetAtomMapNum(0)
        
        return networkx_graph

products = list(map(lambda x: x.rdkit_molecule, set(map(lambda y: WrappedMol(y), products))))

Отфильтурем оставшееся множество соединений по правилу Липински.

In [7]:
filtered_products = []
for product in products:
    cnt = 0
    Chem.GetSSSR(product)
    if rdMolDescriptors.CalcNumHBD(product) <= 5:
        cnt += 1
    if rdMolDescriptors.CalcNumHBA(product) <= 10:
        cnt += 1
    if Descriptors.ExactMolWt(product) <= 500:
        cnt += 1
    if Crippen.MolLogP(product) <= 5:
        cnt += 1
    if cnt >= 3:
        filtered_products.append(product)

Отобразим некоторые из оставшихся соединений.

In [17]:
Draw.MolsToGridImage(filtered_products[:10], molsPerRow=3)
Out[17]: