В этом примере показано, как использовать Global Optimization Toolbox с Bioinformatics Toolbox™, чтобы оптимизировать поиск функций, чтобы классифицировать масс-спектрометрию (SELDI) данные.
Генетические алгоритмы оптимизируют результаты поиска для проблем с большими наборами данных. Можно использовать функцию генетического алгоритма MATLAB®, чтобы решить эти задачи в Биоинформатике. Генетические алгоритмы были применены к филогенетическому древовидному созданию, экспрессии гена и анализу данных масс-спектрометрии и многим другим областям Биоинформатики, которые имеют большие и в вычислительном отношении дорогие проблемы. Этот пример ищет оптимальные функции (peaks) в данных о масс-спектрометрии. Мы будем искать определенный peaks в данных, которые отличают больных раком от пациентов управления.
Сначала ознакомьте себя с Global Optimization Toolbox. Документация описывает, как работает генетический алгоритм и как использовать его в MATLAB. Чтобы получить доступ к документации, используйте команду документа.
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, группа. MZ является вектором массы/заряда, Y является значениями интенсивности для всех 216 пациентов (управление и рак), и группа содержит информацию индекса, относительно которой из этих выборок представляют больных раком и которые представляют нормальных пациентов. Чтобы визуализировать эти данные, смотрите, что пример Идентифицирует Значительные Функции и Классифицирует Профили Белка.
Инициализируйте переменные, используемые в примере.
[numPoints, numSamples] = size(Y); % total number of samples and data points id = grp2idx(grp); % ground truth: Cancer=1, Control=2
Генетический алгоритм требует целевой функции, также известной как функцию фитнеса, которая описывает явление, которое мы хотим оптимизировать. В этом примере машинное оборудование генетического алгоритма тестирует небольшие подмножества значений M/Z с помощью функции фитнеса и затем определяет, которому значения M/Z передаются или удалили из каждой последующей генерации. Функция фитнеса biogafit передается решателю генетического алгоритма с помощью указателя на функцию. В этом примере biogafit максимизирует отделимость двух классов при помощи линейной комбинации 1) по опыту вероятности и 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 поколений. Однако можно установить 'Поколения' на большее значение.
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');
Наблюдайте интересные пиковые приблизительно 8 100 дальтонов., который, кажется, смещен направо на здоровых выборках.
axis([8082 8128 -.5 4])
[1] Conrads, T P, V Fusaro, С Росс, Д Йохан, V Раджапаксе, Б А Хитт, С М Стайнберг, и др. “Сыворотка с высоким разрешением Протеомные Функции Диагностики Рака яичника”. Эндокринно-связанный Рак, июнь 2004, 163–78.
[2] Petricoin, Эмануэль Ф, Али М Ардекэни, Бен А Хитт, Питер Дж Левин, Винсент А Фусаро, Сет М Стайнберг, Гордон Б Миллз, и др. “Использование Протеомных Шаблонов в Сыворотке, чтобы Идентифицировать Рак яичника”. The Lancet 359, № 9306 (февраль 2002): 572–77.