drawPathway

 drawPathway
    Draws a metabolic network to a figure.

   pathway         pathway structure representing the pathway to be drawn
   handle          handle to a figure (opt)
   cutOff          the fluxes are only printed if the absolute value of
                   at least one of the fluxes is above the cutoff value
                   (opt, default 0)
   defaultColor    color in Matlab format to be used as the background
                   color for enzymes if no color is specified in the
                   pathway structure (opt, default [1 1 1])

   NOTE:   At the moment all text sizes and some positions are hard coded.
           This means that this code is not appliable for any map

   Usage: drawPathway(pathway, handle, cutOff, defaultColor)

0001 function drawPathway(pathway, h, cutOff, defaultColor)
0002 % drawPathway
0003 %    Draws a metabolic network to a figure.
0004 %
0005 %   pathway         pathway structure representing the pathway to be drawn
0006 %   handle          handle to a figure (opt)
0007 %   cutOff          the fluxes are only printed if the absolute value of
0008 %                   at least one of the fluxes is above the cutoff value
0009 %                   (opt, default 0)
0010 %   defaultColor    color in Matlab format to be used as the background
0011 %                   color for enzymes if no color is specified in the
0012 %                   pathway structure (opt, default [1 1 1])
0013 %
0014 %   NOTE:   At the moment all text sizes and some positions are hard coded.
0015 %           This means that this code is not appliable for any map
0016 %
0017 %   Usage: drawPathway(pathway, handle, cutOff, defaultColor)
0018 
0019 if nargin<4
0020     defaultColor=[1 1 1];
0021 end
0022 if nargin<3
0023     cutOff=0;
0024 end
0025 if nargin<2
0026     if ishandle(gcf)
0027         h=gcf;
0028     else
0029         h=figure();
0030     end
0031 end
0032 
0033 %Will contain all metabolite info and where to print them
0034 metText={};
0035 metX=[];
0036 metY=[];
0037 metRec=[];
0038 
0039 %Same for reaction text (names, fluxes and maybe other)
0040 rxnText={};
0041 rxnX=[];
0042 rxnY=[];
0043 rxnLinesX=[];
0044 rxnLinesY=[];
0045 rxnRec=[];
0046 rxnFace=[];
0047 rxnEdge=[];
0048 rxnArrows=[];
0049 rxnArrowsColor=[];
0050 
0051 %Expression boxes
0052 expRec=[];
0053 expFace=[];
0054 
0055 %Lines to enzymes
0056 enzLinesX=[];
0057 enzLinesY=[];
0058 
0059 %Same for compartments
0060 compText={};
0061 compX=[];
0062 compY=[];
0063 compRec=[];
0064 
0065 %Check if it should plot expression data
0066 if isfield(pathway.listOfSpecies(1),'orfs')
0067     drawExpression=true;
0068 else
0069     drawExpression=false;
0070 end
0071 
0072 %Get the dimensions of the network
0073 dimension=getPathwayDimensions(pathway);
0074 
0075 %Specify the properties of the figure
0076 set(h, 'PaperOrientation', 'landscape');
0077 set(h, 'Color', [1 1 1]);
0078 xmargin=0.83;
0079 ymargin=1.32;
0080 set(h, 'PaperPosition', [xmargin ymargin 29.67743169791-2*xmargin 20.98404194812-2*ymargin]);
0081 set(h, 'Renderer', 'painters');
0082 set(h, 'Position', [10 10 1000 700]);
0083 
0084 %The height has to be at least the height of the information panel
0085 %(approximately 3000)
0086 height=dimension(2)+dimension(4)+200;
0087 height=max(height,3000);
0088 
0089 axes('ytick',[],'ydir','reverse','xtick',[],'ycolor',[1 1 1], 'xcolor',...
0090     [1 1 1],'position',[0 0 1 1], 'xLim',[dimension(1)-40 dimension(1)+dimension(3)+800],...
0091     'yLim', [dimension(2)-200 height]);
0092 daspect([1,1,1]);
0093 
0094 %Loops through the compartments.
0095 for i=1:length(pathway.listOfCompartments)
0096     position=[pathway.listOfCompartments(1,i).x pathway.listOfCompartments(1,i).y...
0097         pathway.listOfCompartments(1,i).w pathway.listOfCompartments(1,i).h];
0098     compRec=[compRec;position];
0099     compText=[compText;upper(pathway.listOfCompartments(i).name)];
0100     compX=[compX;position(1)+0.5*position(3)];
0101     compY=[compY;position(2)+position(4)-70];
0102 end
0103 
0104 %Get positions of lines
0105 for i=1:length(pathway.listOfReactions)
0106     for j=1:length(pathway.listOfReactions(i).componentList)
0107         %The line is always drawn from the middle to the component
0108         middle=pathway.listOfReactions(i).middlePoint;
0109         x(1)=middle(1);
0110         y(1)=middle(2);
0111         x(2)=pathway.listOfReactions(i).componentList(j).anchor(1);
0112         y(2)=pathway.listOfReactions(i).componentList(j).anchor(2);
0113         
0114         %Check to see if a line or an arrow should be drawn
0115         if strcmpi(pathway.listOfReactions(i).componentList(j).toArrow,'true')
0116             rxnArrows=[rxnArrows;x(1) y(1) x(2) y(2)];
0117             
0118             %Draw a red arrow if it's a base product
0119             if strcmpi(pathway.listOfReactions(i).componentList(j).baseProduct,'true')
0120                 rxnArrowsColor=[rxnArrowsColor;1 0 0];
0121             else
0122                 rxnArrowsColor=[rxnArrowsColor;0 0 0];
0123             end
0124         else
0125             %If the component is an enzyme then a different line should be
0126             %used
0127             if strcmpi(pathway.listOfReactions(i).componentList(j).type,'ENZYME')
0128                 enzLinesX=[enzLinesX;x(1) x(2)];
0129                 enzLinesY=[enzLinesY;y(1) y(2)];
0130             else
0131                 rxnLinesX=[rxnLinesX;x(1) x(2)];
0132                 rxnLinesY=[rxnLinesY;y(1) y(2)];
0133             end
0134         end
0135     end
0136 end
0137 
0138 %Loops through the species
0139 for i=1:length(pathway.listOfSpecies)
0140     %Species should be represented with ellipses and enzymes with
0141     %rectangles
0142     position=[pathway.listOfSpecies(i).x pathway.listOfSpecies(i).y...
0143         pathway.listOfSpecies(i).w pathway.listOfSpecies(i).h];
0144     if strcmpi(pathway.listOfSpecies(i).type, 'SIMPLE_MOLECULE')
0145         metRec=[metRec;position];
0146         metText=[metText;pathway.listOfSpecies(i).name];
0147         metX=[metX;position(1)+0.3*position(3)];
0148         metY=[metY;position(2)+0.5*position(4)];
0149     end
0150     if strcmpi(pathway.listOfSpecies(i).type, 'PROTEIN')
0151         %The color of the enzyme can be specified in the pathway object.
0152         %The default color is used if there is no 'color' field present.
0153         %If neither of the fluxes are above the cutoff value, then use the
0154         %default color
0155         faceColor=defaultColor;
0156         if isfield(pathway.listOfSpecies(i),'color')
0157             if any(pathway.listOfSpecies(i).color)
0158                 %Check that the fluxes are specified
0159                 if isfield(pathway.listOfSpecies(i),'flux') && isfield(pathway.listOfSpecies(i),'referenceFlux')
0160                     %Check that the value of one of the fluxes is above the cutoff value
0161                     if abs(pathway.listOfSpecies(i).referenceFlux)>=cutOff || abs(pathway.listOfSpecies(i).flux)>=cutOff
0162                         faceColor=pathway.listOfSpecies(i).color;
0163                     end
0164                 end
0165             end
0166         end
0167         
0168         %If the reaction is associated with a sign change and either of the
0169         %fluxes are larger than the cutoff value, then use a different frame
0170         edgeColor=[0 0 0];
0171         if isfield(pathway.listOfSpecies(i),'signChange')
0172             if any(pathway.listOfSpecies(i).signChange)
0173                 if (pathway.listOfSpecies(i).signChange==true) && (abs(pathway.listOfSpecies(i).referenceFlux)>=cutOff || abs(pathway.listOfSpecies(i).flux)>=cutOff)
0174                     edgeColor=[1 0.6 0.2];
0175                 end
0176             end
0177         end
0178         
0179         rxnFace=[rxnFace;faceColor];
0180         rxnEdge=[rxnEdge;edgeColor];
0181         
0182         %Draw smaller boxes if expression data should also be printed
0183         if drawExpression==true
0184             rxnRec=[rxnRec;position-[0 0 40 0]];
0185             
0186             %Draw rectangles representing each gene. Max 8 circles are used.
0187             %Their width is hardcoded and so is the width of the area where
0188             %they are printed (40)
0189             if any(pathway.listOfSpecies(i).expA)
0190                 %Calculate the position of the dots. 4 in each column.
0191                 firstFourX=pathway.listOfSpecies(i).x+pathway.listOfSpecies(i).w-34;
0192                 fourY=pathway.listOfSpecies(i).y:pathway.listOfSpecies(i).h/4:pathway.listOfSpecies(i).y+pathway.listOfSpecies(i).h;
0193                 
0194                 nExp=numel(pathway.listOfSpecies(i).expA);
0195                 
0196                 %One log10-fold change is total up or down
0197                 colorCodes=getColorCodes(pathway.listOfSpecies(i).expA, pathway.listOfSpecies(i).expB);
0198                 
0199                 %Plot the circles
0200                 for j=1:min(nExp,4)
0201                     expRec=[expRec;firstFourX fourY(j) 16 pathway.listOfSpecies(i).h/4];
0202                     expFace=[expFace;colorCodes{j}];
0203                 end
0204                 
0205                 %If a second row is needed
0206                 if nExp>4
0207                     secondFourX=pathway.listOfSpecies(i).x+pathway.listOfSpecies(i).w-16;
0208                     for j=5:min(nExp,8)
0209                         expRec=[expRec;secondFourX fourY(j-4) 16 pathway.listOfSpecies(i).h/4];
0210                         expFace=[expFace;colorCodes{j}];
0211                     end
0212                 end
0213             end
0214         else
0215             rxnRec=[rxnRec;position];
0216         end
0217         
0218         %If no fluxes are specified then should only the name be printed.
0219         %Flux values should only be printed if at least one of the fluxes
0220         %is above the cutoff value
0221         textToPrint={};
0222         
0223         %NOTE:  This is done since I know that I have names including '_' which is a reserved
0224         %       character in LaTeX. This should be done more systematically
0225         
0226         %NOTE:  The notation used means that only flux values with with less than
0227         %       nine characters can be used
0228         
0229         textToPrint{1,1}=strrep(pathway.listOfSpecies(i).name,'_','\_');
0230         
0231         %Check that the fluxes are specified
0232         %NOTE: Should you only print if there are two fluxes?
0233         if isfield(pathway.listOfSpecies(i),'flux') && isfield(pathway.listOfSpecies(i),'referenceFlux')
0234             %Check that the value of one of the fluxes is above the cutoff value
0235             if ~isempty(pathway.listOfSpecies(i).flux) && ~isempty(pathway.listOfSpecies(i).referenceFlux)
0236                 if abs(pathway.listOfSpecies(i).referenceFlux)>=cutOff || abs(pathway.listOfSpecies(i).flux)>=cutOff
0237                     textToPrint{2,1}=num2str(pathway.listOfSpecies(i).flux,'%0.4g');
0238                     textToPrint{3,1}=num2str(pathway.listOfSpecies(i).referenceFlux,'%0.4g');
0239                 end
0240             end
0241         end
0242         
0243         xPos=ones(numel(textToPrint),1)*(position(1)+3);
0244         distance=position(4)/(numel(textToPrint)+1);
0245         yPos=position(2)+distance:distance:(position(2)+position(4)-distance);
0246         
0247         rxnText=[rxnText;textToPrint];
0248         rxnX=[rxnX;xPos];
0249         rxnY=[rxnY;yPos'];
0250     end
0251 end
0252 
0253 %*** Print everything ***
0254 
0255 %Compartment rectangles
0256 for i=1:size(compRec,1)
0257     rectangle('edgecolor',[1 0.8 0.2], 'facecolor', [1 0.91 0.55], 'linewidth', 2, 'curvature', [0.05 0.05],...
0258         'position',compRec(i,:));
0259     rectangle('edgecolor',[1 0.8 0.2], 'facecolor', [1 1 1], 'linewidth', 1, 'curvature', [0.05 0.05],...
0260         'position',compRec(i,:)+[12 12 -24 -24]);
0261 end
0262 
0263 %Reaction lines
0264 line(rxnLinesX',rxnLinesY','linewidth',0.1,'color',[0 0 0]);
0265 
0266 %Arrow lines
0267 hold on;
0268 for i=1:size(rxnArrows,1)
0269     plotArrow(rxnArrows(i,1),rxnArrows(i,2),rxnArrows(i,3),rxnArrows(i,4),'linewidth',0.1,'edgecolor',rxnArrowsColor(i,:),'facecolor',rxnArrowsColor(i,:));
0270 end
0271 hold off;
0272 
0273 %Enzyme lines
0274 line(enzLinesX',enzLinesY','marker','o','markeredgecolor',[0 0 0],'markersize',0.5,...
0275     'markerfacecolor',[1 0.969 0.922],'linewidth',0.1,'color',[0 0 0]);
0276 
0277 %Metabolite rectangles
0278 for i=1:size(metRec,1)
0279     rectangle('edgecolor', [0 0 0], 'facecolor', [1 1 0.78], 'curvature', [0.8 0.8],...
0280         'position',metRec(i,:),'linewidth',0.1);
0281 end
0282 
0283 %Reaction rectangles
0284 for i=1:size(rxnRec,1)
0285     rectangle('edgecolor', rxnEdge(i,:), 'facecolor', rxnFace(i,:), 'linewidth',0.1,...
0286         'curvature', [0.1 0.1], 'position',rxnRec(i,:));
0287 end
0288 
0289 %Expression rectangles
0290 for i=1:size(expRec,1)
0291     rectangle('facecolor', expFace(i,:), 'linewidth',0.1, 'curvature', [0.1 0.1],...
0292         'position',expRec(i,:));
0293 end
0294 
0295 %Metabolite text
0296 %Temporary thing to get it right for small networks as well. Should be
0297 %calculated better!
0298 if dimension(3)<5000
0299     textSize=3;
0300 else
0301     textSize=0.5;
0302 end
0303 
0304 text(metX,metY,metText,'fontname','Small Fonts','fontsize',textSize,'interpreter','tex',...
0305     'HorizontalAlignment','center','verticalalignment','middle');
0306 
0307 %Reaction text
0308 text(rxnX,rxnY,rxnText,'fontname','Small Fonts','fontsize',textSize,'interpreter','tex',...
0309     'verticalalignment','middle','HorizontalAlignment','left');
0310 
0311 %Compartment text. The text is centered and at the bottom of the compartment
0312 text(compX, compY,compText,'fontname','Small Fonts','fontsize',7,'interpreter','none',...
0313     'HorizontalAlignment','center','verticalalignment','middle');
0314 end
0315 
0316 function handles = plotArrow( x1,y1,x2,y2,varargin )
0317 %
0318 % plotArrow - plots an arrow to the current plot
0319 %
0320 % format:   handles = plotArrow( x1,y1,x2,y2 [,options...] )
0321 %
0322 % input:    x1,y1   - starting point
0323 %           x2,y2   - end point
0324 %           options - come as pairs of "property","value" as defined for "line" and "patch"
0325 %                     controls, see matlab help for listing of these properties.
0326 %                     note that not all properties where added, one might add them at the end of this file.
0327 %
0328 %                     additional options are:
0329 %                     'headwidth':  relative to complete arrow size, default value is 0.07
0330 %                     'headheight': relative to complete arrow size, default value is 0.15
0331 %                     (encoded are maximal values if pixels, for the case that the arrow is very long)
0332 %
0333 % output:   handles - handles of the graphical elements building the arrow
0334 %
0335 % Example:  plotArrow( -1,-1,15,12,'linewidth',2,'color',[0.5 0.5 0.5],'facecolor',[0.5 0.5 0.5] );
0336 %           plotArrow( 0,0,5,4,'linewidth',2,'headwidth',0.25,'headheight',0.33 );
0337 %           plotArrow;   % will launch demo
0338 
0339 % =============================================
0340 % constants (can be edited)
0341 % =============================================
0342 %NOTE: These should be 100 times smaller
0343 alpha       = 15;   % head length
0344 beta        = 7;   % head width
0345 max_length  = 11;
0346 %max_width   = 7;
0347 
0348 % =============================================
0349 % check if head properties are given
0350 % =============================================
0351 % if ratio is always fixed, this section can be removed!
0352 if ~isempty( varargin )
0353     for c = 1:floor(length(varargin)/2)
0354         switch lower(varargin{c*2-1})
0355             % head properties - do nothing, since handled above already
0356             case 'headheight',alpha = max( min( varargin{c*2},1 ),0.01 );
0357             case 'headwidth', beta = max( min( varargin{c*2},1 ),0.01 );
0358         end
0359     end
0360 end
0361 
0362 % =============================================
0363 % calculate the arrow head coordinates
0364 % =============================================
0365 den         = x2 - x1 + eps;                                % make sure no devision by zero occurs
0366 teta        = atan( (y2-y1)/den ) + pi*(x2<x1) - pi/2;      % angle of arrow
0367 cs          = cos(teta);                                    % rotation matrix
0368 ss          = sin(teta);
0369 R           = [cs -ss;ss cs];
0370 line_length = sqrt( (y2-y1)^2 + (x2-x1)^2 );                % sizes
0371 head_length = min( line_length*alpha,max_length );
0372 head_width  = min( line_length*beta,max_length );
0373 x0          = x2*cs + y2*ss;                                % build head coordinats
0374 y0          = -x2*ss + y2*cs;
0375 coords      = R*[x0 x0+head_width/2 x0-head_width/2; y0 y0-head_length y0-head_length];
0376 
0377 % =============================================
0378 % plot arrow  (= line + patch of a triangle)
0379 % =============================================
0380 h1          = plot( [x1,x2],[y1,y2],'k' );
0381 h2          = patch( coords(1,:),coords(2,:),[0 0 0] );
0382 
0383 % =============================================
0384 % return handles
0385 % =============================================
0386 handles = [h1 h2];
0387 
0388 % =============================================
0389 % check if styling is required
0390 % =============================================
0391 % if no styling, this section can be removed!
0392 if ~isempty( varargin )
0393     for c = 1:floor(length(varargin)/2)
0394         switch lower(varargin{c*2-1})
0395             
0396             % only patch properties
0397             case 'edgecolor',   set( h2,'EdgeColor',varargin{c*2} );
0398             case 'facecolor',   set( h2,'FaceColor',varargin{c*2} );
0399             case 'facelighting',set( h2,'FaceLighting',varargin{c*2} );
0400             case 'edgelighting',set( h2,'EdgeLighting',varargin{c*2} );
0401                 
0402                 % only line properties
0403             case 'color'    , set( h1,'Color',varargin{c*2} );
0404                 
0405                 % shared properties
0406             case 'linestyle', set( handles,'LineStyle',varargin{c*2} );
0407             case 'linewidth', set( handles,'LineWidth',varargin{c*2} );
0408             case 'parent',    set( handles,'parent',varargin{c*2} );
0409                 
0410                 % head properties - do nothing, since handled above already
0411             case 'headwidth',
0412             case 'headheight',
0413                 
0414         end
0415     end
0416 end
0417 end