Поиск генетического алгоритма функций в данных о масс-спектрометрии

В этом примере показано, как использовать Global Optimization Toolbox с Bioinformatics Toolbox™, чтобы оптимизировать поиск функций, чтобы классифицировать масс-спектрометрию (SELDI) данные.

Введение

Генетические алгоритмы оптимизируют результаты поиска для проблем с большими наборами данных. Можно использовать функцию генетического алгоритма MATLAB®, чтобы решить эти задачи в Биоинформатике. Генетические алгоритмы были применены к филогенетическому древовидному созданию, экспрессии гена и анализу данных масс-спектрометрии и многим другим областям Биоинформатики, которые имеют большие и в вычислительном отношении дорогие проблемы. Этот пример ищет оптимальные функции (peaks) в данных о масс-спектрометрии. Мы будем искать определенный peaks в данных, которые отличают больных раком от пациентов управления.

Global Optimization Toolbox

Сначала ознакомьте себя с Global Optimization Toolbox. Документация описывает, как работает генетический алгоритм и как использовать его в MATLAB. Чтобы получить доступ к документации, используйте команду документа.

doc ga

Предварительно обработайте данные о масс-спектрометрии

Данные в этом примере доступны от FDA-NCI Клинический Банк данных Программы Протеомики в http://home.ccr.cancer.gov/ncifdaproteomics/ppatterns.asp. Это - набор выборок от 121 больного раком яичника и 95 пациентов управления.

Этот пример принимает, что у вас уже есть предварительно обработанные данные 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®

Если у вас есть предварительно обработанные данные, можно загрузить их в 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, чтобы найти 20 отличительных функций

Используйте 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])

Для просмотра документации необходимо авторизоваться на сайте