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
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