printFluxes

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

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

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