Этот пример показывает вторичный метод предсказания структуры, который использует feedforward нейронную сеть и функциональность, доступную с Deep Learning Toolbox™.
Это упрощенный пример, предназначенный для иллюстрации этапов создания нейронной сети с целью предсказания вторичной структуры белков. Его строение и методы обучения не предназначены для того, чтобы обязательно быть лучшим решением для рассматриваемой задачи.
Модели нейронных сетей пытаются моделировать обработку информации, которая происходит в мозге и широко используется в различных приложениях, включая автоматическое распознавание шаблонов.
Набор данных Рост-Сандера [1] состоит из белков, структуры которых охватывают относительно широкую область значений типов области, состава и длины. Файл RostSanderDataset.mat содержит подмножество этого набора данных, где структурное назначение каждого остатка сообщается для каждой белковой последовательности.
load RostSanderDataset.mat N = numel(allSeq); id = allSeq(7).Header % annotation of a given protein sequence seq = int2aa(allSeq(7).Sequence) % protein sequence str = allSeq(7).Structure % structural assignment
id =
'1CSE-ICOMPLEX(SERINEPROTEINASE-INHIBITOR)03-JU'
seq =
'KSFPEVVGKTVDQAREYFTLHYPQYNVYFLPEGSPVTLDLRYNRVRVFYNPGTNVVNHVPHVG'
str =
'CCCHHHCCCCHHHHHHHHHHHCCCCEEEEEECCCCEECCCCCCEEEEEEECCCCEECCCCEEC'
В этом примере вы создадите нейронную сеть, чтобы узнать структурное состояние (спираль, лист или катушка) каждого остатка в данном белке, основываясь на структурных шаблонах, наблюдаемых во время фазы обучения. Из-за случайной природы некоторых шагов в следующем подходе числовые результаты могут быть немного различными каждый раз, когда сеть обучена или моделируется предсказание. Чтобы гарантировать воспроизводимость результатов, мы сбрасываем глобальный случайный генератор в сохраненное состояние, включенное в загруженный файл, как показано ниже:
rng(savedState);
Для текущей задачи зададим нейронную сеть с одним входным слоем, одним скрытым слоем и одним выходным слоем. Входной слой кодирует скользящее окно в каждой входной аминокислотной последовательности, и делается предсказание структурного состояния центрального остатка в окне. Мы выбираем окно размера 17 на основе статистической корреляции, обнаруженной между вторичной структурой заданного положения остатка и восемью остатками с каждой стороны точки предсказания [2]. Каждое положение окна кодируется с использованием двоичного массива размера 20, имеющей по одному элементу для каждого типа аминокислот. В каждой группе из 20 входов элемент, относящийся к аминокислотному типу в заданном положении, устанавливается равным 1, в то время как все другие входы установлены равным 0. Таким образом, входной слой состоит из R = 17x20 входных модули, то есть 17 групп по 20 входов каждая.
В следующем коде мы сначала определяем для каждой последовательности белков все возможные подпоследовательности, соответствующие скользящему окну размера W, путем создания матрицы Ханкеля, где i-й столбец представляет подпоследовательность, начиная с i-го положения в исходной последовательности. Затем для каждой позиции в окне мы создадим массив размера 20 и установите j-й элемент равным 1, если остаток в заданной позиции имеет числовое представление, равное j.
W = 17; % sliding window size % === binarization of the inputs for i = 1:N seq = double(allSeq(i).Sequence); % current sequence win = hankel(seq(1:W),seq(W:end)); % all possible sliding windows myP = zeros(20*W,size(win,2)); % input matrix for current sequence for k = 1:size(win, 2) index = 20*(0:W-1)' + win(:,k); % input array for each position k myP(index,k) = 1; end allSeq(i).P = myP; end
Слой выхода нашей нейронной сети состоит из трёх модулей, по одному для каждого из рассматриваемых структурных состояний (или классов), которые кодируются с помощью двоичной схемы. Чтобы создать целевую матрицу для нейронной сети, мы сначала получаем из данных структурные назначения всех возможных подпоследовательностей, соответствующих скользящему окну. Затем рассмотрим центральное положение в каждом окне и преобразуем соответствующее структурное назначение с помощью следующей двоичной кодировки: 1 0 0 для катушки, 0 1 0 для листа, 0 0 1 для спирали.
cr = ceil(W/2); % central residue position % === binarization of the targets for i = 1:N str = double(allSeq(i).Structure); % current structural assignment win = hankel(str(1:W),str(W:end)); % all possible sliding windows myT = false(3,size(win,2)); myT(1,:) = win(cr,:) == double('C'); myT(2,:) = win(cr,:) == double('E'); myT(3,:) = win(cr,:) == double('H'); allSeq(i).T = myT; end
Можно выполнить бинаризацию входной и целевой матрицы, описанной в двух шагах выше, более кратким способом, выполнив следующий эквивалентный код:
% === concise binarization of the inputs and targets for i = 1:N seq = double(allSeq(i).Sequence); win = hankel(seq(1:W),seq(W:end)); % concurrent inputs (sliding windows) % === binarization of the input matrix allSeq(i).P = kron(win,ones(20,1)) == kron(ones(size(win)),(1:20)'); % === binarization of the target matrix allSeq(i).T = allSeq(i).Structure(repmat((W+1)/2:end-(W-1)/2,3,1)) == ... repmat(('CEH')',1,length(allSeq(i).Structure)-W+1); end
Когда мы задаем вход и целевые матрицы для каждой последовательности, мы создадим вход матрицу P, и целевая матрица, T, представляющего кодирование для всех последовательностей, подаваемых в сеть.
% === construct input and target matrices P = double([allSeq.P]); % input matrix T = double([allSeq.T]); % target matrix
Проблема предсказания вторичной структуры может быть рассмотрена как задача распознавания шаблонов, где сеть обучена распознавать структурное состояние центрального остатка, наиболее вероятно, когда наблюдаются определенные остатки в данном скользящем окне. Создадим нейронную сеть распознавания шаблона с помощью входа и целевой матриц, заданных выше и задающих скрытый слой размера 3.
hsize = 3;
net = patternnet(hsize);
net.layers{1} % hidden layer
net.layers{2} % output layer
ans =
Neural Network Layer
name: 'Hidden'
dimensions: 3
distanceFcn: (none)
distanceParam: (none)
distances: []
initFcn: 'initnw'
netInputFcn: 'netsum'
netInputParam: (none)
positions: []
range: [3x2 double]
size: 3
topologyFcn: (none)
transferFcn: 'tansig'
transferParam: (none)
userdata: (your custom info)
ans =
Neural Network Layer
name: 'Output'
dimensions: 0
distanceFcn: (none)
distanceParam: (none)
distances: []
initFcn: 'initnw'
netInputFcn: 'netsum'
netInputParam: (none)
positions: []
range: []
size: 0
topologyFcn: (none)
transferFcn: 'softmax'
transferParam: (none)
userdata: (your custom info)
Сеть распознавания шаблонов использует алгоритм Scaled Conjugate Gradient по умолчанию для обучения, но доступны другие алгоритмы (список доступных функций см. в документации Deep Learning Toolbox). В каждом цикле обучения обучающие последовательности представляются сети через скользящее окно, заданное выше, по одному остатку за раз. Каждый скрытый модуль преобразует сигналы, принятые с входом уровня, используя передаточную функцию logsig чтобы создать сигнал выхода, который находится между и близко к 0 или 1, симулируя запуск нейрона [2]. Веса корректируются так, чтобы ошибка между наблюдаемым выходом от каждого модуля и желаемым выходом, заданным целевой матрицей, была минимизирована.
% === use the log sigmoid as transfer function net.layers{1}.transferFcn = 'logsig'; % === train the network [net,tr] = train(net,P,T);
Во время обучения открывается окно инструмента обучения, в котором отображается прогресс. Показаны детали обучения, такие как алгоритм, критерии эффективности, тип учитываемой ошибки и т.д.
Используйте функцию view чтобы сгенерировать графическое представление нейронной сети.
view(net)
Одной из распространенных проблем, которая возникает во время обучения нейронной сети, является сверхподбор кривой данных, где сеть, как правило, запоминает обучающие примеры, не научившись обобщать новые ситуации. Метод по умолчанию для улучшения обобщения называется ранней остановкой и заключается в делении доступного набора обучающих данных на три подмножества: (i) обучающий набор, который используется для вычисления градиента и обновления весов и смещений сети; (ii) набор валидации, ошибка которого контролируется в процессе обучения, поскольку он имеет тенденцию увеличиваться при избыточной подгонке данных; и iii) тестовый набор, ошибку которого можно использовать для оценки качества деления набора данных.
При использовании функции trainПо умолчанию данные случайным образом делятся так, что 60% выборок присваиваются набору обучающих данных, 20% - набору валидации и 20% - тестовому набору, но другие типы разбиения могут применяться путем определения свойства net.divideFnc (по умолчанию dividerand). Структурный состав остатков в трех подмножествах сопоставим, как видно из следующего исследования:
[i,j] = find(T(:,tr.trainInd)); Ctrain = sum(i == 1)/length(i); Etrain = sum(i == 2)/length(i); Htrain = sum(i == 3)/length(i); [i,j] = find(T(:,tr.valInd)); Cval = sum(i == 1)/length(i); Eval = sum(i == 2)/length(i); Hval = sum(i == 3)/length(i); [i,j] = find(T(:,tr.testInd)); Ctest = sum(i == 1)/length(i); Etest = sum(i == 2)/length(i); Htest = sum(i == 3)/length(i); figure() pie([Ctrain; Etrain; Htrain]); title('Structural assignments in training data set'); legend('C', 'E', 'H') figure() pie([Cval; Eval; Hval]); title('Structural assignments in validation data set'); legend('C', 'E', 'H') figure() pie([Ctest; Etest; Htest]); title('Structural assignments in testing data set '); legend('C', 'E', 'H')
Функция plotperform отображать тренды ошибок обучения, валидации и тестирования по мере прохождения итераций обучения.
figure() plotperform(tr)
Процесс обучения останавливается, когда одно из нескольких условий (см net.trainParam) соблюдается. Для примера в рассматриваемом обучении процесс обучения останавливается, когда ошибка валидации увеличений для заданного количества итераций (6) или достигнуто максимальное количество разрешенных итераций (1000).
% === display training parameters net.trainParam % === plot validation checks and gradient figure() plottrainstate(tr)
ans =
Function Parameters for 'trainscg'
Show Training Window Feedback showWindow: true
Show Command Line Feedback showCommandLine: false
Command Line Frequency show: 25
Maximum Epochs epochs: 1000
Maximum Training Time time: Inf
Performance Goal goal: 0
Minimum Gradient min_grad: 1e-06
Maximum Validation Checks max_fail: 6
Sigma sigma: 5e-05
Lambda lambda: 5e-07
Чтобы проанализировать сетевую реакцию, мы исследуем матрицу неточностей, рассматривая выходы обученной сети и сравнивая их с ожидаемыми результатами (целями).
O = sim(net,P); figure() plotconfusion(T,O);
Диагональные камеры показывают количество положений остатков, которые были правильно классифицированы для каждого структурного класса. Недиагональные камеры показывают количество положений остатков, которые были неправильно классифицированы (например, спиральные положения, предсказанные как спиральные положения). Диагональные камеры соответствуют наблюдениям, которые правильно классифицированы. Как количество наблюдений, так и стойкость общего количества наблюдений показаны в каждой камере. В столбце в крайнем правом углу графика показаны проценты всех предсказанных примеров, принадлежащих каждому классу, которые правильно и неправильно классифицированы. Эти метрики часто называют точностью (или положительным прогностическим значением) и частотой ложных открытий, соответственно. Строка в нижней части графика показывает проценты всех примеров, принадлежащих каждому классу, которые правильно и неправильно классифицированы. Эти метрики часто называют отзывом (или истинной положительной частотой) и ложноотрицательной частотой, соответственно. Камера в правом нижнем углу графика показывает общую точность.
Можно также рассмотреть кривой Приемник рабочей характеристики (ROC), график истинной положительной скорости (чувствительности) от ложноположительной скорости (1 - специфичности).
figure() plotroc(T,O);
Нейронная сеть, которую мы определили, относительно проста. Чтобы добиться некоторых улучшений точности предсказания, мы могли бы попробовать одно из следующего:
Увеличьте количество обучающих векторов. Увеличение количества последовательностей, предназначенных для обучения, требует большей курированной базы данных белковых структур с соответствующим распределением спиральных, спиральных и листовых элементов.
Увеличьте количество значений входа. Увеличение размера окна или добавление более релевантной информации, такой как биохимические свойства аминокислот, являются допустимыми опциями.
Используйте другой алгоритм настройки. Различные алгоритмы различаются требованиями к памяти и скорости. Например, алгоритм Scaled Conjugate Gradient относительно медлен, но эффективен для памяти, в то время как Levenberg-Marquardt быстрее, но более требователен с точки зрения памяти.
Увеличьте количество скрытых нейронов. Добавляя больше скрытые модули, мы обычно получаем более сложную сеть с потенциалом для повышения производительности, но мы должны быть осторожны, чтобы не подгонять данные.
Мы можем задать больше скрытых слоев или увеличить размер скрытого слоя при создании сети распознавания шаблонов, как показано ниже:
hsize = [3 4 2]; net3 = patternnet(hsize); hsize = 20; net20 = patternnet(hsize);
Мы также можем назначить начальные веса сети случайным значениям в области значений от -0,1 до 0,1, как предложено исследованием, сообщенным в [2], путем установки net20.IW и net20.LW свойства следующим образом:
% === assign random values in the range -.1 and .1 to the weights
net20.IW{1} = -.1 + (.1 + .1) .* rand(size(net20.IW{1}));
net20.LW{2} = -.1 + (.1 + .1) .* rand(size(net20.LW{2}));
В целом большие сети (с 20 или более скрытыми модулями) достигают лучшей точности на наборе обучающих данных белка, но худшей точности в точности предсказания. Поскольку 20-скрытая сеть включает почти 7000 весов и смещений, сеть, как правило, способна тесно соответствовать набору обучающих данных, но теряет способность обобщения. Компромисс между интенсивным обучением и точностью предсказания является одним из фундаментальных ограничений нейронных сетей.
net20 = train(net20,P,T); O20 = sim(net20,P); numWeightsAndBiases = length(getx(net20))
numWeightsAndBiases =
6883
Можно отобразить матрицы неточностей для подмножеств обучения, валидации и тестирования, нажав на соответствующую кнопку в окне инструмента обучения.
Можно детально оценить предсказания структуры путем вычисления индексов качества предсказания [3], которые указывают, насколько хорошо предсказано конкретное состояние и произошло ли сверхпредсказание или недопрогнозирование. Зададим индекс pcObs (S) для состояния S (S = {C, E, H}) поскольку количество остатков, правильно предсказанное в состоянии S, разделено на количество остатков, наблюдаемых в состоянии S. Точно так же мы задаем индекс pcPred(S) для состояния S как количество остатков, правильно предсказанных в состоянии S, деленное на количество остатков, предсказанных в состоянии S.
[i,j] = find(compet(O)); [u,v] = find(T); % === compute fraction of correct predictions when a given state is observed pcObs(1) = sum(i == 1 & u == 1)/sum (u == 1); % state C pcObs(2) = sum(i == 2 & u == 2)/sum (u == 2); % state E pcObs(3) = sum(i == 3 & u == 3)/sum (u == 3); % state H % === compute fraction of correct predictions when a given state is predicted pcPred(1) = sum(i == 1 & u == 1)/sum (i == 1); % state C pcPred(2) = sum(i == 2 & u == 2)/sum (i == 2); % state E pcPred(3) = sum(i == 3 & u == 3)/sum (i == 3); % state H % === compare quality indices of prediction figure() bar([pcObs' pcPred'] * 100); ylabel('Correctly predicted positions (%)'); ax = gca; ax.XTickLabel = {'C';'E';'H'}; legend({'Observed','Predicted'});
Эти индексы качества полезны для интерпретации точности предсказания. На самом деле, в случаях, когда метод предсказания имеет тенденцию переопределять/занижать заданное состояние, высокая/низкая точность предсказания может быть просто программным продуктом и не обеспечивает меры качества для самого метода.
Представленный здесь способ предсказывает структурное состояние данного белкового остатка на основе структурного состояния его соседей. Однако существуют дополнительные ограничения при прогнозировании содержимого структурных элементов в белке, такие как минимальная длина каждого структурного элемента. В частности, спираль относят к любой группе из четырех или более смежных остатков, а лист - к любой группе из двух или более смежных остатков. Для включения этого типа информации может быть создана дополнительная сеть, так что первая сеть предсказывает структурное состояние из аминокислотной последовательности, а вторая сеть предсказывает структурный элемент из структурного состояния.
nntraintool('close') % close Neural Network Training Tool
[1] Rost, B. and Sander, C., «Предсказания белковой вторичной структуры лучше, чем 70% точность», Journal of Molecular Biology, 232 (2): 584-99, 1993.
[2] Holley, L.H. and Karplus, M., «Protein secondary structure предсказания with a neural network», PNAS, 86 (1): 152-6, 1989.
[3] Kabsch, W., and Sander, C., «Насколько хороши предсказания белковой вторичной структуры?», FEBS Букв, 155 (2): 179-82, 1983.