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

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-08-01
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         dispEM('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 %Add padding to the reaction annotation to prevent an error in solveLP
0123 padding=cell(numel(model.rev),1);
0124 padding(:)={''};
0125 model.rxns=padding;
0126 model.rxnNames=padding;
0127 model.eccodes=padding;
0128 model.rxnMiriams=padding;
0129 model.grRules=padding;
0130 if isfield(model,'genes')
0131     model.rxnGeneMat=sparse(numel(model.rev),numel(model.genes));
0132 end
0133 model.subSystems=padding;
0134 model.rxnFrom=padding;
0135 model.rxnComps=ones(numel(model.rev),1);
0136 
0137 sol=solveLP(model,1);   
0138 if any(sol.x)
0139    %It could be that several metabolites were produced in order to get the
0140    %best solution.
0141    %The setdiff is to avoid including the last fake metabolite
0142    I=setdiff(find(sol.x(nRxns+1:end)>0.1),size(model.S,1));
0143     
0144    if any(I) %This should always be true
0145         %Change the coefficients so that only the first is
0146         %produced. This is not always possible, but it is tested for since it it
0147         %results in more easily interpretable results
0148         
0149         oldS=model.S;
0150         foundSingle=false;
0151         %Test if any of the metabolites could be produced on their own
0152         for i=1:numel(I)
0153             model.S=oldS;
0154             J=nRxns+1:numel(model.lb)-1;
0155             J(I(i))=[];
0156             model.S(:,J)=0;
0157             sol=solveLP(model);
0158             if any(sol.x)
0159                 foundSingle=true;
0160                 break;
0161             end
0162         end
0163         %This means that there was no metabolite which could be produced on
0164         %its own. Then let all the produceable metabolites be produced.
0165         if foundSingle==false
0166             model.S=oldS;
0167         end
0168         if minNrFluxes==true
0169             %Has to add names for the rxns to prevent an error in
0170             %minNrFluxes
0171             model.rxns=cell(numel(model.lb),1);
0172             model.rxns(:)={'TEMP'};
0173             model.mets=cell(size(model.b,1),1);
0174             model.mets(:)={'TEMP'};
0175             sol=solveLP(model,3,params);
0176         else
0177             sol=solveLP(model,1);
0178         end
0179         solution=sol.x(1:nRxns);
0180         metabolite=find(sol.x(nRxns+1:end-1)>0.1);
0181    end
0182 end
0183 end

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