\documentclass[a4paper]{article}
%\VignetteIndexEntry{Drawing Phylogenies}
%\VignettePackage{ape}
\usepackage{ape}

\author{Emmanuel Paradis}
\title{Drawing Phylogenies in \R: Basic and Advanced Features With \pkg{ape}}

\begin{document}
\DefineVerbatimEnvironment{Sinput}{Verbatim}{formatcom=\color{darkblue}}
\DefineVerbatimEnvironment{Soutput}{Verbatim}{formatcom=\color{black}\vspace{-1.5em}}

\maketitle
\tableofcontents\vspace*{1pc}\hrule

<<echo=false,quiet=true>>=
options(width = 80, prompt = "> ")
@

\vspace{1cm}

\section{Introduction}

Graphical functions have been present in \ape\ since its first version
(0.1, released in August 2002). Over the years, these tools have been
improved to become quite sophisticated although complicated to use
efficiently. This document gives an overview of these functionalities.
Section~\ref{sec:basic} explains the basic concepts
and tools behind graphics in \ape. A figure made with \ape\ usually
starts by calling the function \code{plot.phylo} which is detailed in
Section~\ref{sec:plotphylo}, and further graphical annotations can be done with
functions covered in Section~\ref{sec:annot}. Section~\ref{sec:spec}
shows some specialized functions available in \ape, and finally,
Sections~\ref{sec:geom} and \ref{sec:build} give an overview of some
ideas to help making complicated figures.

\section{Basic Concepts}\label{sec:basic}

The core of \ape's graphical tools is the \code{plot} method for the
class \code{"phylo"}, the function \code{plot.phylo}. This function is
studied in details in Section~\ref{sec:plotphylo}, but first we see
the basic ideas behind it and other functions mentioned in this
document.

\subsection{Graphical Model}

