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, 2013-08-01

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, 2013-08-01
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         dispEM('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',false);
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         dispEM('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 if useModelConstraints==true
0145     newConnected={};
0146     cannotConnect={};
0147     
0148     %Check that the input model isn't solveable without any input
0149     sol=solveLP(model);
0150     if ~isempty(sol.f)
0151         addedRxns={};
0152         newModel=model;
0153         exitFlag=1;
0154         return;
0155     end
0156     
0157     %Then check that the merged model is solveable
0158     sol=solveLP(allModels);
0159     if isempty(sol.f)
0160         dispEM('There are no reactions in the template model(s) that can make the model constraints satisfied'); 
0161     end
0162     
0163     %Remove dead ends for speed reasons. This has to be done here and
0164     %duplicate below because there is otherwise a risk that a reaction
0165     %which is constrained to something relevant is removed
0166     allModels=simplifyModel(allModels,false,false,false,true,false,false,[],true);
0167     allModels.c(:)=0;
0168 else
0169     %Remove dead ends for speed reasons
0170     allModels=simplifyModel(allModels,false,false,false,true,false,false,[],true);
0171     allModels.c(:)=0;
0172 
0173     %If model constraints shouldn't be used, then determine which reactions to
0174     %force to have flux
0175 
0176     %Get the reactions that can carry flux in the original model
0177     originalFlux=haveFlux(model,1);
0178     
0179     %For the ones that can't carry flux, see if they can do so in the merged
0180     %model
0181     toCheck=intersect(allModels.rxns(strcmp(allModels.rxnFrom,model.id)),model.rxns(~originalFlux));
0182 
0183     %Get the ones that still cannot carry flux
0184     I=haveFlux(allModels,1,toCheck);
0185 
0186     %Get the reactions that can't carry flux in the original model, but can
0187     %in the merged one
0188     K=toCheck(I);
0189 
0190     %This is a temporary thing to only look at the non-reversible rxns.
0191     %This is because all reversible rxns can have a flux in the irreversible
0192     %model format that is used by getMinNrFluxes
0193     [crap I]=ismember(K,model.rxns);
0194     K(model.rev(I)~=0)=[];
0195 
0196     %Constrain all reactions in the original model to have a flux
0197     allModels.lb(ismember(allModels.rxns,K))=0.1;    
0198     
0199     %Return stuff
0200     newConnected=K;
0201     cannotConnect=setdiff(model.rxns(~originalFlux ),newConnected);
0202 end
0203 
0204 %Then minimize for the number of fluxes used. The fixed rxns doesn't need
0205 %to participate
0206 templateRxns=find(~strcmp(allModels.rxnFrom,model.id));
0207 [crap J exitFlag]=getMinNrFluxes(allModels,templateRxns,params,allModels.rxnScores(templateRxns));
0208 
0209 %Remove everything except for the added ones
0210 I=true(numel(allModels.rxns),1);
0211 I(templateRxns(J))=false;
0212 addedModel=removeRxns(allModels,I,true);
0213 
0214 newModel=mergeModels({model;addedModel},true);
0215 addedRxns=setdiff(newModel.rxns,model.rxns);
0216 newModel=rmfield(newModel,'rxnScores');
0217 end

Generated on Mon 06-Jan-2014 14:58:12 by m2html © 2005