runSimpleOptKnock

 runSimpleOptKnock
   Simple OptKnock algorithm that checks all gene or reaction deletions
   for growth-coupled metabolite production, by testing all possible
   combinations. This is not defined as MILP, and is therefore slow (but
   simple).

 Input:
    model          a model structure
    targetRxn      identifier of target reaction
    biomassRxn     identifier of biomass reaction
    deletions      cell array with gene or reaction identifiers that
                   should be considered for knockout
                   (optional, default = model.rxns)
    genesOrRxns    string indicating whether deletions parameter is given
                   with 'genes' or 'rxns' identifiers (optional, default
                   'rxns')
    maxNumKO       numeric with maximum number of simulatenous knockout
                   (optional, default 1)
    minGrowth      numeric of minimum growth rate (optional, default 0.05)

 Output:
    out            structure with deletions strategies that result in
                   growth-coupled production
       KO          cell array with gene(s) or reaction(s) to be deleted
       growthRate  vector with growth rates after deletion
       prodRate    vector with production rates after deletion

 Usage: out = runSimpleOptKnock(model, targetRxn, biomassRxn, deletions,...
                   genesOrRxns, maxNumKO, minGrowth)

0001 function out = runSimpleOptKnock(model, targetRxn, biomassRxn, deletions, genesOrRxns, maxNumKO, minGrowth)
0002 % runSimpleOptKnock
0003 %   Simple OptKnock algorithm that checks all gene or reaction deletions
0004 %   for growth-coupled metabolite production, by testing all possible
0005 %   combinations. This is not defined as MILP, and is therefore slow (but
0006 %   simple).
0007 %
0008 % Input:
0009 %    model          a model structure
0010 %    targetRxn      identifier of target reaction
0011 %    biomassRxn     identifier of biomass reaction
0012 %    deletions      cell array with gene or reaction identifiers that
0013 %                   should be considered for knockout
0014 %                   (optional, default = model.rxns)
0015 %    genesOrRxns    string indicating whether deletions parameter is given
0016 %                   with 'genes' or 'rxns' identifiers (optional, default
0017 %                   'rxns')
0018 %    maxNumKO       numeric with maximum number of simulatenous knockout
0019 %                   (optional, default 1)
0020 %    minGrowth      numeric of minimum growth rate (optional, default 0.05)
0021 %
0022 % Output:
0023 %    out            structure with deletions strategies that result in
0024 %                   growth-coupled production
0025 %       KO          cell array with gene(s) or reaction(s) to be deleted
0026 %       growthRate  vector with growth rates after deletion
0027 %       prodRate    vector with production rates after deletion
0028 %
0029 % Usage: out = runSimpleOptKnock(model, targetRxn, biomassRxn, deletions,...
0030 %                   genesOrRxns, maxNumKO, minGrowth)
0031 
0032 if nargin < 4
0033     params.deletions = model.rxns;
0034 else
0035     params.deletions = deletions;
0036 end
0037 if nargin < 5
0038     params.genesOrRxns = 'rxns';
0039 else
0040     params.genesOrRxns = genesOrRxns;
0041 end
0042 if nargin < 6
0043     params.maxNumKO = 1;
0044 else
0045     params.maxNumKO = maxNumKO;
0046 end
0047 if nargin < 7
0048     params.minGrowth = 0.05;
0049 else
0050     params.minGrowth = minGrowth;
0051 end
0052 
0053 % Number of deletions
0054 out.KO         = cell(0,params.maxNumKO); % The KO genes/rxns
0055 out.growthRate = zeros(0);
0056 out.prodRate   = zeros(0);
0057 out.score      = zeros(0);
0058 
0059 params.biomassIdx  = getIndexes(model,biomassRxn,'rxns');
0060 params.targetIdx   = getIndexes(model,targetRxn,'rxns');
0061 
0062 model = setParam(model,'obj',params.biomassIdx,1);
0063 [solWT, hsSol] = solveLP(model);
0064 WT.minScore = solWT.x(params.targetIdx)*solWT.f;
0065 
0066 fprintf('Running simple OptKnock analysis...   0%% complete');
0067 KO=zeros(1,params.maxNumKO);
0068 [~,~,out,~] = knockoutIteration(model,params,WT,out,params.maxNumKO,KO,[],hsSol);
0069 
0070 if size(out.KO,2)>1
0071     singleKO = cellfun(@isempty,out.KO(:,1));
0072     if any(singleKO)
0073         singleKO = out.KO{singleKO,2};
0074         singleKO = strcmp(out.KO(:,1),singleKO);
0075         out.KO(singleKO,:) = [];
0076         out.growthRate(singleKO,:) = [];
0077         out.prodRate(singleKO,:) = [];
0078         out.score(singleKO,:) = [];
0079     end
0080 end
0081 
0082 
0083 fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\bCOMPLETE\n');
0084 end
0085 
0086 function [model,iteration,out,KO,hsSol] = knockoutIteration(model,params,WT,out,iteration,KO,minScore,hsSol)
0087 if nargin<7 || isempty(minScore)
0088     minScore = WT.minScore;
0089 end
0090 iteration = iteration - 1;
0091 for i = 1:numel(params.deletions)
0092     if iteration+1==params.maxNumKO
0093         progress=pad(num2str(floor(i/numel(params.deletions)*100)),3,'left');
0094         fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b%s%% complete',progress);
0095     end
0096     KO(iteration+1)=i;
0097     switch params.genesOrRxns
0098         case 'rxns'
0099             modelKO = setParam(model,'eq',params.deletions{i},0);
0100         case 'genes'
0101             modelKO = removeGenes(model,params.deletions{i},false,false,false);
0102     end
0103     solKO = solveLP(modelKO,0,[],hsSol);
0104     if ~isempty(solKO.f)
0105         growthRate = solKO.f;
0106         prodRate   = solKO.x(params.targetIdx);
0107         prodRate(prodRate<1e-10)=0; % Filter out results from solver tolerance
0108         if growthRate > params.minGrowth
0109             if growthRate*prodRate > minScore*1.01
0110                 out.KO(end+1,find(KO)) = transpose(params.deletions(KO(find(KO))));
0111                 out.growthRate(end+1,1) = growthRate;
0112                 out.prodRate(end+1,1)   = prodRate;
0113                 out.score(end+1,1)      = growthRate*prodRate;
0114             end
0115             if iteration>0
0116                 [~,~,out] = knockoutIteration(modelKO,params,WT,out,iteration,KO,growthRate*prodRate,hsSol);
0117             end
0118         end
0119     end
0120 end
0121 end