Хемоинформатика

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import RDConfig
from rdkit.Chem.Draw import IPythonConsole 
from rdkit.Chem import Draw
import numpy as np
from IPython.display import display,Image

RDKit WARNING: [18:43:07] Enabling RDKit 2019.09.3 jupyter extensions

Сконструируем и визуализируем объект — молекулу ибупрофена.

ibu=Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(C)C(=O)O')
AllChem.Compute2DCoords(ibu)
display(ibu)

Проверим правила Липински.

import rdkit.Chem.Lipinski as Lipinksi
print(Lipinksi.NumHDonors(ibu))
print (Lipinksi.NumHAcceptors(ibu))
print (Lipinksi.rdMolDescriptors.CalcExactMolWt(ibu))
print (Lipinksi.rdMolDescriptors.CalcCrippenDescriptors(ibu)[0])

1
1
206.130679816
3.073200000000001

Правила выполнены. Далее поработаем с pubchempy.

import pubchempy as pcp

Реализуем процедуру для скачивания соединений.

def get_compounds(smiles, listkey_count=10**3):
    result = []
    
    for i in range(200):
        try:
            extension = pcp.get_properties(
                                            properties="CanonicalSMILES", 
                                            identifier=smiles, namespace="smiles", 
                                            searchtype="substructure",
                                            RingsNotEmbedded=True,
                                            listkey_count=listkey_count, listkey_start=i*listkey_count,
                                          )
            print("Retrieved page {} of {} search".format(i+1, smiles))
        except:
            break
        result.extend(extension)
    return result

Потом я их долго скачивал и сохранил в файл (в дальнейшем загружал из файла).

%%time all_smiles = ["N=N=N", "NN#N", "N=[N+]=[N-]"] all_compounds = [] for smiles in all_smiles: all_compounds.extend(get_compounds(smiles))

from json import load 
with open("compounds.json", "r") as infile:
    all_compounds = load(infile)

len(all_compounds)

405970

# from json import dump
# with open("compounds.json", "w") as outfile:
#     dump(all_compounds, outfile)

all_compounds[0]

{'CID': 162624792,
 'CanonicalSMILES': 'CCOC(=O)C(C(OCCOCCN=[N+]=[N-])OC1C(C(C(C(O1)COC(=O)C)OC(=O)C)OC(=O)C)OC(=O)C)Br'}

Попробуем модифицировать ибупрофен, чтобы было похоже на продукт клик-химии. Постараемся расположить три азота так, чтобы можно было заменять азид на эту молекулу в виде подстроки SMILES и мало что ломалось. Сделаем, чтобы запись SMILES с этих азотов начиналась.

template_smiles = 'N2N=NC=C2C(C1=CC=C(C=C1)CC(C)C)C(=O)O'
template = Chem.MolFromSmiles(template_smiles)
AllChem.Compute2DCoords(template)
display(template)

Напишем процедуру для проверки правил Липински.

def check_lipinski(mol):
    rule1 = Lipinksi.NumHDonors(mol) <= 5
    rule2 = Lipinksi.NumHAcceptors(mol) <= 10
    rule3 = Lipinksi.rdMolDescriptors.CalcExactMolWt(mol) <= 500
    rule4 = Lipinksi.rdMolDescriptors.CalcCrippenDescriptors(ibu)[0] <= 5
    return rule1 and rule2 and rule3 and rule4

Нагенерируем продукты клик-химии.

%%time
transformed_compounds = []

for cmp in all_compounds[:1500]:
    transformed_smiles = cmp["CanonicalSMILES"]
    was_transformed = False
    if "N=N=N" in cmp["CanonicalSMILES"]:
        transformed_smiles = transformed_smiles.replace("N=N=N", template_smiles)
        was_transformed = True
    if "NN#N" in transformed_smiles:
        transformed_smiles = transformed_smiles.replace("NN#N", template_smiles)
        was_transformed = True
    if "N=[N+]=[N-]" in transformed_smiles:
        transformed_smiles = transformed_smiles.replace("N=[N+]=[N-]", template_smiles)
        was_transformed = True
    
    try:
        newmol = Chem.MolFromSmiles(transformed_smiles)
        Chem.SanitizeMol(newmol) # to raise errors
    except:
        continue
    if check_lipinski(newmol):
        transformed_compounds.append(newmol)

len(transformed_compounds)

359

Визуализируем несколько первых продуктов.

Draw.MolsToGridImage(transformed_compounds[:21],useSVG=True, molsPerRow=3, subImgSize=(300, 300))

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

from rdkit.Chem.Draw import SimilarityMaps

fp = SimilarityMaps.GetMorganFingerprint(transformed_compounds[4], fpType='bv')

fig, maxweight = SimilarityMaps.GetSimilarityMapForFingerprint(ibu, transformed_compounds[4],
                                                               SimilarityMaps.GetMorganFingerprint)

И визуализируем эту молекулу в 3D.

m3d = Chem.AddHs(transformed_compounds[4])
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )

0

import nglview as nv

nv.show_rdkit(m3d)