Home > RAVEN > getElementalBalance.m

getElementalBalance

PURPOSE ^

getElementalBalance

SYNOPSIS ^

function balanceStructure=getElementalBalance(model,rxns,printUnbalanced,printUnparsable)

DESCRIPTION ^

 getElementalBalance
   Checks a model to see if the reactions are elementally balanced

   model             a model structure
   rxns              either a cell array of reaction IDs, a logical vector 
                     with the same number of elements as reactions in the model,
                     of a vector of indexes. Only these reactions will be
                     checked (opt, default model.rxns)
   printUnbalanced   print warnings about the reactions that were 
                     unbalanced (opt, default false)
   printUnparsable   print warnings about the reactions that cannot be
                     parsed (opt, default false)

   balanceStructure
       balanceStatus    1 if the reaction is balanced, 0 if it's unbalanced, 
                      -1 if it couldn't be balanced due to missing information, 
                      -2 if it couldn't be balanced due to an error
       leftComp        An array with the sum of coefficients for carbon, phosphate, 
                       sulfur, nitrogen, "R-groups", oxygen, hydrogen for the left side
       rightComp       The corresponding array for the right side

   Usage: balanceStructure=getElementalBalance(model,rxns,printUnbalanced,printUnparsable)

   Rasmus Agren, 2013-02-20

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function balanceStructure=getElementalBalance(model,rxns,printUnbalanced,printUnparsable)
0002 % getElementalBalance
0003 %   Checks a model to see if the reactions are elementally balanced
0004 %
0005 %   model             a model structure
0006 %   rxns              either a cell array of reaction IDs, a logical vector
0007 %                     with the same number of elements as reactions in the model,
0008 %                     of a vector of indexes. Only these reactions will be
0009 %                     checked (opt, default model.rxns)
0010 %   printUnbalanced   print warnings about the reactions that were
0011 %                     unbalanced (opt, default false)
0012 %   printUnparsable   print warnings about the reactions that cannot be
0013 %                     parsed (opt, default false)
0014 %
0015 %   balanceStructure
0016 %       balanceStatus    1 if the reaction is balanced, 0 if it's unbalanced,
0017 %                      -1 if it couldn't be balanced due to missing information,
0018 %                      -2 if it couldn't be balanced due to an error
0019 %       leftComp        An array with the sum of coefficients for carbon, phosphate,
0020 %                       sulfur, nitrogen, "R-groups", oxygen, hydrogen for the left side
0021 %       rightComp       The corresponding array for the right side
0022 %
0023 %   Usage: balanceStructure=getElementalBalance(model,rxns,printUnbalanced,printUnparsable)
0024 %
0025 %   Rasmus Agren, 2013-02-20
0026 %
0027 
0028 if nargin<2
0029     rxns=[];
0030 end
0031 
0032 if nargin<3
0033     printUnbalanced=false;
0034 end
0035 
0036 if nargin<4
0037     printUnparsable=false;
0038 end
0039 
0040 if ~isempty(rxns)
0041    indexes=~getIndexes(model,rxns,'rxns',true);
0042    model=removeRxns(model,indexes);
0043 end
0044 
0045 %Get the formulas. Formulas should be parsed from the InChIs when defined,
0046 %otherwise the metFormulas values are used.
0047 if isfield(model,'metFormulas')
0048     %Get the elemental composition of all metabolites
0049     [nC, nP, nS, nN, nR, nO, nH, exitFlag]=parseFormulas(model.metFormulas,true);
0050 else
0051     nC=zeros(numel(model.mets),1);
0052     nP=nC;
0053     nS=nC;
0054     nN=nC;
0055     nR=nC;
0056     nO=nC;
0057     nH=nC;
0058     exitFlag=nC;
0059 end
0060 if isfield(model,'inchis')
0061     [nCI, nPI, nSI, nNI, nRI, nOI, nHI, exitFlagI]=parseFormulas(model.inchis,true,true);
0062     
0063     %Use the values calculated from the InChI codes if successful
0064     nC(exitFlagI==1)=nCI(exitFlagI==1);
0065     nP(exitFlagI==1)=nPI(exitFlagI==1);
0066     nS(exitFlagI==1)=nSI(exitFlagI==1);
0067     nN(exitFlagI==1)=nNI(exitFlagI==1);
0068     nR(exitFlagI==1)=nRI(exitFlagI==1);
0069     nO(exitFlagI==1)=nOI(exitFlagI==1);
0070     nH(exitFlagI==1)=nHI(exitFlagI==1);
0071     exitFlag(exitFlagI==1)=exitFlagI(exitFlagI==1);
0072 end
0073 
0074 balanceStructure.balanceStatus=zeros(numel(model.rxns),1);
0075 balanceStructure.leftComp=zeros(numel(model.rxns),7);
0076 balanceStructure.rightComp=zeros(numel(model.rxns),7);
0077 
0078 %Loop through each reaction and calculate the total amount of each compound
0079 for i=1:numel(model.rxns)
0080     leftMets=find(model.S(:,i)<0);
0081     rightMets=find(model.S(:,i)>0);
0082     
0083     %Calculate the composition of left/right side even if the reaction as a
0084     %whole cannot be balanced
0085     balanceStructure.leftComp(i,1)=abs(sum(model.S(leftMets,i).*nC(leftMets)));
0086     balanceStructure.leftComp(i,2)=abs(sum(model.S(leftMets,i).*nP(leftMets)));
0087     balanceStructure.leftComp(i,3)=abs(sum(model.S(leftMets,i).*nS(leftMets)));
0088     balanceStructure.leftComp(i,4)=abs(sum(model.S(leftMets,i).*nN(leftMets)));
0089     balanceStructure.leftComp(i,5)=abs(sum(model.S(leftMets,i).*nR(leftMets)));
0090     balanceStructure.leftComp(i,6)=abs(sum(model.S(leftMets,i).*nO(leftMets)));
0091     balanceStructure.leftComp(i,7)=abs(sum(model.S(leftMets,i).*nH(leftMets)));
0092 
0093     balanceStructure.rightComp(i,1)=abs(sum(model.S(rightMets,i).*nC(rightMets)));
0094     balanceStructure.rightComp(i,2)=abs(sum(model.S(rightMets,i).*nP(rightMets)));
0095     balanceStructure.rightComp(i,3)=abs(sum(model.S(rightMets,i).*nS(rightMets)));
0096     balanceStructure.rightComp(i,4)=abs(sum(model.S(rightMets,i).*nN(rightMets)));
0097     balanceStructure.rightComp(i,5)=abs(sum(model.S(rightMets,i).*nR(rightMets)));
0098     balanceStructure.rightComp(i,6)=abs(sum(model.S(rightMets,i).*nO(rightMets)));
0099     balanceStructure.rightComp(i,7)=abs(sum(model.S(rightMets,i).*nH(rightMets)));
0100     
0101     %Are they ok?
0102     if all(exitFlag([leftMets; rightMets]))
0103         if all(exitFlag([leftMets; rightMets])==1)             
0104              %Check if it's balanced (use a cut off for fractional
0105              %coefficients)
0106              if sum(abs(balanceStructure.leftComp(i,:)-balanceStructure.rightComp(i,:)))<10^-6
0107                  balanceStructure.balanceStatus(i)=1;
0108              end
0109         else
0110             balanceStructure.balanceStatus(i)=-2;
0111         end
0112     else
0113         balanceStructure.balanceStatus(i)=-1;
0114     end
0115 end
0116 
0117 %Print warnings
0118 toPrint=[];
0119 if printUnbalanced==true
0120    toPrint=[toPrint;find(balanceStructure.balanceStatus==0)];
0121 end
0122 if printUnparsable==true
0123    toPrint=[toPrint;find(balanceStructure.balanceStatus<0)];
0124 end
0125 
0126 compounds={'carbon' 'phosphate' 'sulfate' 'nitrogen' 'R-groups' 'oxygen' 'hydrogen'};
0127 toPrint=sort(toPrint);
0128 for i=1:numel(toPrint)
0129     if balanceStructure.balanceStatus(toPrint(i))<0
0130         fprintf(['WARNING: The reaction ' model.rxns{toPrint(i)} ' could not be parsed\n']);
0131     else
0132        %Find the compounds that it's not balanced for
0133        notBalanced=find(abs(balanceStructure.leftComp(toPrint(i),:)-balanceStructure.rightComp(toPrint(i),:))>10^-6);
0134        for j=1:numel(notBalanced)
0135             fprintf(['WARNING: The reaction ' model.rxns{toPrint(i)} ' is not balanced with respect to ' compounds{notBalanced(j)} '\n']);
0136        end
0137     end
0138 end
0139 
0140 end

Generated on Tue 23-Apr-2013 15:18:37 by m2html © 2005