analyzeGCdesign

PURPOSE ^

analyzeGCdesign Analyze results with replacement knockouts

SYNOPSIS ^

function [improvedRxns,intermediateSlns] = analyzeGCdesign(modelRed,selectedRxns,target,deletions,maxKOs,objFunction,delPenalty,intermediateSlns)

DESCRIPTION ^

analyzeGCdesign Analyze results with replacement knockouts
should get closer to local maxima.  must have num KOs > 1

 [improvedRxns,intermediateSlns] = analyzeGCdesign(modelRed,selectedRxns,target,deletions,maxKOs,objFunction,delPenalty,intermediateSlns)

INPUTS
 modelRed          reduced model
 selectedRxns      selected reaction list from the reduced model
 target            exchange rxn to optimize
 deletions         initial set of KO rxns (must have at least 1 rxn)

OPTIONAL INPUTS
 maxKOs            maximum number of rxn KOs to allow (Default = 10)
 objFunction       pick an objective function to use (Default = 1):
                      1: obj = maxRate (yield)
                      2: obj = growth*maxRate (SSP)
                      3: obj = maxRate*(delPenalty^numDels) (yield with KO penalty)
                      4: obj = growth*maxRate*(delPenalty^numDels)  (SSP with KO penalty)
                      5: obj = maxRate*(slope^(-1))  (GC_yield)
                      6: obj = growth*maxRate*(slope^(-1))  (GC_SSP)
                      7: obj = maxRate*(delPenalty^numDels)*(slope^(-1)) (GC_yield with KO penalty)
                      8: obj = growth*maxRate*(delPenalty^numDels)*(slope^(-1))  (GC_SSP with KO penalty)
 delPenalty        penalty on extra rxn deletions (Default = .99)
 intermediateSlns  Previous set of solutions (Default = deletions)

