Введение

Сывороточная протеомная диагностика может быть использована для дифференциации образцов от пациентов с заболеванием и без него. Профили формируются с использованием масс-спектрометрии белка с улучшенной лазерной десорбцией и ионизацией (SELDI). Эта технология может улучшить клинические диагностические тесты на патологии рака. Цель состоит в том, чтобы выбрать уменьшенный набор измерений или «особенностей», которые можно использовать для различия между раком и контрольными пациентами. Этими характеристиками будут уровни интенсивности ионов при определенных значениях массы/заряда.

Данные предварительной обработки

Набор данных о раке яичников в этом примере взят из банка данных программы клинической протеомики FDA-NCI. Набор данных был создан с использованием массива WCX2 белков. Набор данных включает 95 контролей и 121 рак яичников. Подробное описание этого набора данных см. в разделах [1] и [4].

В этом примере предполагается, что предварительно обработанные данные уже имеются. OvarianCancerQAQCdataset.mat. Однако, если файл данных отсутствует, его можно создать заново, выполнив шаги в примере Пакетная обработка спектров с использованием последовательных и параллельных вычислений.

Кроме того, можно запустить сценарий msseqprocessing.m.

addpath(fullfile(matlabroot,'examples','bioinfo','main')) % Make sure the supporting files are on the search path
type msseqprocessing

% MSSEQPROCESSING Script to create OvarianCancerQAQCdataset.mat (used in
% CANCERDETECTDEMO). Before running this file initialize the variable
% "repository" to the full path where you placed you mass-spectrometry
% files. For Example:
%
%   repository = 'F:/MassSpecRepository/OvarianCD_PostQAQC/';
% 
% or
%
%   repository = '/home/username/MassSpecRepository/OvarianCD_PostQAQC/';
%
% The approximate time of execution is 18 minutes (Pentium 4, 4GHz). If you
% have the Parallel Computing Toolbox refer to BIODISTCOMPDEMO to see
% how you can speed this analysis up.

%   Copyright 2003-2008 The MathWorks, Inc.


repositoryC = [repository 'Cancer/'];
repositoryN = [repository 'Normal/'];

filesCancer = dir([repositoryC '*.txt']);
NumberCancerDatasets = numel(filesCancer);
fprintf('Found %d Cancer mass-spectrograms.\n',NumberCancerDatasets)
filesNormal = dir([repositoryN '*.txt']);
NumberNormalDatasets = numel(filesNormal);
fprintf('Found %d Control mass-spectrograms.\n',NumberNormalDatasets)

files = [ strcat('Cancer/',{filesCancer.name}) ...
          strcat('Normal/',{filesNormal.name})];
N = numel(files);   % total number of files

fprintf('Total %d mass-spectrograms to process...\n',N)

[MZ,Y] = msbatchprocessing(repository,files);

disp('Finished; normalizing and saving to OvarianCancerQAQCdataset.mat.')
Y = msnorm(MZ,Y,'QUANTILE',0.5,'LIMITS',[3500 11000],'MAX',50);

grp = [repmat({'Cancer'},size(filesCancer));...
       repmat({'Normal'},size(filesNormal))];
   
save OvarianCancerQAQCdataset.mat Y MZ grp

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

Загрузить данные

Предварительно обработанные данные можно загрузить в MATLAB.

load OvarianCancerQAQCdataset
whos

  Name          Size                Bytes  Class     Attributes

  MZ        15000x1                120000  double              
  Y         15000x216            25920000  double              
  grp         216x1                 25056  cell                

Существует три переменные: MZ, Y, grp. MZ - вектор масса/заряд, Y - значения интенсивности для всех 216 пациентов (контроль и рак), и grp содержит индексную информацию о том, какие из этих образцов представляют раковых пациентов, а какие - нормальных пациентов.

Инициализируйте некоторые переменные, которые будут использоваться в примере.

N = numel(grp);                         % Number of samples
Cidx = strcmp('Cancer',grp);            % Logical index vector for Cancer samples
Nidx = strcmp('Normal',grp);            % Logical index vector for Normal samples
Cvec = find(Cidx);                      % Index vector for Cancer samples
Nvec = find(Nidx);                      % Index vector for Normal samples
xAxisLabel = 'Mass/Charge (M/Z)';       % x label for plots
yAxisLabel = 'Ion Intensity';           % y label for plots

