printFluxes

 printFluxes
   Prints reactions and fluxes to the screen or to a file

 Input:
   model           a model structure
   fluxes          a vector with fluxes
   onlyExchange    only print exchange fluxes (optional, default true)
   cutOffFlux      only print fluxes with absolute values above or equal
                   to this value (optional, default 10^-8)
   outputFile      a file to save the print-out to (optional, default is
                   output to the command window)
   outputString    a string that specifies the output of each reaction
                   (optional, default '%rxnID\t(%rxnName):\t%flux\n')
   metaboliteList  cell array of metabolite names. Only reactions
                   involving any of these metabolites will be
                   printed (optional)

 The following codes are available for user-defined output strings:

 %rxnID      reaction ID
 %rxnName    reaction name
 %lower      lower bound
 %upper      upper bound
 %obj        objective coefficient
 %eqn        equation
 %flux       flux
 %element    equation using the metabolite formulas rather than
             metabolite names
 %unbalanced "(*)" if the reaction is unbalanced and "(-)" if it could
             not be parsed
 %lumped     equation where the elemental compositions for the left/right
             hand sides are lumped

 Usage: printFluxes(model, fluxes, onlyExchange, cutOffFlux, outputFile,...
                    outputString, metaboliteList)

0001 function printFluxes(model, fluxes, onlyExchange, cutOffFlux, outputFile,outputString,metaboliteList)
0002 % printFluxes
0003 %   Prints reactions and fluxes to the screen or to a file
0004 %
0005 % Input:
0006 %   model           a model structure
0007 %   fluxes          a vector with fluxes
0008 %   onlyExchange    only print exchange fluxes (optional, default true)
0009 %   cutOffFlux      only print fluxes with absolute values above or equal
0010 %                   to this value (optional, default 10^-8)
0011 %   outputFile      a file to save the print-out to (optional, default is
0012 %                   output to the command window)
0013 %   outputString    a string that specifies the output of each reaction
0014 %                   (optional, default '%rxnID\t(%rxnName):\t%flux\n')
0015 %   metaboliteList  cell array of metabolite names. Only reactions
0016 %                   involving any of these metabolites will be
0017 %                   printed (optional)
0018 %
0019 % The following codes are available for user-defined output strings:
0020 %
0021 % %rxnID      reaction ID
0022 % %rxnName    reaction name
0023 % %lower      lower bound
0024 % %upper      upper bound
0025 % %obj        objective coefficient
0026 % %eqn        equation
0027 % %flux       flux
0028 % %element    equation using the metabolite formulas rather than
0029 %             metabolite names
0030 % %unbalanced "(*)" if the reaction is unbalanced and "(-)" if it could
0031 %             not be parsed
0032 % %lumped     equation where the elemental compositions for the left/right
0033 %             hand sides are lumped
0034 %
0035 % Usage: printFluxes(model, fluxes, onlyExchange, cutOffFlux, outputFile,...
0036 %                    outputString, metaboliteList)
0037 
0038 if nargin<3
0039     onlyExchange=true;
0040 end
0041 if nargin<4
0042     cutOffFlux=10^-8;
0043 end
0044 if isempty(cutOffFlux)
0045     cutOffFlux=10^-8;
0046 end
0047 if nargin<5
0048     fid=1;
0049 else
0050     if ~isempty(outputFile)
0051         outputFile=char(outputFile);
0052         fid=fopen(outputFile,'w');
0053     else
0054         fid=1;
0055     end
0056 end
0057 if nargin<6 || isempty(outputString)
0058     outputString='%rxnID\t(%rxnName):\t%flux\n';
0059 else
0060     outputString=char(outputString);
0061 end
0062 if nargin<7
0063     metaboliteList={};
0064 else
0065     metaboliteList=convertCharArray(metaboliteList);
0066 end
0067 if isempty(fluxes)
0068     EM='Empty vector of fluxes, solveLP possibly returned infeasible';
0069     dispEM(EM);
0070 elseif size(fluxes,1)~=numel(model.rxns)
0071     EM='The number of fluxes and the number of reactions must be the same';
0072     dispEM(EM);
0073 end
0074 
0075 %Only keep reactions involving the defined metabolites
0076 if ~isempty(metaboliteList)
0077     I=ismember(upper(model.metNames),upper(metaboliteList));
0078     [~, K]=find(model.S(I,:));
0079     
0080     %Delete all other reactions
0081     toDelete=true(numel(model.rxns),1);
0082     toDelete(K)=false;
0083     model=removeReactions(model,toDelete);
0084     fluxes(toDelete,:)=[];
0085 end
0086 
0087 if onlyExchange==true
0088     fprintf(fid,'EXCHANGE FLUXES:\n');
0089 else
0090     fprintf(fid,'FLUXES:\n');
0091 end
0092 
0093 %Remove reactions which are below the cut off
0094 toDelete=abs(fluxes)<cutOffFlux;
0095 toDelete=all(toDelete,2);
0096 model=removeReactions(model,toDelete,true,true);
0097 fluxes(toDelete,:)=[];
0098 
0099 if any(strfind(outputString,'%eqn'))
0100     %Construct the equations
0101     eqn=constructEquations(model);
0102 else
0103     eqn=cell(numel(model.rxns),1);
0104     eqn(:)={''};
0105 end
0106 if any(strfind(outputString,'%element'))
0107     %For printing equations using the composition
0108     cModel=model;
0109     cModel.metNames=cModel.metFormulas;
0110     cModel.metNames(cellfun(@isempty,cModel.metNames))={'?'};
0111     element=constructEquations(cModel);
0112 else
0113     element=cell(numel(model.rxns),1);
0114     element(:)={''};
0115 end
0116 
0117 if any(strfind(outputString,'%unbalanced')) || any(strfind(outputString,'%lumped'))
0118     balanceStructure=getElementalBalance(model);
0119 end
0120 
0121 unbalanced=cell(numel(model.rxns),1);
0122 unbalanced(:)={''};
0123 if any(strfind(outputString,'%unbalanced'))
0124     unbalanced(balanceStructure.balanceStatus==0)={'(*)'};
0125     unbalanced(balanceStructure.balanceStatus<0)={'(-)'};
0126 end
0127 
0128 lumped=cell(numel(model.rxns),1);
0129 lumped(:)={''};
0130 if any(strfind(outputString,'%lumped'))
0131     for i=1:numel(model.rxns)
0132         leftGroup='';
0133         rightGroup='';
0134         for j=1:numel(balanceStructure.elements.names)
0135             I=balanceStructure.leftComp(i,j);
0136             if I~=0
0137                 if I==1
0138                     leftGroup=[leftGroup balanceStructure.elements.abbrevs{j}];
0139                 else
0140                     leftGroup=[leftGroup balanceStructure.elements.abbrevs{j} num2str(I)];
0141                 end
0142             end
0143             I=balanceStructure.rightComp(i,j);
0144             if I~=0
0145                 if I==1
0146                     rightGroup=[rightGroup balanceStructure.elements.abbrevs{j}];
0147                 else
0148                     rightGroup=[rightGroup balanceStructure.elements.abbrevs{j} num2str(I)];
0149                 end
0150             end
0151         end
0152         if model.rev(i)
0153             lumped{i}=[leftGroup ' <=> ' rightGroup];
0154         else
0155             lumped{i}=[leftGroup ' => ' rightGroup];
0156         end
0157     end
0158 end
0159 
0160 for i=1:numel(model.rxns)
0161     %Only print if it's an exchange reaction or if all reactions should be
0162     %printed. Exchange reactions only have reactants or only products.
0163     reactants=model.S(:,i)<0;
0164     products=model.S(:,i)>0;
0165     
0166     %Only print if the absolute value is >= cutOffFlux
0167     if (onlyExchange==false || (~any(reactants) || ~any(products)))
0168         printString=outputString;
0169         
0170         %Produce the final string
0171         printString=strrep(printString,'%rxnID',model.rxns{i});
0172         printString=strrep(printString,'%eqn',eqn{i});
0173         printString=strrep(printString,'%rxnName',model.rxnNames{i});
0174         printString=strrep(printString,'%lower',num2str(model.lb(i)));
0175         printString=strrep(printString,'%upper',num2str(model.ub(i)));
0176         printString=strrep(printString,'%obj',num2str(model.c(i)));
0177         printString=strrep(printString,'%flux',num2str(fluxes(i,:)));
0178         printString=strrep(printString,'%element',element{i});
0179         printString=strrep(printString,'%unbalanced',unbalanced{i});
0180         printString=strrep(printString,'%lumped',lumped{i});
0181         fprintf(fid,printString);
0182     end
0183 end
0184 
0185 if fid~=1
0186     fprintf('File successfully saved.\n');
0187     fclose(fid);
0188 end
0189 end