OUTPUTS
 improvedRxns      the KO rxns for an improved strain
 intermediateSlns  all the sets of best KO rxns that are picked before the
                   final set is reached

 Jeff Orth  7/25/07
 Richard Que 1/19/10       Replaced try/catch blocks

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [improvedRxns,intermediateSlns] = analyzeGCdesign(modelRed,selectedRxns,target,deletions,maxKOs,objFunction,delPenalty,intermediateSlns)
0002 %analyzeGCdesign Analyze results with replacement knockouts
0003 %should get closer to local maxima.  must have num KOs > 1
0004 %
0005 % [improvedRxns,intermediateSlns] = analyzeGCdesign(modelRed,selectedRxns,target,deletions,maxKOs,objFunction,delPenalty,intermediateSlns)
0006 %
0007 %INPUTS
0008 % modelRed          reduced model
0009 % selectedRxns      selected reaction list from the reduced model
0010 % target            exchange rxn to optimize
0011 % deletions         initial set of KO rxns (must have at least 1 rxn)
0012 %
0013 %OPTIONAL INPUTS
0014 % maxKOs            maximum number of rxn KOs to allow (Default = 10)
0015 % objFunction       pick an objective function to use (Default = 1):
0016 %                      1: obj = maxRate (yield)
0017 %                      2: obj = growth*maxRate (SSP)
0018 %                      3: obj = maxRate*(delPenalty^numDels) (yield with KO penalty)
0019 %                      4: obj = growth*maxRate*(delPenalty^numDels)  (SSP with KO penalty)
0020 %                      5: obj = maxRate*(slope^(-1))  (GC_yield)
0021 %                      6: obj = growth*maxRate*(slope^(-1))  (GC_SSP)
0022 %                      7: obj = maxRate*(delPenalty^numDels)*(slope^(-1)) (GC_yield with KO penalty)
0023 %                      8: obj = growth*maxRate*(delPenalty^numDels)*(slope^(-1))  (GC_SSP with KO penalty)
0024 % delPenalty        penalty on extra rxn deletions (Default = .99)
0025 % intermediateSlns  Previous set of solutions (Default = deletions)
0026 %
0027 %OUTPUTS
0028 % improvedRxns      the KO rxns for an improved strain
0029 % intermediateSlns  all the sets of best KO rxns that are picked before the
0030 %                   final set is reached
0031 %
0032 % Jeff Orth  7/25/07
0033 % Richard Que 1/19/10       Replaced try/catch blocks
0034 
0035 if (nargin < 5)
0036     maxKOs = 10;
0037 end
0038 if (nargin < 6)
0039     objFunction = 1;
0040 end
0041 if (nargin < 7)
0042     delPenalty = .99;
0043 end
0044 if (nargin < 8)
0045     intermediateSlns = {deletions};
0046 end
0047 
0048 %set the objective function
0049 switch objFunction
0050     case 1
0051         objectiveFunction = 'maxRate';
0052         hasSlope = false;
0053     case 2
0054         objectiveFunction = 'growth*maxRate';
0055         hasSlope = false;
0056     case 3
0057         objectiveFunction = 'maxRate*(delPenalty^numDels)';
0058         hasSlope = false;
0059     case 4
0060         objectiveFunction = 'growth*maxRate*(delPenalty^numDels)';
0061         hasSlope = false;
0062     case 5
0063         objectiveFunction = 'maxRate*(slope^(-1))';
0064         hasSlope = true;
0065     case 6
0066         objectiveFunction = 'growth*maxRate*(slope^(-1))';
0067         hasSlope = true;
0068     case 7
0069         objectiveFunction = 'maxRate*(delPenalty^numDels)*(slope^(-1))';
0070         hasSlope = true;
0071     case 8
0072         objectiveFunction = 'growth*maxRate*(delPenalty^numDels)*(slope^(-1))';
0073         hasSlope = true;
0074 end
0075 
0076 if isempty(deletions)
0077     error('no knockout reactions defined')
0078 end
0079 
0080 delArraySize = size(deletions); %make sure deletions list is horizontal
0081 if delArraySize(1) > 1
0082     rxns = deletions';
0083 else
0084     rxns = deletions;
0085 end
0086 
0087 BOF = modelRed.rxns(modelRed.c==1); %get biomass objective function
0088 
0089 
0090 modelKO = changeRxnBounds(modelRed,rxns,0,'b');
0091 FBAsol1 = optimizeCbModel(modelKO,'max',true); %find max growth rate of strain
0092 if FBAsol1.stat>0
0093     modelKOfixed = changeRxnBounds(modelKO,BOF,FBAsol1.f-1e-6,'l'); %fix the growth rate
0094     modelKOfixed = changeObjective(modelKOfixed,target); %set target as the objective
0095     FBAsol2 = optimizeCbModel(modelKOfixed,'min',true); %find minimum target rate at this growth rate
0096     growth = FBAsol1.f;
0097     maxRate = FBAsol2.f;
0098     numDels = length(rxns);
0099     
0100     if hasSlope %only calculate these if the obj function includes slope
0101         modelTarget = changeObjective(modelKO,target); %set target as the objective
0102         FBAsol4 = optimizeCbModel(modelTarget,'min',true); %find min production rate
0103         modelTargetFixed = changeRxnBounds(modelKO,target,FBAsol4.f,'b'); %fix production to minimum
0104         FBAsol5 = optimizeCbModel(modelTargetFixed,'max',true); %find max growth at min production
0105         minProdRate = FBAsol4.f;
0106         maxGrowthMinRate = FBAsol5.f;
0107         
0108         if growth ~= maxGrowthMinRate
0109             slope = (maxRate-minProdRate)/(growth-maxGrowthMinRate);
0110         else
0111             slope = 1; %don't consider slope if div by 0
0112         end
0113     end
0114     
0115     objective = eval(objectiveFunction);
0116     
0117     bestObjective = objective
0118     bestRxns = rxns;
0119     % if the initial reactions are lethal
0120 else
0121     bestObjective = 0
0122     bestRxns = rxns;
0123 end
0124 
0125 % loop through each KO rxn and replace with every rxn from selectedRxns to
0126 % search for a possible improvement
0127 h = waitbar(0,'improving knockout design');
0128 for i = 1:length(rxns)+1
0129     bestObjective2 = bestObjective;
0130     bestRxns2 = bestRxns;
0131     for j = 1:length(selectedRxns)+1
0132         waitbar((j+(i-1)*length(selectedRxns))/((length(rxns)+1)*(length(selectedRxns)+1)),h);
0133         newRxns = rxns;
0134         if (i==length(rxns)+1)&&(j==length(selectedRxns)+1)
0135             %don't do anything at the very end
0136         elseif j ~= length(selectedRxns)+1
0137             newRxns{i} = selectedRxns{j}; %replace rxn with different one
0138         elseif i == 1 %or else remove one of the rxns
0139             newRxns = rxns(2:length(rxns));
0140         elseif i == length(rxns)
0141             newRxns = rxns(1:length(rxns)-1);
0142         else
0143             newRxns = cat(2,rxns(1:i-1),rxns(i+1:length(rxns)));
0144         end
0145         
0146         if length(newRxns) <= maxKOs %limit the total number of knockouts
0147             
0148             modelKO = changeRxnBounds(modelRed,newRxns,0,'b');
0149             FBAsol1 = optimizeCbModel(modelKO,'max',true); %find max growth rate of strain
0150             if FBAsol1.stat>0
0151                 modelKOfixed = changeRxnBounds(modelKO,BOF,FBAsol1.f-1e-6,'l'); %fix the growth rate
0152                 modelKOfixed = changeObjective(modelKOfixed,target); %set target as the objective
0153                 FBAsol2 = optimizeCbModel(modelKOfixed,'min',true); %find minimum target rate at this growth rate
0154                 FBAsol3 = optimizeCbModel(modelKOfixed,'max',true); %find maximum target rate at this growth rate
0155                 growth = FBAsol1.f;
0156                 maxRate = FBAsol2.f;
0157                 numDels = length(newRxns);
0158                 
0159                 if hasSlope %only calculate these if the obj function includes slope
0160                     modelTarget = changeObjective(modelKO,target); %set target as the objective
0161                     FBAsol4 = optimizeCbModel(modelTarget,'min',true); %find min production rate
0162                     modelTargetFixed = changeRxnBounds(modelKO,target,FBAsol4.f,'b'); %fix production to minimum
0163                     FBAsol5 = optimizeCbModel(modelTargetFixed,'max',true); %find max growth at min production
0164                     minProdRate = FBAsol4.f;
0165                     maxGrowthMinRate = FBAsol5.f;
0166                     
0167                     if growth ~= maxGrowthMinRate
0168                         slope = (maxRate-minProdRate)/(growth-maxGrowthMinRate);
0169                     else
0170                         slope = 1; %don't consider slope if div by 0
0171                     end
0172                 end
0173                 
0174                 newObjective = eval(objectiveFunction);
0175                 
0176                 %see if objective is increased by this new gene
0177                 if newObjective > bestObjective2
0178                     bestObjective2 = newObjective
0179                     bestRxns2 = newRxns
0180                     intermediateSlns{length(intermediateSlns)+1} = bestRxns2; %add new intermediateSln to the list
0181                 end
0182             end
0183         end
0184     end
0185     
0186     if bestObjective2 > bestObjective
0187         bestObjective = bestObjective2
0188         bestRxns = bestRxns2
0189     end
0190 end
0191 close(h);
0192 
0193 bestObjective
0194 bestRxns
0195 
0196 % recursively call analyzeGCdesign again until no improvement is found
0197 if length(bestRxns) ~= length(rxns)
0198     [bestRxns,intermediateSlns] = analyzeGCdesign(modelRed,selectedRxns,target,bestRxns,maxKOs,objFunction,delPenalty,intermediateSlns);
0199 elseif length(find(strcmp(bestRxns,rxns)))~=length(rxns)
0200     [bestRxns,intermediateSlns] = analyzeGCdesign(modelRed,selectedRxns,target,bestRxns,maxKOs,objFunction,delPenalty,intermediateSlns);
0201 end
0202 
0203 % print final results
0204 improvedRxns = sort(bestRxns)
0205 
0206 
0207 
0208

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