checkMassChargeBalance

PURPOSE ^

checkMassChargeBalance tests for a list of reactions if these reactions are

SYNOPSIS ^

function [massImbalance,imBalancedMass,imBalancedCharge,imBalancedBool,Elements] = checkMassChargeBalance(model,rxnBool,printLevel)

DESCRIPTION ^

checkMassChargeBalance tests for a list of reactions if these reactions are
mass-balanced by adding all elements on left hand side and comparing them
with the sums of elements on the right hand side of the reaction.

 [UnbalancedRxns] = checkMassChargeBalance(model,RxnList)

INPUT
 model                         COBRA model structure

OPTIONAL INPUT
 rxnBool       Boolean vector corresponding to reactions in model to be
               tested. If empty, then all tested.
               Alternatively, can be the indices of reactions to test:
               i.e. rxnBool(indixes)=1;
 printLevel    {-1,(0),1} 
               -1 = print out diagnostics on problem reactions to a file 
                0 = silent
                1 = print out diagnostics on problem reactions to screen

OUTPUTS
 massImbalance                 nRxn x nElement matrix with mass imblance
                               for each element checked. 0 if balanced.
 imBalancedMass                nRxn x 1 cell with charge imbalance
                               e.g. -3 H means three hydrogens disappear
                               in the reaction.
 imBalancedCharge              nRxn x 1 vector with charge imbalance,
                               empty if no imbalanced reactions

 imbalancedBool                boolean vector indicating imbalanced reactions
       
 Elements                      nElement x 1 cell array of element
                               abbreviations checked 
 Ines Thiele 12/09
 IT, 06/10, Corrected some bugs and improved speed.
 RF, 09/09/10, Support for very large models and printing to file.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [massImbalance,imBalancedMass,imBalancedCharge,imBalancedBool,Elements] = checkMassChargeBalance(model,rxnBool,printLevel)
