constructS

 constructS
   Constructs a stoichiometric matrix from a cell array of equations

   equations   cell array of equations on the form 'A + 2 B <=> 3 C',
               where <=> indicates reversible and => irreversible reactions
   mets        cell array of metabolites. All metabolites in the equations
               must be present in the list (optional, 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 (optional,
               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)

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 <=> 3 C',
0006 %               where <=> indicates reversible and => irreversible reactions
0007 %   mets        cell array of metabolites. All metabolites in the equations
0008 %               must be present in the list (optional, default generated from
0009 %               the equations)
0010 %   rxns        cell array of reaction ids. This is only used for printing
0011 %               reaction ids instead of equations in warnings/errors (optional,
0012 %               default [])
0013 %
0014 %   S           the resulting stoichiometric matrix mets cell array with
0015 %               metabolites that corresponds to the order in 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
0018 %               be the phosphotransferase ATP + ADP <=> ADP + ATP. In the
0019 %               stoichiometric matrix this equals to an empty reaction
0020 %               which 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 equations=convertCharArray(equations);
0026 switch nargin
0027     case 2
0028         mets=convertCharArray(mets);
0029     case 3
0030         rxns=convertCharArray(rxns);
0031 end
0032 
0033 badRxns=false(numel(equations),1);
0034 
0035 %Check that no equations are too short to have reversibility data
0036 I=cellfun(@numel,equations);
0037 I=find(I<4,1);
0038 if any(I)
0039     if isempty(rxns)
0040         EM=['The following equation does not have reversibility data: ' equations{I} ];
0041         dispEM(EM);
0042     else
0043         EM=['The reaction ' rxns{I} ' does not have reversibility data'];
0044         dispEM(EM);
0045     end
0046 end
0047 
0048 %Makes life a little easier
0049 equations=strtrim(equations);
0050 equations=fixEquations(equations);
0051 
0052 if nargin<2
0053     mets=parseRxnEqu(equations);
0054 end
0055 if nargin<3
0056     rxns=[];
0057 end
0058 
0059 %Get which reactions are reversible
0060 reversible=cellfun(@any,strfind(equations,' <=> '));
0061 
0062 %Make them all reversible. This is not all that neat, but nevermind
0063 equations=strrep(equations,' => ',' <=> ');
0064 
0065 %Replace the the plus signs with some weird character that will be used for
0066 %parsing
0067 equations=strrep(equations,' + ', '€');
0068 
0069 %Generate the stoichiometric matrix
0070 S=zeros(numel(mets),numel(equations));
0071 
0072 %Keep track of coefficients to be added to S-matrix
0073 metsToS = cell(100000,1);
0074 rxnsToS = zeros(100000,1);
0075 coefToS = zeros(100000,1);
0076 newEntry = 1;
0077 %Loop through the equations and add the info to the S matrix
0078 for i=1:numel(equations)
0079     %Start by finding the position of the (=> or <=>)
0080     arrowIndex=strfind(equations{i},' <=> ');
0081     
0082     if numel(arrowIndex)~=1
0083         if isempty(rxns)
0084             EM=['The following equation does not have reversibility data: ' equations{i} ];
0085             dispEM(EM);
0086         else
0087             EM=['The reaction ' rxns{i} ' does not have reversibility data'];
0088             dispEM(EM);
0089         end
0090     end
0091     
0092     reactants=regexp(equations{i}(1:arrowIndex-1),'€','split');
0093     products=regexp(equations{i}(arrowIndex+5:end),'€','split');
0094     
0095     %If the splitting character is at the end (if exchange rxns), then an
0096     %empty string will exist together with the real ones. Remove it
0097     reactants(cellfun(@isempty,reactants))=[];
0098     products(cellfun(@isempty,products))=[];
0099     
0100     %A vector where an element is -1 is the corresponding metabolite is a
0101     %reactant and 1 if it's a product
0102     multiplyWith=[ones(numel(reactants),1)*-1; ones(numel(products),1)];
0103     
0104     metabolites=strtrim([reactants products]);
0105     
0106     %Now loop through the reactants and see if the metabolite has a
0107     %coefficient (it will look as 'number name')
0108     for j=1:numel(metabolites)
0109         space=strfind(metabolites{j},' ');
0110         
0111         if isempty(space)
0112             %No coefficient
0113             coeff=1;
0114             name=metabolites{j};
0115         else
0116             coeff=str2double(metabolites{j}(1:space(1)));
0117             
0118             %If it was not a coefficiant
0119             if isnan(coeff)
0120                 coeff=1;
0121                 name=strtrim(metabolites{j});
0122             else
0123                 name=strtrim(metabolites{j}(space(1)+1:end));
0124             end
0125         end
0126         
0127         %Find the name in the mets list [a b]=ismember(name,mets);
0128         metsToS{newEntry}=name;
0129         rxnsToS(newEntry)=i;
0130         coefToS(newEntry)=coeff*multiplyWith(j);
0131         newEntry=newEntry+1;
0132     end
0133 end
0134 %Remove unused fields
0135 metsToS(newEntry:end)=[];
0136 rxnsToS(newEntry:end)=[];
0137 coefToS(newEntry:end)=[];
0138 
0139 %Match to mets array
0140 [metsPresent,metsLoc]=ismember(metsToS,mets);
0141 
0142 %Find badRxns
0143 [~,I]=unique([rxnsToS,metsLoc],'rows','stable');
0144 x=1:length(rxnsToS);
0145 x(I)=[];
0146 x=unique(rxnsToS(x));
0147 badRxns(x)=true;
0148 
0149 if any(~metsPresent)
0150     if isempty(rxns)
0151         error(['Could not find the following metabolites in the metabolite list: ',...
0152         strjoin(unique(metsToS(~metsPresent)),', ')],'')
0153     else
0154         missingMet = find(~metsPresent);
0155         missingMet = strcat(metsToS(missingMet),' (reaction:',rxns(rxnsToS(missingMet)),')\n');
0156         missingMet = strjoin(missingMet,'');
0157         error(['Could not find the following metabolites (reaction indicated) in the metabolite list: \n' ...
0158             missingMet '%s'],'');
0159     end
0160 end
0161 linearIndices=sub2ind(size(S),metsLoc,rxnsToS);
0162 S(linearIndices)=coefToS;
0163 S=sparse(S);
0164 end
0165 
0166 function equ=fixEquations(equ)
0167 %If the equation starts with "=>" or "<=>" then add a space again. This is
0168 %an alternative way to represent uptake reactions. The opposite way for
0169 %producing reactions
0170 equ=equ(:);
0171 for i=1:numel(equ)
0172     if strcmp(equ{i}(1:2),'=>') || strcmp(equ{i}(1:3),'<=>')
0173         equ{i}=[' ' equ{i}];
0174     else
0175         if strcmp(equ{i}(end-1:end),'=>') || strcmp(equ{i}(end-2:end),'<=>')
0176             equ{i}=[equ{i} ' '];
0177         end
0178     end
0179 end
0180 end