Home > RAVEN > constructS.m

constructS

PURPOSE ^

constructS

SYNOPSIS ^

function [S mets badRxns reversible]=constructS(equations,mets,rxns)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Mon 06-Jan-2014 14:58:12 by m2html © 2005