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

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 (optional, 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 (optional, 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 (optional, default
0021 %                   false)
0022 %   allowExcretion  allow for excretion of all other metabolites (optional,
0023 %                   default true)
0024 %   params          parameter structure as used by getMILPParams (optional)
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 (optional,
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 if nargin<2
0046     ignoreMets=[];
0047 elseif ~islogical(ignoreMets) && ~isnumeric(ignoreMets)
0048     ignoreMets=convertCharArray(ignoreMets);
0049 end
0050 if nargin<3
0051     isNames=false;
0052 end
0053 if nargin<4
0054     minNrFluxes=false;
0055 end
0056 if nargin<5
0057     allowExcretion=true;
0058 end
0059 if nargin<6
0060     params.relGap=0.8;
0061 end
0062 if nargin<7
0063     ignoreIntBounds=false;
0064 end
0065 
0066 if isNames==true && ~isempty(ignoreMets)
0067     %Check that metsToRemove is a cell array
0068     if iscellstr(ignoreMets)==false
0069         EM='Must supply a cell array of strings if isNames=true';
0070         dispEM(EM);
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     [~, 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=1:numel(model.rev);
0124 padding=num2cell(padding)';
0125 padding=cellfun(@num2str,padding,'uniformoutput',false);
0126 model.rxns=padding;
0127 model.rxnNames=padding;
0128 model.eccodes=padding;
0129 model.rxnMiriams=padding;
0130 model.grRules=padding;
0131 if isfield(model,'genes')
0132     model.rxnGeneMat=sparse(numel(model.rev),numel(model.genes));
0133 end
0134 model.subSystems=padding;
0135 model.rxnFrom=padding;
0136 model.rxnComps=ones(numel(model.rev),1);
0137 model.rxnNotes=padding;
0138 model.rxnReferences=padding;
0139 model.rxnConfidenceScores=NaN(numel(model.rev),1);
0140 
0141 sol=solveLP(model,1);
0142 if any(sol.x)
0143     %It could be that several metabolites were produced in order to get the
0144     %best solution. The setdiff is to avoid including the last fake
0145     %metabolite
0146     I=setdiff(find(sol.x(nRxns+1:end)>0.1),size(model.S,1));
0147     
0148     if any(I) %This should always be true
0149         %Change the coefficients so that only the first is produced. This
0150         %is not always possible, but it is tested for since it it results
0151         %in more easily interpretable results
0152         
0153         oldS=model.S;
0154         foundSingle=false;
0155         %Test if any of the metabolites could be produced on their own
0156         for i=1:numel(I)
0157             model.S=oldS;
0158             J=nRxns+1:numel(model.lb)-1;
0159             J(I(i))=[];
0160             model.S(:,J)=0;
0161             sol=solveLP(model);
0162             if any(sol.x)
0163                 foundSingle=true;
0164                 break;
0165             end
0166         end
0167         %This means that there was no metabolite which could be produced on
0168         %its own. Then let all the produceable metabolites be produced.
0169         if foundSingle==false
0170             model.S=oldS;
0171         end
0172         if minNrFluxes==true
0173             %Has to add names for the rxns to prevent an error in
0174             %minNrFluxes
0175             model.rxns=cell(numel(model.lb),1);
0176             model.rxns(:)=sprintfc('tmp_%d',1:numel(model.lb));
0177             model.mets=cell(size(model.b,1),1);
0178             model.mets(:)={'TEMP'};
0179             sol=solveLP(model,3,params);
0180         else
0181             sol=solveLP(model,1);
0182         end
0183         solution=sol.x(1:nRxns);
0184         metabolite=find(sol.x(nRxns+1:end-1)>0.1);
0185     end
0186 end
0187 end

Generated by m2html © 2005