0002 %checkMassChargeBalance tests for a list of reactions if these reactions are
0003 %mass-balanced by adding all elements on left hand side and comparing them
0004 %with the sums of elements on the right hand side of the reaction.
0005 %
0006 % [UnbalancedRxns] = checkMassChargeBalance(model,RxnList)
0007 %
0008 %INPUT
0009 % model                         COBRA model structure
0010 %
0011 %OPTIONAL INPUT
0012 % rxnBool       Boolean vector corresponding to reactions in model to be
0013 %               tested. If empty, then all tested.
0014 %               Alternatively, can be the indices of reactions to test:
0015 %               i.e. rxnBool(indixes)=1;
0016 % printLevel    {-1,(0),1}
0017 %               -1 = print out diagnostics on problem reactions to a file
0018 %                0 = silent
0019 %                1 = print out diagnostics on problem reactions to screen
0020 %
0021 %OUTPUTS
0022 % massImbalance                 nRxn x nElement matrix with mass imblance
0023 %                               for each element checked. 0 if balanced.
0024 % imBalancedMass                nRxn x 1 cell with charge imbalance
0025 %                               e.g. -3 H means three hydrogens disappear
0026 %                               in the reaction.
0027 % imBalancedCharge              nRxn x 1 vector with charge imbalance,
0028 %                               empty if no imbalanced reactions
0029 %
0030 % imbalancedBool                boolean vector indicating imbalanced reactions
0031 %
0032 % Elements                      nElement x 1 cell array of element
0033 %                               abbreviations checked
0034 % Ines Thiele 12/09
0035 % IT, 06/10, Corrected some bugs and improved speed.
0036 % RF, 09/09/10, Support for very large models and printing to file.
0037 
0038 [nMet,nRxn]=size(model.S);
0039 if exist('rxnBool','var')
0040     if ~isempty(rxnBool)
0041         if length(rxnBool)~=nRxn
0042             rxnBool2=false(nRxn,1);
0043             rxnBool2(rxnBool)=1;
0044             rxnBool=rxnBool2;
0045         end
0046     else
0047         model=findSExRxnInd(model);
0048         %only check mass balance of internal reactions
0049         rxnBool=model.SIntRxnBool;
0050     end
0051 else
0052     model=findSExRxnInd(model);
0053     %only check mass balance of internal reactions
0054     rxnBool=model.SIntRxnBool;
0055 end
0056 if ~exist('printLevel','var')
0057     printLevel=0;
0058 end
0059 
0060 % List of Elements
0061 Elements = {'H','C', 'O', 'P', 'S', 'N', 'Mg','X','Fe','Zn','Co','R'};
0062 
0063 E=sparse(nMet,length(Elements));
0064 massImbalance=sparse(nRxn,length(Elements));
0065 for j = 1 : length(Elements)
0066     if j==1
0067         [dE,E_el]=checkBalance(model,Elements{j},printLevel);
0068         massImbalance(:,j)=dE;
0069         E(:,j)=E_el;
0070         fprintf('%s\n',['Checked element ' Elements{j}]);  
0071     else
0072         %no need to print out for each element which metabolites have no
0073         %formula
0074         [massImbalance(:,j),E(:,j)]=checkBalance(model,Elements{j},0);
0075         fprintf('%s\n',['Checking element ' Elements{j}]);
0076     end
0077 end
0078 massImbalance(~rxnBool,:)=0;
0079 imBalancedBool=sum(abs(massImbalance'))'~=0;
0080 
0081 imBalancedBool=rxnBool & imBalancedBool;
0082 
0083 imBalancedMass=cell(nRxn,1);
0084 for i = 1 : nRxn
0085     imBalancedMass{i,1}='';   
0086     if imBalancedBool(i)
0087         for j = 1 : length(Elements)
0088             if massImbalance(i,j)~=0
0089                 if ~strcmp(imBalancedMass{i,1},'')
0090                     imBalancedMass{i,1} = [imBalancedMass{i,1} ', ' int2str(massImbalance(i,j)) ' ' Elements{j}];
0091                 else
0092                     imBalancedMass{i,1} = [int2str(massImbalance(i,j)) ' ' Elements{j}];
0093                 end
0094             end
0095             
0096         end
0097         if strfind(imBalancedMass{i,1},'NaN')
0098             imBalancedMass{i,1}='NaN';
0099         end
0100     end
0101     if mod(i,1000)==0
0102         fprintf('%n\t%s\n',i,['reactions checked for ' Elements{j} ' balance']);
0103     end
0104 end
0105 if printLevel==-1
0106     firstMissing=0;
0107     for p=1:nRxn
0108         if ~strcmp(imBalancedMass{p,1},'')
0109             %at the moment, ignore reactions with a metabolite that have
0110             %no formula
0111             if ~strcmp(imBalancedMass{p,1},'NaN')
0112                 if ~firstMissing
0113                     fid=fopen('mass_imbalanced_reactions.txt','w');
0114                     fprintf(fid,'%s;%s;%s;%s\n','#Rxn','rxnAbbr','imbalance','equation');
0115 
0116                     warning('There are mass imbalanced reactions, see mass_imbalanced_reactions.txt')
0117                     firstMissing=1;
0118                 end
0119                 equation=printRxnFormula(model,model.rxns(p),0);
0120                 fprintf(fid,'%s;%s;%s;%s\n',int2str(p),model.rxns{p},imBalancedMass{p,1},equation{1});
0121                 for m=1:size(model.S,1)
0122                     if model.S(m,p)~=0
0123                         fprintf(fid,'%s\t%s\t%s\t%s\t%s\n',int2str(m),model.mets{m},int2str(model.S(m,p)),int2str(E(m)),model.metFormulas{m});
0124                     end
0125                 end
0126             end
0127         end
0128     end
0129     if firstMissing
0130         fclose(fid);
0131     end
0132 end
0133 if printLevel==1
0134     for p=1:nRxn
0135         if ~strcmp(imBalancedMass{p,1},'')
0136             %at the moment, ignore reactions with a metabolite that have
0137             %no formula
0138             if ~strcmp(imBalancedMass{p,1},'NaN')
0139                 equation=printRxnFormula(model,model.rxns(p),0);
0140                 fprintf('%6s\t%30s\t%10s\t%s\n',int2str(p),model.rxns{p},imBalancedMass{p,1},equation{1});
0141                 if 0
0142                 for m=1:size(model.S,1)
0143                     if model.S(m,p)~=0
0144                         fprintf(fid,'%s\t%s\t%s\t%s\t%s\n',int2str(m),model.mets{m},int2str(model.S(m,p)),int2str(E(m)),model.metFormulas{m});
0145                     end
0146                 end
0147                 end
0148             end
0149         end
0150     end
0151 end
0152 
0153 %
0154 if nnz(strcmp('',imBalancedMass))==nRxn
0155     imBalancedMass=[];
0156 end
0157 
0158 % Check for charge balance
0159 imBalancedCharge=[];
0160 firstMissing=0;
0161 if isfield(model, 'metCharges')
0162     for m=1:nMet
0163         if isnan(model.metCharges(m)) && ~isempty(model.metFormulas{m})
0164             if printLevel==1
0165                 fprintf('%s\t%s\n',int2str(m),[model.mets{m} ' has no charge but has formula.'])
0166                 if ~firstMissing
0167                     warning('model structure must contain model.metCharges field for each metabolite');
0168                 end
0169                 firstMissing=1;
0170             end
0171             if printLevel==-1
0172                 if ~firstMissing
0173                     fid=fopen('metabolites_without_charge.txt','w');
0174                 end
0175                 firstMissing=1;
0176                 fprintf(fid,'%s\t%s\n',int2str(m),model.mets{m})
0177             end
0178         else
0179             dC=model.S'*model.metCharges;
0180         end
0181     end
0182     if any(dC(rxnBool))~=0
0183         imBalancedCharge=dC;
0184         imBalancedCharge(~rxnBool)=0;
0185     else
0186         imBalancedCharge=[];
0187     end
0188 end
0189 
0190 if printLevel==-1
0191     firstMissing=0;
0192     if ~isempty(imBalancedCharge)
0193         for q=1:nRxn
0194             if model.SIntRxnBool(q) && dC(q)~=0 && strcmp(imBalancedMass{p,1},'')
0195                 if ~firstMissing
0196                     fid=fopen('charge_imbalanced_reactions.txt','w');
0197                     warning('There are charged imbalanced reactions (that are mass balanced), see charge_imbalanced_reactions.txt')
0198                     firstMissing=1;
0199                 end
0200                 equation=printRxnFormula(model,model.rxns(q),0);
0201                 fprintf(fid,'%s\t%s\t%s\n',int2str(q),model.rxns{q},equation{1});
0202                 if 0
0203                     for m=1:size(model.S,1)
0204                         if model.S(m,q)~=0
0205                             fprintf(fid,'%s\t%15s\t%3s\t%3s\t%s\n',int2str(m),model.mets{m},int2str(model.S(m,q)),int2str(model.metCharges(m)),model.metFormulas{m});
0206                         end
0207                     end
0208                 end
0209             end
0210         end
0211         if firstMissing
0212             fclose(fid);
0213         end
0214     end
0215 end
0216 
0217 if printLevel==1
0218     if ~isempty(imBalancedCharge)
0219         fprintf('%s\n','Mass balanced, but charged imbalanced reactions:')
0220         for q=1:nRxn
0221             if model.SIntRxnBool(q) && dC(q)~=0 && strcmp(imBalancedMass{p,1},'')
0222                 equation=printRxnFormula(model,model.rxns(q),0);
0223                 fprintf('%s\t%s\t%s\n',int2str(q),model.rxns{q},equation{1});
0224                 if 1
0225                     for m=1:size(model.S,1)
0226                         if model.S(m,q)~=0
0227                             fprintf('%s\t%15s\t%3s\t%3s\t%s\n',int2str(m),model.mets{m},int2str(model.S(m,q)),int2str(model.metCharges(m)),model.metFormulas{m});
0228                         end
0229                     end
0230                 end
0231             end
0232         end
0233     end
0234 end
0235 
0236 if ~isempty(imBalancedCharge)
0237     imBalancedBool = imBalancedBool |  imBalancedCharge~=0;
0238 end
0239 
0240 
0241 
0242 
0243

Generated on Sun 26-Dec-2010 12:06:07 by m2html © 2003