drawMap.m

function notMapped=drawMap(title,pathway,modelA,conditionA,conditionB,modelB,filename,cutOff,supressOpen)
% drawMap
%   Imports a previously drawn map of the metabolic network and plots
%   the fluxes on that map. If the pathway contains expression data the log-fold
%   changes are plotted as well
%
%   title           title to be printed on the map
%   pathway         map structure as generated by
%                   constructPathwayFromCelldesigner
%   modelA        	model structure for condition A
%   conditionA      flux vector for condition A
%   conditionB    	flux vector for condition B (optional)
%   modelB        	model structure for condition B (optional, default modelA)
%   filename        exports the map as a pdf (optional)
%   cutOff          only print fluxes where one of the fluxes are above this
%                   value (optional, default 10^-7)
%   supressOpen     true if the pdf file should not be opened (optional, default
%                   false)
%
%   notMapped     	the reaction ids that carried a flux in either of the
%                   conditions but that were not present in the metabolic map
%
%   The objective function in modelA is assumed to be the objective function
%   in both cases. If only one flux is supplied then that flux is printed two
%   times in each box.
%
% Usage: notMapped=drawMap(title,pathway,modelA,conditionA,...
%           conditionB,modelB,filename,cutOff,supressOpen)

if nargin<4
    conditionB=conditionA;
end

if nargin<5
    modelB=modelA;
end

if nargin<8
    cutOff=10^-7;
end

if nargin<9
    supressOpen=false;
end

if isempty(conditionB)
    conditionB=conditionA;
end

if isempty(modelB)
    modelB=modelA;
end

objectiveNames=modelA.rxns(modelA.c~=0);
objectiveValues=modelA.c(modelA.c~=0);

%Replace all values which are below the cut-off with 0.0
conditionA(abs(conditionA)<cutOff)=0;
conditionB(abs(conditionB)<cutOff)=0;

%Get the reactions that could not be mapped
carriesFlux=union(modelA.rxns,modelB.rxns);

rxnsInModel={};
%Loop through the pathway structure and get the reaction that can be mapped
for i=1:numel(pathway.listOfSpecies)
    if strcmp(pathway.listOfSpecies(i).type,'PROTEIN')
        if ~isempty(pathway.listOfSpecies(i).note)
            rxnsInModel=[rxnsInModel;pathway.listOfSpecies(i).note(1)];
        end
    end
end

%Find all reaction bounds that are "non-standard". Since the default bounds
%are not indicated in the model structure it is assumed that the maximal
%bound value is the default one.
rxnBoundsA=find(modelA.ub~=max(modelA.ub) | (modelA.lb~=-max(modelA.ub) & modelA.rev==1) | (modelA.lb~=0 & modelA.rev==0));
rxnBoundsB=find(modelB.ub~=max(modelB.ub) | (modelB.lb~=-max(modelB.ub) & modelB.rev==1) | (modelB.lb~=0 & modelB.rev==0));

notMapped=setdiff(carriesFlux,rxnsInModel);

pathway=colorPathway(pathway,modelA.rxns,conditionA,conditionB,cutOff);
pathway=markPathwayWithFluxes(pathway,modelA.rxns,conditionA,conditionB);

h=figure('Visible','Off');
drawPathway(pathway,h,cutOff);

%Plot labels and title at pathway map
plotLabels(h,pathway);
setTitle(h,pathway,title);

%Plot example box, color example and additional information
exampleBoxText(1)={'Reaction ID'};
exampleBoxText(2)={'Condition A'};
exampleBoxText(3)={'Condition B'};
if nargin>7
    exampleBoxText(4)={'Other data'};
end

additionalText={};
additionalText=[additionalText;{'Objective function'}];

additionalText=[additionalText;{getObjectiveString(true, objectiveNames, objectiveValues)}];
additionalText=[additionalText;{''}];
additionalText=[additionalText;{'Constraints (condition A)'}];

counter=numel(additionalText)+1;
for i=1:length(rxnBoundsA)
    additionalText{counter}=[modelA.rxns{rxnBoundsA(i)} '  [' ...
        num2str([modelA.lb(rxnBoundsA(i)) modelA.ub(rxnBoundsA(i))]) ']'];
    counter=counter+1;
end

additionalText=[additionalText;{''}];
additionalText=[additionalText;{'Constraints (condition B)'}];

counter=numel(additionalText)+1;
for i=1:length(rxnBoundsB)
    additionalText{counter}=[modelB.rxns{rxnBoundsB(i)} '  [' ...
        num2str([modelB.lb(rxnBoundsB(i)) modelB.ub(rxnBoundsB(i))]) ']'];
    counter=counter+1;
end

plotAdditionalInfo(h, pathway, additionalText, exampleBoxText, 1);

if nargin>5
    if any(filename)
        set(gcf, 'PaperType', 'A4');
        orient landscape;
        print('-dpdf', '-r600', filename);
        fprintf('File saved as %s\n',filename);
    end
end

if supressOpen==false
    open(filename);
end
end