В этом примере показан метод прогнозирования вторичной структуры, который использует нейронную сеть прямой связи и функциональные возможности, доступные в 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) обучающий набор, который используется для вычисления градиента и обновления весов и смещений сети; набор проверки достоверности, ошибка которого отслеживается в процессе обучения, поскольку она имеет тенденцию увеличиваться при переполнении данных; и 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 является относительно медленным, но эффективным с точки зрения памяти, в то время как алгоритм Левенберга-Марквардта является более быстрым, но более требовательным с точки зрения памяти.
Увеличьте количество скрытых нейронов. Добавляя больше скрытых единиц, мы, как правило, получаем более сложную сеть с потенциалом для повышения производительности, но мы должны быть осторожны, чтобы не перегружать данные.
При создании сети распознавания образов можно указать больше скрытых слоев или увеличить размер скрытых слоев, как показано ниже:
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], которые указывают, насколько хорошо предсказано конкретное состояние и имело ли место чрезмерное или недостаточное предсказание. Мы определяем индекс pcOb (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
Рост, В. и Сандер, С., «Предсказание вторичной структуры белка с точностью выше 70%», Journal of Molecular Biology, 232 (2): 584-99, 1993.
[2] Холли, Л.Х. и Карплюс, М., «Предсказание вторичной структуры белка с помощью нейронной сети», PNAS, 86 (1): 152-6, 1989.
[3] Кабш, В. и Сандер, С., «Насколько хороши прогнозы вторичной структуры белка?», FEBS Letters, 155 (2): 179-82, 1983.