Home > RAVEN > makeSomething.m

makeSomething

PURPOSE ^

makeSomething

SYNOPSIS ^

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

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)
   ignoreIntBounds    true if internal bounds (including reversibility)
                   should be ignored. Exchange reactions are not affected.
                   This can be used to find unbalanced solutions which are
                   not possible using the default constraints (opt,
                   default false)

   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,ignoreIntBounds)

   Rasmus Agren, 2013-05-16

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Tue 16-Jul-2013 21:50:02 by m2html © 2005