Визуализация некоторых образцов

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

figure; hold on;
hC = plot(MZ,Y(:,Cvec(1:5)),'b');
hN = plot(MZ,Y(:,Nvec(1:5)),'g');
xlabel(xAxisLabel); ylabel(yAxisLabel);
axis([2000 12000 -5 60])
legend([hN(1),hC(1)],{'Control','Ovarian Cancer'})
title('Multiple Sample Spectrograms')

Увеличение в области от 8500 до 8700 M/Z показывает некоторые пики, которые могут быть полезны для классификации данных.

axis([8450,8700,-1,7])

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

mean_N = mean(Y(:,Nidx),2);  % group average for control samples
max_N = max(Y(:,Nidx),[],2); % top envelopes of the control samples
min_N = min(Y(:,Nidx),[],2); % bottom envelopes of the control samples
mean_C = mean(Y(:,Cidx),2);  % group average for cancer samples
max_C = max(Y(:,Cidx),[],2); % top envelopes of the control samples
min_C = min(Y(:,Cidx),[],2); % bottom envelopes of the control samples

figure; hold on;
hC = plot(MZ,mean_C,'b');
hN = plot(MZ,mean_N,'g');
gC = plot(MZ,[max_C min_C],'b--');
gN = plot(MZ,[max_N min_N],'g--');
xlabel(xAxisLabel); ylabel(yAxisLabel);
axis([8450,8700,-1,7])
legend([hN,hC,gN(1),gC(1)],{'Control Group Avg.','Ovarian Cancer Group Avg',...
                            'Control Envelope','Ovarian Cancer Envelope'},...
                            'Location', 'NorthWest')
title('Group Average and Group Envelopes')

Заметьте, что, по-видимому, нет единой функции, которая могла бы отлично различать обе группы.

Ранжирование ключевых функций

Простой подход для нахождения значимых признаков заключается в предположении, что каждое значение M/Z является независимым и вычисляет двусторонний t-тест. rankfeatures возвращает индекс к наиболее значимым значениям M/Z, например, 100 индексов, ранжированных абсолютным значением тестовой статистики. Этот метод выбора функций также известен как метод фильтрации, в котором алгоритм обучения не участвует в выборе функций.

[feat,stat] = rankfeatures(Y,grp,'CRITERION','ttest','NUMBER',100);

Первый выход rankfeatures может использоваться для извлечения значений M/Z значимых признаков.

sig_Masses = MZ(feat);
sig_Masses(1:7)' %display the first seven

ans =

   1.0e+03 *

    8.1009    8.1016    8.1024    8.1001    8.1032    7.7366    7.7359

Второй выход rankfeatures - вектор с абсолютным значением тестовой статистики. Вы можете нарисовать его по спектрам, используя yyaxis.

figure;

yyaxis left
plot(MZ, [mean_N mean_C]);
ylim([-1,20])
xlim([7950,8300])
title('Significant M/Z Values')
xlabel(xAxisLabel);
ylabel(yAxisLabel);

yyaxis right
plot(MZ,stat);
ylim([-1,22])
ylabel('Test Statistic');

legend({'Control Group Avg.','Ovarian Cancer Group Avg.', 'Test Statistics'})

Обратите внимание, что существуют значительные области с высокими значениями M/Z, но низкой интенсивностью (~ 8100 Да.). Другие подходы к измерению разделяемости классов доступны вrankfeatures, такие как на основе энтропии, бхаттачарья или область под кривой эмпирической рабочей характеристики приемника (ROC).

Слепая классификация с использованием линейного дискриминантного анализа (LDA)

Теперь, когда вы определили некоторые важные особенности, вы можете использовать эту информацию для классификации рака и нормальных образцов. Из-за небольшого количества выборок можно выполнить перекрестную проверку с использованием 20% holdout, чтобы лучше оценить производительность классификатора. cvpartition позволяет установить показатели обучения и тестирования для различных типов методов оценки системы, таких как удержание, K-fold и Leave-M-Out.

per_eval = 0.20;          % training size for cross-validation
rng('default');             % initialize random generator to the same state
                          % used to generate the published example
