Home > RAVEN > makeSomething.m

makeSomething

PURPOSE ^

makeSomething

SYNOPSIS ^

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

DESCRIPTION ^

 makeSomething
   Tries to excrete any metabolite using as few reactions as possible.
   The intended use is when you want to make sure that you model cannot
   synthesize anything from nothing. It is then a faster way than to use
   checkProduction or canProduce

   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)
   allowExcretion  allow for excretion of all other metabolites (opt,
                   default true)
   params          parameter structure as used by getMILPParams (opt)

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

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

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

   Rasmus Agren, 2013-02-06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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