В этом примере показано, как классифицировать данные о масс-спектрометрии и использовать некоторые статистические инструменты, чтобы искать потенциальные маркеры болезни и протеомную диагностику шаблона.
Сыворотка протеомная диагностика шаблона может использоваться, чтобы дифференцировать выборки от пациентов с и без болезни. Шаблоны профиля сгенерированы с помощью улучшенной поверхностью лазерной десорбции и ионизации (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).
Теперь, когда вы идентифицировали некоторые значительные функции, можно использовать эту информацию, чтобы классифицировать рак и нормальные выборки. Из-за небольшого количества выборок, можно запустить перекрестную проверку с помощью 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
Лилиэн и др. представил в [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.
classperf
| msnorm
| rankfeatures