findSExRxnInd

PURPOSE ^

Returns a model with boolean vectors indicating internal vs exchange/demand/sink reactions.

SYNOPSIS ^

function model=findSExRxnInd(model,nRealMet)

DESCRIPTION ^

Returns a model with boolean vectors indicating internal vs exchange/demand/sink reactions.

finds the reactions in the model which export/import from the model
boundary
e.g. Exchange reactions
     Demand reactions
     Sink reactions

INPUT
 model
 model.biomassRxnAbbr      abbreviation of biomass reaction 

OPTIONAL INPUT
 nRealMet                  specified in case extra rows in S which dont
                           correspond to metabolties
OUTPUT
 model.SIntRxnBool         Boolean of internal (mass balanced) reactions.

 OPTIONAL OUTPUT
 model.DMRxnBool           Boolean of demand reactions. Prefix 'DM_'
 model.SinkRxnBool         Boolean of sink reactions. Prefix 'sink_'

 Ronan M.T. Fleming

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function model=findSExRxnInd(model,nRealMet)
0002 %Returns a model with boolean vectors indicating internal vs exchange/demand/sink reactions.
0003 %
0004 %finds the reactions in the model which export/import from the model
0005 %boundary
0006 %e.g. Exchange reactions
0007 %     Demand reactions
0008 %     Sink reactions
0009 %
0010 %INPUT
0011 % model
0012 % model.biomassRxnAbbr      abbreviation of biomass reaction
0013 %
0014 %OPTIONAL INPUT
0015 % nRealMet                  specified in case extra rows in S which dont
0016 %                           correspond to metabolties
0017 %OUTPUT
0018 % model.SIntRxnBool         Boolean of internal (mass balanced) reactions.
0019 %
0020 % OPTIONAL OUTPUT
0021 % model.DMRxnBool           Boolean of demand reactions. Prefix 'DM_'
0022 % model.SinkRxnBool         Boolean of sink reactions. Prefix 'sink_'
0023 %
0024 % Ronan M.T. Fleming
0025 
0026 [nMet,nRxn]=size(model.S);
0027 
0028 if ~exist('nRealMet','var')
0029     nRealMet=length(model.mets);
0030     if nMet~=nRealMet
0031         fprintf('%s\n','Detected extra rows of S without corresponding metabolite abbreviations.')
0032     end
0033 end
0034 
0035 biomassBool=false(nRxn,1);
0036 
0037 %locate biomass reaction if there is one
0038 if ~isfield(model,'biomassRxnAbbr')
0039     if 0
0040         fprintf('%s\n','No model.biomassRxnAbbr ? Give abbreviation of biomass reaction if there is one.');
0041     else
0042         bool=model.c~=0;
0043         if nnz(bool)==1
0044             model.biomassRxnAbbr=model.rxns{model.c~=0};
0045             fprintf('%s%s\n','Assuming biomass reaction is: ', model.biomassRxnAbbr);
0046             biomassBool(bool)=1;
0047         else
0048             if nnz(bool)==0
0049                 fprintf('%s\n','No model.biomassRxnAbbr ? Give abbreviation of biomass reaction if there is one.');
0050             else
0051                 error('More than one biomass reaction?');
0052             end
0053         end
0054     end
0055 else
0056     bool=strcmp(model.biomassRxnAbbr,model.rxns);
0057     if nnz(bool)==1
0058         fprintf('%s%s\n','Found biomass reaction: ', model.biomassRxnAbbr);
0059         biomassBool(bool)=1;
0060     else
0061         if nnz(bool)==0
0062             fprintf('%s\n','No model.biomassRxnAbbr ? Give abbreviation of biomass reaction if there is one.');
0063         else
0064             error('More than one biomass reaction?');
0065         end
0066     end
0067 end
0068 
0069 if nMet > 2000
0070 
0071 
0072     % Human model or E.coli merged matrix
0073     model.ExchRxnBool=strncmp('Exch_', model.rxns, 5)==1;
0074     model.EXRxnBool=strncmp('EX_', model.rxns, 3)==1;
0075     %demand reactions going out of model
0076     model.DMRxnBool=strncmp('DM_', model.rxns, 3)==1;
0077     
0078     bool=strcmp('DM_atp(c)',model.rxns);
0079     if any(bool)
0080         fprintf('%s\n','ATP demand reaction is not considered an exchange reaction by default.')
0081         model.DMRxnBool(bool)=0;
0082     end
0083 
0084     %sink reactions going into or out of model
0085     model.SinkRxnBool=strncmp('sink_', model.rxns, 5)==1;
0086     
0087     %input/output
0088     SExRxnBool = model.ExchRxnBool | model.EXRxnBool | model.DMRxnBool | model.SinkRxnBool | biomassBool;
0089     
0090     %double check now by identifying reactions with only one metabolite
0091     SExRxnBool2=false(nRxn,1);
0092     for n=1:nRxn
0093         %find reactions with only one coefficient
0094         if nnz(model.S(1:nRealMet,n))==1
0095             SExRxnBool2(n)=1;
0096         end
0097     end
0098     SExRxnBool2 = SExRxnBool2 | biomassBool;
0099     
0100     diffBool= ~SExRxnBool & SExRxnBool2;
0101     if any(diffBool)
0102         fprintf('%s\n','Missed Exchanges:')
0103         fprintf('%s\t%s\t%s\t\t%s\t\t%s\n','Coefficient','Metabolite','#','Reaction','#')
0104         for n=1:nRxn
0105             if diffBool(n)
0106                 objMetInd=find(model.S(:,n));
0107                 for m=1:length(objMetInd)
0108                     Sij=model.S(objMetInd(m),n);
0109                     if length(model.mets{objMetInd(m)})<4
0110                         fprintf('%g\t\t\t%s\t\t\t%i\t%s\t\t%i\n',Sij,model.mets{objMetInd(m)},objMetInd(m),model.rxns{n},n)
0111                     else
0112                         if length(model.mets{objMetInd(m)})<8
0113                             fprintf('%g\t\t\t%s\t\t%i\t%s\t\t%i\n',Sij,model.mets{objMetInd(m)},objMetInd(m),model.rxns{n},n)
0114                         else
0115                             if length(model.mets{objMetInd(m)})<12
0116                                 fprintf('%g\t\t\t%s\t%i\t%s\t\t%i\n',Sij,model.mets{objMetInd(m)},objMetInd(m),model.rxns{n},n)
0117                             end
0118                         end
0119                     end
0120                 end
0121             end
0122         end
0123     end
0124     
0125     %dont check if there are coupling constraints
0126     %(E. coli E matrix specific)
0127     if ~isfield(model,'A')
0128         diffBool= SExRxnBool & ~SExRxnBool2;
0129         if any(diffBool)
0130             fprintf('%s\n','Exchanges (by prefix) with more than one coefficient:')
0131             fprintf('%s\t%s\n','#', 'Exchange')
0132             for n=1:length(diffBool)
0133                 if diffBool(n)
0134                     equation=printRxnFormula(model,model.rxns(n),0);
0135                     fprintf('%i\t%s\t%s\n',n,model.rxns{n},equation{1});
0136                 end
0137             end
0138         end
0139     end
0140     
0141 else
0142     SExRxnBool=false(nRxn,1);
0143     
0144     for n=1:nRxn
0145         %find reactions with only one coefficient
0146         if nnz(model.S(1:nRealMet,n))==1
0147             SExRxnBool(n,1)=1;
0148             if 0
0149                 if nonzeros(model.S(1:nRealMet,n))>0
0150                     fprintf('%s\t%s\n','Positive coefficient:',model.rxns{n});
0151                 else
0152                     fprintf('%s\t%s\n','Negative coefficient:',model.rxns{n});
0153                     %                 fprintf('%s%s%s%s%s\n','''',model.rxns{n},''',0 ,0 ''',model.mets{find(model.S(1:nRealMet,n)~=0)},''',0 ,0 ;');
0154                 end
0155             end
0156         end
0157     end
0158     if ~isempty(strcmp('ATPM',model.rxns))
0159         fprintf('%s\n','ATP maintenance reaction is not considered an exchange reaction by default.')
0160         ATPM_Ind=find(strcmp('ATPM',model.rxns)==1);
0161         SExRxnBool(ATPM_Ind,1)=0;
0162     end
0163     SExRxnBool(biomassBool)=1;
0164 end
0165 model.SIntRxnBool=~SExRxnBool;

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003