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 (optional, 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)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 (optional, 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 if nargin<2 0037 printResults=true; 0038 end 0039 0040 %The metabolites for which there is no elemental composition 0041 originalMissing=unique(model.metNames(cellfun(@isempty,model.metFormulas))); 0042 guessedFor={}; 0043 0044 %This is to keep track of if new metabolite compositions were predicted 0045 predicted=true; 0046 while predicted==true 0047 predicted=false; 0048 0049 %Get the unique names (composition should be independent of 0050 %compartment) for the metabolites without composition 0051 metNames=unique(model.metNames(cellfun(@isempty,model.metFormulas))); 0052 0053 %Parse the formulas in the model 0054 [elements, useMat, exitFlag]=parseFormulas(model.metFormulas, true,false); 0055 0056 for i=1:numel(metNames) 0057 %Get the metabolites with this name. Not so neat, but this is a 0058 %fast function anyways 0059 mets=find(ismember(model.metNames,metNames(i))); 0060 0061 currentComp=[]; 0062 0063 %Loop through the metabolites -1: Could not assign due to missing 0064 %info. -2: Could not assign due to contradiction 1: Composition 0065 %assigned 0066 metStatus=-1; 0067 for j=1:numel(mets) 0068 %Get the reactions that the metabolite participates in 0069 [~, I]=find(model.S(mets(j),:)); 0070 if any(I) 0071 for k=1:numel(I) 0072 %Loop through the reactions and check if all other mets 0073 %in them have known composition 0074 eqn=model.S(:,I(k)); 0075 eqn(mets(j))=0; 0076 if all(exitFlag(eqn~=0)==1) 0077 %This means that all other mets had composition. 0078 %Calculate the resulting composition for the 0079 %unknown one 0080 comp=useMat'*eqn; 0081 0082 %This can result in round off errors if there are 0083 %stoichiometries with many decimals. Ignore values 0084 %below 10^-12 0085 comp(abs(comp)<10^-12)=0; 0086 0087 %Check if the composition consist of both negative 0088 %and positive values. If so, throw an error 0089 if all(comp<=0) || all(comp>=0) 0090 comp=abs(comp); 0091 if isempty(currentComp) 0092 currentComp=comp; 0093 end 0094 %Also to deal with round off errors 0095 if all(abs(currentComp-comp)<10^-10) 0096 metStatus=1; 0097 else 0098 metStatus=-2; 0099 break; 0100 0101 %%Check if there is an inconcistency 0102 %if any(currentComp~=comp) 0103 % dispEM(['Could not predict composition 0104 % of ' model.metNames{mets(i)} ],false); 0105 %end 0106 end 0107 else 0108 %Check if there is an inconcistency if 0109 %any(currentComp~=comp) 0110 % dispEM(['Could not predict composition of 0111 % ' 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 0119 %for 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 EM=['Could not predict composition for "' metNames{i} '" due to inconsistencies']; 0133 dispEM(EM,false); 0134 case 1 0135 %Calculate and add the composition 0136 str=getCompString(elements,comp); 0137 model.metFormulas(mets)={str}; 0138 if printResults==true 0139 fprintf(['Predicted composition for "' metNames{i} '" to be ' str '\n']); 0140 end 0141 0142 %Keep track 0143 guessedFor=[guessedFor;metNames(i)]; 0144 0145 predicted=true; %To loop again 0146 end 0147 end 0148 end 0149 0150 couldNotGuess=setdiff(originalMissing,guessedFor); 0151 end 0152 0153 %Helper function for getting the composition string 0154 function str=getCompString(elements,comp) 0155 str=''; 0156 0157 for i=1:numel(comp) 0158 if comp(i)~=0 0159 if comp(i)==1 0160 str=[str elements.abbrevs{i}]; 0161 else 0162 str=[str elements.abbrevs{i} num2str(comp(i))]; 0163 end 0164 end 0165 end 0166 end