optGeneFitness

PURPOSE ^

optGeneFitness GeneOptFitness the fitness function

SYNOPSIS ^

function [val] = optGeneFitness(rxn_vector_matrix, model, targetRxn, rxnListInput, isGeneList)

DESCRIPTION ^

optGeneFitness GeneOptFitness the fitness function

 [val] = optGeneFitness(rxn_vector_matrix, model, targetRxn, rxnListInput, isGeneList)

INPUTS
 rxn_vector_matrix
 model
 targetRxn
 rxnListInput
 isGeneList

OUTPUT
 val

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [val] = optGeneFitness(rxn_vector_matrix, model, targetRxn, rxnListInput, isGeneList)
0002 %optGeneFitness GeneOptFitness the fitness function
0003 %
0004 % [val] = optGeneFitness(rxn_vector_matrix, model, targetRxn, rxnListInput, isGeneList)
0005 %
0006 %INPUTS
0007 % rxn_vector_matrix
0008 % model
0009 % targetRxn
0010 % rxnListInput
0011 % isGeneList
0012 %
0013 %OUTPUT
0014 % val
0015 %
0016 %
0017 global MaxKnockOuts
0018 %size(rxn_vector_matrix)
0019 
0020 popsize = size(rxn_vector_matrix,1);
0021 val = zeros(1,popsize);
0022 
0023 for i = 1:popsize
0024     rxn_vector = rxn_vector_matrix(i,:);
0025     rxnList = rxnListInput(logical(rxn_vector));
0026 
0027     
0028     %see if we've done this before
0029     val_temp = memoize(rxn_vector);
0030     if ~ isempty(val_temp)
0031         val(i) = val_temp;
0032         continue;
0033     end
0034    
0035     % check to see if mutations is above the max number allowed
0036     nummutations = sum(rxn_vector);
0037     if nummutations > MaxKnockOuts
0038         continue;
0039     end
0040     
0041     % generate knockout.
0042     if isGeneList
0043         modelKO = deleteModelGenes(model, rxnList);
0044     else % is reaction list
0045         [isValidRxn,removeInd] = ismember(rxnList,model.rxns);
0046         removeInd = removeInd(isValidRxn);
0047         modelKO = model;
0048         modelKO.ub(removeInd) = 0;
0049         modelKO.lb(removeInd) = 0;
0050     end
0051     
0052     % find growthrate;
0053 %     slnKO = optimizeCbModel(modelKO);
0054 %     growthrate1 = slnKO.f; %max growth rate.
0055     if exist('LPBasis', 'var')
0056         modelKO.LPBasis = LPBasis;
0057     end
0058 
0059     [slnKO, LPOUT] = solveCobraLPCPLEX(modelKO, 0,1);
0060     LPBasis = LPOUT.LPBasis;
0061     growthrate = slnKO.obj;
0062     
0063     
0064     % check to ensure that GR is above a certain value
0065     if growthrate < .10
0066         continue;
0067     end
0068     
0069 %    display('second optimization');
0070     % find the lowesest possible production rate (a hopefully high number)
0071     % at the max growth rate minus some set factor gamma (a growth rate slightly
0072     % smaller than the max). A positive value will eliminate solutions where the
0073     % production envelope has a vertical line at the max GR, a "non-unique"
0074     % solution. Set value to zero if "non-unique" solutions are not an issue.
0075     gamma = 0.01; % proportional to Grwoth Rate (hr-1), a value around 0.5 max.
0076 
0077     %find indicies of important vectors
0078     indBOF = find(modelKO.c);
0079     indTar = findRxnIDs(modelKO, targetRxn);
0080     % generate a model with a fixed max KO growth rate
0081     modelKOsetGR = modelKO;
0082     modelKOsetGR.lb(indBOF) = growthrate - gamma; % this growth rate is required as lb.
0083     modelKOsetGR.c = zeros(size(modelKO.c));
0084     modelKOsetGR.c(indTar) = -1; % minimize for this variable b/c we want to look at the very minimum production.
0085 
0086     % find the minimum production rate for the targeted reaction.
0087 
0088 %     slnKOsetGR = optimizeCbModel(modelKOsetGR);
0089 %     minProdAtSetGR1 = -slnKOsetGR.f;  % This should be a negative value b/c of the minimization setup, so -1 is necessary.
0090 
0091     if exist('LPBasis2', 'var')
0092         modelKOsetGR.LPBasis = LPBasis2;
0093     end
0094 
0095     [slnKOsetGR, LPOUT2] = solveCobraLPCPLEX(modelKOsetGR, 0,1);
0096     LPBasis2 = LPOUT2.LPBasis;
0097     minProdAtSetGR = -slnKOsetGR.obj;
0098 
0099     
0100     % objective function for optGene algorithm = val (needs to be a negative value, since it is
0101     % a minimization)
0102     val(i) = -minProdAtSetGR; 
0103     % penalty for a greater number of mutations
0104     % val(i) = -minProdAtSetGR * (.98^nummutations);
0105     % select best substrate-specific productivity
0106     % val(i) = -minProdAtSetGR * (.98^nummutations) * growthrate;
0107 
0108     % check to prevent very small values from being considerered improvments
0109     if val(i) > -1e-3
0110         val(i) = 0;
0111     end
0112 
0113     memoize(rxn_vector, val(i));
0114 end
0115 
0116 return;
0117 
0118 
0119 
0120 %% Memoize
0121 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0122 %%%%% MEMOIZE %%%%%%%%%%%%%%%%%%%%%
0123 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0124 % internal function to speed things up.
0125 
0126 function [value] =  memoize(gene_vector, value)
0127 global HTABLE
0128 hashkey = num2str(gene_vector);
0129 hashkey = strrep(hashkey,' ',''); % cut out white space from string (more space efficient).
0130 
0131 if nargin == 1
0132     value = HTABLE.get(hashkey);
0133     return;
0134 else
0135     if HTABLE.size() > 10000
0136         HTABLE = java.util.Hashtable;  %reset the hashtable if more than 10,000 entries.
0137     end
0138     HTABLE.put(hashkey, value);
0139     value = [];
0140     return;
0141 end
0142 return

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