The graphical functions in \ape\ use the package \pkg{graphics}.
Overall, the conventions of this package are followed quite
closely (see Murrell's book \cite{Murrell2006}), so users familiar
with graphics in \R\ are expected to find their way relatively easily when
plotting phylogenies with \ape.

\ape\ has several functions to perform computations before drawing a
tree, so that they may be used to implement the same graphical
functionalities with other graphical engines such as the \pkg{grid}
package. These functions are detailed in the next section.

To start simply, we build a small tree with three genera of primates
which we will use in several examples in this document:

<<>>=
library(ape)
mytr <- read.tree(text = "((Pan:5,Homo:5):2,Gorilla:7);")
@

\noindent Now let's build a small function to show the frame around
the plot with dots, and the $x$- and $y$-axes in green:

<<>>=
foo <- function() {
    col <- "green"
    for (i in 1:2)
        axis(i, col = col, col.ticks = col, col.axis = col, las = 1)
    box(lty = "19")
}
@

\noindent We then plot the tree in four different ways (see below for
explanations about the options) and call for each of them the previous
small function:

<<fig=true>>=
layout(matrix(1:4, 2, 2, byrow = TRUE))
plot(mytr); foo()
plot(mytr, "c", FALSE); foo()
plot(mytr, "u"); foo()
par(xpd = TRUE)
plot(mytr, "f"); foo()
box("outer")
@

\noindent The last command (\code{box("outer")}) makes visible the
most outer frame of the figure showing more clearly the margins around
each tree (more on this in Sect.~\ref{sec:geom}). We note also the command \code{par(xpd
  = TRUE)}: by default this parameter is \code{FALSE} so that graphical
elements (points, lines, text, \dots) outside the plotting region (i.e., in
the margins or beyond) are cut (clipped).\footnote{\code{par(xpd = TRUE)} is
  used in several examples in this document mainly because of the
  small size of the trees drawn here. However, in practice, this is rarely
  needed.} These small figures illustrate the way trees are drawn with
\ape. This can be summarised with the following (pseudo-)algorithm:

\bigskip\hrule height 1pt\relax
%\renewcommand{\theenumi}{\alph{enumi}}
\renewcommand{\labelenumi}{\textbf{\theenumi.}}
\begin{enumerate}\small
\item Compute the node coordinates depending on the type of tree plot,
  the branch lengths, and other parameters.
\item Evaluate the space required for printing the tip labels.
\item Depending on the options, do some rotations and/or translations.
\item Set the limits of the $x$- and $y$-axes.
\item Open a graphical device (or reset it if already open) and draw
  an empty plot with the limits found at the previous step.
\item Call \code{segments()} to draw the branches.
\item Call \code{text()} to draw the labels.
\end{enumerate}
\hrule height 1pt\relax\bigskip

There are a lot of ways to control these steps. The main variations
along these steps are given below.

\textbf{Step 1. } The option \code{type} specifies the shape
of the tree plot: five values are possible, \code{"phylogram"},
\code{"cladogram"}, \code{"fan"}, \code{"unrooted"}, and
\code{"radial"} (the last one is not considered in this document). The
first three types are valid representations for rooted trees, while
the fourth one should be selected for unrooted trees.

The node coordinates depend also on whether the tree has branch
lengths or not, and on the options \code{node.pos} and
\code{node.depth}. This is illustrated below using a tree with eight
tips and all branch length equal to one (these options have little
effect if the tree has only three tips):

<<>>=
tr <- compute.brlen(stree(8, "l"), 0.1)
tr$tip.label[] <- ""
@

\noindent We now draw this tree using the option \code{type =
  "phylogram"} (first column of plots) or \code{type = "cladogram"}
(second column) and different options:

<<echo=false,quiet=true>>=
foo <- function() {
    col <- "green"
    axis(1, col = col, col.ticks = col, col.axis = col)
    axis(2, col = col, col.ticks = col, col.axis = col, at = 1:Ntip(tr), las = 1)
    box(lty = "19")
}
@
<<>>=
@

<<fig=true>>=
layout(matrix(1:12, 6, 2))
par(mar = c(2, 2, 0.3, 0))
for (type in c("p", "c")) {
   plot(tr, type); foo()
   plot(tr, type, node.pos = 2); foo()
   plot(tr, type, FALSE); foo()
   plot(tr, type, FALSE, node.pos = 1, node.depth = 2); foo()
   plot(tr, type, FALSE, node.pos = 2); foo()
   plot(tr, type, FALSE, node.pos = 2, node.depth = 2); foo()
}
@

\noindent Some combinations of options may result in the same
tree shape as shown by the last two rows of trees. For unrooted and
circular trees, only the option \code{use.edge.length} has an effect
on the layout and/or the scales of the axes:

<<echo=false,quiet=true>>=
foo <- function() {
    col <- "green"
    for (i in 1:2) axis(i, col = col, col.ticks = col, col.axis = col, las = 1)
    box(lty = "19")
}
@

<<fig=true>>=
layout(matrix(1:4, 2, 2))
par(las = 1)
plot(tr, "u"); foo()
plot(tr, "u", FALSE); foo()
plot(tr, "f"); foo()
plot(tr, "f", FALSE); foo()
@

\textbf{Step 2.} In the \pkg{graphics} package, text are
printed with a fixed size, which means that whether you draw a small
tree or a large tree, on a small or large device, the labels will
have the same size. However, before anything is plotted or
drawn on the device it is difficult to find the correspondence between
this size (in inches) and the user coordinates used for the node
coordinates. Therefore, the following steps are implemented to
determine the limits on the $x$-axis:

\renewcommand{\labelenumi}{\theenumi.}
\begin{enumerate}
\item Find the width of the device in inches (see
  Sect.~\ref{sec:overlay}).
\item Find the widths of all labels in inches: if at least one of them is
wider than the device, assign two thirds of the device for the
branches and one third to the tip labels. (This makes sure that by
default the tree is visible in the case there are very long tip
labels.)
\item Otherwise, the space allocated to the tip labels is increased
incrementally until all labels are visible on the device.
\end{enumerate}

The limits on the $y$-axis are easier to determine since it depends
only on the number of branches in the tree. The limits on both axes
can be changed manually with the options \code{x.lim} and \code{y.lim}
which take one or two values: if only one value is given this will
set the rightmost or uppermost limit, respectively; if two values are given these
will set both limits on the respecive axis.\footnote{These two options differ
from their standard counterparts \code{xlim} and
\code{ylim} which always require two values.}

By default, there is no space between the tip labels and the tips of
the terminal branches; however, text strings are printed with a
bounding box around them making sure there is actually a small
space (besides, the default font is italics making this space more
visible). The option \code{label.offset} (which is 0 by default) makes
possible to add an explicit space between them (this must be in user
coordinates).

\textbf{Step 3.} For rooted trees, only 90\textdegree\
rotations are supported using the option \code{direction}.\footnote{To have
full control of the tree rotation, the option `rotate' in
\LaTeX\ does the job very well.} For
unrooted (\code{type = "u"}) and circular (\code{type = "fan"}) trees, full rotation is
supported with the option \code{rotate.tree}. If these options are
used, the tip labels are not rotated. Label rotation is controlled by
other options: \code{srt}\footnote{\code{srt} is for \textit{string
    rotation}, not to be confused with the function \code{str} to
  print the \textit{structure} of an object.} for all trees, and
\code{lab4ut} for unrooted trees.

\textbf{Step 4.} These can be fully controlled with the options
\code{x.lim} and \code{y.lim}. Note that the options \code{xlim} and
\code{ylim} \emph{cannot} be used from \code{plot.phylo}.

\textbf{Step 5.} If the options \code{plot = FALSE} is used, then
steps 6 and 7 are not performed.

\subsection{Computations}\label{sec:comput}

As we can see from the previous section, a lot of computations are done
before a tree is plotted. Some of these computations are performed
by special functions accessible to all users, particularly the three
functions used to calculate the node coordinates. First, two functions
calculate ``node depths'' which are the coordinates of the nodes on the
$x$-axis for rooted trees:

<<>>=
args(node.depth.edgelength)
args(node.depth)
@

\noindent Here, \code{phy} is an object of class \code{"phylo"}. The
first function uses edge lengths to calculate these coordinates, while
the second one calculates these coordinates proportional to the number
of tips descending from each node (if \code{method = 1}), or evenly
spaced (if \code{method = 2}).

The third function is \code{node.height} and is used to calculate ``node
heights'', the coordinates of the nodes on the $y$-axis:

<<>>=
args(node.height)
@

\noindent If \code{clado.style = TRUE}, the node heights are
calculated for a ``triangular cladogram'' (see figure
above). Otherwise, by default they are calculated to fall in the
middle of the vertical segments with the default \code{type =
  "phylogram"}.\footnote{It may be good to remind here than these segments,
  vertical since \code{direction = "rightwards"} is the default, are
  not part of the edges of the tree.}

For unrooted trees, the node coordinates are calculated with the
``equal angle'' algorithm described by Felsenstein
\cite{Felsenstein2004}. This is done by an internal function which
arguments are:

<<>>=
args(unrooted.xy)
@

\noindent There are three other internal functions used to plot the
segments of the tree after the above calculations have been performed
(steps 1--4 in the previous section):

<<>>=
args(phylogram.plot)
args(cladogram.plot)
args(circular.plot)
@

\noindent Although these four functions are not formally documented,
they are anyway exported because they are used by several packages
outside \ape.

\section{The \code{plot.phylo} Function}\label{sec:plotphylo}

The \code{plot} method for \code{"phylo"} objects follows quite
closely the \R\ standard practice. It has a relatively
large number of arguments: the first one (\code{x}) is mandatory
and is the tree to be drawn. It is thus not required to name it, so in
practice the tree \code{tr} can be plotted with the command
\code{plot(tr)}. All other arguments have default values:

<<>>=
args(plot.phylo)
@

\noindent The second and third arguments are the two commonly
used in practice, so they can be modified without explicitly naming them like
in the above examples. Besides, \code{"cladogram"} can be abbreviated
with \code{"c"}, \code{"unrooted"} with \code{"u"}, and so on.
For the other arguments, it is better to name them if they are used or
modified (e.g., \code{lab4ut = "a"}).

\subsection{Overview on the Options}

The logic of this long list of options is double: the user can modify
the aspect of the tree plot, and/or use some of these options to display
some data in association with the tree. Therefore, the table below
group these options into three categories. The following two sections
show how data can be displayed in connection to the tips or to the branches of the tree.

\begin{center}
\begin{tabular}{lll}
\toprule
Aspect of the tree & Attributes of the labels & Attributes of the edges\\
\midrule
\code{type} & \code{show.tip.label} & \code{edge.color}\\
\code{use.edge.length} & \code{show.node.label} & \code{edge.width}\\
\code{node.pos} & \code{font} & \code{edge.lty}\\
\code{x.lim} & \code{tip.color}\\
\code{y.lim} & \code{cex}\\
\code{direction} & \code{adj}\\
\code{no.margin} & \code{underscore}\\
\code{root.edge} & \code{srt}\\
\code{rotate.tree} & \code{lab4ut}\\
\code{open.angle} & \code{label.offset}\\
\code{node.depth} & \code{align.tip.label}\\
\bottomrule
\end{tabular}
\end{center}

\subsection{Connecting Data to the Tips}

It is common that some data are associated with the tips of a tree:
body mass, population, treatment, \dots\ The options \code{font},
\code{tip.color}, and \code{cex} make possible to show this kind of
information by changing the font (normal, bold, italics, or
bold-italics), the colour, or the size of the tip labels, or any
combination of these. These three arguments work in the usual \R\ way:
they can a vector of any length whose values are eventually recycled if
this length is less than the number of tips. This makes possible to
change all tips if a single value is given.

For instance, consider the small primate tree where we want to show
the geographic distributions stored in a factor:

<<>>=
geo <- factor(c("Africa", "World", "Africa"))
@

\noindent We can define a color for each region and use the above
factor as a numeric index vector and pass it to \code{tip.color}:

\begin{center}
\setkeys{Gin}{width=.5\textwidth}
<<fig=true,width=4,height=3.5>>=
(mycol <- c("blue", "red")[geo])
plot(mytr, tip.color = mycol)
@
\end{center}
The values must be in the same order than in the
vector of tip labels, here \code{mytr\$tip.label}. Reordering can be
done in the usual \R\ way (e.g., with \code{names} or with
\code{row.names} if the data are in a data frame). This can be
combined with another argument, for instance to show (relative) body
size:

\begin{center}
\setkeys{Gin}{width=.5\textwidth}
<<fig=true,width=4,height=3.5>>=
par(xpd = TRUE)
plot(mytr, tip.color = mycol, cex = c(1, 1, 1.5))
@
\end{center}

The function \code{def} gives another way to define the above
arguments given a vector of labels (\code{x}):

<<>>=
args(def)
@

\noindent The `\code{...}' are arguments separated by commas of the
form \code{\textsl{<label> = <value>}} with \code{\textsl{<label>}} being one of
the labels in \code{x} and \code{\textsl{<value>}} the value which will be
given to this element, whereas the value \code{default} will be given to
the others. This default value is either 1 if \code{\textsl{<value>}} is numeric,
or \code{"black"} if \code{\textsl{<value>}} is character. The above set of colours could
have thus be defined with:

<<>>=
mycol2 <- def(mytr$tip.label, Homo = "red", default = "blue")
identical(mycol, mycol2)
@

The function \code{tiplabels}, presented below, gives more
possibilities to display data on the tips of a tree.

\subsection{Connecting Data to the Branches}

The three options in the third column of the above table control the
aspect of the branches of the tree. Like the options for the tips,
the arguments are recycled. The values
given as argument(s) must be in the same order than in the \code{edge}
matrix of the tree. There are several ways to find these numbers.

The function \code{which.edge} returns the row indices of the
\code{edge} matrix given a list of tip labels so that the edges define
the clade of the tips. The returned indices can then be used to change
the aspect of these branches:

\begin{center}
\setkeys{Gin}{width=.5\textwidth}
<<fig=true,width=4,height=3.5>>=
(i <- which.edge(mytr, c("Homo", "Pan")))
co <- rep("black", Nedge(mytr))
co[i] <- "blue"
par(xpd = TRUE)
plot(mytr, edge.col = co)
edgelabels()
@
\end{center}
These indices can be displayed on the tree simply by calling
\code{edgelabels()} without option as illustrated on the previous plot
(this function is further explained in the next section).

It is possible for a function to
return data or values that are directly indexed with respect to the
rows of the \code{edge} matrix. One example is \code{chronos()}
which returns an ultrametric tree with an additional numeric vector
giving the estimated substitution rate for each branch of the
chronogram. This vector can be used directly to define the aspect of
the branches when plotting the chronogram.

\section{Annotating Trees}\label{sec:annot}

Once a tree has been plotted, it is possible to add graphical elements
with the low-level plotting functions\footnote{These functions add
  graphical elements to an existing plot, by contrast to high-level
  plotting functions which create a new plot.} in \R: \code{text},
\code{segments}, \code{points}, \code{arrows}, \code{rect}, and
\code{polygon} (all in the package \pkg{graphics}). These functions
require to give the coordinates where to draw these elements which may
be a bit difficult when plotting a tree since the axes are not drawn
by default. In practice, the function \code{locator} may be helpful here.

\ape\ has five specialized functions that also facilitate this task:
three of them add elements on the tree, and the two others do so beside the
tree. Additionally, the functions \code{axisPhylo} and \code{add.scale.bar} display scale information of the tree branches.

\subsection{Scales}

The function \code{axisPhylo} draws a time axis assuming the tips
represent the present time and the time scales backward towards the
root. For non-ultrametric trees, the most distant tip from the root is
taken to represent present time. There are a few options:

<<>>=
args(axisPhylo)
@

\noindent The option \code{backward} and \code{root.time} change the
settings of the scale, \code{side} has the same meaning than in
\code{axis()}, and the extra arguments (`\code{...}') are passed to
\code{text()}.

\code{add.scale.bar} draws a scale bar of the branch length, by
default on the lower left corner of the plot. Several options control
the position and/or the formatting:

<<>>=
args(add.scale.bar)
@

\noindent If \code{ask = TRUE}, the function calls \code{locator()}
and the user is asked to click on the plot to specify where to draw the
scale bar.

The example below illustrates both functions and draws the standard
$x$-axis above the tree:

\begin{center}
\setkeys{Gin}{width=.5\textwidth}
<<fig=true,width=4,height=3.5>>=
plot(mytr)
add.scale.bar()
axisPhylo()
axis(3)
@
\end{center}

\subsection{Tips, Nodes, and Edges}

The three functions \code{nodelabels}, \code{tiplabels}, and
\code{edgelabels} have almost identical options to add annotations on
the nodes, tips, or edges of a plotted tree:

<<>>=
args(nodelabels)
args(tiplabels)
args(edgelabels)
@

\noindent These functions differ in the second argument which specifies the
positions where the annotations should be printed or drawn. If these
functions are used without any argument, they print the numbers of the
nodes, tips, or edges, respectively, used in the \code{"phylo"}
object:\footnote{The structure of objects of class \code{"phylo"} is
  described on \ape's web site: \texttt{https://emmanuelparadis.github.io/ape\_development.html}}

\begin{center}
  \setkeys{Gin}{width=.5\textwidth}
<<fig=true,width=4,height=3.5>>=
par(xpd = TRUE)
plot(mytr)
nodelabels()
tiplabels()
edgelabels()
@
\end{center}

There are a lot of possibilities offered by these functions: a range
of examples are given in the help page which can be displayed with
\code{example(nodelabels)} and are not repeated here. A common
application of these function, and in particular \code{nodelabels}, is
to diplay pie charts of ancestral values. Several functions in \ape\
(e.g., \code{ace}), or in other packages, return their ancestral
reconstructions indexed along the nodes of the tree so that
using \code{nodelabels} is usually straightforward.

Like for \code{plot.phylo}, the options can be used for either showing
some data, or formatting the annotations. The option \code{text}
works in the same way than the standard function of the same name;
formatting of the character strings is done with the option \code{col}
and others (e.g., \code{font} or \code{cex} can be passed with
`\code{...}'). A frame is drawn around the character strings unless
\code{frame = "none"}. The option \code{pch} is used to plot
standard plotting symbols, this can take the values from 1 to 25, as
shown below:

<<fig=true>>=
plot(rep(1:5, 5), rep(1:5, each = 5), pch = 1:25, xlim = c(1, 5.2),
     col = "blue", bg = "yellow", cex = 2)
text(rep(1:5, 5) + 0.2, rep(1:5, each = 5), 1:25)
@

\noindent It should be kept in mind that these three functions are
low-level plotting commands, so they can be called as many times as
the user wishes with possibly different formatting options which could
be easier than using a single call (see the first example in
\code{?nodelabels}).

The option \code{adj} is used to set the adjustment of the
character strings passed to \code{text} (and also the plotting symbols
drawn with \code{pch}, unlike the function \code{points}). This can
used to print several series of numbers around nodes (see
\code{?nodelabels}). The adjustment is made with respect to the size of each
string, so useful values are 0 (left- or bottom-aligned), 0.5
(centred), and 1 (right- or top-aligned):

<<fig=true>>=
v <- c(0, 0.5, 1)
layout(matrix(1:9, 3, 3))
par(mar = c(3, 3, 3, 0), las = 1)
for (i in v) {
    for (j in v) {
        plot(0, 0, "n", main = paste0("adj = c(", i, ", ", j, ")"))
        abline(v = 0, h = 0, lty = 3)
        text(0, 0, "Gorilla", adj = c(i, j))
    }
}
@

\noindent Other values are accepted but they will likely result in
strings badly aligned if they have different (printed) widths:

\begin{center}
\setkeys{Gin}{width=.5\textwidth}
<<fig=true,width=4,height=3.5>>=
plot(rep(0, 3), 0:2, "n", las = 1)
abline(v = 0, h = 0:2, lty = 3)
text(0, 0:2, mytr$tip.label, adj = -1)
text(0, 0:2, mytr$tip.label, adj = 2, font = 2)
@
\end{center}

The options \code{pie} and \code{thermo}, by contrast to the previous
ones, take a matrix whose values are taken as proportions which
are displayed as piecharts or as thermometres. The colours of the
``slices'' of the pies, or the ``levels''  of the thermometers, are
specified with \code{piecol} (by default,
the function \code{rainbow} is used to define these colours). The
options \code{horiz}, \code{width}, and \code{height} are used to
format the aspect of the thermometres.

\subsection{Data Plots Besides Trees}

The two following functions add a graph beside and in connection to a
plotted tree (or around it if it is circular):

<<>>=
args(phydataplot)
args(ring)
@

\noindent Their usage is a bit more complicated than the previous
functions because it is usually needed to leave some space when
plotting the tree using
\code{x.lim} and/or \code{y.lim} (see above). It may not be
straightforward to find how much space must be left: in practice the
simplest solution is to try different values, possibly using
\code{axis} as shown above, or using the (invisibly) returned value of
\code{plot.phylo}. Once some values of these limits are found for a
given data set, they are expected to work similarly accross different
graphical devices (\code{pdf()}, \code{png()}, \dots) For more
sophisticated or automated solutions, it is possible to calculate
these limits a priori (see below).

These two functions use the labels of the data (given by their
\code{names} if the data are in a vector, their \code{rownames} if in
a matrix, or their \code{row.names} if in a data frame) to match with
the tip labels of the tree (argument \code{phy}). Thus, it is not
needed to reorder the data if they have such (row)names.

The option \code{style} can take seven different values giving a
lot of possibilities. Furthermore, further arguments can be given
which are then passed to a function of the \pkg{graphics} package giving
flexibility to customize the appearance of the plot. The table
below gives the correspondence between the values of \code{style} and
the \pkg{graphics} functions called by \code{phydataplot}:\footnote{\code{fancyarrows()} is in
  \ape: it has the same arguments than \code{arrows()} plus the option
  \code{type} which can be either \code{"triangle"} (the default) or
  \code{"harpoon"}.}

\begin{center}
  \begin{tabular}{cc}
  \toprule
  \code{style} & Extra arguments (`\code{...}') passed to:\\
  \midrule
  \code{"bars"} & \code{barplot}\\
  \code{"segments"} & \code{segments}\\
  \code{"image"} & \code{image}\\
  \code{"arrows"} & \code{fancyarrows}\\
  \code{"boxplot"} & \code{bxp}\\
  \code{"dotchart"} & \code{points}\\
  \code{"mosaic"} & \code{rect}\\
  \bottomrule
\end{tabular}
\end{center}

A simple example is given here with the primate tree to illustrate
what can be done with \code{type = "mosaic"}. We build a matrix with
ten columns filled with values 1 to 6, and give names to its rows and columns:

<<>>=
x <- matrix(1:6, 3, 10)
dimnames(x) <- list(c("Homo", "Gorilla", "Pan"), LETTERS[1:ncol(x)])
@

\noindent Note that these data are not ordered like the tip labels of
the tree:

<<>>=
x
mytr$tip.label
@

\noindent It is possible to treat these values as categories or as
continuous variables. With \code{phydataplot}, this is controlled by the
option \code{continuous} which is either a logical value (\code{FALSE}
by default, so the variables are considered as categorical), or an
integer specifying the number of classes used to represent the data (if
\code{continuous = TRUE}, then ten classes are used). Because we
call \code{phydataplot} twice, we use the option
\code{offset} which sets the space between the tree and the graph:

<<fig=true,width=11,height=6>>=
par(mar = c(10, 2, 5, 5))
plot(mytr, x.lim = 30)
phydataplot(x, mytr, "m", border = "white", offset = 2, width = 1)
phydataplot(x, mytr, "m", border = "white", offset = 15,
            width = 1, continuous = 2)

@

\noindent The options \code{border = "white"} and \code{width = 1}
help to have a better looking figure. The colour scheme used for the
first mosaic makes easy to check that the data were correctly matched
between the rows of the matrix and tips of the tree.

For circular trees, \code{offset} is the key to combine different
`rings' of data, like in the next example with a random coalescent
tree with 100 tips and three variables from a uniform distribution:

\setkeys{Gin}{width=\textwidth}
<<fig=true>>=
n <- 100
p <- 3
set.seed(3)
tr <- rcoal(n)
x <- matrix(runif(p * n), n, p)
rownames(x) <- tr$tip.label
COL <- c("red", "green", "blue")
par(xpd = TRUE)
plot(tr, "f", x.lim = c(-5, 5), show.tip.label = FALSE, no.margin = TRUE)
for (j in 1:p) ring(x[, j], tr, offset = j - 1 + 0.1, col = COL[j])
@

Like for the functions discussed in the previous section,
\code{phydataplot} and \code{ring} offer a very wide range of
possibilities which are (partially) illustrated on their help page
(\code{example(phydataplot)}) and are not repeated here.

\section{Specialized Functions}\label{sec:spec}

We see in this section several functions, built from
\code{plot.phylo}, that each have a specific purpose.

\subsection{Function \code{plotBreakLongEdges}}

It happens sometimes that, in a tree, one (or several) branch is much
longer than the others, for instance the branch between the ingroup
and a very distant outgroup. When plotting such a tree, the contrast
between the shortest branches will likely not be
visible. \code{plotBreakLongEdges()}
shortens the longest branch of the tree and represents it with a
broken segment. We illustrate this with a random tree where the first
branch is given a length of 1000 (all others have random lengths
between 0 and 1):

\setkeys{Gin}{width=.8\textwidth}
<<fig=true>>=
tree <- rtree(10)
tree$edge.length[1] <- 1000
layout(matrix(1:2, 1))
plot(tree); axisPhylo()
plotBreakLongEdges(tree); axisPhylo()
@

\noindent This function has the option \code{n = 1} which specified
the number of branches to be broken. This can be automated in some
way, for instance if we want to find how many branches in the above
tree are longer than the mean plus twice the standard-deviation:

<<>>=
el <- tree$edge.length
sum(el > mean(el) + 2 * sd(el))
@

\noindent We test with approach after changing another branch length:

<<fig=true>>=
tree$edge.length[8] <- 1000
el <- tree$edge.length
sum(el > mean(el) + 2 * sd(el))
layout(matrix(1:2, 1))
plotBreakLongEdges(tree)
plotBreakLongEdges(tree, n = 2)
@

\subsection{Function \code{drawSupportOnEdges}}

The phylogenetic bootstrap assessed the uncertainty of a
phylogenetic tree inferred from a method (parsimony, maximum likehood,
neighbour-joining, \dots) If the inference method outputs an
unrooted tree (e.g., ML, NJ), then the bootstrap values relate to the
internal branches of the tree which define bipartitions (also known as
splits). Following the review by Czech et al.\ \cite{Czech2017}, two
modifications were done in \ape. First, the option \code{edgelabel}
was added to the function \code{root} so that bootstrap values,
normally attached to the nodes (see \code{?boot.phylo}) are considered
attached to the edges. Second, the function
\code{drawSupportOnEdges}, derived from \code{edgelabels()}, was
written so that bootstrap values are correctly printed on the internal
branches of a tree.

\subsection{Function \code{kronoviz}}

\code{kronoviz()} plots a series of ultrametric trees so that all tips
are aligned on the same line and the scale is the same for all trees:

<<fig=true>>=
TR <- replicate(10, rcoal(sample(11:20, size = 1)), simplify = FALSE)
kronoviz(TR, type = "c", show.tip.label = FALSE)
@

\subsection{Function \code{plotTreeTime}}

If some dates are associated with the tips of a tree (e.g., from an
epidemiological study), the function \code{plotTreeTime} can plot both
information where the branch lengths are in their own units and the
$x$-axis is drawn as a time axis. For example, with a simulated random
coalescent tree and dates taken from the leap seconds in \R:

<<fig=true>>=
dates <- as.Date(.leap.seconds)
tr <- rcoal(length(dates))
plotTreeTime(tr, dates)
@

\section{Geometry and Composite Figures}\label{sec:geom}

\subsection{Single Plotting Region}

We now see some details of the geometry of a figure when plotting
a tree. We first set the outer margins (graphical parameter
\code{oma}) which are by default zero and then display the inner and
outer margins by calling \code{box()} three times with the appropriate
options:

<<fig=true>>=
par(oma = rep(2, 4))
plot(mytr)
box(lty = 3)
box("figure", col = "blue")
for (i in 1:4)
    mtext("Outer margin", i, 0.5, outer = TRUE, font = 2, col = "blue")
box("outer")
@

\noindent It is, of course, also possible to changer the inner margins
with \code{par(mar = ....)}. Note that doing \code{plot(mytr,
  no.margin = TRUE)} is identical to calling \code{par(mar = rep(0,
  4))} before doing \code{plot(mytr)}.

\subsection{Multiple Layout}

It is possible to draw several trees on the same figure by first
splitting the graphical device with \code{layout}. In that case, each
plot will have its own set of inner margins (which can be different
like in the next example) whereas the outer margins are the ``global''
margins of the figure:

<<fig=true>>=
layout(matrix(1:4, 2, 2))
par(oma = rep(2, 4))
for (i in 1:4) {
    par(mar = rep(i, 4))
    plot(mytr)
    box(lty = 3)
    box("figure", col = "blue")
}
for (i in 1:4)
    mtext("Outer margin", i, 0.5, outer = TRUE, font = 2, col = "blue")
box("outer")
@

An alternative to split a graphical device is to call \code{par(mfcol
  = c(2, 2))}. However, I recommend to use \code{layout} which is much
more flexible and powerful:

\begin{itemize}
\item \code{layout} has the options \code{width} and \code{height}
  which give the (relative) sizes of the rows and columns of plots.
\item The matrix given as main argument to \code{layout} is a symbolic
  representation of the plot layout: a single plot can span on several
  cells of the matrix; for instance, \code{layout(matrix(c(1, 2, 1,
    3), 2, 2))} will specify three plots on the device with the first
  one spanning on the first row, and the second row split into two
  plots:
\end{itemize}

<<fig=true>>=
layout(matrix(c(1, 2, 1, 3), 2, 2))
for (i in 1:3) plot(rtree(5))
@

\section{Building Your Own Code and Functions}\label{sec:build}

In this section, we examine two ways to build code or functions using
\ape's graphical functionalities.

\subsection{Getting the Recorded Parameters}

Every time \code{plot.phylo} is called, some data and parameters are
saved in a specific \R\ environment\footnote{Here ``\R\ environment''
  means a part of the memory of the computer used by \R\ with a set of
  objects.} which are then used by functions such as
\code{nodelabels}. This environment is called \code{.PlotPhyloEnv} and
is accessible to all users so that the list \code{last\_plot.phylo} can
be obtained with:


<<>>=
(pp <- get("last_plot.phylo", envir = .PlotPhyloEnv))
@

\noindent Most of the recorded parameters correspond to the arguments
of \code{plot.phylo}. The last two vectors stored in this list,
\code{xx} and \code{yy}, are the coordinates of the tips and nodes of
the tree.

\subsection{Overlaying Layouts}\label{sec:overlay}

The \pkg{graphics} package does not have an easy way to interact with
the plots on a multiple layout. For instance, after calling
\code{layout} and after the second plot has been made, it not
possible to add elements (points, arrows, \dots) to the first
plot which limits the possibility of connecting these plots in a
graphical way. However, the command \code{par(new = TRUE)} offers a way
to do that as described here.

The idea is to keep track of the coordinates of the
nodes (and other graphical elements if needed) on the whole graphical
device. However, because the plots are very likely to use different
coordinate systems (e.g., due to different scales of the branch
lengths) these coordinates must be converted into physical units on the graphical device
(usually inches). This is possible thanks to \code{par()} which
returns, among many parameters, the physical dimensions of the
plotting region, the margins, and the device in inches with the
parameters \code{"pin"}, \code{"mai"}, and \code{"din"}, respectively
(it also has \code{"fin"} for the size of the figure, which
excludes the outer margins, and will be needed if these outer
margins are not zero). An additional graphical parameter that
we will use here is given by \code{par("usr")} giving the extremes of the
coordinates of the plotting region.

The procedure follows these steps:

\begin{enumerate}
\item Rescale the node coordinates of each tree plot so they are
  between 0 and 1.
\item Convert these coordinates in inches by multiplying them with
  the relevant value in \code{"fin"} and adding the size (also in
  inches) of the relevant margin. This gives the coordinates in inches
  in each figure (i.e., including the margins).
\item Use the geometry eventually defined by \code{layout} to add to
  the coordinates output at the previous step.
\end{enumerate}

After these three steps have been performed, we obtain the coordinates
of the nodes in inches over the whole device independently of the
various scales used in each plot. Then, the call to \code{par(new =
  TRUE)} makes like nothing has been plotted and \R\ will not delete what
is already drawn on the device. We then call \code{plot(NA)} with the
options \code{type = "n", ann = FALSE, axes = FALSE} which results in a
``bare'' plot. The important options here are \code{[x|y]lim = 0:1}
which set the scale on each axis, and \code{[x|y]axs = "i"} which set
the axes to extend exactly to the data range or the specified limits.\footnote{By default,
  \code{xaxs} and \code{yaxs} are equal to \code{"r"} which adds 4\% on each side
  of the data range or limits.}

Unfortunately, in a vignette like the present document, each call to
\code{plot} opens a new device so that it is not possible to execute
the above steps interactively. Instead, the next bloc of commands
includes some commentaries and prints some intermediate results. In
this example, two random trees are plotted and an arrow connects the
first tip of the first tree to the root of the second tree:

<<fig=true>>=
## split the device into 2:
layout(matrix(1:2, 2))
## plot the 1st tree
plot(rtree(4))
## get the parameters of this first plot:
pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
## keep the coordinates of the first tip:
x <- pp$xx[1]
y <- pp$yy[1]
## extremes of the coordinates on both axes:
(pu <- par("usr"))
## fraction of the plot region:
(x - pu[1]) / (pu[2] - pu[1])
(y - pu[3]) / (pu[4] - pu[3])
## get the dimensions of plotting region and margins in inches:
pi <- par("pin")
mi <- par("mai")
## convert the coordinates into inches:
xi1 <- (x - pu[1]) / (pu[2] - pu[1]) * pi[1] + mi[2]
yi1 <- (y - pu[3]) / (pu[4] - pu[3]) * pi[2] + mi[1]

## xi1 and yi1 are the final coordinates of this tip in inches

## plot the 2nd tree:
plot(rtree(4))
## same as above:
pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
## ... except we take the coordinates of the root:
x <- pp$xx[5]
y <- pp$yy[5]
pu <- par("usr")
pi <- par("pin")
mi <- par("mai")
xi2 <- (x - pu[1]) / (pu[2] - pu[1]) * pi[1] + mi[2]
yi2 <- (y - pu[3]) / (pu[4] - pu[3]) * pi[2] + mi[1]

## xi2 and yi2 are the final coordinates of this root in inches

## we add the height of this second plot to the 'y' coordinate
## of the first tip of the first tree which is above the second
## one according to layout()
yi1 <- yi1 + par("fin")[2]
## => this operation depends on the specific layout of plots

## get the dimension of the device in inches:
di <- par("din")

## reset the layout
layout(1)
## set new = TRUE and the margins to zero:
par(new = TRUE, mai = rep(0, 4))
## set the scales to be [0,1] on both axes (in user coordinates):
plot(NA, type = "n", ann = FALSE, axes = FALSE, xlim = 0:1,
     ylim = 0:1, xaxs = "i", yaxs = "i")
## graphical elements can now be added:
fancyarrows(xi1/di[1], yi1/di[2], xi2/di[1], yi2/di[2], 1,
            lwd = 10, col = rgb(1, .5, 0, .5), type = "h")
@

\noindent This code must be adapted to the exact layout in order to
shift appropriately the coordinates in the $x$- and/or
$y$-direction(s). The next bloc defines a function which uses the above
code to plot two trees and add an arrow from the tip specified with
the argument \code{from} of the first tree to the root of the second
tree:

<<>>=
foo <- function(phy1, phy2, from)
{
    layout(matrix(1:2, 2))
    plot(phy1, font = 1)
    pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
    from <- which(phy1$tip.label == from)
    x <- pp$xx[from]
    y <- pp$yy[from]
    ## fraction of the plot region:
    pu <- par("usr")
    ## convert into inches:
    pi <- par("pin")
    mi <- par("mai")
    xi1 <- (x - pu[1]) / (pu[2] - pu[1]) * pi[1] + mi[2]
    yi1 <- (y - pu[3]) / (pu[4] - pu[3]) * pi[2] + mi[1]
    plot(phy2)
    pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
    to <- Ntip(phy2) + 1
    x <- pp$xx[to]
    y <- pp$yy[to]
    ## same as above:
    pu <- par("usr")
    pi <- par("pin")
    mi <- par("mai")
    xi2 <- (x - pu[1]) / (pu[2] - pu[1]) * pi[1] + mi[2]
    yi2 <- (y - pu[3]) / (pu[4] - pu[3]) * pi[2] + mi[1]
    yi1 <- yi1 + par("fin")[2]
    di <- par("din")
    layout(1)
    par(new = TRUE, mai = rep(0, 4))
    plot(NA, type = "n", ann = FALSE, axes = FALSE, xlim = 0:1,
         ylim = 0:1, xaxs = "i", yaxs = "i")
    fancyarrows(xi1/di[1], yi1/di[2], xi2/di[1], yi2/di[2], 1,
                lwd = 10, col = rgb(1, .5, 0, .5), type = "h")
}
@

\noindent We try this function with a tree of mammalian orders that we
want to connect with our small tree \code{mytr}:

<<fig=true>>=
trb <- read.tree(text = "((Primates,Carnivora),Monotremata);")
par(xpd = TRUE)
foo(trb, mytr, "Primates")
@

\noindent Note that the first tree has no branch length, so both
trees have obviously different scales.

Another use for these functionalities is to draw coloured rectangles
to delimit clades. Again we can use \code{par(new = TRUE)} this
time with the option \code{plot = FALSE} of \code{plot.phylo} which
opens and set the graphical device with the coordinates taken from the
tree as detailed in Section~\ref{sec:comput}, but nothing is drawn. The user can
then call any low-level plotting command, then use \code{plot()} after
\code{par(new = TRUE)}. For instance, if we want to draw a rectangle
to show the clade made by humans and chimpanzees:

\begin{center}
\setkeys{Gin}{width=.5\textwidth}
<<fig=true,width=4,height=3.5>>=
plot(mytr, plot = FALSE)
pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
rect(pp$xx[5] - 0.1, pp$yy[1] - 0.1, pp$xx[1] + 2, pp$yy[2] + 0.1,
     col = "yellow", border = NA)
par(new = TRUE)
plot(mytr)
@
\end{center}
There are several ways to find the limits of the rectangle: trying
different values empirically, using \code{locator()} on a first draft
plot with the tree, or with calculations similar to the previous
example. The examples in \code{?phydataplot} gives other examples of
using \code{par(new = TRUE)}.

\bibliographystyle{plain}
\bibliography{ape}
%\setlength{\bibsep}{0pt}
\addcontentsline{toc}{section}{References}

\end{document}