Домашнее задание 3¶
Выполнил: Ким Адамейко, группа мАДБМ16
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(color_codes=True)
import pubchempy as pcp
import py3Dmol
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 rdkit.Chem.Lipinski as Lip
from IPython.display import display,Image
def show_compound_info(name, from_smiles=False):
if from_smiles:
smiles = name
else:
pc_compound = pcp.get_compounds(name, 'name')[0]
smiles = pc_compound.canonical_smiles.encode("ascii")
print('Compound: %s\n IUPAC: %s\n SMILES: %s' % (name, pc_compound.iupac_name, smiles))
rd_compound = Chem.MolFromSmiles(smiles)
AllChem.Compute2DCoords(rd_compound)
return rd_compound
Часть 1. Модификация ибупрофена¶
Получим с помощью библиотеки pubchempy smiles-представление ибупрофена и визуализируем его с помощью библиотеки rdkit
ibuprofen_smiles = pcp.get_compounds('ibuprofen', 'name')[0].canonical_smiles.encode("ascii")
ibuprofen = show_compound_info('ibuprofen')
ibuprofen
Получим также информацию об этине, на который нужно заменить изопропил
show_compound_info('ethyne')
Самого изопропила в базе нет, но есть изопропиловый спирт
show_compound_info('isopropyl alcohol')
Отсюда, а также IUPAC-имени можно сделать вывод о требуемой замене:
mod_ibuprofen_smiles = ibuprofen_smiles.replace('CC(C)', 'C#C')
print(mod_ibuprofen_smiles)
mod_ibuprofen = Chem.MolFromSmiles(mod_ibuprofen_smiles)
AllChem.Compute2DCoords(mod_ibuprofen)
display(mod_ibuprofen)
Часть 2. Поиск радикалов с азидом¶
Википедия сообщает, что азид имеет несколько мезомеров. Я подготовил SMILES-нотации для этих трёх видов с учётом возможного написания в двух "направлениях".
azide_smiles = ['N=[N+]=[N-]', '[N-][N+]#N', '[N-]N=[N+]', '[N-]=[N+]=N', 'N#[N+][N-]', '[N+]=N[N-]']
for a_s in azide_smiles[:3]: # показано только одно "направление" записи, второе -- симметрично
display(show_compound_info(a_s, from_smiles=True))
Реакция идёт по схеме: $\require{mhchem} \ce{R1-N=N+=N- + R2-C≡C->[Cu]R1-N^{(1)}=N-N-C=C^{(1)}-R2}$
где в качестве $\ce{R2-C≡C}$ будет выступать приведенный выше модифицированный ибупрофен
Что будет присоединено к $R_1$ после реакции (вместо радикала на рисунке автоматически добавлен водород):
ibuprofen_res_smiles = mod_ibuprofen_smiles.replace('C#C', 'N2C=C(N=N2)')
print(ibuprofen_res_smiles)
ibuprofen_res = Chem.MolFromSmiles(ibuprofen_res_smiles)
AllChem.Compute2DCoords(ibuprofen_res)
display(ibuprofen_res)
Для поиска по базе PubChem уберем из образцов информацию о заряде, получим следующие:
azide_smiles_for_search = set(map(lambda x: re.sub(r'[\[\]\+\-]', '', x), azide_smiles))
azide_smiles_for_search
dfs_azides = []
for a_s in azide_smiles_for_search:
page_no = 0
results_per_page = 100000
while True:
try:
df_tmp = pcp.get_properties(
properties="CanonicalSMILES", namespace="smiles", identifier=a_s,
searchtype="substructure", RingsNotEmbedded=True, as_dataframe=True,
listkey_count=results_per_page, listkey_start=page_no*results_per_page)
dfs_azides.append(df_tmp)
page_no += 1
except:
break
df_azides = pd.concat(dfs_azides)
del dfs_azides
print('Found %d items' % df_azides.shape[0])
df_azides.head()
Сразу отбросим вещества, в SMILES нотации которых есть точка, обозначающая нековалентную связь:
df_azides.drop(
df_azides[df_azides.CanonicalSMILES.str.contains('.', regex=False)].index,
inplace=True)
print('% d items left' % df_azides.shape[0])
Посчитаем, сколько разных представлений азида встречается в каждой молекуле (будем считать два "симметричных представления одним). Затем сразу удалим строки, где все счётчики нулевые.
for i in range(3):
df_azides.loc[:,'cnt_a' + str(i)] = \
df_azides.CanonicalSMILES.str.count(re.escape(azide_smiles[i]) + '|' + re.escape(azide_smiles[i+3]))
df_azides.head()
df_azides.drop(
df_azides[df_azides.iloc[:, -3:].sum(axis=1) == 0].index,
inplace=True)
print('% d items left' % df_azides.shape[0])
Посчитаем ещё примерное количество атомов путём подсчёта только алфавитных символов в строках
df_azides.loc[:,'atoms_cnt'] = df_azides.CanonicalSMILES.str.count('\w')
df_azides.head()
Можно построить несколько графиков и таблиц, чтобы оценить разнообразие найденных результатов. Распределение примерного количества атомов:
ax = sns.kdeplot(df_azides.atoms_cnt)
ax.set_xlabel('approx. atoms count')
plt.show()
Суммарные значения счётчиков трёх представлений азидов:
df_azides.iloc[:, -4:-1].sum(axis=0)
Видим, что можно пренебречь последним представлением. Как распределяются представления в комбинациях? Видим, что нетривиальные комбинации разных типов встречаются крайне редко, и более 80% комбинаций это единственное в строке появление канонической комбинации
N=[N+]=[N-].
df_counts = pd.DataFrame(df_azides.groupby(list(df_azides.columns[-4:-2])).size(), columns=['count'])
plt.scatter(*zip(*df_counts.index.tolist()), s=3*np.sqrt(df_counts['count']), alpha=0.66)
plt.xlabel(azide_smiles[0]); plt.ylabel(azide_smiles[1]);
plt.show()
df_counts
Наконец, применим "правила пяти" Липински
def Lipinksi_rules(compound):
return ((Lip.NumHDonors(compound) <= 5) # <= 5 hydrogen bond donors
and (Lip.NumHAcceptors(compound) <= 10) # <= 10 hydrogen bond acceptors
and (Lip.rdMolDescriptors.CalcExactMolWt(compound)<500) # mol.mass < 500 daltons
and (Lip.rdMolDescriptors.CalcCrippenDescriptors(compound)[0]<=5)) # octanol-water partit.coeff. logP < 5
azide_full_re = re.compile('|'.join(map(re.escape, azide_smiles)))
df_azides.loc[:, 'resSMILES'] = df_azides.CanonicalSMILES.str.replace(azide_full_re, ibuprofen_res_smiles)
from rdkit import rdBase
rdBase.DisableLog('rdApp.error')
results = []
for res_smile in df_azides.resSMILES:
try:
compound = Chem.MolFromSmiles(res_smile)
if Lipinksi_rules(compound):
results.append(res_smile)
except:
pass
print('%d items left after filtering by Lipinski rules' % len(results))
Сохраним результаты:
import pickle
import gzip
with gzip.open('results.pklz', 'wb') as f:
pickle.dump(results, f)
Отобразим некоторые из них в 2D и в 3D:
indexes = random.sample(range(len(results)), 12)
images = [Chem.MolFromSmiles(results[i]) for i in indexes]
display(Chem.Draw.MolsToGridImage(images, molsPerRow=3, subImgSize=(300,200)))
К сожалению, у меня пока не получилось заставить работать nglviewer, зато я разобрался, как использовать другую библиотеку -- py3Dmol, которая является обёрткой над известной библиотекой 3Dmol.js. Она легко устанавливается с помощью команды pip install py3Dmol
def show_compound_3d(smiles, view, grid=(0, 0)):
m2d = Chem.MolFromSmiles(smiles)
m3d = Chem.AddHs(m2d)
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
mb = Chem.MolToMolBlock(m3d)
view.addModel(mb,'sdf', viewer=grid)
view.setStyle({'model':0},{'stick': {}}, viewer=grid)
view.zoomTo(viewer=grid)
return view
# view = py3Dmol.view(width=330, height=250) # образец кода для просмотра одной молекулы
# show_compound_3d(results[0], view)
view = py3Dmol.view(width=320*3, height=250*4, linked=False, viewergrid=(4,3))
for i in range(4):
for j in range(3):
show_compound_3d(results[indexes[i*3 + j]], view, grid=(i, j))
view.render()
Источники¶
- Azide - Wikipedia
- Примеры использования py3Dmol в Jupyter
- Реализация метода IPythonConsole.drawMol3D(), которая почему-то отказывалась работать. Там можно увидеть, как загрузить молекулу в просмотрщик из файла, причём двумя способами -- для больших и малых молекул.
- Полезный ресурс для отладки регулярных выражений