Этот пример показывает, как классифицировать данные о масс-спектрометрии и показывает некоторые статистические инструменты, которые могут использоваться, чтобы искать потенциальные маркеры болезни и протеомную диагностику шаблона.
Сыворотка протеомная диагностика шаблона может использоваться, чтобы дифференцировать выборки от пациентов с и без болезни. Шаблоны профиля сгенерированы с помощью улучшенной поверхностью лазерной десорбции и ионизации (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 26784 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, таковы как базирующаяся энтропия, Bhattacharyya или область под эмпирической кривой рабочей характеристики получателя (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-nearest по рандомизированным подмножествам функций.
Примечание: следующий пример в вычислительном отношении интенсивен, таким образом, он был отключен от примера. Кроме того, для лучших результатов необходимо увеличить размер пула и строгость классификатора от значений по умолчанию в 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-nearest соседи (fitcknn), нейронные сети (Deep Learning Toolbox™) и деревья классификации (fitctree). Для выбора функции можно также использовать последовательный выбор функции подмножества (sequentialfs) или оптимизировать рандомизированные методы поиска при помощи генетического алгоритма (Global Optimization Toolbox). Например, смотрите Поиск Генетического алгоритма Функций в Данных о Масс-спектрометрии.
[1] Conrads, T.P., и др., "Сыворотка с высоким разрешением протеомные функции яичникового обнаружения", Эндокринно-связанный Рак, 11 (2):163-78, 2004.
[2] Лилиэн, R.H., и др., "Вероятностная Классификация Болезней Зависимых Выражением Протеомных Данных из Масс-спектрометрии Человеческой Сыворотки", Журнал Вычислительной Биологии, 10 (6):925-46, 2003.
[3] Литий, L., и др., "Применение метода GA/KNN к данным о протеомике SELDI", Биоинформатика, 20 (10):1638-40, 2004.
[4] Petricoin, E.F., и др., "Использование протеомных шаблонов в сыворотке, чтобы идентифицировать рак яичника", Ланцет, 359 (9306):572-7, 2002.