В этом примере показано, как использовать Global Optimization Toolbox с Bioinformatics Toolbox™ для оптимизации поиска функций для классификации данных масс-спектрометрии (SELDI).
Генетические алгоритмы оптимизируют результаты поиска для задач с большими наборами данных. Вы можете использовать функцию генетического алгоритма MATLAB ®, чтобы решить эти проблемы в биоинформатике. Генетические алгоритмы были применены к филогенетическим созданиям дерева, экспрессии генов и масс-спектрометрии данных анализу и многим другим областям биоинформатики, которые имеют большие и вычислительно дорогие проблемы. Этот пример ищет оптимальные функции ( peaks) в данных масс-спектрометрии. Мы будем искать конкретный peaks в данных, которые отличают больных раком от контрольных пациентов.
Сначала ознакомьтесь с Global Optimization Toolbox. Документация описывает, как работает генетический алгоритм и как его использовать в MATLAB. Для доступа к документации используйте команду doc.
doc ga
Исходные данные в этом примере получены из Банка данных программы клинической протеомики FDA-NCI. Это набор выборок от 121 пациента с раком яичников и 95 пациентов с контролем. Подробное описание этого набора данных см. в разделах [1] и [2].
Этот пример предполагает, что у вас уже есть предварительно обработанные данные 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 содержит индекс информацию о том, какой из этих выборок представляет онкологических пациентов и какие таковые представляют нормальных пациентов. Чтобы визуализировать эти данные, смотрите пример «Идентификация значимых функций и классификация профилей белка».
Инициализируйте переменные, используемые в примере.
[numPoints, numSamples] = size(Y); % total number of samples and data points id = grp2idx(grp); % ground truth: Cancer=1, Control=2
Генетический алгоритм требует целевой функции, также известной как функция соответствия, которая описывает феномен, который мы хотим оптимизировать. В этом примере машина генетического алгоритма проверяет небольшие подмножества значений M/Z, используя функцию соответствия, и затем определяет, какие значения M/Z передаются или удаляются из каждой последующей генерации. Биогафит функции соответствия передается решателю генетического алгоритма с помощью указателя на функцию. В этом примере биогафит максимизирует разделяемость двух классов с помощью линейной комбинации 1) a-апостериорная вероятность и 2) эмпирическая вероятность ошибки линейного классификатора (классификация). Можно создать свою собственную функцию соответствия, чтобы попробовать различные классификаторы или альтернативные методы для оценки эффективности классификаторов.
type biogafit
function classPerformance = biogafit(thePopulation,Y,id)
%BIOGAFIT The fitness function for BIOGAMSDEMO
%
% This function uses the classify function to measure how well mass
% spectrometry data is grouped using certain masses. The input argument
% thePopulation is a vector of row indices from the mass spectrometry
% data Y. Classification performance is a linear combination of the error
% rate and the posteriori probability of the classifier.
% Copyright 2003-2013 The MathWorks, Inc.
thePopulation = round(thePopulation);
try
[c,~,p] = classify(Y(thePopulation,:)',Y(thePopulation,:)',double(id),'linear');
cp = classperf(id,c);
classPerformance = 100*cp.ErrorRate + 1 - mean(max(p,[],2));
catch
% In case pooled covariance matrix is not positive definite we try a
% naive-Bayes classifier:
try
[c,~,p] = classify(Y(thePopulation,:)',Y(thePopulation,:)',double(id),'diaglinear');
cp = classperf(id,c);
classPerformance = 100*cp.ErrorRate + 1 - mean(max(p,[],2));
catch
classPerformance = Inf;
end
end
Пользователи могут изменить, как оптимизация выполняется генетическим алгоритмом, создав пользовательские функции для кроссовера, масштабирования соответствия, мутации, выбора и создания населения. В этом примере вы будете использовать функцию biogacreate, написанную для этого примера, чтобы создать начальные случайные точки данных из данных масс-спектрометрии. Заголовок функции требует определенных входных параметров, заданных документацией GA. В тулбоксе есть функция создания по умолчанию для создания начальных генеральных совокупностей точек данных.
type biogacreate
function pop = biogacreate(GenomeLength,~,options,Y,id) %BIOGACREATE Population creation function for MSGADEMO % % This function creates a population matrix with dimensions of % options.PopulationSize rows by the number of independent variables % (GenomeLength) columns. These values are integers that correspond to % randomly selected rows of the mass spectrometry data Y. Each row of the % population matrix is a random sample of row indices of the mass spec % data. % Note: This function's input arguments are required by the GA function. % See GA documentation for further detail. % Copyright 2005-2013 The MathWorks, Inc. pop = zeros(options.PopulationSize,GenomeLength); npop = numel(pop); ranked_features = rankfeatures(Y,id,'NumberOfIndices',npop,'NWeighting',.5); pop(:) = randsample(ranked_features,npop);
Функция GA использует структуру опций, чтобы удерживать параметры алгоритма, которые она использует при выполнении минимизации с помощью генетического алгоритма. Функция optimoptions создаст эту структуру опций. В целях этого примера генетический алгоритм будет запускаться только в течение 50 генерации. Однако вы можете задать 'Generations' большее значение.
options = optimoptions('ga','CreationFcn',{@biogacreate,Y,id},... 'PopulationSize',120,... 'Generations',50,... 'Display','iter')
options =
ga options:
Set properties:
CreationFcn: {1x3 cell}
Display: 'iter'
MaxGenerations: 50
PopulationSize: 120
Default properties:
ConstraintTolerance: 1.0000e-03
CrossoverFcn: @crossoverscattered
CrossoverFraction: 0.8000
EliteCount: '0.05*PopulationSize'
FitnessLimit: -Inf
FitnessScalingFcn: @fitscalingrank
FunctionTolerance: 1.0000e-06
HybridFcn: []
InitialPopulationMatrix: []
InitialPopulationRange: []
InitialScoresMatrix: []
MaxStallGenerations: 50
MaxStallTime: Inf
MaxTime: Inf
MutationFcn: {@mutationgaussian [1] [1]}
NonlinearConstraintAlgorithm: 'auglag'
OutputFcn: []
PlotFcn: []
PopulationType: 'doubleVector'
SelectionFcn: @selectionstochunif
UseParallel: 0
UseVectorized: 0
Используйте ga, чтобы запустить функцию генетического алгоритма. 100 групп по 20 точек данных каждая будет развиваться более 50 генерация. События отбора, кроссовера и мутации генерируют новое население в каждой генерации.
% initialize the random generators to the same state used to generate the % published example rng('default') nVars = 12; % set the number of desired features FitnessFcn = {@biogafit,Y,id}; % set the fitness function feat = ga(FitnessFcn,nVars,options); % call the Genetic Algorithm feat = round(feat); Significant_Masses = MZ(feat) cp = classperf(classify(Y(feat,:)',Y(feat,:)',id),id); cp.CorrectRate
Best Mean Stall
Generation Func-count f(x) f(x) Generations
1 240 2.827 8.928 0
2 354 2.827 8.718 1
3 468 0.9663 8.001 0
4 582 0.9516 7.249 0
5 696 0.9516 6.903 1
6 810 0.4926 6.804 0
7 924 0.4926 6.301 1
8 1038 0.02443 5.215 0
9 1152 0.02443 4.77 1
10 1266 0.02101 4.084 0
11 1380 0.02101 3.792 1
12 1494 0.01854 3.437 0
13 1608 0.01606 3.44 0
14 1722 0.01372 2.768 0
15 1836 0.01218 2.74 0
16 1950 0.01204 2.471 0
17 2064 0.01204 2.649 1
18 2178 0.01189 2.326 0
19 2292 0.01189 2.003 0
20 2406 0.0118 2.341 0
21 2520 0.01099 1.714 0
22 2634 0.01094 1.828 0
23 2748 0.01094 1.94 1
24 2862 0.01094 2.285 2
25 2976 0.009843 2.026 0
26 3090 0.009843 1.899 1
27 3204 0.009183 1.802 0
28 3318 0.007877 1.5 0
29 3432 0.007788 1.793 0
30 3546 0.007788 1.756 1
Best Mean Stall
Generation Func-count f(x) f(x) Generations
31 3660 0.007091 1.719 0
32 3774 0.006982 1.598 0
33 3888 0.006982 1.269 1
34 4002 0.006732 1.279 0
35 4116 0.005008 1.229 0
36 4230 0.004325 1.179 0
37 4344 0.004325 1.534 1
38 4458 0.003982 1.15 0
39 4572 0.003982 0.9602 1
40 4686 0.003982 0.8547 2
41 4800 0.003891 0.9083 0
42 4914 0.003683 0.7409 0
43 5028 0.003683 0.516 1
44 5142 0.003364 0.5153 0
45 5256 0.003172 0.4218 0
46 5370 0.003172 0.3783 1
47 5484 0.002997 0.1883 0
48 5598 0.002675 0.1297 0
49 5712 0.002611 0.04382 0
50 5826 0.002519 0.007859 0
Optimization terminated: maximum number of generations exceeded.
Significant_Masses =
1.0e+03 *
7.6861
7.9234
8.9834
8.6171
7.1808
7.3057
8.1132
8.5241
7.0527
7.7600
7.7442
7.7245
ans =
1
Чтобы визуализировать, какие функции были выбраны генетическим алгоритмом, данные строятся с пиковыми положениями, отмеченными красными вертикальными линиями.
xAxisLabel = 'Mass/Charge (M/Z)'; % x label for plots yAxisLabel = 'Relative Ion Intensity'; % y label for plots figure; hold on; hC = plot(MZ,Y(:,1:15) ,'b'); hN = plot(MZ,Y(:,141:155),'g'); hG = plot(MZ(feat(ceil((1:3*nVars )/3))), repmat([0 100 NaN],1,nVars),'r'); xlabel(xAxisLabel); ylabel(yAxisLabel); axis([1900 12000 -1 40]); legend([hN(1),hC(1),hG(1)],{'Control','Ovarian Cancer', 'Discriminative Features'}, ... 'Location', 'NorthWest'); title('MS Data with Discriminative Features found by Genetic Algorithm');
Наблюдайте интересный пик около 8100 Да., который, по-видимому, смещен вправо на здоровых выборках.
axis([8082 8128 -.5 4])
[1] Conrads, T P, V A Fusaro, S Ross, D Johann, V Rajapakse, B A Hitt, S M Steinberg, et al. «Сывороточные протеомные функции высокого разрешения для выявления рака яичников». Рак, связанный с эндокринной системой, июнь 2004, 163-78.
[2] Petricoin, Emanuel F, Ali M Ardekani, Ben A Hitt, Peter J Levine, Vincent A Fusaro, Seth M Steinberg, Gordon B Mills, et al. «Использование протеомных шаблонов в сыворотке для идентификации рака яичников». Lancet 359, № 9306 (февраль 2002): 572-77.