cv  = cvpartition(grp,'holdout',per_eval)

cv = 

Hold-out cross validation partition
   NumObservations: 216
       NumTestSets: 1
         TrainSize: 173
          TestSize: 43

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

cp_lda1 = classperf(grp); % initializes the CP object
for k=1:10 % run cross-validation 10 times
    cv = repartition(cv);
    feat = rankfeatures(Y(:,training(cv)),grp(training(cv)),'NUMBER',100);
    c = classify(Y(feat,test(cv))',Y(feat,training(cv))',grp(training(cv)));
    classperf(cp_lda1,c,test(cv)); % updates the CP object with current validation
end

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

cp_lda1

                        Label: ''
                  Description: ''
                  ClassLabels: {2x1 cell}
                  GroundTruth: [216x1 double]
         NumberOfObservations: 216
               ControlClasses: 2
                TargetClasses: 1
            ValidationCounter: 10
           SampleDistribution: [216x1 double]
            ErrorDistribution: [216x1 double]
    SampleDistributionByClass: [2x1 double]
     ErrorDistributionByClass: [2x1 double]
               CountingMatrix: [3x2 double]
                  CorrectRate: 0.8488
                    ErrorRate: 0.1512
              LastCorrectRate: 0.8837
                LastErrorRate: 0.1163
             InconclusiveRate: 0
               ClassifiedRate: 1
                  Sensitivity: 0.8208
                  Specificity: 0.8842
      PositivePredictiveValue: 0.8995
      NegativePredictiveValue: 0.7962
           PositiveLikelihood: 7.0890
           NegativeLikelihood: 0.2026
                   Prevalence: 0.5581
              DiagnosticTable: [2x2 double]

Этот наивный подход к выбору признаков может быть улучшен путем исключения некоторых признаков на основе региональной информации. Например, «NWEIGHT» в rankfeatures перевешивает тестовую статистику соседних признаков M/Z, так что другие значимые значения M/Z могут быть включены в подмножество выбранных признаков

cp_lda2 = classperf(grp); % initializes the CP object
for k=1:10 % run cross-validation 10 times
    cv = repartition(cv);
    feat = rankfeatures(Y(:,training(cv)),grp(training(cv)),'NUMBER',100,'NWEIGHT',5);
    c = classify(Y(feat,test(cv))',Y(feat,training(cv))',grp(training(cv)));
    classperf(cp_lda2,c,test(cv)); % updates the CP object with current validation
end
cp_lda2.CorrectRate % average correct classification rate

ans =

    0.9023

Уменьшение размерности данных PCA/LDA

Lilien et al. представлен в [2] алгоритм уменьшения размерности данных, использующий анализ основных компонентов (PCA), затем LDA используется для классификации групп. В этом примере 2000 наиболее значимых элементов в пространстве M/Z отображаются на 150 главных компонентов.

