Home > RAVEN > consumeSomething.m

consumeSomething

PURPOSE ^

consumeSomething

SYNOPSIS ^

function [solution metabolite]=consumeSomething(model,ignoreMets,isNames,minNrFluxes,params)

DESCRIPTION ^

 consumeSomething
   Tries to consume any metabolite using as few reactions as possible.
   The intended use is when you want to make sure that you model cannot
   consume anything without producing something. It is intended to be used
   with no active exchange reactions.

   model           a model structure
   ignoreMets      either a cell array of metabolite IDs, a logical vector 
                   with the same number of elements as metabolites in the model,
                   of a vector of indexes for metabolites to exclude from
                   this analysis (opt, default [])
   isNames         true if the supplied mets represent metabolite names
                   (as opposed to IDs). This is a way to delete
                   metabolites in several compartments at once without
                   knowing the exact IDs. This only works if ignoreMets
                   is a cell array (opt, default false)
   minNrFluxes     solves the MILP problem of minimizing the number of
                   fluxes instead of the sum. Slower, but can be
                   used if the sum gives too many fluxes (opt, default
                   false)
   params          parameter structure as used by getMILPParams (opt)

   solution        flux vector for the solution
   metabolite      the index of the metabolite(s) which was consumed. If
                   possible only one metabolite is reported, but there are
                   situations where metabolites can only be consumed in
                   pairs (or more)

   NOTE: This works by forcing at least 1 unit of "any metabolites" to be
   consumed and then minimize for the sum of fluxes. If more than one
   metabolite is consumed, it picks one of them to be consumed and then 
   minimizes for the sum of fluxes.

   Usage: [solution metabolite]=consumeSomething(model,ignoreMets,isNames,...
               minNrFluxes,params)

   Rasmus Agren, 2013-02-06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [solution metabolite]=consumeSomething(model,ignoreMets,isNames,minNrFluxes,params)
0002 % consumeSomething
0003 %   Tries to consume any metabolite using as few reactions as possible.
0004 %   The intended use is when you want to make sure that you model cannot
0005 %   consume anything without producing something. It is intended to be used
0006 %   with no active exchange reactions.
0007 %
0008 %   model           a model structure
0009 %   ignoreMets      either a cell array of metabolite IDs, a logical vector
0010 %                   with the same number of elements as metabolites in the model,
0011 %                   of a vector of indexes for metabolites to exclude from
0012 %                   this analysis (opt, default [])
0013 %   isNames         true if the supplied mets represent metabolite names
0014 %                   (as opposed to IDs). This is a way to delete
0015 %                   metabolites in several compartments at once without
0016 %                   knowing the exact IDs. This only works if ignoreMets
0017 %                   is a cell array (opt, default false)
0018 %   minNrFluxes     solves the MILP problem of minimizing the number of
0019 %                   fluxes instead of the sum. Slower, but can be
0020 %                   used if the sum gives too many fluxes (opt, default
0021 %                   false)
0022 %   params          parameter structure as used by getMILPParams (opt)
0023 %
0024 %   solution        flux vector for the solution
0025 %   metabolite      the index of the metabolite(s) which was consumed. If
0026 %                   possible only one metabolite is reported, but there are
0027 %                   situations where metabolites can only be consumed in
0028 %                   pairs (or more)
0029 %
0030 %   NOTE: This works by forcing at least 1 unit of "any metabolites" to be
0031 %   consumed and then minimize for the sum of fluxes. If more than one
0032 %   metabolite is consumed, it picks one of them to be consumed and then
0033 %   minimizes for the sum of fluxes.
0034 %
0035 %   Usage: [solution metabolite]=consumeSomething(model,ignoreMets,isNames,...
0036 %               minNrFluxes,params)
0037 %
0038 %   Rasmus Agren, 2013-02-06
0039 %
0040 
0041 if nargin<2
0042     ignoreMets=[];
0043 end
0044 
0045 if nargin<3
0046     isNames=false;
0047 end
0048 
0049 if nargin<4
0050     minNrFluxes=false;
0051 end
0052 
0053 if nargin<5
0054     params.relGap=0.8;
0055 end
0056 
0057 if isNames==true && ~isempty(ignoreMets)
0058    %Check that metsToRemove is a cell array
0059    if iscellstr(ignoreMets)==false
0060         throw(MException('','Must supply a cell array of strings if isNames=true'));
0061    end 
0062 end
0063 
0064 if isNames==false
0065     indexesToIgnore=getIndexes(model,ignoreMets,'mets');
0066 else
0067     indexesToIgnore=[];
0068     for i=1:numel(ignoreMets)
0069        indexesToIgnore=[indexesToIgnore;find(strcmp(ignoreMets(i),model.metNames))]; 
0070     end
0071 end
0072 
0073 solution=[];
0074 metabolite=[];
0075 
0076 nRxns=numel(model.rxns);
0077 nMets=numel(model.mets);
0078 
0079 %Add uptake reactions for all metabolites
0080 model.S=[model.S speye(nMets)];
0081 
0082 %Add so that they all consume a fake metabolite
0083 model.S=[model.S;[sparse(1,nRxns) ones(1,nMets)*-1]];
0084 
0085 %Change so that the ignoreMets have a coefficient 0 in their reactions.
0086 %Does not remove the actual reaction to make mapping easier later
0087 model.S(:,indexesToIgnore+nRxns)=0;
0088 
0089 %Add an uptake reaction for this last fake metabolite
0090 model.S(size(model.S,1),size(model.S,2)+1)=1;
0091 model.b=[model.b;zeros(1,size(model.b,2))];
0092 model.lb=[model.lb;zeros(nMets,1);1];
0093 model.ub=[model.ub;inf(nMets+1,1)];
0094 model.rev=[model.rev;zeros(nMets+1,1)];
0095 model.c=zeros(size(model.S,2),1);
0096 
0097 sol=solveLP(model,1);   
0098 if any(sol.x)
0099    %It could be that several metabolites were consumed in order to get the
0100    %best solution.
0101    %The setdiff is to avoid including the last fake metabolite
0102    I=setdiff(find(sol.x(nRxns+1:end)>0.1),size(model.S,1));
0103     
0104    if any(I) %This should always be true
0105         %Change the coefficients so that only the first is
0106         %consumed. This is not always possible, but it is tested for since it it
0107         %results in more easily interpretable results
0108         
0109         oldS=model.S;
0110         foundSingle=false;
0111         %Test if any of the metabolites could be consumed on their own
0112         for i=1:numel(I)
0113             model.S=oldS;
0114             J=nRxns+1:numel(model.lb)-1;
0115             J(I(i))=[];
0116             model.S(:,J)=0;
0117             sol=solveLP(model);
0118             if any(sol.x)
0119                 foundSingle=true;
0120                 break;
0121             end
0122         end
0123         %This means that there was no metabolite which could be consumed on
0124         %its own. Then let all the consumeable metabolites be consumed.
0125         if foundSingle==false
0126             model.S=oldS;
0127         end
0128         if minNrFluxes==true
0129             %Has to add names for the rxns to prevent an error in
0130             %minNrFluxes
0131             model.rxns=cell(numel(model.lb),1);
0132             model.rxns(:)={'TEMP'};
0133             model.mets=cell(size(model.b,1),1);
0134             model.mets(:)={'TEMP'};
0135             sol=solveLP(model,3,params);
0136         else
0137             sol=solveLP(model,1);
0138         end
0139         solution=sol.x(1:nRxns);
0140         metabolite=find(sol.x(nRxns+1:end-1)>0.1);
0141    end
0142 end
0143 end

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