Home > RAVEN > fillGaps.m

fillGaps

PURPOSE ^

fillGaps

SYNOPSIS ^

function [newConnected cannotConnect addedRxns newModel exitFlag]=fillGaps(model,models,allowNetProduction,useModelConstraints,supressWarnings,rxnScores,params)

DESCRIPTION ^

 fillGaps
   Uses template model(s) to fill gaps in a model

   model               a model structure that may contains gaps to be filled
   models              a cell array of reference models or a model structure. 
                       The gaps will be filled using reactions from these models
   allowNetProduction  true if net production of all metabolites is
                       allowed. A reaction can be unable to carry flux because one of
                       the reactants is unavailable or because one of the
                       products can't be further processed. If this
                       parameter is true, only the first type of
                       unconnectivity is considered (opt, default false)
   useModelConstraints true if the constraints specified in the model
                       structure should be used. If false then reactions
                       included from the template model(s) so that as many
                       reactions as possible in model can carry flux
                       (opt, default false)
   supressWarnings     false if warnings should be displayed (opt, default
                       false)
   rxnScores           array with scores for each of the reactions in the 
                       reference model(s). If more than one model is supplied,
                       then rxnScores should be a cell array of vectors.
                       The solver will try to maximize the sum of the
                       scores for the included reactions (opt, default
                       is -1 for all reactions)
   params              parameter structure as used by getMILPParams (opt)

   newConnected        cell array with the reactions that could be
                       connected. This is not calulated if
                       useModelConstraints is true
   cannotConnect       cell array with reactions that could not be
                       connected. This is not calculated if
                       useModelConstraints is true
   addedRxns           cell array with the reactions that were added from
                       "models"
   newModel            the model with reactions added to fill gaps
   exitFlag            1: optimal solution found
                      -1: no feasible solution found
                      -2: optimization time out

   This method works by merging the model to the reference model(s) and 
   checking which reactions can carry flux. All reactions that can't 
   carry flux are removed (cannotConnect).
   If useModelConstraints is false it then solves the MILP problem of 
   minimizing the number of active reactions from the reference models 
   that are required to have flux in all the reactions in model. This 
   requires that the input model has exchange reactions present for the 
   nutrients that are needed for its metabolism. If useModelConstraints is 
   true then the problem is to include as few reactions as possible from 
   the reference models in order to satisfy the model constraints. 
   The intended use is that the user can attempt a general gap-filling using
   useModelConstraint=false or a more targeted gap-filling by setting
   constraints in the model structure and then use
   useModelConstraints=true. Say that the user want to include reactions
   so that all biomass components can be synthesized. He/she could then
   define a biomass equation and set the lower bound to >0. Running this
   function with useModelConstraints=true would then give the smallest set
   of reactions that have to be included in order for the model to produce
   biomass.
   
   Usage: [newConnected cannotConnect addedRxns newModel exitFlag]=...
           fillGaps(model,models,allowNetProduction,useModelConstraints,...
           supressWarnings,rxnScores,params)

   Rasmus Agren, 2012-10-24

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [newConnected cannotConnect addedRxns newModel exitFlag]=fillGaps(model,models,allowNetProduction,useModelConstraints,supressWarnings,rxnScores,params)
0002 % fillGaps
0003 %   Uses template model(s) to fill gaps in a model
0004 %
0005 %   model               a model structure that may contains gaps to be filled
0006 %   models              a cell array of reference models or a model structure.
0007 %                       The gaps will be filled using reactions from these models
0008 %   allowNetProduction  true if net production of all metabolites is
0009 %                       allowed. A reaction can be unable to carry flux because one of
0010 %                       the reactants is unavailable or because one of the
0011 %                       products can't be further processed. If this
0012 %                       parameter is true, only the first type of
0013 %                       unconnectivity is considered (opt, default false)
0014 %   useModelConstraints true if the constraints specified in the model
0015 %                       structure should be used. If false then reactions
0016 %                       included from the template model(s) so that as many
0017 %                       reactions as possible in model can carry flux
0018 %                       (opt, default false)
0019 %   supressWarnings     false if warnings should be displayed (opt, default
0020 %                       false)
0021 %   rxnScores           array with scores for each of the reactions in the
0022 %                       reference model(s). If more than one model is supplied,
0023 %                       then rxnScores should be a cell array of vectors.
0024 %                       The solver will try to maximize the sum of the
0025 %                       scores for the included reactions (opt, default
0026 %                       is -1 for all reactions)
0027 %   params              parameter structure as used by getMILPParams (opt)
0028 %
0029 %   newConnected        cell array with the reactions that could be
0030 %                       connected. This is not calulated if
0031 %                       useModelConstraints is true
0032 %   cannotConnect       cell array with reactions that could not be
0033 %                       connected. This is not calculated if
0034 %                       useModelConstraints is true
0035 %   addedRxns           cell array with the reactions that were added from
0036 %                       "models"
0037 %   newModel            the model with reactions added to fill gaps
0038 %   exitFlag            1: optimal solution found
0039 %                      -1: no feasible solution found
0040 %                      -2: optimization time out
0041 %
0042 %   This method works by merging the model to the reference model(s) and
0043 %   checking which reactions can carry flux. All reactions that can't
0044 %   carry flux are removed (cannotConnect).
0045 %   If useModelConstraints is false it then solves the MILP problem of
0046 %   minimizing the number of active reactions from the reference models
0047 %   that are required to have flux in all the reactions in model. This
0048 %   requires that the input model has exchange reactions present for the
0049 %   nutrients that are needed for its metabolism. If useModelConstraints is
0050 %   true then the problem is to include as few reactions as possible from
0051 %   the reference models in order to satisfy the model constraints.
0052 %   The intended use is that the user can attempt a general gap-filling using
0053 %   useModelConstraint=false or a more targeted gap-filling by setting
0054 %   constraints in the model structure and then use
0055 %   useModelConstraints=true. Say that the user want to include reactions
0056 %   so that all biomass components can be synthesized. He/she could then
0057 %   define a biomass equation and set the lower bound to >0. Running this
0058 %   function with useModelConstraints=true would then give the smallest set
0059 %   of reactions that have to be included in order for the model to produce
0060 %   biomass.
0061 %
0062 %   Usage: [newConnected cannotConnect addedRxns newModel exitFlag]=...
0063 %           fillGaps(model,models,allowNetProduction,useModelConstraints,...
0064 %           supressWarnings,rxnScores,params)
0065 %
0066 %   Rasmus Agren, 2012-10-24
0067 %
0068 
0069 %If the user only supplied a single template model
0070 if ~iscell(models)
0071     models={models};
0072 end
0073 
0074 if nargin<3
0075     allowNetProduction=false;
0076 end
0077 if nargin<4
0078     useModelConstraints=false;
0079 end
0080 if nargin<5
0081     supressWarnings=false;
0082 end
0083 if nargin<6
0084     rxnScores=cell(numel(models),1);
0085     for i=1:numel(models)
0086        rxnScores{i}=ones(numel(models{i}.rxns),1)*-1; 
0087     end
0088 end
0089 if isempty(rxnScores)
0090     rxnScores=cell(numel(models),1);
0091     for i=1:numel(models)
0092        rxnScores{i}=ones(numel(models{i}.rxns),1)*-1; 
0093     end
0094 end
0095 if nargin<7
0096     params=[];
0097 end
0098 
0099 if ~iscell(rxnScores)
0100     rxnScores={rxnScores};
0101 end
0102 
0103 models=models(:);
0104 rxnScores=rxnScores(:);
0105 
0106 %Check if the original model has an unconstrained field. If so, give a warning
0107 if supressWarnings==false
0108     if isfield(model,'unconstrained');
0109         fprintf('WARNING: This algorithm is meant to function on a model with exchange reactions for uptake and excretion of metabolites. The current model still has the "unconstrained" field.\n');
0110     else
0111         if isempty(getExchangeRxns(model,'both'))
0112             fprintf('NOTE: This algorithm is meant to function on a model with exchange reactions for uptake and excretion of metabolites. The current model does not seem to contain any such reactions.\n');
0113         end
0114     end
0115 end
0116 
0117 %Simplify the template models to remove constrained rxns. At the same time,
0118 %check that the id of the template models isn't the same as the model.
0119 %That would cause an error further down
0120 for i=1:numel(models)
0121    models{i}.rxnScores=rxnScores{i};
0122    models{i}=simplifyModel(models{i},false,false,true);
0123    if strcmpi(models{i}.id,model.id)
0124         throw(MException('','The reference model(s) cannot have the same id as the model')); 
0125    end
0126 end
0127 
0128 %This is a rather ugly solution to the issue that it's a bit tricky to keep
0129 %track of which scores belong to which reactions. This requires that
0130 %removeRxns and mergeModels are modified to check for the new field.
0131 model.rxnScores=zeros(numel(model.rxns),1); 
0132 
0133 %First merge all models into one big one
0134 allModels=mergeModels([{model};models],true);
0135 
0136 %Add that net production is ok
0137 if allowNetProduction==true
0138    %A second column in model.b means that the b field is lower and upper
0139    %bound on the RHS.
0140    model.b=[model.b(:,1) inf(numel(model.mets),1)];
0141    allModels.b=[allModels.b(:,1) inf(numel(allModels.mets),1)];
0142 end
0143 
0144 %Remove dead ends for speed reasons
0145 allModels=simplifyModel(allModels,false,false,false,true,false,false,[],true);
0146 allModels.c(:)=0;
0147 
0148 if useModelConstraints==true
0149     newConnected={};
0150     cannotConnect={};
0151     
0152     %Check that the input model isn't solveable without any input
0153     sol=solveLP(model);
0154     if ~isempty(sol.f)
0155         addedRxns={};
0156         newModel=model;
0157         return;
0158     end
0159     
0160     %Then check that the merged model is solveable
0161     sol=solveLP(allModels);
0162     if isempty(sol.f)
0163         throw(MException('','There are no reactions in the template model(s) that can make the model constraints satisfied')); 
0164     end
0165 else
0166     %If model constraints shouldn't be used, then determine which reactions to
0167     %force to have flux
0168 
0169     %Get the reactions that can carry flux in the original model
0170     originalFlux=haveFlux(model,1);
0171     
0172     %For the ones that can't carry flux, see if they can do so in the merged
0173     %model
0174     toCheck=intersect(allModels.rxns(strcmp(allModels.rxnFrom,model.id)),model.rxns(~originalFlux));
0175 
0176     %Get the ones that still cannot carry flux
0177     I=haveFlux(allModels,1,toCheck);
0178 
0179     %Get the reactions that can't carry flux in the original model, but can
0180     %in the merged one
0181     K=toCheck(I);
0182 
0183     %This is a temporary thing to only look at the non-reversible rxns.
0184     %This is because all reversible rxns can have a flux in the irreversible
0185     %model format that is used by getMinNrFluxes
0186     [crap I]=ismember(K,model.rxns);
0187     K(model.rev(I)~=0)=[];
0188 
0189     %Constrain all reactions in the original model to have a flux
0190     allModels.lb(ismember(allModels.rxns,K))=0.1;    
0191     
0192     %Return stuff
0193     newConnected=K;
0194     cannotConnect=setdiff(model.rxns(~originalFlux ),newConnected);
0195 end
0196 
0197 %Then minimize for the number of fluxes used. The fixed rxns doesn't need
0198 %to participate
0199 templateRxns=find(~strcmp(allModels.rxnFrom,model.id));
0200 [crap J exitFlag]=getMinNrFluxes(allModels,templateRxns,params,allModels.rxnScores(templateRxns));
0201 
0202 %Remove everything except for the added ones
0203 I=true(numel(allModels.rxns),1);
0204 I(templateRxns(J))=false;
0205 addedModel=removeRxns(allModels,I,true);
0206 
0207 newModel=mergeModels({model;addedModel},true);
0208 addedRxns=setdiff(newModel.rxns,model.rxns);
0209 newModel=rmfield(newModel,'rxnScores');
0210 end

Generated on Tue 23-Apr-2013 15:18:37 by m2html © 2005