Что мы делаем здесь:
- Поднимаем среду.
- Choice A: Подключаемся к shadbox'у, где всё уже собрано, поднимаем там jupyter-notebook и работаем там.
- Плюс: работает, даже если у нас нет своей конды с пакетами и jupyter'ом.
- Минус: если оные есть, то на сервере не будет привычных кастомных свистелок.
- При отладке обнаруживаем, что 127.0.0.1 ≠ localhost. Ну ок.
- Choice B: Разворачиваем на локальной машине правильный набор пакетов.
- Плюсы: нет разрывов, знакомая среда.
- Минусы: таки нужно собрать коллекцию нужных пакетов.
- Опять же, научиться ходить на удалённый jupyter ssh'ем — полезное упражнение. $$\require{mhchem}$$
- Choice A: Подключаемся к shadbox'у, где всё уже собрано, поднимаем там jupyter-notebook и работаем там.
- Вытаскиваем из PubChem'а (а именно, из Structure Search) все соединения с азидной группой: $\ce{-N=N+=N-},$ или же $\ce{-N^{-}-N+#N}.$
- Каждое такое соединение присоединяем к ибупрофену по одной из схем click chemistry, а именно по схеме CuAAC: $\ce{R1-N=N+=N- + R2-C#C ->[\ce{Cu(I)}] R1-N^1=N-N-C=C^1-R2}$
- (верхние индексы сиволизируют циклическую связь)
- Каждую такую склейку проверяем на соответствие правилу пяти Липински.
- Молекулярные свойства (молярную массу и пр.) вычисляем с помощью RDKit.
- Чего не пишут в подсказках к заданию (но пишут в лекции): к PubChem'у можно обращаться прямо из питона. Пакетом pubchempy.
- pubchempy нет в скрипте подготовки среды, но ничего не стоит его туда добавить.
In [3]:
import pubchempy as pcp
Вытаскиваем из PubChem'а все соединения с азидной группой.
- Стоит иметь в виду, какая нам нужна информация.
- В PubChem'е есть много всего: и молекулярная масса, и числа доноров-акцепторов водородных связей, и так далее.
- Но нам нужна информация не об оригинальных, а об изменённых соединениях! А от оригинальных нам достаточно SMILES-представления.
- А всё необходимое про изменённые соединения (которых в PubChem'е может и не быть) мы посчитаем сами, с помощью rdkit'а.
- Откуда мы знаем, что есть rdkit, который такое умеет? Из подсказок, вестимо.
С какими проблемами мы столкнулись и как их решали.
- Во-первых, скорость загрузки. Добрая минута на 200 соединений. Таким темпом 150000 найденных азидов мы будем качать целый месяц.
- Приходит в голову скачивать только SMILES, а не полную информации о веществе. Но как? В документации get_compounds таких тонкостей не предусмотрено.
- Случайным блужданием по мануалам, документации и гитхабу натыкаемся на get_properties, который делает именно то, что нам нужно. Ура!
- Это не решает проблему медленности.
- Возможно, дело в том, что у каждого запроса гигансткий оверхед, и стоит выставить listkey_count побольше, чтобы «отображать» побольше результатов на одной «странице»? Скажем, 1000. Или вовсе без разделения по страницам. Что нам, полторы сотни тысяч строк по сети перебросить сложно?
- Оказывается, не так уж просто: запросы с $listkey\_count\ge10^6$ возвращают 504 Gateway Timeout. По счастью, $10^5$ работает, и три-пять страниц результатов выгружаются достаточно быстро.
- Заодно обнаруживаем, что синтаксис в подсказке из лекции обманчив: listkey_start обозначает не номер страницы, с которой начать, а номер записи, с которой начать. И записи нумеруются с нуля.
- Это поведение описано в документации PubChem Power User Gateway REST API, на которую можно выйти из документации pubchempy про расширенный поиск.
- Во-вторых, PubChem по умолчанию игнорирует указанные заряды.
- Да и с кратностью связей обращается довольно-таки вольно.
- И это важно, поскольку нас интересуют конструкции наподобие $\ce{-N=N+=N^-}$ и $\ce{-N^{-}-N^+#N}$, и не интересуют $\ce{=N-N+#N}$, $\ce{-N^1-N+#N}$, $\ce{-N(N+#N)\ce-}$,</nobr> атомы азота, находящиеся в циклах, и так далее.
- В конце концов решаем проблему так:
- Делаем несколько поисковых запросов.
- В первом ищем соединения по содержанию $\ce{N=N=N}$.
- Во втором ищем соединения по содержанию $\ce{N-N#N}$.
- При поиске устанавливаем параметр RingsNotEmbedded=True, что исключает из поиска не интересующие нас атомы азота в кольцах.
- Определяем три группы подстрок:
- соответствующие настоящим азидным группам (
azides), - соответствующие не-азидным группам, похожим на азидные (
not_azides), - и соответствующие группам, с которыми мы не хотим работать ни при каких условиях (
nope). - К запрещённым относятся:
., точка, обозначающая в SMILES не-ковалентную связь,- $\ce{N=N+=N}$, в которой непонятно, где конец,
- $\ce{=N-N+#N}$, которую очень легко перепутать с $\ce{-N-N+#N}$.</nobr>
- (Непонятно и очень легко — т.е. с точки зрения работы со строками.)
- (Мораль: не надо работать с химическими соединениями как со строками, хуже будет.)
- соответствующие настоящим азидным группам (
- Проверяем себя: смотрим на найденные соединения, в которых нет ни одной подстроки из вышеперечисленных групп. Если таких почти нет, то мы не упустили ни одной существенно представленной формы записи азидной группы в SMILES.
- Фильтруем результаты поиска: оставляем только соединения, имеющие хотя бы одну истинную азидную группу, и не имеющие не одной запрещённой группы.
- Делаем несколько поисковых запросов.
Здесь происходит выгрузка результатов поиска:
In [4]:
compounds = []
per_page = 10**5
for smiles in ["N=N=N", "NN#N"]:
for i in range(20):
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 KeyboardInterrupt:
print("Interrupted")
raise
except:
print("Page {} returns 500, which is "
"our signal to stop.".format(i+1))
break
print("Retrieved page {} of {} search".format(i+1, smiles))
compounds.extend(a)
# Remove duplicates
compounds = {c["CID"]:c for c in compounds}.values()
print("Retrieved {} SMILES of (supposedly) "
"azide-containing compounds".format(len(compounds)))
Здесь происходит проверка фильтра на адекватность:
In [5]:
azides = ["N=[N+]=[N-]", "N[NH+]=[N-]",
"[N-]=[N+]=N", "[N-]=[NH+]N",
"[N-][N+]#N", "N[N+]#N", "[N][N+]#N",
"N#[N+][N-]", "N#[N+]N", "N#[N+][N]",
]
not_azides = ["=N[N+]#N", "=N[N+]=[N-]",
"([N+]#N)", ")[N+]#N", "(=[N-])", ")[NH+]=[N-]", "NN(N"
] + ["N{}[N+]#N".format(i) for i in range(1, 4+1)] \
+ ["N{}NN{}".format(i, i) for i in range(1, 4+1)]
nope = [".", "=N[N+]#N", "N=[N+]=N"]
In [6]:
nonmatched = sum(1 for compound in compounds if
all(substr not in compound["CanonicalSMILES"]
for substr in azides+not_azides+nope))
print(nonmatched)
for compound in compounds:
if all((substr not in compound["CanonicalSMILES"]
for substr in azides+not_azides+nope)):
# print(compound)
pass
Just for fun: здесь можно посмотреть, какие представления азидов и не-азидов наиболее популярны в результатах нашего поиска:
In [ ]:
azidelike_counts = {y:0 for y in azides+not_azides+nope}
for x in compounds:
for y in azides+not_azides+nope:
if y in x["CanonicalSMILES"]:
azidelike_counts[y] += 1
print(sum(azidelike_counts.itervalues()) - len(compounds))
for key, value in azidelike_counts.iteritems():
print(key, value)
Здесь, наконец, происходит фильтрация.
In [7]:
compounds = [compound for compound in compounds if
any(substr in compound["CanonicalSMILES"]
for substr in azides) and
all(substr not in compound["CanonicalSMILES"]
for substr in nope)]
print(len(compounds))
Теперь мы можем добавлять ибупрофен.
Для начала подключим rdkit.
In [8]:
# Those imports are copy-pasted;
# I have not verified that they're all necessary.
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 Lipinski
from IPython.display import display,Image
Извлечём из того же PubChem'а SMILES ибупрофена и для наглядности нарисуем его.
In [9]:
ibuprofen = Chem.MolFromSmiles("CC(C)CC1=CC=C(C=C1)C(C)C(=O)O")
display(ibuprofen)
Для CuAAC нам нужна ацетиленовая группа. Объявляем, что заменить в ибупрофене изопропил на неё легко и просто, и это и будет делаться в живых реакциях.
In [10]:
ibuprofen_2 = Chem.MolFromSmiles("C#CCC1=CC=C(C=C1)C(C)C(=O)O")
display(ibuprofen_2)
- Проходимся по отобранным веществам, проводим над каждым реакцию.
- Пффф!
- Не верьте!
- На самом деле мы заменяем в SMILES фрагменты, которые наверное обозначают азидную группу, на фрагменты, обозначающие присоединённый ибупрофен.
- В большей части случаев мы действительно получим SMILES вещества, в котором к азидной группе присоединён ибупрофен, но не всегда.
- По-хорошему, здесь нужно найти все азидные группы (их может быть много больше одной), и перебрать все варианты присоединения ибупрофенов к ним.
- Но это сложно.
- Тем сложнее, если мы пытаемся работать с веществами как со строками.
- Поэтому мы заменяем все азидные группы, которые нашлись.
- Вообще-то следует иметь в виду, записан ли азид в направлении от центра молекулы к концу или в обратном: в первом случае связь у него будет «слева», а во втором — «справа».
- Но можно подсединнный ибупрофен записать в SMILES так, чтобы он одинаково хорошо присоединялся как к связи слева, так и к связи справа.
- Иногда SMILES, который мы получили, не является корректным и RDKit отказывается строить по нему молекулу. Например, в нескольких десятках молекул хлор имеет неканоничную валентность 3.
- Эти молекулы мы пропускаем.
In [11]:
ibu = "C(=C9)(N=NN9CC8=CC=C(C=C8)C(C)C(=O)O)"
compounds2 = {}
for i, compound in enumerate(compounds):
s = compound["CanonicalSMILES"]
for azide in azides:
s = s.replace(azide, ibu)
mol = Chem.MolFromSmiles(s)
if i%10000 == 0:
print(i)
#display(mol)
if mol is None:
continue
if (Lipinski.NumHDonors(mol) <= 5 and
Lipinski.NumHAcceptors(mol) <= 10 and
Lipinski.rdMolDescriptors.CalcExactMolWt(mol) < 500 and
Lipinski.rdMolDescriptors.CalcCrippenDescriptors(mol)[0] <= 5):
compounds2[compound['CID']] = mol
print("Done: {} molecules Lipinski-approved".format(len(compounds2)))
- Собственно, всё?
- Можно нарисовать какие-нибудь из полученных молекул.
- Было бы круто визуализировать их в 3д, но nglview не взлетает.
In [12]:
ibuprofen_3 = Chem.MolFromSmiles(ibu)
AllChem.Compute2DCoords(ibuprofen_3)
to_draw = {key:compounds2[key] for key in compounds2.keys()[::5000]}
for mol in to_draw.values():
AllChem.GenerateDepictionMatching2DStructure(mol, ibuprofen_3)
AllChem.Compute2DCoords(mol)
#display(mol)
img = Draw.MolsToGridImage(to_draw.values(),
molsPerRow=2,
subImgSize=(400,200),
#legends=[str(k) for k in to_draw.keys()],
)
display(img)
#for mol in compounds2.values()[::10000]:
# display(mol)
#nv.show_rdkit(mol3d)
- Update: nglview починили.
In [13]:
import nglview as nv
In [17]:
m3d=Chem.AddHs(to_draw.values()[7])
Chem.AllChem.EmbedMolecule(m3d)
AllChem.MMFFOptimizeMolecule(m3d,maxIters=500,nonBondedThresh=200 )
nv.show_rdkit(m3d)
In [ ]: