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