Идентификация значительных функций и классификация профилей белка

Открыть скрипт

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

Введение

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

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

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')

Увеличивание масштаб области от 8 500 до 8 700 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/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%-й затяжки, чтобы иметь лучшую оценку эффективности классификатора. cvpartition позволяет вам устанавливать обучение и тестовые индексы для различных типов системных методов оценки, таких как затяжка, K-сгиб и 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 Размерности Данных

Лилиэн и др. представил в [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

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

Выбор признаков может также быть укреплен классификацией, этот подход обычно упоминается как метод выбора обертки. Рандомизированный поиск выбора признаков генерирует случайные подмножества функций и оценивает их качество независимо с алгоритмом обучения. Позже, это выбирает пул самых частых хороших функций. Литий и др. в [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));

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

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). Например, смотрите Поиск Генетического алгоритма Функций в Данных о Масс-спектрометрии.

Ссылки

[1] Conrads, T P, V Fusaro, С Росс, Д Йохан, V Раджапаксе, Б А Хитт, С М Стайнберг, и др. “Сыворотка с высоким разрешением Протеомные Функции Диагностики Рака яичника”. Эндокринно-связанный Рак, июнь 2004, 163–78.

[2] Лилиэн, Райан Х., Хэни Фарид и Брюс Р. Дональд. “Вероятностная Классификация Болезней Зависимых Выражением Протеомных Данных из Масс-спектрометрии Человеческой Сыворотки”. Журнал Вычислительной Биологии 10, № 6 (декабрь 2003): 925–46.

[3] Литий, L., Д. М. Умбах, П. Терри и Дж. А. Тейлор. “Применение Метода GA/KNN к Данным о Протеомике SELDI”. Биоинформатика 20, № 10 (1 июля 2004): 1638–40.

[4] Petricoin, Эмануэль Ф, Али М Ардекэни, Бен А Хитт, Питер Дж Левин, Винсент А Фусаро, Сет М Стайнберг, Гордон Б Миллз, и др. “Использование Протеомных Шаблонов в Сыворотке, чтобы Идентифицировать Рак яичника”. The Lancet 359, № 9306 (февраль 2002): 572–77.

Смотрите также

msnorm | rankfeatures | classperf

Документация