Практикум 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

RDKit WARNING: [23:47:38] 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 Lipinksy
print(f"#H-donors is {Lipinksy.NumHDonors(ibu)}")
print(f"#H-acceptors is {Lipinksy.NumHAcceptors(ibu)}")
print(f"Molecular weight is {Lipinksy.rdMolDescriptors.CalcExactMolWt(ibu)}")
print(f"LogP is {Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(ibu)[0]}")

#H-donors is 1
#H-acceptors is 1
Molecular weight is 206.130679816
LogP is 3.073200000000001

import pubchempy as pcp

Скачаем SMILES некоторого количества азидов из PubChem

smiles = []
per_page = 10
n_pages = 10
for azide_smiles in ["[N-][N+]#N", 'N=[N+]=[N-]']:
    for i in range(n_pages):     
        try:
            query = pcp.get_properties(
            properties="CanonicalSMILES", 
            identifier=azide_smiles, 
            namespace="smiles", 
            searchtype="substructure",
            listkey_count=per_page, listkey_start=i*per_page),
            smiles.extend(query[0])
        except:
            print("Failed successfully")
        print(f"Retrieved page {i+1} of {azide_smiles} search, total azides found: {len(smiles)}")

Retrieved page 1 of [N-][N+]#N search, total azides found: 10
Retrieved page 2 of [N-][N+]#N search, total azides found: 20
Retrieved page 3 of [N-][N+]#N search, total azides found: 30
Retrieved page 4 of [N-][N+]#N search, total azides found: 40
Retrieved page 5 of [N-][N+]#N search, total azides found: 50
Retrieved page 6 of [N-][N+]#N search, total azides found: 60
Retrieved page 7 of [N-][N+]#N search, total azides found: 70
Retrieved page 8 of [N-][N+]#N search, total azides found: 80
Retrieved page 9 of [N-][N+]#N search, total azides found: 90
Retrieved page 10 of [N-][N+]#N search, total azides found: 100
Retrieved page 1 of N=[N+]=[N-] search, total azides found: 110
Retrieved page 2 of N=[N+]=[N-] search, total azides found: 120
Retrieved page 3 of N=[N+]=[N-] search, total azides found: 130
Retrieved page 4 of N=[N+]=[N-] search, total azides found: 140
Retrieved page 5 of N=[N+]=[N-] search, total azides found: 150
Retrieved page 6 of N=[N+]=[N-] search, total azides found: 160
Retrieved page 7 of N=[N+]=[N-] search, total azides found: 170
Retrieved page 8 of N=[N+]=[N-] search, total azides found: 180
Retrieved page 9 of N=[N+]=[N-] search, total azides found: 190
Retrieved page 10 of N=[N+]=[N-] search, total azides found: 200

smiles_only = [x['CanonicalSMILES'] for x in smiles]

print(f"Total duplicates: {len(smiles_only) - len(set(smiles_only))}")
print(f"Duplicates in [N-][N+]#N search : {100 - len(set(smiles_only[:100]))}")
print(f"Duplicates in N=[N+]=[N-] search : {100 - len(set(smiles_only[100:]))}")

Total duplicates: 4
Duplicates in [N-][N+]#N search : 4
Duplicates in N=[N+]=[N-] search : 0

Почему-то нашлись дубликаты. Скорее всего, дело в том, что canonical SMILES, которые мы получили, могут соответствовать нескольким стереоизомерам.

from collections import Counter
counts = Counter(smiles_only)
counts.most_common(4)

[('C1C(C(C(CO1)OC(C2=CC=CC=C2)(C3=CC=CC=C3)C4=CC=CC=C4)O)N=[N+]=[N-]', 3),
 ('CC(=O)OCC1C(C(C(C(O1)OC2C(C(OC3C2OC(OC3)C4=CC=CC=C4)OC5=CC=C(C=C5)OC)N=[N+]=[N-])OC(=O)C)OC(=O)C)OC(=O)C',
  2),
 ('CC(=O)OC1C(COCC1OC(C2=CC=CC=C2)(C3=CC=CC=C3)C4=CC=CC=C4)N=[N+]=[N-]', 2),
 ('CC(=O)OCC(=O)NC1C(CC(OC1C(C(CN=[N+]=[N-])OC(=O)C)OC(=O)C)(C(=O)OC)SC2=CC=CC=C2)OC(=O)C',
  1)]

Пытался по Canonical SMILES найти соединение, которое мы нашли 3 раза, на сайте PubChem. PubChem выдал 0 Identity находок. Загадка.
По 7 находкам в Substructure стало понятно, что этот smiles соотвествует модифицированному сахару. А сахара славятся тем, что у них много стереоизомеров.

Построим продукты реакции азид-алкильного циклоприсоединения

В качестве алкина возьмем вот такую молекулу

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

Перейдем к азидам. Внезапно, есть SMILES, которые не содержат азидной группы вообще. Выглядят они как соединения с тремя последовательно связанными азотами.

i = 0
j = 0
for smi in smiles_only:
    if "N=[N+]=[N-]" in smi:
        i += 1
    if "[N-][N+]#N" in smi:
        j += 1
print(f"Mols with N=[N+]=[N-]: {i},\nMols with [N-][N+]#N: {j}")
print(f"Total unique smiles: {len(set(smiles_only))}")

Mols with N=[N+]=[N-]: 99,
Mols with [N-][N+]#N: 0
Total unique smiles: 196

Пример молекулы, которая нашлась, но азидной группы не содержит:

Chem.MolFromSmiles(smiles_only[160])

Будем проверять наличие азидной группы перед тем как создавать продукт клик-реакции

products = []
for smi in set(smiles_only):
    if "N=[N+]=[N-]" in smi:
        newsmi=smi.replace('N=[N+]=[N-]',
                           'N8-N=NC(=C8)c9c(C(C)C(=O)O)ccc(CC(C)C)c9')
        
        try:
            newmol=Chem.MolFromSmiles(newsmi)

            HD_count = Lipinksy.NumHDonors(newmol)
            HAcc_count = Lipinksy.NumHAcceptors(newmol)
            MW = Lipinksy.rdMolDescriptors.CalcExactMolWt(newmol)
            LogP = Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(newmol)[0]

            if (HD_count <= 5) and \
                (HAcc_count <= 10) and \
                (MW <= 500) and LogP <= 5:
                products.append(newmol)
        except:
            print("Failed successfully")
    else:
        continue

print(f"Удовлетворяют правилу Лепински {len(products)} соединения из 99")

Удовлетворяют правилу Лепински 24 соединения из 99

Посмотрим на наши продукты

Draw.MolsToGridImage(products, 
                         useSVG=True,
                         molsPerRow=2,
                         subImgSize=(400,200))

Строим карту сходства ибупрофена с пятым соединением, которое мы получили

from rdkit.Chem.Draw import SimilarityMaps
fp = SimilarityMaps.GetMorganFingerprint(products[4], fpType='bv')
fig, maxweight = SimilarityMaps.GetSimilarityMapForFingerprint(ibu, products[4], SimilarityMaps.GetMorganFingerprint)

Перейдем от 2Д к 3Д. Для этого надо добавить водороды и оптимизировать геометрию. Для оптимизации используем поле MMFF.

mol = products[4]
mol3d=Chem.AddHs(mol)
Chem.AllChem.EmbedMolecule(mol3d)
AllChem.MMFFOptimizeMolecule(mol3d, maxIters=500, 
                             nonBondedThresh=200 )

0

Посмотрим на то, что мы натворили

import nglview as nv
nv.show_rdkit(mol3d)