Home > RAVEN > consumeSomething.m

consumeSomething

PURPOSE ^

consumeSomething

SYNOPSIS ^

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

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

   Rasmus Agren, 2013-08-01

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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