Хемоинформатика¶
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit import RDConfig
from rdkit.Chem.Draw import IPythonConsole, SimilarityMaps
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 Lipinski
print (Lipinski.NumHDonors(ibu))
print(Lipinski.NumHAcceptors(ibu))
print(Lipinski.rdMolDescriptors.CalcExactMolWt(ibu))
print(Lipinski.rdMolDescriptors.CalcCrippenDescriptors(ibu)[0])
Правила Липински:
- No more than 5 hydrogen bond donors
- No more than 10 hydrogen bond acceptors
- A molecular mass less than 500 daltons
- An octanol-water partition coefficient log P not greater than 5
Модифицируем ибупрофен для клик-химии
modified_ibu = 'N1N=NC(=C1)C(C)CC1=CC=C(C=C1)C(C)C(=O)O'
diene_ibu = Chem.MolFromSmiles('{}'.format(modified_ibu))
AllChem.Compute2DCoords(diene_ibu)
display(diene_ibu)
Будем присоединять к нему азиды
Для этого скачаем структуры азидных-соединений из базы PubChem (я делала это вручную)
PubChem -> Substructure/Superstructure -> CID, SMILES/SMARTS, InChI Вписываем в поле поиска (CID азида: 33558) Options: Substructure И поставим фильтры а массу и доноров/акцепторов. Исходя из параметров нашей затравки, я решил ставить ограничение по массе <=300 Да, по донорам <=10 и по акцепторам <= 15 (на всякий случай с запасом).
ограничение по массе <=300 Да (т.к. mw фрагмента, который будем присоединять ~200, а в сумме надо получить <500), доноры водородной связи <=10 и акцепторы <= 15 (с запасом).
Всего получилось найти 558093 соединений
Скачаем найденные соединения Actions on your result -> download structures -> smiles
#Отфильтруем соединения
# "not '.' in line" == убираем соединения с нековалентными связями
strings=np.genfromtxt('Molsim/my_azides.txt',dtype=np.str)
smiles = []
for line in strings:
if (len(line[1]) < 30) and (len(line[1]) > 11) and not ( '.' in line[1]):
smiles.append(line[1])
len (smiles) #кол-во отфильтрованных соединений
Теперь построим новые молекулы. Сразу будем проверять, соответствуют ли они правилу Липински
def check_Lipinski(mol):
is_Lipinsky = bool((Lipinski.NumHDonors(mol) <= 5) \
and (Lipinski.NumHAcceptors(mol) <= 10) \
and (Lipinski.rdMolDescriptors.CalcExactMolWt(mol) < 500) \
and(Lipinski.rdMolDescriptors.CalcCrippenDescriptors(mol)[0]) <= 5)
return is_Lipinsky
products = []
for smi in smiles[:len(smiles)]:
if "N=[N+]=[N-]" in smi:
newsmi = smi.replace("N=[N+]=[N-]", modified_ibu)
else:
continue
try:
newmol = Chem.MolFromSmiles(newsmi)
if check_Lipinski(newmol):
products.append(newmol)
except Exception:
pass
Есть какие-то ошибки, но вообще, что-то отфильтровалось:
print(len(products)) #сколько получилось соединений
Посмотрим, как выглядят наши соединения
import random
molecules = []
for pics in random.sample(products, 25):
AllChem.Compute2DCoords(pics)
molecules.append(pics)
Chem.Draw.MolsToGridImage(molecules, molsPerRow = 5, subImgSize = (500, 500))
Построим Similiraty Map ибупрофена с пятым веществом из полученного массива.
Chem.MolToSmiles(molecules[0])
ibu = Chem.MolFromSmiles('CC(C)CC1=CC=C(C=C1)C(C)C(=O)O')
try_mol = Chem.MolFromSmiles('COCCC(=O)n1cc(C(C)Cc2ccc(C(C)C(=O)O)cc2)nn1')
#fig, maxweight = SimilarityMaps.GetSimilarityMapForFingerprint(ibu,try_mol, SimilarityMaps.GetMorganFingerprint)
Я пыталась, но у меня не получилось постоить карту сходства ;(
Построим 3D структуру
m3d = Chem.AddHs(molecules[0])
AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d, maxIters = 500, nonBondedThresh = 200)
m3d
import nglview as nv
nv.show_rdkit(m3d)
# не ругается, но и не работает
#записываем в виде smiles
smiles_to_pymol = Chem.MolToSmiles(m3d)
with open('new_mol.smi', 'w') as file:
file.write(smiles_to_pymol)
# с помощью open-babel переводим в 3d формат
import subprocess
subprocess.run('obgen new_mol.smi > new_mol.mol', shell = True)
Картиночку сделала в Pymol
Image("Molsim/my_mol.png")
