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