Home > RAVEN > constructS.m

constructS

PURPOSE ^

constructS

SYNOPSIS ^

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

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)

   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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Tue 16-Jul-2013 21:50:02 by m2html © 2005