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

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

Вспомним параметры для правила Лепински: соединение, чтобы быть похожим на лекарство, должно:

• Иметь менее 5 атомов-доноров водородной связи;

• Обладать молекулярным весом менее 500;

• Иметь липофильность (log P — коэффициент распределения вещества на границе раздела вода-октанол) менее 5;

• Иметь суммарно не более 10 атомов азота и кислорода (грубая оценка количества акцепторов водородной связи)

Теперь модифицируем ибупрофен так, чтобы к нему можно было присоединить азид, то есть сделаем ему концевую алкильную группу. А также рассмотрим то, что должно аолучиться в результате клик-реакции.

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

ibuzid=Chem.MolFromSmiles('N1C=C(N=N1)C(C)CC1=CC=C(C=C1)C(C)C(=O)O')
AllChem.Compute2DCoords(ibuzid)
display(ibuzid)

Скачаем как можно больше различных азидов.

import pubchempy as pcp

compounds = []
per_page = 1000
for smiles in ["N=[N+]=[N-]", "N=N=N", "NN#N"]:
    for i in range(5):
        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:
            print ('oh')
        print ("Retrieved page {} of {} search".format(i+1, smiles))
        compounds.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 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 1 of NN#N search
Retrieved page 2 of NN#N search
Retrieved page 3 of NN#N search
Retrieved page 4 of NN#N search
Retrieved page 5 of NN#N search

smiles_list = [x["CanonicalSMILES"] for x in compounds]
print(len(smiles_list))
smiles_list = list(filter(lambda x: len(x) < 30 and not '.' in x,
                          smiles_list))
print(len(smiles_list))

15000
315

from random import shuffle, seed

groups = ['N=[N+]=[N-]', '[N-]=[N+]=[N-]', '[N-][N+]#[N]', '[N]#[N+][N-]', '[N-][N]#[N+]', '[N+]#[N][N-]']
template = 'N1C=C(N=N1)C(C)CC1=CC=C(C=C1)C(C)C(=O)O'

result = []
num_to_check = 2000

# Новую молекулу лучше создавать в try из-за битых Smiles
#При первом запуске ячейки убрать решетки, чтобы перемешать
#shuffle(smiles_list)
for smi in smiles_list[:num_to_check]:
    for group in groups:
        if group in smi:
            newsmi = smi.replace(group, template)
            break
    else:
        continue
# Если новая молекула удолтворяет правилу пяти, сохраним в массив и покажем
    try:
        new_mol=Chem.MolFromSmiles(newsmi)
        if (Lipinksy.NumHDonors(new_mol) <= 5 and 
            Lipinksy.NumHAcceptors(new_mol) <= 10 and 
            Lipinksy.rdMolDescriptors.CalcExactMolWt(new_mol) <=500 and 
            Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(new_mol)[0] <= 5):
                result.append((new_mol, Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(new_mol)[0]))
    except Exception:
        continue

print(len(result))

result.sort(key=lambda x: x[1])
mols = [x[0] for x in result]

140

Посмотрим на топ-12 гидрофильных молекул. Вот они:

Draw.MolsToGridImage(mols[:12],useSVG=True, molsPerRow=3, subImgSize=(400, 320))

Посмотрим на 5-го представителя из наиболее гидрофильных молекул.

from rdkit.Chem.Draw import SimilarityMaps
somemol = mols[4]
SimilarityMaps.GetMorganFingerprint(somemol, fpType='bv')
SimilarityMaps.GetSimilarityMapForFingerprint(ibu,
                    somemol, SimilarityMaps.GetMorganFingerprint)

(<Figure size 250x250 with 1 Axes>, 0.19272727272727275)

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

Теперь получим 3D структуру этой молекулы и покрутим ее.

m3d=Chem.AddHs(somemol)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )

0

import nglview as nv
nv.show_rdkit(m3d)

NGLWidget()

К сожалению, модуль выдает ошибку при импорте, из-за чего покрутить модель в html невозможно.

Поэтому просто отобразим снимок этой молекулы.

display(m3d)