Введение

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

Набор данных Рост-Сандера [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);

Уточнение нейронной сети для получения более точных результатов

Нейронная сеть, которую мы определили, относительно проста. Чтобы добиться некоторых улучшений точности предсказания, мы могли бы попробовать одно из следующего:

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

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.