constructS Constructs a stoichiometric matrix from a cell array of equations equations cell array of equations on the form 'A + 2 B <=> or => 3 C' mets cell array of metabolites. All metabolites in the equations must be present in the list (opt, default generated from the equations) S the resulting stoichiometric matrix mets cell array with metabolites that corresponds to the order in the S matrix badRxns boolean vector with the reactions that have one or more metabolites as both substrate and product. An example would be the phosphotransferase ATP + ADP <=> ADP + ATP. In the stoichiometric matrix this equals to an empty reaction which can be problematic reversible boolean vector with true if the equation is reversible Usage: [S mets badRxns reversible]=constructS(equations,mets) Rasmus Agren, 2012-10-31
0001 function [S mets badRxns reversible]=constructS(equations,mets) 0002 % constructS 0003 % Constructs a stoichiometric matrix from a cell array of equations 0004 % 0005 % equations cell array of equations on the form 'A + 2 B <=> or => 3 C' 0006 % mets cell array of metabolites. All metabolites in the equations 0007 % must be present in the list (opt, default generated from the equations) 0008 % 0009 % S the resulting stoichiometric matrix 0010 % mets cell array with metabolites that corresponds to the order in 0011 % the S matrix 0012 % badRxns boolean vector with the reactions that have one or more 0013 % metabolites as both substrate and product. An example would be 0014 % the phosphotransferase ATP + ADP <=> ADP + ATP. In the 0015 % stoichiometric matrix this equals to an empty reaction which 0016 % can be problematic 0017 % reversible boolean vector with true if the equation is reversible 0018 % 0019 % Usage: [S mets badRxns reversible]=constructS(equations,mets) 0020 % 0021 % Rasmus Agren, 2012-10-31 0022 % 0023 0024 badRxns=false(numel(equations),1); 0025 0026 %Makes life a little easier 0027 equations=strtrim(equations); 0028 equations=fixEquations(equations); 0029 0030 if nargin<2 0031 mets=parseRxnEqu(equations); 0032 end 0033 0034 %Get which reactions are reversible 0035 reversible=cellfun(@any,strfind(equations,' <=> ')); 0036 0037 %Make them all reversible. This is not all that neat, but nevermind 0038 equations=strrep(equations,' => ',' <=> '); 0039 0040 %Replace the the plus signs with a weird character that will be used for 0041 %parsing 0042 equations=strrep(equations,' + ', '¤'); 0043 0044 %Generate the stoichiometric matrix 0045 S=sparse(numel(mets),numel(equations)); 0046 0047 %Loop through the equations and add the info to the S matrix 0048 for i=1:numel(equations) 0049 %Start by finding the position of the (=> or <=>) 0050 arrowIndex=strfind(equations{i},' <=> '); 0051 0052 if numel(arrowIndex)~=1 0053 throw(MException('',['Error in equation: ' equations{i}])); 0054 end 0055 0056 reactants=regexp(equations{i}(1:arrowIndex-1),'¤','split'); 0057 products=regexp(equations{i}(arrowIndex+5:end),'¤','split'); 0058 0059 %If the splitting character is at the end (if exchange rxns), then an 0060 %empty string will exist together with the real ones. Remove it 0061 reactants(cellfun(@isempty,reactants))=[]; 0062 products(cellfun(@isempty,products))=[]; 0063 0064 %A vector where an element is -1 is the corresponding metabolite is a 0065 %reactant and 1 if it's a product 0066 multiplyWith=[ones(numel(reactants),1)*-1; ones(numel(products),1)]; 0067 0068 metabolites=[reactants products]; 0069 0070 %Now loop through the reactants and see if the metabolite has a coefficient 0071 %(it will look as 'number name') 0072 for j=1:numel(metabolites) 0073 space=strfind(metabolites{j},' '); 0074 0075 if isempty(space) 0076 %No coefficient 0077 coeff=1; 0078 name=metabolites{j}; 0079 else 0080 coeff=str2double(metabolites{j}(1:space(1))); 0081 0082 %If it was not a coefficiant 0083 if isnan(coeff) 0084 coeff=1; 0085 name=metabolites{j}; 0086 else 0087 name=metabolites{j}(space+1:end); 0088 end 0089 end 0090 0091 %Find the name in the mets list 0092 [a b]=ismember(name,mets); 0093 0094 if any(a) 0095 %Check if the metabolite already participates in this reaction 0096 if S(b,i)~=0 0097 badRxns(i)=true; 0098 end 0099 S(b,i)=S(b,i)+coeff*multiplyWith(j); 0100 else 0101 throw(MException('',['Could not find metabolite ' name ' in metabolite list'])); 0102 end 0103 end 0104 end 0105 0106 end 0107 0108 function equ=fixEquations(equ) 0109 %If the equation starts with "=>" or "<=>" then add a space again. This is 0110 %an alternative way to represent uptake reactions. The opposite way for 0111 %producing reactions 0112 equ=equ(:); 0113 for i=1:numel(equ) 0114 if strcmp(equ{i}(1:2),'=>') || strcmp(equ{i}(1:3),'<=>') 0115 equ{i}=[' ' equ{i}]; 0116 else 0117 if strcmp(equ{i}(end-1:end),'=>') || strcmp(equ{i}(end-2:end),'<=>') 0118 equ{i}=[equ{i} ' ']; 0119 end 0120 end 0121 end 0122 end