guessComposition Attempts to guess the composition of metabolites without information about elemental composition model a model structure printResults true if the output should be printed (opt, default true) model a model structure with information about elemental composition added guessedFor indexes for the metabolites for which a composition could be guessed couldNotGuess indexes for the metabolites for which no composition could be assigned This function works in a rather straight forward manner: 1. Get the metabolites which lack composition and participates in at least one reaction where all other metabolites have composition information 2. Loop through them and calculate their composition based on the rest of the involved metabolites. If there are any inconsistencies, so that a given metabolite should have different composition in different equations, then throw an error 3. Go to 1 This simple approach requires that the rest of the metabolites have correct composition information, and that the involved reactions are correct. The function will exit with an error on any inconsistencies, which means that it could also be used as a way of checking the model for errors. Note that just because this exits sucessfully, the calculated compositions could still be wrong (in case that the existing compositions were wrong) Usage: [newModel, guessedFor, couldNotGuess]=guessComposition(model, printResults) Rasmus Agren, 2013-11-06
0001 function [model, guessedFor, couldNotGuess]=guessComposition(model, printResults) 0002 % guessComposition 0003 % Attempts to guess the composition of metabolites without information 0004 % about elemental composition 0005 % 0006 % model a model structure 0007 % printResults true if the output should be printed (opt, default true) 0008 % 0009 % model a model structure with information about elemental 0010 % composition added 0011 % guessedFor indexes for the metabolites for which a composition 0012 % could be guessed 0013 % couldNotGuess indexes for the metabolites for which no 0014 % composition could be assigned 0015 % 0016 % This function works in a rather straight forward manner: 0017 % 0018 % 1. Get the metabolites which lack composition and participates in 0019 % at least one reaction where all other metabolites have composition information 0020 % 2. Loop through them and calculate their composition based on the rest 0021 % of the involved metabolites. If there are any inconsistencies, so that 0022 % a given metabolite should have different composition in different 0023 % equations, then throw an error 0024 % 3. Go to 1 0025 % 0026 % This simple approach requires that the rest of the metabolites have 0027 % correct composition information, and that the involved reactions are 0028 % correct. The function will exit with an error on any inconsistencies, 0029 % which means that it could also be used as a way of checking the model 0030 % for errors. Note that just because this exits sucessfully, the 0031 % calculated compositions could still be wrong (in case that the existing 0032 % compositions were wrong) 0033 % 0034 % Usage: [newModel, guessedFor, couldNotGuess]=guessComposition(model, printResults) 0035 % 0036 % Rasmus Agren, 2013-11-06 0037 % 0038 0039 if nargin<2 0040 printResults=true; 0041 end 0042 0043 %The metabolites for which there is no elemental composition 0044 originalMissing=unique(model.metNames(cellfun(@isempty,model.metFormulas))); 0045 guessedFor={}; 0046 0047 %This is to keep track of if new metabolite compositions were predicted 0048 predicted=true; 0049 while predicted==true 0050 predicted=false; 0051 0052 %Get the unique names (composition should be independent of compartment) 0053 %for the metabolites without composition 0054 metNames=unique(model.metNames(cellfun(@isempty,model.metFormulas))); 0055 0056 %Parse the formulas in the model 0057 [elements, useMat, exitFlag]=parseFormulas(model.metFormulas, true,false); 0058 0059 for i=1:numel(metNames) 0060 %Get the metabolites with this name. Not so neat, but this is a 0061 %fast function anyways 0062 mets=find(ismember(model.metNames,metNames(i))); 0063 0064 currentComp=[]; 0065 0066 %Loop through the metabolites 0067 %-1: Could not assign due to missing info. -2: Could not assign due to contradiction 0068 %1: Composition assigned 0069 metStatus=-1; 0070 for j=1:numel(mets) 0071 %Get the reactions that the metabolite participates in 0072 [crap, I]=find(model.S(mets(j),:)); 0073 if any(I) 0074 for k=1:numel(I) 0075 %Loop through the reactions and check if all other mets in them 0076 %have known composition 0077 eqn=model.S(:,I(k)); 0078 eqn(mets(j))=0; 0079 if all(exitFlag(eqn~=0)==1) 0080 %This means that all other mets had composition. Calculate 0081 %the resulting composition for the unknown one 0082 comp=useMat'*eqn; 0083 0084 %This can result in round off errors if there are 0085 %stoichiometries with many decimals. Ignore values 0086 %below 10^-12 0087 comp(abs(comp)<10^-12)=0; 0088 0089 %Check if the composition consist of both negative and 0090 %positive values. If so, throw an error 0091 if all(comp<=0) || all(comp>=0) 0092 comp=abs(comp); 0093 if isempty(currentComp) 0094 currentComp=comp; 0095 end 0096 %Also to deal with round off errors 0097 if all(abs(currentComp-comp)<10^-10) 0098 metStatus=1; 0099 else 0100 metStatus=-2; 0101 break; 0102 0103 %%Check if there is an inconcistency 0104 %if any(currentComp~=comp) 0105 % dispEM(['Could not predict composition of ' model.metNames{mets(i)} ],false); 0106 %end 0107 end 0108 else 0109 %Check if there is an inconcistency 0110 %if any(currentComp~=comp) 0111 % dispEM(['Could not predict composition of ' model.metNames{loopThrough(i)} ],false); 0112 %end 0113 metStatus=-2; 0114 break; 0115 end 0116 end 0117 end 0118 %If there was contradictions in one compartment, then abort for 0119 %all compartments 0120 if metStatus==-2 0121 break; 0122 end 0123 else 0124 %The metabolite doesn't participate, no composition can be 0125 %calculated 0126 metStatus=-1; 0127 end 0128 end 0129 %Check status of the metabolite 0130 switch metStatus 0131 case -2 0132 dispEM(['Could not predict composition for "' metNames{i} '" due to inconsistencies'],false); 0133 case 1 0134 %Calculate and add the composition 0135 str=getCompString(elements,comp); 0136 model.metFormulas(mets)={str}; 0137 if printResults==true 0138 fprintf(['Predicted composition for "' metNames{i} '" to be ' str '\n']); 0139 end 0140 0141 %Keep track 0142 guessedFor=[guessedFor;metNames(i)]; 0143 0144 predicted=true; %To loop again 0145 end 0146 end 0147 end 0148 0149 couldNotGuess=setdiff(originalMissing,guessedFor); 0150 end 0151 0152 %Helper function for getting the composition string 0153 function str=getCompString(elements,comp) 0154 str=''; 0155 0156 for i=1:numel(comp) 0157 if comp(i)~=0 0158 if comp(i)==1 0159 str=[str elements.abbrevs{i}]; 0160 else 0161 str=[str elements.abbrevs{i} num2str(comp(i))]; 0162 end 0163 end 0164 end 0165 end