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