GDLS

PURPOSE ^

GDLS (Genetic Design through Local Search) attempts to find genetic

SYNOPSIS ^

function [gdlsSolution, bilevelMILPProblem, gdlsSolutionStructs] = GDLS(model, targetRxns, varargin)

DESCRIPTION ^

GDLS (Genetic Design through Local Search) attempts to find genetic
 designs with greater in silico production of desired metabolites.

 [gdlsSolution, bilevelMILPProblem, gdlsSolutionStructs] = GDLS(model, varargin)

INPUTS
 model             Cobra model structure
 targetRxn         Reaction(s) to be maximized (Cell array of strings)

OPTIONAL INPUTS
 varargin          parameters entered using either a structure or list of
                   parameter, parameter value
   List of optional parameters
   'nbhdsz'            Neighborhood size (default: 1)
   'M'                 Number of search paths (default: 1)
   'maxKO'             Maximum number of knockouts (default: 50)
   'koCost'            Cost for knocking out a reaction, gene set, or gene
                       A different cost can be set for each knockout.
                       (default: 1 for each knockout)
   'selectedRxns'      List of reactions/geneSets that can be knocked out
   'koType'            What to knockout: reactions, gene sets, or genes
                       {('rxns'), 'geneSets', 'genes'}
   'iterationLimit'    Maximum number of iterations (default: 70)
   'timeLimit'         Maximum run time in seconds (default: 252000)
   'minGrowth'         Minimum growth rate

OUTPUTS
 gdlsSolution          GDLS solution structure (similar to OptKnock sol struct)
 bilevelMILPProblem    Problem structure used in computation

 Adapted from Desmond S Lun's gdls scripts.
 Richard Que (1/28/2010)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [gdlsSolution, bilevelMILPProblem, gdlsSolutionStructs] = GDLS(model, targetRxns, varargin)
