optGene

PURPOSE ^

OPTGENE implements the optgene algorithm.

SYNOPSIS ^

function [x, population, scores, optGeneSol] = optGene(model, targetRxn, substrateRxn, generxnList, MaxKOs, population)

DESCRIPTION ^

OPTGENE implements the optgene algorithm.

 [x, population, scores, optGeneSol] = optGene(model, targetRxn, substrateRxn, generxnList, MaxKOs, population)

INPTUS
 mode          Model of reconstruction
 targetRxn     String name of reaction which is to be maximized.
 generxnList   List of genes or rxns which can be knocked out.  The
               program will guess which of the two it is, based on the 
               content in model.

OPTIONAL INPUT
 population    population matrix (binary matrix).  Use this
               parameter to interrupt simulation and resume afterwards.

OUTPUTS
   x           best optimized value found
   population  Population of individuals.  Pass this back into optgene to
               continue simulating where you left off.
   scores      An array of scores
   optGeneSol  optGene solution strcture

 Jan Schellenberger and Adam Feist 04/08/08

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [x, population, scores, optGeneSol] = optGene(model, targetRxn, substrateRxn, generxnList, MaxKOs, population)
0002 %OPTGENE implements the optgene algorithm.
0003 %
0004 % [x, population, scores, optGeneSol] = optGene(model, targetRxn, substrateRxn, generxnList, MaxKOs, population)
0005 %
0006 %INPTUS
0007 % mode          Model of reconstruction
0008 % targetRxn     String name of reaction which is to be maximized.
0009 % generxnList   List of genes or rxns which can be knocked out.  The
0010 %               program will guess which of the two it is, based on the
0011 %               content in model.
0012 %
0013 %OPTIONAL INPUT
0014 % population    population matrix (binary matrix).  Use this
0015 %               parameter to interrupt simulation and resume afterwards.
0016 %
0017 %OUTPUTS
0018 %   x           best optimized value found
0019 %   population  Population of individuals.  Pass this back into optgene to
0020 %               continue simulating where you left off.
0021 %   scores      An array of scores
0022 %   optGeneSol  optGene solution strcture
0023 %
0024 % Jan Schellenberger and Adam Feist 04/08/08
0025 
0026 % hash table for hashing results... faster than not using it.
0027 global HTABLE
0028 HTABLE = java.util.Hashtable;
0029 global MaxKnockOuts
0030 
0031 ngenes = length(generxnList);
0032 
0033 % figure out if list is genes or reactions
0034 rxnok = 1;
0035 geneok = 1;
0036 for i = 1:length(generxnList)
0037     if(~ ismember(generxnList{i}, model.rxns)),  rxnok = 0; end
0038     if(~ ismember(generxnList{i}, model.genes)),geneok = 0; end
0039 end
0040 if geneok
0041     display('assuming list is genes');
0042 elseif rxnok
0043     display('assuming list is reactions');
0044 else
0045     display('list appears to be neither genes nor reactions:  aborting');
0046     return;
0047 end
0048 
0049 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0050 %%% PARAMETERS - set parameters here %%%%%%
0051 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0052 if nargin < 5
0053     MaxKnockOuts = 10;
0054 else
0055     MaxKnockOuts = MaxKOs;
0056 end
0057 mutationRate = 1/ngenes; % paper: a mutation rate of 1/(genome size) was found to be optimal for both representations.
0058 crossovermutationRate = mutationRate*.2;  % the rate of mutation after a crossover.  This value should probably be fairly low.  It is only there to ensure that not every member of the population ends up with the same genotype.
0059 CrossoverFraction = .80;    % Percentage of offspring created by crossing over (as opposed to mutation). 0.7 - 0.8 were found to generate the highest mean, but this can be adjusted.
0060 PopulationSize = [125 125 125 125]; % paper: it was found that an increase beyond 125 individuals did not improve the results significantly.
0061 Generations = 10000;    % paper: 5000.  maximum number of generations to perform
0062 TimeLimit =  3600*24*2;  % global time limit in seconds
0063 StallTimeLimit = 3600*24*1;   % Stall time limit (terminate after this much time of not finding an improvement in fitness)
0064 StallGenLimit =  10000;       % terminate after this many generations of not finding an improvement
0065 PlotFcns =  {@gaplotscores, @gaplotbestf, @gaplotscorediversity, @gaplotstopping, @gaplotmutationdiversity}; % what to plot.
0066 crossfun = @(a,b,c,d,e,f) crossoverCustom(a,b,c,d,e,f,crossovermutationRate);
0067 MigrationFraction = .1;   % how many individuals migrate (.1 * 125 ~ 12 individuals).
0068 MigrationInterval = 100;  % how often individuals migrate from one population to another.
0069 
0070 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0071 %%% END PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%
0072 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0073 
0074 InitialPopulation = [];
0075 if nargin > 5
0076     InitialPopulation = double(population);
0077 end
0078 options = gaoptimset(                                   ...
0079     'PopulationType', 'bitstring',                          ...
0080     'CreationFcn', @lowmutationcreation,                    ...
0081     'MutationFcn', {@mutationUniformEqual, mutationRate},   ...
0082     'PopulationSize', PopulationSize,                       ...
0083     'StallTimeLimit', StallTimeLimit,                       ...
0084     'TimeLimit', TimeLimit,                                 ...
0085     'PlotFcns', PlotFcns,                                   ...
0086     'InitialPopulation', InitialPopulation,                 ...
0087     'CrossoverFraction', CrossoverFraction,                 ...
0088     'CrossoverFcn', crossfun,                               ...
0089     'StallGenLimit', StallGenLimit,                         ...
0090     'Generations', Generations,                             ...
0091     'TolFun', 1e-10,                                        ...
0092     'Vectorize', 'on', ...
0093     'MigrationFraction', MigrationFraction, ...
0094     'MigrationInterval', MigrationInterval ...
0095     ... % 'SelectionFcn',  @selectionstochunif ...
0096     );
0097 
0098 %     options
0099 %     pause;
0100 
0101 %finess function call
0102 %FitnessFunction = @(x) optGeneFitness(x,model,targetRxn, generxnList, geneok);
0103 FitnessFunction = @(x) optGeneFitnessTilt(x,model,targetRxn, generxnList, geneok);
0104 
0105 gap.fitnessfcn  = FitnessFunction;
0106 gap.nvars = ngenes;
0107 gap.options = options;
0108 
0109 [x,FVAL,REASON,OUTPUT,population, scores] = ga(gap);
0110 
0111 % save the solution
0112 [optGeneSol] = GetOptGeneSol(model, targetRxn, substrateRxn, generxnList, population, x, scores, geneok); % in case of genes
0113 
0114 return;
0115 
0116 
0117 %% Creation Function
0118 % generates initial warmup with much lower number of mutations (on average
0119 % one mutation per
0120 function [Population] = lowmutationcreation(GenomeLength,FitnessFcn,options)
0121 totalPopulation = sum(options.PopulationSize);
0122 initPopProvided = size(options.InitialPopulation,1);
0123 individualsToCreate = totalPopulation - initPopProvided;
0124 
0125 % Initialize Population to be created
0126 Population = true(totalPopulation,GenomeLength);
0127 % Use initial population provided already
0128 if initPopProvided > 0
0129     Population(1:initPopProvided,:) = options.InitialPopulation;
0130 end
0131 % Create remaining population
0132 Population(initPopProvided+1:end,:) = logical(1/GenomeLength > rand(individualsToCreate,GenomeLength));
0133 return;
0134 
0135 %% Mutation Function
0136 %%%%%%%%%%%%%%%%%%%%%%%%%%%
0137 %%% mutation function %%%%%
0138 %%%%%%%%%%%%%%%%%%%%%%%%%%%
0139 function mutationChildren = mutationUniformEqual(parents,options,GenomeLength,FitnessFcn,state,thisScore,thisPopulation,mutationRate)
0140 global MaxKnockOuts
0141 if(nargin < 8)
0142     mutationRate = 0.01; % default mutation rate
0143 end
0144 mutationChildren = zeros(length(parents),GenomeLength);
0145 for i=1:length(parents)
0146     child = thisPopulation(parents(i),:);
0147     kos = sum(child);
0148     mutationPoints = find(rand(1,length(child)) < mutationRate);
0149     child(mutationPoints) = ~child(mutationPoints);
0150 
0151     if MaxKnockOuts > 0
0152         while(sum(child(:))> MaxKnockOuts)
0153             ind2 = find(child);
0154             removeindex = ind2(randint(1,1,length(ind2))+1);
0155             child(removeindex) = 0;
0156         end
0157     end
0158 
0159     % with 50% chance, you will have fewer knockouts after mutation
0160     % than before.  This is to stop aquiring so many mutations.
0161     if rand > .5 && kos > 1
0162         while(sum(child(:))>= kos)
0163             ind2 = find(child);
0164             removeindex = ind2(randint(1,1,length(ind2))+1);
0165             child(removeindex) = 0;
0166         end
0167     end
0168     mutationChildren(i,:) = child;
0169 end
0170 return;
0171 
0172 %% Crossover Function
0173 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0174 %%%% crossover function %%%%%%%%
0175 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0176 function xoverKids  = crossoverCustom(parents,options,GenomeLength,FitnessFcn,unused,thisPopulation, mutationRate)
0177 nKids = length(parents)/2;
0178 % Extract information about linear constraints
0179 % Allocate space for the kids
0180 xoverKids = zeros(nKids,GenomeLength);
0181 
0182 global MaxKnockOuts;
0183 
0184 % To move through the parents twice as fast as thekids are
0185 % being produced, a separate index for the parents is needed
0186 index = 1;
0187 
0188 % for each kid...
0189 for i=1:nKids
0190     % get parents
0191     r1 = parents(index);
0192     index = index + 1;
0193     r2 = parents(index);
0194     index = index + 1;
0195 
0196     % Randomly select half of the genes from each parent
0197     % This loop may seem like brute force, but it is twice as fast as the
0198     % vectorized version, because it does no allocation.
0199     for j = 1:GenomeLength
0200         if(rand > 0.5)
0201             xoverKids(i,j) = thisPopulation(r1,j);
0202         else
0203             xoverKids(i,j) = thisPopulation(r2,j);
0204         end
0205     end
0206     if MaxKnockOuts>0
0207         while(sum(xoverKids(i,:))> MaxKnockOuts)
0208             ind2 = find(xoverKids(i,:));
0209             removeindex = ind2(randint(1,1,length(ind2))+1);
0210             xoverKids(i,removeindex) = 0;
0211         end
0212     end
0213 end
0214 % also apply mutations to crossover kids...
0215 xoverKids = mutationUniformEqual(1:size(xoverKids,1) ,[],GenomeLength,[],[],[],xoverKids,mutationRate);
0216 return;
0217 
0218 
0219 function state = gaplotmutationdiversity(options,state,flag,p1)
0220 %GAPLOTSCOREDIVERSITY Plots a histogram of this generation's scores.
0221 %   STATE = GAPLOTSCOREDIVERSITY(OPTIONS,STATE,FLAG) plots a histogram of current
0222 %   generation's scores.
0223 %
0224 %   Example:
0225 %   Create an options structure that uses GAPLOTSCOREDIVERSITY
0226 %   as the plot function
0227 %     options = gaoptimset('PlotFcns',@gaplotscorediversity);
0228 
0229 %   Copyright 2003-2007 The MathWorks, Inc.
0230 %   $Revision: 1.6.4.3 $  $Date: 2007/05/23 18:49:53 $
0231 global MaxKnockOuts
0232 if nargin < 4
0233     p1 = 10;
0234 end
0235 
0236 p1 = MaxKnockOuts+1;
0237 p1 = [0:(MaxKnockOuts)];
0238 switch flag
0239     case 'init'
0240         title('Mutation Histogram','interp','none')
0241         xlabel('number of mutations');
0242         ylabel('Number of individuals');
0243     case 'iter'
0244         % Check if Rank is a field and there are more than one objectives, then plot for Rank == 1
0245         if size(state.Score,2) > 1 && isfield(state,'Rank') 
0246             index = (state.Rank == 1);
0247             % When there is one point hist will treat it like a vector
0248             % instead of matrix; we need to add one more duplicate row
0249             if nnz(index) > 1
0250                 set(gca,'ylimmode','auto');
0251                 hist(sum(state.Population(index,:)),p1);
0252             else
0253                 set(gca,'ylim',[0 1]);
0254                 hist([sum(state.Population(index,:)); sum(state.Population(index,:))],p1);
0255             end
0256             % Legend for each function <min max> values on the Pareto front
0257             nObj = size(state.Score,2);
0258             fminval = min(state.Score(index,:),[],1);
0259             fmaxval = max(state.Score(index,:),[],1);
0260             legendText = cell(1,nObj);
0261             for i = 1:nObj
0262                legendText{i} = ['fun',num2str(i),' [',sprintf('%g  ',fminval(i)), ...
0263                    sprintf('%g',fmaxval(i)),']'];
0264             end
0265             legend(legendText);
0266         else % else plot all score
0267             hist(sum(state.Population,2),p1);
0268         end
0269 end
0270

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003