optGeneFitnessTilt

PURPOSE ^

optGeneFitnessTilt GeneOptFitness the fitness function

SYNOPSIS ^

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

DESCRIPTION ^

optGeneFitnessTilt GeneOptFitness the fitness function

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

INPUTS
 rxn_vector_matrix     
 model                 
 targetRxn             
 rxnListInput          
 isGeneList            

OUTPUT
 val                   fitness value

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [val] = optGeneFitnessTilt(rxn_vector_matrix, model, targetRxn, rxnListInput, isGeneList)
0002 %optGeneFitnessTilt GeneOptFitness the fitness function
0003 %
0004 % [val] = optGeneFitnessTilt(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                   fitness value
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  % augment BOF (tilt)
0053  [modelKO] = augmentBOF(modelKO, targetRxn, .001);
0054  
0055  % find growthrate
0056     if exist('LPBasis', 'var')
0057         modelKO.LPBasis = LPBasis;
0058     end
0059     
0060     [slnKO, LPOUT] = solveCobraLPCPLEX(modelKO, 0,1);
0061     LPBasis = LPOUT.LPBasis;
0062     growthrate = slnKO.obj;
0063     [tmp,tar_loc] = ismember(targetRxn,modelKO.rxns);
0064     minProdAtSetGR = slnKO.full(tar_loc);
0065     
0066     % check to ensure that GR is above a certain value
0067     if growthrate < .10
0068         continue;
0069     end
0070     
0071 % %    display('second optimization');
0072 %     % find the lowesest possible production rate (a hopefully high number)
0073 %     % at the max growth rate minus some set factor gamma (a growth rate slightly
0074 %     % smaller than the max). A positive value will eliminate solutions where the
0075 %     % production envelope has a vertical line at the max GR, a "non-unique"
0076 %     % solution. Set value to zero if "non-unique" solutions are not an issue.
0077 %     gamma = 0.01; % proportional to Grwoth Rate (hr-1), a value around 0.5 max.
0078 %
0079 %     %find indicies of important vectors
0080 %     indBOF = find(modelKO.c);
0081 %     indTar = findRxnIDs(modelKO, targetRxn);
0082 %     % generate a model with a fixed max KO growth rate
0083 %     modelKOsetGR = modelKO;
0084 %     modelKOsetGR.lb(indBOF) = growthrate - gamma; % this growth rate is required as lb.
0085 %     modelKOsetGR.c = zeros(size(modelKO.c));
0086 %     modelKOsetGR.c(indTar) = -1; % minimize for this variable b/c we want to look at the very minimum production.
0087 %
0088 %     % find the minimum production rate for the targeted reaction.
0089 %
0090 % %     slnKOsetGR = optimizeCbModel(modelKOsetGR);
0091 % %     minProdAtSetGR1 = -slnKOsetGR.f;  % This should be a negative value b/c of the minimization setup, so -1 is necessary.
0092 %
0093 %     if exist('LPBasis2', 'var')
0094 %         modelKOsetGR.LPBasis = LPBasis2;
0095 %     end
0096 %
0097 %     [slnKOsetGR, LPOUT2] = solveCobraLPCPLEX(modelKOsetGR, 0,1);
0098 %     LPBasis2 = LPOUT2.LPBasis;
0099 %     minProdAtSetGR = -slnKOsetGR.obj;
0100 
0101     
0102     % objective function for optGene algorithm = val (needs to be a negative value, since it is
0103     % a minimization)
0104     val(i) = -minProdAtSetGR; 
0105     % penalty for a greater number of mutations
0106     
0107     %val(i) = -minProdAtSetGR * (.98^nummutations);
0108     
0109     % select best substrate-specific productivity
0110     % val(i) = -minProdAtSetGR * (.98^nummutations) * growthrate;
0111 
0112     % check to prevent very small values from being considerered improvments
0113     if val(i) > -1e-3
0114         val(i) = 0;
0115     end
0116 
0117     memoize(rxn_vector, val(i));
0118 end
0119 
0120 return;
0121 
0122 
0123 
0124 %% Memoize
0125 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0126 %%%%% MEMOIZE %%%%%%%%%%%%%%%%%%%%%
0127 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0128 % internal function to speed things up.
0129 
0130 function [value] =  memoize(gene_vector, value)
0131 global HTABLE
0132 hashkey = num2str(gene_vector);
0133 hashkey = strrep(hashkey,' ',''); % cut out white space from string (more space efficient).
0134 
0135 if nargin == 1
0136     value = HTABLE.get(hashkey);
0137     return;
0138 else
0139     if HTABLE.size() > 50000
0140         HTABLE = java.util.Hashtable;  %reset the hashtable if more than 50,000 entries.
0141     end
0142     HTABLE.put(hashkey, value);
0143     value = [];
0144     return;
0145 end
0146 return

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