** Код, использованный в этой статье, доступен здесь на моем GitHub.
Проблема
Процесс разработки и выпуска нового лекарственного препарата - это грандиозный подвиг. Часто первым шагом на пути к открытию лекарств является тестирование сотен тысяч лекарств-кандидатов in vitro с использованием высокопроизводительного скрининга - автоматизированных процессов, которые позволяют проводить простые фармакологические эксперименты параллельно с целевым лекарством. Как и следовало ожидать, этот процесс смехотворно дорог, требует много времени и ресурсов. Поэтому вычислительные подходы, позволяющие проводить скрининг кандидатов в лекарственные препараты in silico, очень востребованы; Такой тип системы может сократить количество необходимых физических экспериментов на порядки, резко ускорив открытие лекарств и снизив затраты. Эти типы подходов уже существуют, и вы можете увидеть мой любимый пример этой концепции в этой статье, опубликованной в Cell ранее в этом году.
В этой статье я собираюсь продемонстрировать, как мы можем создать очень простой скрипт Python и Torch, который принимает молекулу в качестве входных данных (представленную строкой SMILES) и классифицирует ее как нерастворимую, малорастворимую или растворимую в воды через глубокую нейронную сеть, обученную на наборе данных Harvard AqSolDB, который содержит список из 10000 соединений и их различных свойств. Как и следовало ожидать, растворимость - чрезвычайно важная характеристика потенциального лекарства. Лекарство должно быть способным распространяться в целевой (ые) области (ах) тела, и он не может этого сделать, если при пероральном приеме остается в основном в виде твердого вещества. Более того, было бы совершенно катастрофически, если бы соединение выскочило из раствора при внутривенном введении пациенту. Вы можете узнать больше о важности растворимости лекарств здесь.
Как мы вообще подходим к этому с помощью глубокого обучения?
Большинство нейронных сетей обычно принимают входные данные фиксированной размерности или, по крайней мере, входные данные, которые имеют фиксированные ограничения на их размерность. Например, сверточная нейронная сеть, обученная на изображениях размером 256x256 пикселей, не может даже обработать изображение размером 255x257 пикселей исключительно с архитектурной точки зрения (конечно, изображение можно обрезать, изменить размер или преобразовать таким образом, чтобы оно могло использоваться как вход). Другие сети, такие как рекуррентные нейронные сети (RNN), могут принимать входные данные произвольной длины, но требуют, чтобы все измерения, кроме длины, были фиксированными. Например, RNN, обученная на данных на уровне символов, может принимать последовательности размерности (Nx26), где N - длина последовательности, а 26 - вектор горячих признаков для каждой буквы алфавита. По сути, гибкость входных данных нейронной сети варьируется от модели к модели, но почти всегда существуют ограничения на входное измерение.
С другой стороны, молекулы сильно различаются по структуре. Даже молекулы с одинаковым числом и типом атомов, такие как октан (C8H20), можно расположить 18 различными способами. Другими словами, состав, длина и связность могут варьироваться между молекулами. Какую нейронную архитектуру мы можем использовать для обработки такого рода диких данных?
Во-первых, это помогает думать о молекулах как о неориентированных графах. Все молекулы могут быть представлены набором узлов и ребер. Тогда каждый узел в таком графе будет иметь набор функций, таких как элемент, который представляет узел, стереохимическая информация об этом атоме, состояние гибридизации, местоположение в физическом пространстве и т. Д. Итак, для начала давайте напишем некоторый код Python, который примет строку SMILES и преобразует ее в структуру графа NetworkX с помощью Pysmiles. Это позволит нам, по крайней мере, обрабатывать данные значимым образом.
import numpy as np
from pysmiles import read_smiles
G = read_smiles("CN(C)C(=N)N=C(N)N", explicit_hydrogen=True)
print(G.nodes(data=’element’))
print(G.edges)
Это так просто - теперь у нас есть граф NetworkX, G, который содержит структуру нашей молекулы (в данном случае метформина). Результат этого фрагмента:
[(0, ‘C’), (1, ‘N’), (2, ‘C’), (3, ‘C’), (4, ‘N’), (5, ‘N’), (6, ‘C’), (7, ‘N’), (8, ‘N’), (9, ‘H’), (10, ‘H’), (11, ‘H’), (12, ‘H’), (13, ‘H’), (14, ‘H’), (15, ‘H’), (16, ‘H’), (17, ‘H’), (18, ‘H’), (19, ‘H’)] [(0, 1), (0, 9), (0, 10), (0, 11), (1, 2), (1, 3), (2, 12), (2, 13), (2, 14), (3, 4), (3, 5), (4, 15), (5, 6), (6, 7), (6, 8), (7, 16), (7, 17), (8, 18), (8, 19)]
Обратите внимание, как каждый узел от 0 до 19 имеет соответствующий элемент и как каждое ребро описывается упорядоченной парой (родительский, дочерний).
Веселая часть
Пришло время заняться глубоким обучением. Оказывается, существует множество исследований, посвященных обработке графов с помощью глубоких нейронных сетей, поскольку эта проблема имеет широкое применение. В этом проекте мы будем использовать сверточную сеть графов (GCN). GCN очень элегантно обрабатывают произвольную и разнообразную структуру графов. По сути, каждый узел на нашем графике имеет собственный вектор одномерных признаков фиксированной размерности. На каждом этапе обработки в нейронной сети узлы обмениваются информацией друг с другом, отправляя преобразованную версию своего вектора признаков своим соседним узлам. Он также обрабатывает свою собственную информацию, отправляя преобразованный вектор признаков самому себе. Преобразование, применяемое к вектору признаков, может быть представлено как обучаемая матрица весов, которую мы можем обучать с помощью градиентного спуска. Каждый узел берет все преобразованные векторы признаков от своих соседей и объединяет их с помощью операции сохранения размерности, такой как поэлементное среднее, максимальный пул или что-то еще. В нашем частном случае операция сохранения размерности немного сложна, но очень красива. Предлагаю прочитать эту отличную статью о GCN. Наконец, применяется функция активации нелинейности, такая как RELU, Tanh, Sigmoid или что-то еще. Наконец, мы можем использовать эти обработанные функции на каждом узле, чтобы вычислить что-то о свойствах каждого узла или, возможно, графа в целом. Например, если мы пытаемся предсказать возраст пользователя в социальной сети, мы могли бы использовать регрессию для функций в каждом узле, чтобы сделать прогноз. В нашем случае мы хотим провести классификацию графов, т.е. мы хотим классифицировать граф в целом как растворимый, нерастворимый или слабо растворимый. Для этого мы можем выполнить одну последнюю агрегацию с сохранением размерности всех узлов в графе и запустить классификационную сеть на этом последнем объекте.
Код
Мы будем использовать Python и PyTorch Geometric для обучения нашей модели. Pytorch Geometric построен на основе PyTorch и реализует множество графических сетевых слоев и алгоритмов, которые значительно упростят нам этот процесс. Во-первых, давайте уберем все наши операторы импорта, чтобы они не загромождали остальную часть нашего кода:
import numpy as np import random import matplotlib.pyplot as plt from pysmiles import read_smiles import pandas as pd import logging from tqdm import tqdm import torch from torch.nn import Sequential as Seq, Linear, ReLU, CrossEntropyLoss import torch.nn.functional as F from torch_geometric.nn import MessagePassing, GCNConv from torch_geometric.utils import remove_self_loops, add_self_loops, degree from torch_geometric.data import Data logging.getLogger(‘pysmiles’).setLevel(logging.CRITICAL) # Anything higher than warning
Загрузка данных
Во-первых, давайте загрузим наш набор данных с помощью Pandas и отформатируем его в структуру данных, которую использует PyTorch Geometric.
df = pd.read_csv(‘dataset.csv’) #read dataset X_smiles = list(df[‘SMILES’]) #get smiles strings Y = np.asarray(df[‘Solubility’]) #get solubility values #list of all elements in the dataset, which I've precomputed elements = [‘K’, ‘Y’, ‘V’, ‘Sm’, ‘Dy’, ‘In’, ‘Lu’, ‘Hg’, ‘Co’, ‘Mg’, ‘Cu’, ‘Rh’, ‘Hf’, ‘O’, ‘As’, ‘Ge’, ‘Au’, ‘Mo’, ‘Br’, ‘Ce’, ‘Zr’, ‘Ag’, ‘Ba’, ’N’, ‘Cr’, ‘Sr’, ‘Fe’, ‘Gd’, ‘I’, ‘Al’, ‘B’, ‘Se’, ‘Pr’, ‘Te’, ‘Cd’, ‘Pd’, ‘Si’, ‘Zn’, ‘Pb’, ‘Sn’, ‘Cl’, ‘Mn’, ‘Cs’, ‘Na’, ‘S’, ‘Ti’, ‘Ni’, ‘Ru’, ‘Ca’, ‘Nd’, ‘W’, ‘H’, ‘Li’, ‘Sb’, ‘Bi’, ‘La’, ‘Pt’, ‘Nb’, ‘P’, ‘F’, ‘C’] #convert element to a one-hot vector of dimension len(elements) def element_to_onehot(element): out = [] for i in range(0, len(element)): v = np.zeros(len(elements)) v[elements.index(element[i])] = 1.0 out.append(v) return np.asarray(out)
Выше мы извлекли строки SMILES и информацию о растворимости из файла .csv, определили список всех элементов в строках и определили функцию с именем element_to_onehot. Эта функция преобразует строку элементов («C», «Ru» и т. Д.) В 61-мерный одномерный вектор, который будет стартовыми характеристиками наших узлов.
def val_to_class(val): if val < -3.65: #insoluble return [1, 0, 0] elif val < -1.69: #slightly soluble return [0, 1, 0] else: #soluble return [0, 0, 1]
Далее мы определяем критерии растворимости в этой функции. По сути, если насыщенный раствор может содержать не более 10 ^ -3,65 моль / л, мы считаем это растворенное вещество нерастворимым. Если он не растворим, но может удерживать в растворе не более 10-1,69 моль / л, мы считаем его малорастворимым. В противном случае результат очень растворим. Конечно, это не лучшее определение растворимости, и я просто сделал это так, чтобы набор данных разделился на четные трети. Однако это подойдет для нашего примера.
#process SMILES strings into graphs nodes = [] edge_index = [] for smiles in tqdm(X_smiles): try: G = read_smiles(smiles, explicit_hydrogen=True) feature = element_to_onehot(np.asarray(G.nodes(data=’element’))[:, 1]) edges = np.asarray(G.edges) index = np.asarray([edges[:,0], edges[:,1]]) nodes.append(feature) edge_index.append(index) except: pass
Здесь мы преобразуем строки SMILES в графики, как делали раньше, и сохраняем их в списках. У нас есть блок try / except, чтобы игнорировать любые ошибки при обработке строки SMILES. TQDM используется как индикатор выполнения.
#Generate Data objects data = list() #process graphs into torch_geometric Data objects for i in tqdm(range(0, len(nodes))): x = torch.tensor(nodes[i], dtype=torch.float) edges = torch.tensor(edge_index[i], dtype=torch.long) y = torch.tensor([val_to_class(Y[i])], dtype=torch.float) data.append(Data(x=x,edge_index=edges, y=y)) random.shuffle(data) train = data[:int(len(data)*0.8)] #train set test = data[int(len(data)*0.8):] #val set train = data
Затем мы генерируем объекты PyTorch Geometric Data (). Для этого просто требуется список характеристик узла (x) и индекс ребра (ребра). Мы также добавляем нашу метку данных (y) к объекту в целях обучения. Наконец, мы разделяем наш набор данных на обучающие и проверочные наборы с разделением 80:20.
Определение сети
Наконец-то мы готовы определить нашу сверточную сеть на основе графов. Мы настроили его так же, как и любую другую сеть PyTorch, поэтому я вставлю блок кода ниже, и мы пройдемся по частям кода, которые являются уникальными и важными.
class Net(torch.nn.Module): def __init__(self): super(Net, self).__init__() self.conv1 = GCNConv(61, 32) self.conv2 = GCNConv(32, 32) self.conv3 = GCNConv(32, 32) self.conv4 = GCNConv(32, 32) self.lin1 = Linear(32, 16) self.lin2 = Linear(16, 3) def forward(self, data): x, edge_index= data.x, data.edge_index x = self.conv1(x, edge_index) x = F.relu(x) x = F.dropout(x, p=0.25, training=self.training) x = self.conv2(x, edge_index) x = F.relu(x) x = F.dropout(x, p=0.25, training=self.training) x = self.conv3(x, edge_index) x = F.relu(x) x = F.dropout(x, p=0.25, training=self.training) x = self.conv4(x, edge_index) x = F.relu(x) x = torch.sum(x, dim=0) x = self.lin1(x) x = F.relu(x) x = self.lin2(x) return x
Во-первых, мы используем слои GCNConv, которые принимают входное и выходное измерения. Обратите внимание, как наш первый слой принимает входной размер 61, что является нашим горячим измерением, и выводит размер 32. После этого слоя у нас будет набор узлов, каждый из которых имеет вектор признаков размером 32. Перед нашим первым На линейном слое я суммирую значения всех характеристик узла поэлементно, чтобы получить один 32-мерный вектор, по которому мы будем выполнять классификацию через обычную полносвязную нейронную сеть. Обратите внимание, что мы не применяем SoftMax или какую-либо нормализацию к окончательному результату размера 3 (представляющему наши 3 класса) - это потому, что функция потерь, которую мы используем позже, сделает это автоматически.
Обучение
Пришло время приступить к обучению модели. Давайте сначала инициализируем нашу среду:
#set up device and create model device = torch.device(‘cuda’ if torch.cuda.is_available() else ‘cpu’) #use CUDA if available model = Net().to(device) #create network and send to the device memory optimizer = torch.optim.Adam(model.parameters(), lr=0.001, weight_decay=1e-4) #use Adam optimizer CSE = CrossEntropyLoss()
Здесь мы определяем устройство, которое используем, объявляем нашу модель и переносим ее в память устройства, объявляем наш оптимизатор (Адам) и объявляем нашу функцию потерь (кросс-энтропия).
Теперь мы определим наш основной цикл обучения:
#train model
model.train() #set model to training mode
for epoch in range(32): #run for epochs of training
sum_loss = 0 #used to compute average loss in an epoch
num_correct = 0
random.shuffle(train) #shuffle the training data each epoch
for d in tqdm(train): #go over each training point
data = d.to(device) #send data to device
optimizer.zero_grad() #zero gradients
out = model(data) #evaluate data point
if torch.argmax(out) == torch.argmax(data.y): #if prediction is correct, increment counter for accuracy calculation
num_correct += 1
loss = CSE(torch.reshape(out, [1, 3]), torch.reshape(torch.argmax(data.y),[1])) #compute loss
sum_loss += float(loss) #add loss value to aggregate loss
loss.backward() #compute gradients
optimizer.step() #apply optimization
print(‘Epoch: {:03d}, Average loss: {:.5f}, Accuracy: {:.5f}’.format(epoch, sum_loss/len(train), num_correct/len(train)))
На самом деле это ничем не отличается от обычного цикла обучения PyTorch, поэтому не нужно много говорить об этом фрагменте кода. После 32 эпох у нас есть модель, точность которой составляет около 43% по нашему набору проверки. По сравнению со случайным угадыванием (точность 33%) мы видим повышение точности примерно на 10%. Конечно, это не очень хорошо, но растворимость - невероятно сложная проблема, и удивительно, что мы видим какое-либо повышение производительности при такой упрощенной настройке.
Тестируем это
Давайте немного поиграем с моделью, чтобы увидеть, чему она научилась.
Во-первых, давайте скармливаем ему строку улыбки для декана (C10H22), углеводорода, который в основном нерастворим в воде.
>>>evaluate_smiles(‘CCCCCCCCCC’) insoluble
Здорово! Теперь давайте добавим функциональные группы карбоновых кислот к концам нашей молекулы и посмотрим, как это влияет на растворимость.
>>>evaluate_smiles(‘C(=O)(O)CCCCCCCCC(=O)(O)’) slightly soluble
Действительно, молекула мало растворима в воде. Замечательно, похоже, модель работает на этот случай!
Хорошо, давай попробуем что-нибудь посложнее. Давайте дадим ему структуру УЛЫБКИ аспирина:
>>>evaluate_smiles(‘O=C(C1=C(OC(C)=O)C=CC=C1)O’) soluble
Я очень на это надеюсь! Аспирин лучше растворим, иначе я последние несколько лет принимал плацебо от головных болей ...
Однако модель не идеальна… она сообщает, что бензол «слабо растворим» в воде, что, конечно же, неверно:
>>>evaluate_smiles(‘c1ccccc1’) slightly soluble
Когда мы добавляем к нашему бензолу пару трет-бутильных групп, модель сообщает, что он полностью нерастворим, и это здорово!
>>>evaluate_smiles('CC(C)(C)c1ccc(C(C)(C)C)cc1')
insoluble
Хорошо, для последнего теста давайте попробуем протеин-инсулин и посмотрим, как он с ним справится ...
>>>evaluate_smiles(‘CCC(C)C1C(=O)NC2CSSCC(C(=O)NC(CSSCC(C(=O)NCC(=O)NC(C(=O)NC(C(=O)NC(C(=O)NC(C(=O)NC(C(=O)NC(C(=O)NC(C(=O)NC(C(=O)NC(C(=O)NC(C(=O)NC(CSSCC(NC(=O)C(NC(=O)C(NC(=O)C(NC(=O)C(NC(=O)C(NC(=O)C(NC(=O)C(NC(=O)C(NC2=O)CO)CC(C)C)CC3=CC=C(C=C3)O)CCC(=O)N)CC(C)C)CCC(=O)O)CC(=O)N)CC4=CC=C(C=C4)O)C(=O)NC(CC(=O)N)C(=O)O)C(=O)NCC(=O)NC(CCC(=O)O)C(=O)NC(CCCNC(=N)N)C(=O)NCC(=O)NC(CC5=CC=CC=C5)C(=O)NC(CC6=CC=CC=C6)C(=O)NC(CC7=CC=C(C=C7)O)C(=O)NC(C(C)O)C(=O)N8CCCC8C(=O)NC(CCCCN)C(=O)NC(C)C(=O)O)C(C)C)CC(C)C)CC9=CC=C(C=C9)O)CC(C)C)C)CCC(=O)O)C(C)C)CC(C)C)CC2=CN=CN2)CO)NC(=O)C(CC(C)C)NC(=O)C(CC2=CN=CN2)NC(=O)C(CCC(=O)N)NC(=O)C(CC(=O)N)NC(=O)C(C(C)C)NC(=O)C(CC2=CC=CC=C2)N)C(=O)NC(C(=O)NC(C(=O)N1)CO)C(C)O)NC(=O)C(CCC(=O)N)NC(=O)C(CCC(=O)O)NC(=O)C(C(C)C)NC(=O)C(C(C)CC)NC(=O)CN’) insoluble
Хм, это же неожиданно, не правда ли? Мы можем вводить инсулин прямо в наши вены. На самом деле инсулин нерастворим в чистой воде, но легко растворяется в плазме крови. Это отличный пример ограничений такой модели. Соединения, считающиеся растворимыми / нерастворимыми в воде, могут вести себя иначе в физиологической среде! Кроме того, такие белки, как инсулин, образуют особую трехмерную структуру, которую эта модель явно не принимает во внимание. Таким образом, для любых больших молекул эта модель, вероятно, будет очень неточной.
Наконец, давайте попробуем суперрастворимое органическое соединение, витамин С:
>>>evaluate_smiles(‘OC([C@H](O1)[C@H](CO)O)=C(O)C1=O’) soluble
Фух! Было бы неловко, если бы он ошибся.
Заворачивать
Используя простую сверточную сеть графов в PyTorch Geometric, мы смогли создать систему классификации, которая могла бы взять произвольное графическое представление молекулы и предсказать ее растворимость с точностью на 10% лучше, чем случайное предположение. Не очень хорошо, но потрясающе, учитывая сложность задачи! Мы также немного повеселились, играя с разными молекулами в нашем программном обеспечении.