0002 %GDLS (Genetic Design through Local Search) attempts to find genetic
0003 % designs with greater in silico production of desired metabolites.
0004 %
0005 % [gdlsSolution, bilevelMILPProblem, gdlsSolutionStructs] = GDLS(model, varargin)
0006 %
0007 %INPUTS
0008 % model             Cobra model structure
0009 % targetRxn         Reaction(s) to be maximized (Cell array of strings)
0010 %
0011 %OPTIONAL INPUTS
0012 % varargin          parameters entered using either a structure or list of
0013 %                   parameter, parameter value
0014 %   List of optional parameters
0015 %   'nbhdsz'            Neighborhood size (default: 1)
0016 %   'M'                 Number of search paths (default: 1)
0017 %   'maxKO'             Maximum number of knockouts (default: 50)
0018 %   'koCost'            Cost for knocking out a reaction, gene set, or gene
0019 %                       A different cost can be set for each knockout.
0020 %                       (default: 1 for each knockout)
0021 %   'selectedRxns'      List of reactions/geneSets that can be knocked out
0022 %   'koType'            What to knockout: reactions, gene sets, or genes
0023 %                       {('rxns'), 'geneSets', 'genes'}
0024 %   'iterationLimit'    Maximum number of iterations (default: 70)
0025 %   'timeLimit'         Maximum run time in seconds (default: 252000)
0026 %   'minGrowth'         Minimum growth rate
0027 %
0028 %OUTPUTS
0029 % gdlsSolution          GDLS solution structure (similar to OptKnock sol struct)
0030 % bilevelMILPProblem    Problem structure used in computation
0031 %
0032 % Adapted from Desmond S Lun's gdls scripts.
0033 % Richard Que (1/28/2010)
0034 
0035 MAXFLUX = 1000;
0036 MAXDUAL = 1000;
0037 EPS = 1e-4;
0038 gdlsSolutionStructs = [];
0039 
0040 if nargin < 2
0041     error('Model and target reaction(s) must be specified')
0042 elseif mod((nargin-2),2)==0 %manual entry
0043     options.targetRxns = targetRxns;
0044     for i=1:2:(nargin-2)
0045         if ismember(varargin{i},{'nbhdsz','M','maxKO','selectedRxns','koType','koCost','minGrowth', ...
0046                                  'timeLimit','iterationLimit'})
0047             options.(varargin{i}) = varargin{i+1};
0048         else
0049             display(['Unknown option ' varargin{i}]);
0050         end
0051     end
0052 elseif isstruct(targetRxns) %options structure
0053     options = varargin{1};
0054 else
0055     error('Invalid number of entries')
0056 end
0057 
0058 %Default Values
0059 if ~isfield(options,'koType'), options.koType = 'rxns'; end
0060 if ~isfield(options,'nbhdsz'), options.nbhdsz=1; end
0061 if ~isfield(options,'M'), options.M=1; end
0062 if ~isfield(options,'maxKO'), options.maxKO=50; end
0063 if ~isfield(options,'iterationLimit'), options.iterationLimit=70; end
0064 if ~isfield(options,'timeLimit'), options.timeLimit=252000 ; end
0065 if isfield(options,'targetRxns')
0066     selTargetRxns = logical(ismember(model.rxns,options.targetRxns));
0067     if ~any(selTargetRxns)
0068         error([options.targetRxns ' not found. Double check spelling.']);
0069     end
0070 else
0071     error('No target reaction specified')
0072 end
0073 if isfield(options,'selectedRxns')
0074     selSelectedRxns = logical(ismember(model.rxns,options.selectedRxns));
0075 else
0076     selSelectedRxns = true(length(model.rxns),1);
0077 end
0078 
0079 
0080 switch lower(options.koType)
0081     case 'rxns'
0082         %% Generate selection reaction matrix
0083         model.selRxnMatrix = selMatrix(selSelectedRxns)';
0084         possibleKOList = model.rxns(selSelectedRxns);
0085         
0086     case 'genesets'
0087         %% Generate reaction gene set mapping matrix
0088         %remove biomass reaction from grRules and generate unique gene set list
0089         possibleKOList = unique(model.grRules(selSelectedRxns));
0090         if isempty(possibleKOList{1}), possibleKOList = possibleKOList(2:end); end
0091         for i = 1:length(possibleKOList)
0092             model.selRxnMatrix(:,i) =  double(strcmp(possibleKOList{i},model.grRules));
0093         end
0094         
0095     case 'genes'
0096         %% Use rxnGeneMat as selRxnMatrix
0097         model.selRxnMatrix = model.rxnGeneMat;
0098         possibleKOList = model.genes;
0099     otherwise
0100         error('Unrecognized KO type')
0101 end
0102 
0103 %% Generate koCost if not present
0104 if ~isfield(options,'koCost')
0105     if isfield(model,'koCost')
0106         if length(model.koCost) == 1
0107             options.koCost = ones(length(possibleKOList,1)) * model.koCost;
0108         else
0109             options.koCost = model.koCost;
0110         end
0111     else
0112         options.koCost = ones(length(possibleKOList),1);
0113     end
0114 elseif length(model.koCost) == 1
0115     options.koCost = ones(length(possibleKOList,1)) * model.koCost;
0116 else
0117     options.koCost = model.koCost;
0118 end
0119 
0120 %index exchange reactions
0121 [selExc] = findExcRxns(model,true,true);
0122 
0123 %% Setup model
0124 model.ub(isinf(model.ub)) = MAXFLUX;
0125 model.ub(model.ub>MAXFLUX) = MAXFLUX;
0126 model.lb(isinf(model.lb)) = -MAXFLUX;
0127 model.lb(model.lb<-MAXFLUX) = -MAXFLUX;
0128 model.rxnsPresent = ones(length(model.rxns),1);
0129 
0130 %% Create Baseline
0131 solution0 = fluxBalance(model,selTargetRxns);
0132 
0133 %% Create bi-level MILP problem
0134 
0135 [nMets nRxns] = size(model.S);
0136 nInt = size(model.selRxnMatrix,2);
0137 
0138 jMu = find(~isinf(model.lb));
0139 nMu = length(jMu);
0140 jNu = find(~isinf(model.ub));
0141 nNu = length(jNu);
0142 
0143 y0 = (model.selRxnMatrix' * ~model.rxnsPresent) ~= 0;
0144 y0 = repmat(y0,1,options.M);
0145 y0p = y0;
0146 iiter = 1;
0147 change = true;
0148 t1 = clock;
0149 while change
0150     change = false;
0151     
0152     y = false(nInt, 0);
0153     fmax = zeros(1, 0);
0154     for istart = 1:size(y0, 2)
0155         for irun = 1:options.M
0156             
0157             c = [selTargetRxns;
0158                 zeros(nMets, 1);
0159                 zeros(nMu, 1);
0160                 zeros(nNu, 1);
0161                 zeros(nRxns, 1);
0162                 zeros(nInt, 1)];
0163             %x= [ v                             lambda                      mu                                              nu                                              xi                              y                                           ]
0164             A = [ model.S                       sparse(nMets, nMets)        sparse(nMets, nMu)                              sparse(nMets, nNu)                              sparse(nMets, nRxns)            sparse(nMets, nInt);
0165                 sparse(nRxns, nRxns)            model.S'                    -sparse(jMu, 1:nMu, ones(nMu, 1), nRxns, nMu)   sparse(jNu, 1:nNu, ones(nNu, 1), nRxns, nNu)    speye(nRxns)                    sparse(nRxns, nInt);
0166                 speye(nRxns)                    sparse(nRxns, nMets)        sparse(nRxns, nMu)                              sparse(nRxns, nNu)                              sparse(nRxns, nRxns)            model.selRxnMatrix .* repmat(model.lb, 1, nInt);
0167                 speye(nRxns)                    sparse(nRxns, nMets)        sparse(nRxns, nMu)                              sparse(nRxns, nNu)                              sparse(nRxns, nRxns)            model.selRxnMatrix .* repmat(model.ub, 1, nInt);
0168                 sparse(nRxns, nRxns)            sparse(nRxns, nMets)        sparse(nRxns, nMu)                              sparse(nRxns, nNu)                              speye(nRxns)                    model.selRxnMatrix * MAXDUAL;
0169                 sparse(nRxns, nRxns)            sparse(nRxns, nMets)        sparse(nRxns, nMu)                              sparse(nRxns, nNu)                              speye(nRxns)                    -model.selRxnMatrix * MAXDUAL;
0170                 model.c'                        sparse(1, nMets)            model.lb(jMu)'                                  -model.ub(jNu)'                                 sparse(1, nRxns)                sparse(1, nInt);
0171                 sparse(1, nRxns)                sparse(1, nMets)            sparse(1, nMu)                                  sparse(1, nNu)                                  sparse(1, nRxns)                ((y0(:, istart) == 0) - (y0(:, istart) == 1))';
0172                 sparse(size(y, 2), nRxns)       sparse(size(y, 2), nMets)   sparse(size(y, 2), nMu)                         sparse(size(y, 2), nNu)                         sparse(size(y, 2), nRxns)       ((y == 0) - (y == 1))';
0173                 sparse(1, nRxns)                sparse(1, nMets)            sparse(1, nMu)                                  sparse(1, nNu)                                  sparse(1, nRxns)                options.koCost'; ];
0174             b = [ zeros(nMets, 1);
0175                 model.c;
0176                 model.lb;
0177                 model.ub;
0178                 zeros(nRxns, 1);
0179                 zeros(nRxns, 1);
0180                 0;
0181                 options.nbhdsz - nnz(y0(:, istart));
0182                 ones(size(y, 2), 1) - sum((y ~= 0), 1)';
0183                 options.maxKO; ];
0184        csense = char(['E' * ones(nMets, 1);
0185                 'E' * ones(nRxns, 1);
0186                 'G' * ones(nRxns, 1);
0187                 'L' * ones(nRxns, 1);
0188                 'G' * ones(nRxns, 1);
0189                 'L' * ones(nRxns, 1);
0190                 'E';
0191                 'L';
0192                 'G' * ones(size(y, 2), 1);
0193                 'L'; ]);
0194            lb = [ model.lb;
0195                 -Inf * ones(nMets, 1);
0196                 zeros(nMu, 1);
0197                 zeros(nNu, 1);
0198                 -Inf * ones(nRxns, 1);
0199                 zeros(nInt, 1) ];
0200            ub = [ model.ub;
0201                 Inf * ones(nMets, 1);
0202                 Inf * ones(nMu, 1);
0203                 Inf * ones(nNu, 1);
0204                 Inf * ones(nRxns, 1);
0205                 ones(nInt, 1) ];
0206       vartype = char(['C' * ones(nRxns, 1);
0207                 'C' * ones(nMets, 1);
0208                 'C' * ones(nMu, 1);
0209                 'C' * ones(nNu, 1);
0210                 'C' * ones(nRxns, 1);
0211                 'B' * ones(nInt, 1); ]);
0212        osense = -1; %maximize
0213             
0214             if isfield(options,'minGrowth')
0215                 A = [A; model.c' sparse(1, nMets + nMu + nNu + nRxns + nInt)];
0216                 b = [b; options.minGrowth];
0217                 csense = [csense; 'G'];
0218             end     
0219        
0220             [bilevelMILPProblem.c, bilevelMILPProblem.A,...
0221                 bilevelMILPProblem.b, bilevelMILPProblem.lb,...
0222                 bilevelMILPProblem.ub, bilevelMILPProblem.csense,...
0223                 bilevelMILPProblem.vartype, bilevelMILPProblem.osense,...
0224                 bilevelMILPProblem.x0] = ...
0225                 deal(c, A, b, lb, ub, csense, vartype, osense, []);
0226             
0227             %solve
0228             solution1 = solveCobraMILP(bilevelMILPProblem);
0229 
0230             %check solver status
0231             if solution1.stat~=1
0232                 continue; %non optimal solution
0233             end
0234             yt = solution1.full((end - nInt + 1):end) > EPS;
0235             
0236             model.rxnsPresent = ~(model.selRxnMatrix * yt);
0237             solution2 = fluxBalance(model,selTargetRxns,false);
0238             if abs(solution2.obj - solution1.obj) > EPS
0239                 continue; %inconsistent
0240             end
0241             
0242             fmax(:, end + 1) = solution1.obj;
0243             y(:, end + 1) = yt;
0244         end
0245     end
0246     
0247     if size(y, 2) == 0
0248         continue;
0249     end
0250     
0251     [fmaxsort, ifmaxsort] = sort(fmax);
0252     y = y(:, ifmaxsort);
0253     y = y(:, max([1 size(y, 2) - options.M + 1]):end);
0254     y = shrinkKnockouts(y, model, selTargetRxns);
0255     if size(y, 2) ~= size(y0, 2) || any(any(y~=y0))&&any(any(y~=y0p))
0256         y0p=y0;
0257         y0 = y;
0258         change = true;
0259     end
0260     
0261     fprintf('Iteration %d\n', iiter);
0262     fprintf('----------%s\n', char('-' * ones(1, floor(log10(iiter)) + 1)));
0263     for iend = 1:size(y0, 2)
0264         model.rxnsPresent = ~(model.selRxnMatrix * y0(:, iend));
0265         [solSynMax solSynMin solBiomass] = fluxBalance(model,selTargetRxns);
0266         fprintf('Knockout cost:   %d\n', options.koCost' * y0(:, iend));
0267         if nnz(y0) > 0
0268             fprintf('Knockouts:\n%s', sprintf('\t%s\n', possibleKOList{y0(:, iend)}));
0269         end
0270         printLabeledData(model.rxns(selExc),solSynMax.full(selExc),true);
0271         fprintf('\n');
0272         
0273         %Save Solutions
0274         gdlsSolutionStructs.(sprintf('Iteration_%d',iiter)).(sprintf('solution_%2',i)).solBiomass = solBiomass;
0275         gdlsSolutionStructs.(sprintf('Iteration_%d',iiter)).(sprintf('solution_%2',i)).solSynMin = solSynMin;
0276         gdlsSolutionStructs.(sprintf('Iteration_%d',iiter)).(sprintf('solution_%2',i)).solSynMax = solSynMax;
0277     end
0278     elapsed_time = etime(clock,t1)
0279     iiter = iiter + 1;
0280     if (elapsed_time >= options.timeLimit) || (iiter >= options.iterationLimit)
0281         break;
0282     end
0283 end
0284 
0285 %Generate Solution Structure
0286 fprintf('\nGenerating Output\n');
0287 gdlsSolution.int = y0;
0288 gdlsSolutions.KOs = cell(max(sum(y0,2)),size(y0,2));
0289 for i = 1:size(y0,2)
0290 gdlsSolution.KOs(1:nnz(y0(:,i)),i) = possibleKOList(y0(:,i));
0291 [solSynMax solSynMin solBiomass] = fluxBalance(model,selTargetRxns,false);
0292 gdlsSolution.biomass(1,i) = solBiomass.obj;
0293 gdlsSolution.minTargetProd(1,i) = solSynMin.obj;
0294 gdlsSolution.maxTargetProd(1,i) = solSynMax.obj;
0295 end
0296 
0297 end
0298 function y = shrinkKnockouts(y, model, selTargetRxns)
0299 
0300 for iycol = 1:size(y, 2)
0301     model.rxnsPresent = ~(model.selRxnMatrix * y(:, iycol));
0302     solution1 = fluxBalance(model, selTargetRxns, false);              
0303     for i = find(y(:, iycol))'
0304         yt = y(:, iycol);
0305         yt(i) = 0;
0306         
0307         model.rxnsPresent = ~(model.selRxnMatrix * yt);
0308         solution2 = fluxBalance(model, selTargetRxns, false);   
0309         
0310         if solution2.obj >= solution1.obj
0311             y(:, iycol) = yt;
0312             
0313             y(:, iycol) = shrinkKnockouts(y(:, iycol), model, selTargetRxns);
0314         end
0315     end
0316 end
0317 
0318 y = unique(y', 'rows')';
0319 end
0320 
0321 function [solSynMax solSynMin solBiomass]  = fluxBalance(model,selTargetRxns,verbFlag)
0322 if nargin < 3
0323     verbFlag = true;
0324 end
0325 model.x0 = [];
0326 modelb = model;
0327 model_syn = model;
0328 
0329 [nMets nRxns] = size(model.S);
0330 
0331 yt = model.rxnsPresent;
0332 
0333 modelb.A = [ model.S; 
0334       sparse(1:nnz(~yt), find(~yt), ones(nnz(~yt), 1), nnz(~yt), nRxns) ];
0335 modelb.b = [ zeros(nMets, 1);
0336       zeros(nnz(~yt), 1) ];
0337 modelb.csense = char('E' * ones(1, nMets + nnz(~yt)));
0338 modelb.vartype = char('C' * ones(1, nRxns));
0339 modelb.osense = -1;
0340 solBiomass = solveCobraMILP(modelb);
0341 
0342 
0343 model_syn.A = [ model.S;
0344       sparse(1:nnz(~yt), find(~yt), ones(nnz(~yt), 1), nnz(~yt), nRxns);
0345       model.c' ];
0346 model_syn.b = [ zeros(nMets, 1);
0347       zeros(nnz(~yt), 1);
0348       solBiomass.obj ];
0349 model_syn.c = selTargetRxns;
0350 model_syn.csense = char('E' * ones(1, nMets + nnz(~yt) + 1));
0351 model_syn.vartype = char('C' * ones(1, nRxns));
0352 model_syn.osense = -1;
0353 solSynMax = solveCobraMILP(model_syn);
0354 model_syn.osense = 1;
0355 solSynMin = solveCobraMILP(model_syn);
0356 
0357 if verbFlag
0358 fprintf('Biomass flux:    %f\n', solBiomass.obj);
0359 fprintf('Synthetic flux:  [%f, %f]\n', solSynMin.obj, solSynMax.obj);
0360 end
0361 end

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