Введение

Диагностика протеомного шаблона сыворотки может использоваться, чтобы дифференцировать выборки от пациентов с заболеваниями и без них. Шаблоны профиля генерируют с помощью масс-спектрометрии белка десорбции и ионизации с усиленной поверхностью (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 показывает некоторый peaks, которые могут быть полезны для классификации данных.

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 / значимых функций.

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% -ного удержания, чтобы получить лучшую оценку эффективности классификатора. cvpartition позволяет вам задать индексы обучения и тестирования для различных типов методов оценки системы, таких как hold-out, 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] примените эту концепцию к анализу шаблонов экспрессии белка. The 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), к-ближайшие соседи (fitcknn), нейронных сетей (Deep Learning Toolbox™) и классификационных деревьев (fitctree). Для выбора признаков можно также использовать последовательный выбор признаков подмножества (sequentialfs) или оптимизировать рандомизированные методы поиска с помощью генетического алгоритма (Global Optimization Toolbox). Для примера смотрите Поиск Генетического Алгоритма для Функций в Данных Масс-Спектрометрии.