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

0. Импорты

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
from IPython.display import Image
import rdkit.Chem.Lipinski as Lipinksy
import nglview as nv
from pymol import cmd
import os, shutil, time, imageio

1. Клик-химия и ибупрофен

На википедии написано, что одна из наиболее часто используемых реакций в клик-химии это азид-алкиновое циклоприсоединение:

display(Image('https://upload.wikimedia.org/wikipedia/commons/thumb/6/66/CuAAC.svg/'
              '512px-CuAAC.svg.png'))

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

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

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

display(Image('https://lh4.googleusercontent.com/'
              'KF73lQNu7Qb6vJ818z0DSoH7vueyxp1jqZm7I2N1fvg=w617-h164-no',
              format='png'))

А уже к загалогенированному ибупрофену мы сможем нацепить ацетилен по следующей реакции:

display(Image('acet.png'))

Итак, получается, что наш модифицированный (для клик-химии) ибупрофен должен иметь ацетилен у второго атома углерода от кислоты. Тогда, в результате азид-алкильной реакции, этот ацетилен превратится в цикл. Тогда получаем следующую затравку (которой мы потом будем заменять азиды):

template = Chem.MolFromSmiles('N1N=NC(=C1)C(C)(C(=O)O)C1=CC=C(C=C1)CC(C)C')
display(template)

2. Соединения с азидом из PubChem

Для начала посмотрим на некоторые параметры Лепински для нашей затравки:

print(Lipinksy.NumHDonors(template))
print(Lipinksy.NumHAcceptors(template))
print(Lipinksy.rdMolDescriptors.CalcExactMolWt(template))

2
3
273.147726848

Теперь поищем соединения с азидом на пабмеде.
Вкладка Substructure/Superstructure. Вбиваем CID азида: 33558. В Options выбираем Substructure и ставим флажок на Strip hydrogens (просто потому что все остальное нам не нужно).
И, наконец, ставим фильтры на массу и доноров/акцепторов. Исходя из параметров нашей затравки, я решил ставить ограничение по массе <=300 Да, по донорам <=10 и по акцепторам <= 15 (на всякий случай с запасом). Октанольный коэффициент я не стал трогать, потому что не представляю как он высчитывается.

Всего нашлось 66869 соединений. Загрузим скаченные данные и отфильтруем соединения с нековалентными связями:

comps = []
with open('compounds_smiles.txt') as infile:
    for line in infile:
        smi = line.strip().split('\t')[1]
        if '.' not in smi:
            comps.append(smi)

3. Замещаем азид на наш ибупрофен

# Two forms of azides in smiles.
azides = [Chem.MolFromSmiles('N=[N+]=[N-]'),
          Chem.MolFromSmiles('N=[N+]=[NH]')]

newcomps_smi = set()
for comp_smi in comps:
    comp = Chem.MolFromSmiles(comp_smi)
    if comp is None:  # If smiles is bad.
        continue

    # Replace all azides in the compound.
    for azide in azides:
        comp = Chem.rdmolops.ReplaceSubstructs(
            comp, azide, template, replaceAll=True)[0]

    newcomp_smi = Chem.MolToSmiles(comp)
    if '.' not in newcomp_smi:  # If azide was not terminal.
        newcomps_smi.add(newcomp_smi)

Тут все так сложно по четырем причинам:

• в соединениях встречаются две формы записи азида;
• некоторые соединения очень странные (например вот это) и их смайлсы не удается распарсить;
• в одном соединении может быть несколько азидов -- нам нужно заменить их все, поскольку в клик-химии выход должен быть ~100%;
• есть соединения, в которых азид не терминальный, и в них после замещения образуется нековалентная связь.

4. Отфильтруем получившиеся соединения по пяти Лепински

final_comps = []
for comp_smi in newcomps_smi:
    comp = Chem.MolFromSmiles(comp_smi)

    if all((Lipinksy.NumHDonors(comp) <= 5,
            Lipinksy.NumHAcceptors(comp) <= 10,
            Lipinksy.rdMolDescriptors.CalcExactMolWt(comp) < 500,
            Lipinksy.rdMolDescriptors.CalcCrippenDescriptors(comp)[0] <= 5)):
        final_comps.append(comp)

Всего их получилось 42555. А теперь посмотрим на первые 25 из них. У меня не получилось нормально вставить SVG, поэтому тут будет PNG, но с хорошим разрешением.

display(Draw.MolsToGridImage(final_comps, 25,
                             molsPerRow=5, subImgSize=(600, 600)))

5. Строим Similarity Map ибупрофена с пятым соединением

fig, maxweight = SimilarityMaps.GetSimilarityMapForFingerprint(ibu, final_comps[4],
                                                               SimilarityMaps.GetMorganFingerprint)
maxweight

0.15536135386322653

6. 3D

m3d = Chem.AddHs(final_comps[4])
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d, maxIters=500, nonBondedThresh=200)

display(m3d)

import ipywidgets
nv.show_rdkit(m3d)

А вот это не сработало 😞. Посмотрим на нее в пимоле:

Chem.rdmolfiles.MolToPDBFile(m3d, 'xxx.pdb')

import pymol
pymol.finish_launching(['pymol'])

cmd.reinitialize()
cmd.load('xxx.pdb')
cmd.mset('x120')
cmd.mview('store', power=1)
for x in range(3):
    cmd.turn('y', 90)
    cmd.mview('store', 30 * (x + 1), power=1)

cmd.mplay()

cmd.mstop()

dirpath = 'images/'
if os.path.exists(dirpath):
    shutil.rmtree(dirpath)
os.mkdir(dirpath)

cmd.set('ray_trace_mode', 3)
cmd.set('ray_trace_gain', 0)
cmd.set('antialias', 2)
cmd.mpng(dirpath, mode=2, width=1920//3, height=1080//3)

while not os.path.exists(f'{dirpath}0120.png'):
    time.sleep(0.1)

with imageio.get_writer('movie.gif', mode='I', fps=30) as writer:
    for x in range(0, 120):
        filename = os.path.join(dirpath, f'{x+1:04}.png')
        image = imageio.imread(filename)
        writer.append_data(image)

display(Image('movie.gif', format='png'))