cp_pcalda = classperf(grp); % initializes the CP object
for k=1:10 % run cross-validation 10 times
    cv = repartition(cv);
    % select the 2000 most significant features.
    feat = rankfeatures(Y(:,training(cv)),grp(training(cv)),'NUMBER',2000);
    % PCA to reduce dimensionality
    P = pca(Y(feat,training(cv))');
    % Project into PCA space
    x = Y(feat,:)' * P(:,1:150);
    % Use LDA
    c = classify(x(test(cv),:),x(training(cv),:),grp(training(cv)));
    classperf(cp_pcalda,c,test(cv));
end
cp_pcalda.CorrectRate % average correct classification rate

ans =

    0.9814

Рандомизированный поиск выбора элементов подмножества

Выбор признаков также может быть усилен классификацией, этот подход обычно называют методом выбора обертки. Рандомизированный поиск выбора признаков генерирует случайные подмножества признаков и оценивает их качество независимо с помощью алгоритма обучения. Позже он выбирает пул наиболее частых хороших функций. Li et al. в [3] применяют эту концепцию к анализу паттернов экспрессии белка. randfeatures функция позволяет искать подмножество элементов с помощью LDA или k-ближайшего соседнего классификатора по рандомизированным подмножествам элементов.

Примечание.Следующий пример является вычислительно интенсивным, поэтому он был отключен из примера. Кроме того, для получения лучших результатов следует увеличить размер пула и жесткость классификатора по сравнению со значениями по умолчанию в randfeatures. Напечатать help randfeatures для получения дополнительной информации.

if 0  % <== change to 1 to enable. This may use extensive time to complete.
   cv = repartition(cv);
   [feat,fCount] = randfeatures(Y(:,training(cv)),grp(training(cv)),...
                         'CLASSIFIER','da','PerformanceThreshold',0.90);
else
   load randFeatCancerDetect
end

Оценка качества выбранных элементов с помощью оценочного набора

Первый выход из randfeatures - упорядоченный список индексов значений MZ. Первый элемент чаще всего встречается в подмножествах, где была достигнута хорошая классификация. Второй вывод - это фактическое число выборок каждого значения. Вы можете использовать hist чтобы посмотреть на это распределение.

figure;
hist(fCount,max(fCount)+1);

Большинство значений отображается в выбранном подмножестве не более одного раза. Увеличение дает лучшее представление о деталях для более часто выбираемых значений.

axis([0 80 0 100])

Только несколько значений были выбраны более 10 раз. Их можно визуализировать с помощью штативного графика для отображения наиболее часто выбираемых элементов.

figure; hold on;
sigFeats = fCount;
sigFeats(sigFeats<=10) = 0;
plot(MZ,[mean_N mean_C]);
stem(MZ(sigFeats>0),sigFeats(sigFeats>0),'r');
axis([2000,12000,-1,80])
legend({'Control Group Avg.','Ovarian Cancer Group Avg.','Significant Features'}, ...
    'Location', 'NorthWest')
xlabel(xAxisLabel); ylabel(yAxisLabel);

Эти функции объединяются в несколько групп. Чтобы выяснить, сколько функций имеет значение, выполните следующий эксперимент. Наиболее часто выбираемый элемент используется для классификации данных, затем используются два наиболее часто выбираемых элемента и так далее до тех пор, пока не будут использованы все элементы, которые были выбраны более 10 раз. Затем можно увидеть, улучшает ли добавление дополнительных функций классификатор.

nSig = sum(fCount>10);
cp_rndfeat = zeros(20,nSig);
for i = 1:nSig
    for j = 1:20
        cv = repartition(cv);
        P = pca(Y(feat(1:i),training(cv))');
        x = Y(feat(1:i),:)' * P;
        c = classify(x(test(cv),:),x(training(cv),:),grp(training(cv)));
        cp = classperf(grp,c,test(cv));
        cp_rndfeat(j,i) = cp.CorrectRate; % average correct classification rate
    end
end
figure
plot(1:nSig, [max(cp_rndfeat);mean(cp_rndfeat)]);
legend({'Best CorrectRate','Mean CorrectRate'},'Location', 'SouthEast')

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

[bestAverageCR, bestNumFeatures] = max(mean(cp_rndfeat));

Теперь можно визуализировать элементы, которые дают наилучшую среднюю классификацию. Вы можете видеть, что они на самом деле соответствуют только трем пикам в данных.

figure; hold on;
sigFeats = fCount;
sigFeats(sigFeats<=10) = 0;
ax_handle = plot(MZ,[mean_N mean_C]);
stem(MZ(feat(1:bestNumFeatures)),sigFeats(feat(1:bestNumFeatures)),'r');
axis([7650,8850,-1,80])
legend({'Control Group Avg.','Ovarian Cancer Group Avg.','Significant Features'})
xlabel(xAxisLabel); ylabel(yAxisLabel);

Альтернативные алгоритмы статистического обучения

В MATLAB ® существует множество инструментов классификации, которые также можно использовать для анализа протеомных данных. Среди них опорные векторные машины (fitcsvm), k-ближайшие соседи (fitcknn), нейронные сети (Deep Learning Toolbox™) и деревья классификации (fitctree). Для выбора элемента можно также использовать выбор последовательного подмножества элементов (sequentialfs) или оптимизировать методы рандомизированного поиска с помощью генетического алгоритма (Global Optimization Toolbox). Например, см. раздел Поиск элементов генетического алгоритма в данных масс-спектрометрии.