checkProduction Checks which metabolites that can be produced from a model using the specified constraints. model a model structure checkNeededForProduction for each of the metabolites that could not be produced, include an artificial production reaction and calculate which new metabolites that could be produced as en effect of this (opt, default false) excretionFromCompartments cell array with compartment ids from which metabolites can be excreted (opt, default model.comps) printDetails print details to the screen (opt, default true) notProduced cell array with metabolites that could not be produced notProducedNames cell array with names and compartments for metabolites that could not be produced neededForProductionMat matrix where n x m is true if metabolite n allows for production of metabolite m minToConnect structure with the minimal number of metabolites that need to be connected in order to be able to produce all other metabolites and which metabolites each of them connects model updated model structure with excretion reactions added The function is intended to be used to identify which metabolites must be connected in order to have a fully connected network. It does so by first identifying which metabolites could have a net production in the network. Then it calculates which other metabolites must be able to have net production in order to have production of all metabolites in the network. So, if a network contains the equations A[external]->B, C->D, and D->E it will identify that production of C will connect the metabolites D and E. Usage: [notProduced, notProducedNames,neededForProductionMat,minToConnect,model]=... checkProduction(model,checkNeededForProduction,... excretionFromCompartments,printDetails) Rasmus Agren, 2013-02-06
0001 function [notProduced, notProducedNames, neededForProductionMat,minToConnect,model]=checkProduction(model,checkNeededForProduction,excretionFromCompartments,printDetails) 0002 % checkProduction 0003 % Checks which metabolites that can be produced from a model using the 0004 % specified constraints. 0005 % 0006 % model a model structure 0007 % checkNeededForProduction for each of the metabolites that could not 0008 % be produced, include an artificial 0009 % production reaction and calculate which new 0010 % metabolites that could be produced as en 0011 % effect of this (opt, default false) 0012 % excretionFromCompartments cell array with compartment ids from which 0013 % metabolites can be excreted (opt, default 0014 % model.comps) 0015 % printDetails print details to the screen (opt, default 0016 % true) 0017 % 0018 % notProduced cell array with metabolites that could not 0019 % be produced 0020 % notProducedNames cell array with names and compartments for 0021 % metabolites that could not be produced 0022 % neededForProductionMat matrix where n x m is true if metabolite n 0023 % allows for production of metabolite m 0024 % minToConnect structure with the minimal number of 0025 % metabolites that need to be connected in 0026 % order to be able to produce all other 0027 % metabolites and which metabolites each of 0028 % them connects 0029 % model updated model structure with excretion 0030 % reactions added 0031 % 0032 % The function is intended to be used to identify which metabolites must 0033 % be connected in order to have a fully connected network. It does so by 0034 % first identifying which metabolites could have a net production in the 0035 % network. Then it calculates which other metabolites must be able to 0036 % have net production in order to have production of all metabolites in 0037 % the network. So, if a network contains the equations A[external]->B, 0038 % C->D, and D->E it will identify that production of C will connect 0039 % the metabolites D and E. 0040 % 0041 % Usage: [notProduced, notProducedNames,neededForProductionMat,minToConnect,model]=... 0042 % checkProduction(model,checkNeededForProduction,... 0043 % excretionFromCompartments,printDetails) 0044 % 0045 % Rasmus Agren, 2013-02-06 0046 % 0047 0048 if nargin<2 0049 checkNeededForProduction=false; 0050 end 0051 0052 if nargin<3 0053 excretionFromCompartments=model.comps; 0054 end 0055 0056 if nargin<4 0057 printDetails=true; 0058 end 0059 0060 %Add an exchange reaction for each metabolite in the allowed compartments 0061 %and see if it can carry a flux 0062 allowedMetIds=ismember(model.comps(model.metComps),excretionFromCompartments); 0063 allowedMetIndexes=find(allowedMetIds); 0064 [model addedRxns]=addExchangeRxns(model,'out',allowedMetIndexes); 0065 0066 canProduce=haveFlux(model,10^-5,addedRxns); 0067 0068 notProduced=find(~canProduce); 0069 minToConnect={}; 0070 if checkNeededForProduction==true 0071 %For each of the metabolites that couldn't be produced allow uptake and check 0072 %which of the other metabolites that couldn't be produced that can be 0073 %produced 0074 neededForProductionMat=false(numel(notProduced)); 0075 for i=1:numel(notProduced) 0076 %Add uptake for this metabolite 0077 if i>1 0078 %Reset last iteration 0079 model.S(:,numel(model.rxns)-numel(addedRxns)+notProduced(i-1))=model.S(:,numel(model.rxns)-numel(addedRxns)+notProduced(i-1))*-1; 0080 end 0081 %Change the production reaction to an uptake reaction 0082 model.S(:,numel(model.rxns)-numel(addedRxns)+notProduced(i))=model.S(:,numel(model.rxns)-numel(addedRxns)+notProduced(i))*-1; 0083 0084 %Test which of the metabolites that couldn't be produced that can 0085 %be produced now 0086 neededForProductionMat(i,:)=haveFlux(model,10^-5,addedRxns(notProduced)); 0087 end 0088 %Calculate the smallest number of metabolites that must be connected to 0089 %make everything connected and return their names 0090 0091 %The algorithm is relatively straight forward. It finds the metabolite 0092 %that connects the most unconnected metabolites (iteratively), adds it 0093 %and removes the now connected metabolites until all are connected. 0094 %This is not guaranteed to find the global minimum 0095 neededForProdTemp=neededForProductionMat; 0096 while 1==1 0097 %Get all metabolites a metabolite makes connected 0098 totalConnected=false(size(neededForProdTemp)); 0099 for i=1:numel(notProduced) 0100 totalConnected(i,:)=neededForProdTemp(i,:); 0101 0102 lastIter=0; 0103 while 1==1 0104 [crap a crap]=find(neededForProdTemp(totalConnected(i,:),:)); 0105 totalConnected(i,a)=true; 0106 if numel(a)==lastIter 0107 break; %No more connections were possible 0108 else 0109 lastIter=numel(a); 0110 end 0111 end 0112 end 0113 [connections mostConnected]=max(sum(totalConnected,2)); 0114 0115 if connections>0 0116 %Add the most connected metabolite to the list and remove all 0117 %metabolites that it's connected to 0118 metID=allowedMetIndexes(notProduced(mostConnected)); 0119 entry=[model.metNames{metID},'[',model.comps{model.metComps(metID)},'] (connects ' num2str(connections) ' metabolites)']; 0120 minToConnect=[minToConnect;entry]; 0121 neededForProdTemp(totalConnected(mostConnected,:),:)=false; 0122 neededForProdTemp(:,totalConnected(mostConnected,:))=false; 0123 else 0124 break; 0125 end 0126 end 0127 else 0128 neededForProductionMat=[]; 0129 end 0130 0131 notProducedNames=strcat(model.metNames(allowedMetIndexes(notProduced)),'[',model.comps(model.metComps(allowedMetIndexes(notProduced))),']'); 0132 0133 if printDetails==true 0134 fprintf('The following metabolites could not be produced:\n'); 0135 [notProducedNamesTemp,perm]=sort(notProducedNames); 0136 0137 if checkNeededForProduction==true 0138 neededForProdTemp=neededForProductionMat(:,perm); 0139 neededForProdTemp=neededForProdTemp(perm,:); 0140 fprintf('\tIf the production of a metabolite is dependent on some other metabolites then those are printed under the name\n\n'); 0141 end 0142 for i=1:numel(notProducedNamesTemp) 0143 fprintf([notProducedNamesTemp{i} '\n']); 0144 neededForProdTemp(i,i)=false; %Not neat to do this here. Prevent printing itself 0145 if checkNeededForProduction==true 0146 enablesProduction=find(neededForProdTemp(:,i)); 0147 if any(enablesProduction) 0148 for j=1:numel(enablesProduction) 0149 fprintf(['\t' notProducedNamesTemp{enablesProduction(j)} '\n']); 0150 end 0151 end 0152 end 0153 end 0154 end 0155 end