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

Нарисуем ибупрофен:

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 Lipinksy
print(Lipinksy.NumHDonors(ibu))
print(Lipinksy.NumHAcceptors(ibu))
print(Lipinksy.rdMolDescriptors.CalcExactMolWt(ibu))
print(Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(ibu)[0])

1
1
206.130679816
3.073200000000001

Рассмотрим некоторые примеры, куда в ибупрофен можно ввести тройную связь для проведения клик-реакции:

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

Example 1

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

Example 2

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

Example 3

Нарисуем азид и то, как будут выглядеть соединения, представленные выше, после клик-реакции:

azide=Chem.MolFromSmiles('N=[N+]=[N-]')
AllChem.Compute2DCoords(azide)
display(azide)

ibu=Chem.MolFromSmiles('N1C=C(N=N1)C(C)CC1=CC=C(C=C1)C(C)C(=O)O')
AllChem.Compute2DCoords(ibu)
display(ibu)
print('Example 1')

Example 1

ibu=Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(C(C(N=N1)=CN1))C(=O)O')
AllChem.Compute2DCoords(ibu)
display(ibu)
print('Example 2')

Example 2

ibu=Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1(C(N=N1)=CN1))C(C)C(=O)O')
AllChem.Compute2DCoords(ibu)
display(ibu)
print('Example 3')

Example 3

Соберем наши шаблоны:

templ1 = 'N1C=C(N=N1)C(C)CC1=CC=C(C=C1)C(C)C(=O)O'
templ2 = 'CC(C)CC1=CC=C(C=C1)C(C(C(N=N1)=CN1))C(=O)O'
templ3 = 'CC(C)CC1=CC=C(C=C1(C(N=N1)=CN1))C(C)C(=O)O'

Проверим функцию внесения замены в исходное вещество (ибупрофен) на примере одного азида:

def replace(prev, azide, who):
    return azide.replace(who, prev) # 'N=[N+]=[N-]'

azide = 'c1cc(ccc1[N-][N+]#N)[N+](=O)[O-]'
azide=Chem.MolFromSmiles(azide)
AllChem.Compute2DCoords(azide)
display(azide)
azide = 'c1cc(ccc1[N-][N+]#N)[N+](=O)[O-]'
print('4-Nitrophenyl azide')

4-Nitrophenyl azide

mol = replace('N1C=C(N=N1)C(C)CC1=CC=C(C=C1)C(C)C(=O)O', azide, '[N-][N+]#N')
azide = mol
azide=Chem.MolFromSmiles(azide)
AllChem.Compute2DCoords(azide)
display(azide)

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

import pubchempy

import pubchempy as pcp

structures = []
per_page = 10**5

for smiles in ["N=N=N", "NN#N", "N#NN"]:
    for i in range(200):
        try:
            a = pcp.get_properties(
              properties="CanonicalSMILES", 
              identifier=smiles, namespace="smiles", 
              searchtype="substructure",
              RingsNotEmbedded=True,
              listkey_count=per_page, listkey_start=i*per_page
            )
        except:
            break
        print("Retrieved page {} of {} search".format(i+1, smiles))
        structures.extend(a)

Retrieved page 1 of N=N=N search
Retrieved page 2 of N=N=N search
Retrieved page 3 of N=N=N search
Retrieved page 4 of N=N=N search
Retrieved page 5 of N=N=N search
Retrieved page 6 of N=N=N search
Retrieved page 7 of N=N=N search
Retrieved page 8 of N=N=N search
Retrieved page 9 of N=N=N search
Retrieved page 10 of N=N=N search
Retrieved page 11 of N=N=N search
Retrieved page 12 of N=N=N search
Retrieved page 13 of N=N=N search
Retrieved page 14 of N=N=N search
Retrieved page 15 of N=N=N search
Retrieved page 16 of N=N=N search
Retrieved page 17 of N=N=N search
Retrieved page 18 of N=N=N search
Retrieved page 19 of N=N=N search
Retrieved page 1 of NN#N search
Retrieved page 1 of N#NN search

Создадим список smiles для более удобной работы:

compounds = {str(i['CID']).strip():i['CanonicalSMILES'].strip() for i in structures}
smiles = list(compounds.values())

Сохраним список в файл, на всякий случай:

import pickle
with open('/home/julybel/study/biopol/smiles', 'wb') as fp:
    pickle.dump(smiles, fp)

len(smiles)

1804790

Немного профильтруем наш список и возьмем небольшие соединения, smiles которых содержит группу N=[N+]=[N-]:

nsmiles = [elem for elem in smiles if len(elem) < 30 and '.' not in elem and 'N=[N+]=[N-]' in elem]

len(nsmiles)

231054

with open('/home/julybel/study/biopol/filtered-smiles', 'wb') as fp:
    pickle.dump(nsmiles, fp)

Теперь будем искать такие структуры, которые после присоединения к ним ибупрофена удовлетворяют 5 правилам Лепински.

def replace(prev, azide, who):
    return azide.replace(who, prev) # 'N=[N+]=[N-]'

template = 'N1C=C(N=N1)C(C)CC1=CC=C(C=C1)C(C)C(=O)O'
sample = []

for elem in nsmiles:
    
    mol = replace(template, elem, 'N=[N+]=[N-]')
    try:
        obj = Chem.MolFromSmiles(mol)
        if Lipinksy.NumHDonors(obj) <= 5 and Lipinksy.NumHAcceptors(obj) <= 10 and\
        Lipinksy.rdMolDescriptors.CalcExactMolWt(obj) <= 500 and \
        Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(obj)[0] <= 5 :
            sample.append(obj)
    except:
        pass

with open('/home/julybel/study/biopol/drugs', 'wb') as fp:
    pickle.dump(sample, fp)

Для дальнейшей работы отберем 40 структур из получившейся выборки

step = len(sample) // 40
sample40 = []
i = 0

for _ in range(40):
    sample40.append(sample[i])
    i += step

print(len(sample40))

40

Теперь нарисуем большую картинку с нашими 40 веществами, собранными по одному темплейту:

from IPython.display import SVG
Chem.Draw.MolsToGridImage(sample40, molsPerRow=4, subImgSize=(300,150))

Построим Similiraty Map ибупрофена с пятым веществом из нашего массива sample40:

from rdkit.Chem.Draw import SimilarityMaps

ibu = Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(C)C(=O)O')
fp = SimilarityMaps.GetMorganFingerprint(sample40[5], fpType='bv')
fig, maxweight = SimilarityMaps.GetSimilarityMapForFingerprint(ibu, sample40[5], SimilarityMaps.GetMorganFingerprint)

Мы видим, что добавился достаточно крупный довесок, который может, при случае, образовывать дополнительные водородные связи.

Вспомним, как выглядит ибупрофен:

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

И 5 вещество из нашего массива sample40:

sample40[5]

Теперь построим 3D структуру нашего лиганда:

import nglview as nv
import ipywidgets

m3d = Chem.AddHs(sample40[5])
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
display(m3d)

nv.show_rdkit(m3d)

..ничего нам не нарисовали