followChanged

 followChanged
    Prints fluxes and reactions for each of the reactions that results in
   different fluxes compared to the reference case.

   model           a model structure
   fluxesA         flux vector for the test case
   fluxesB         flux vector for the reference test
   cutOffChange    reactions where the fluxes differ by less than
                   this many percent won't be printed (optional, default 10^-8)
   cutOffFlux      reactions where the absolute value of both fluxes
                   are below this value won't be printed (optional,
                   default 10^-8)
   cutOffDiff      reactions where the fluxes differ by less than
                   cutOffDiff won't be printed (optional, default 10^-8)
   metaboliteList  cell array of metabolite names. Only reactions
                   involving any of these metabolites will be
                   printed (optional)

 Usage: followChanged(model,fluxesA,fluxesB, cutOffChange, cutOffFlux,
           cutOffDiff, metaboliteList)

0001 function followChanged(model,fluxesA,fluxesB, cutOffChange, cutOffFlux, cutOffDiff, metaboliteList)
0002 % followChanged
0003 %    Prints fluxes and reactions for each of the reactions that results in
0004 %   different fluxes compared to the reference case.
0005 %
0006 %   model           a model structure
0007 %   fluxesA         flux vector for the test case
0008 %   fluxesB         flux vector for the reference test
0009 %   cutOffChange    reactions where the fluxes differ by less than
0010 %                   this many percent won't be printed (optional, default 10^-8)
0011 %   cutOffFlux      reactions where the absolute value of both fluxes
0012 %                   are below this value won't be printed (optional,
0013 %                   default 10^-8)
0014 %   cutOffDiff      reactions where the fluxes differ by less than
0015 %                   cutOffDiff won't be printed (optional, default 10^-8)
0016 %   metaboliteList  cell array of metabolite names. Only reactions
0017 %                   involving any of these metabolites will be
0018 %                   printed (optional)
0019 %
0020 % Usage: followChanged(model,fluxesA,fluxesB, cutOffChange, cutOffFlux,
0021 %           cutOffDiff, metaboliteList)
0022 
0023 %Checks if a cut off flux has been set
0024 if nargin<4
0025     cutOffChange=10^-8;
0026 end
0027 if nargin<5
0028     cutOffFlux=10^-8;
0029 end
0030 if nargin<6
0031     cutOffDiff=10^-8;
0032 end
0033 if nargin<7
0034     metaboliteList=[];
0035 else
0036     metaboliteList=convertCharArray(metaboliteList);
0037 end
0038 
0039 %If a metabolite list is to be used, then find all the reactions involving
0040 %any of those metabolites Finds the metabolites
0041 if nargin>6
0042     reactionIndexes=[];
0043     for i=1:length(metaboliteList)
0044         metaboliteIndex=find(strcmpi(metaboliteList(i),model.metNames)); %Should use id maybe, setting
0045         if ~isempty(metaboliteIndex)
0046             [~, b]=find(model.S(metaboliteIndex,:));
0047             reactionIndexes=[reactionIndexes; b(:)];
0048         else
0049             fprintf('Could not find any reactions with the metabolite %s\n\n',char(metaboliteList(i)))
0050         end
0051     end
0052     reactionIndexes=unique(reactionIndexes);
0053 else
0054     reactionIndexes=(1:length(fluxesA))';
0055 end
0056 
0057 %Finds the reactions where either flux is at or above the cutOffFlux value
0058 in1=find(abs(fluxesA(reactionIndexes))>=cutOffFlux)';
0059 in2=find(abs(fluxesB(reactionIndexes))>=cutOffFlux)';
0060 ineither=reactionIndexes(unique([in1 in2]));
0061 
0062 %Keep only those solutions where the difference is larger than or equal to
0063 %cutOffDiff
0064 ineither=ineither(find(abs(fluxesA(ineither)-fluxesB(ineither))>=cutOffDiff));
0065 
0066 %Finds the reactions where the fluxes differ more than cutOffChange percent
0067 %First check those fluxes that are non-zero in solution1.x
0068 nonZeroFluxes=ineither(find(fluxesA(ineither)));
0069 quota=1+cutOffChange/100;
0070 larger=nonZeroFluxes(find((fluxesB(nonZeroFluxes)./fluxesA(nonZeroFluxes))>=(quota)))';
0071 smaller=nonZeroFluxes(find((fluxesB(nonZeroFluxes)./fluxesA(nonZeroFluxes))<(1/quota)))';
0072 fluxIndexes=[larger smaller];
0073 
0074 %Then add those where solution1 has a zero flux
0075 zeroFluxes=ineither(find(fluxesA(ineither)==0));
0076 fluxIndexes=unique([fluxIndexes zeroFluxes(find(abs(fluxesB(zeroFluxes))>=cutOffFlux))']);
0077 
0078 formulas=constructEquations(model,model.rxns(fluxIndexes));
0079 
0080 if nargin>4
0081     if nargin>5
0082         fprintf('These reactions have flux values that differ by more than %s percent, absolute values above %s, and a total difference above %s (%s reactions)\n\n',num2str(cutOffChange),num2str(cutOffFlux),num2str(cutOffDiff),num2str(length(formulas)));
0083     else
0084         fprintf('These reactions have flux values that differ by more than %s percent and absolute values above %s (%s reactions)\n\n',num2str(cutOffChange),num2str(cutOffFlux),num2str(length(formulas)));
0085     end
0086 else
0087     fprintf('These reactions have flux values that differ by more than %s percent (%s reactions)\n\n',num2str(cutOffChange),num2str(length(formulas)));
0088 end
0089 
0090 metaboliteNames=[];
0091 for i=1:length(metaboliteList)
0092     metaboliteNames=[metaboliteNames char(metaboliteList(i)) ' '];
0093 end
0094 
0095 if ~isempty(metaboliteNames)
0096     fprintf('Only prints reactions involving one or more of the following metabolites:\n%s\n\n',metaboliteNames)
0097 end
0098 
0099 for i=1:length(formulas)
0100     fluxText=['Flux: ' num2str(fluxesA(fluxIndexes(i))) ' Reference flux: ' num2str(fluxesB(fluxIndexes(i))) ' Difference: ' num2str(fluxesA(fluxIndexes(i))-fluxesB(fluxIndexes(i)))];
0101     fprintf('%s: %s\n\t%s\n\t%s\n\n', char(model.rxns(fluxIndexes(i))), char(formulas(i)),...
0102         char(model.rxnNames(fluxIndexes(i))),fluxText);
0103 end
0104 end