%\VignetteIndexEntry{Advanced features}
%\VignetteKeywords{Documentation}
%\VignettePackage{phangorn}
\documentclass[12pt]{article}
% setwd("/home/kschliep/Desktop/phangorn/inst/doc")
% Sweave("phangorn-specials.Rnw")
% tools::texi2dvi("phangorn-specials.tex", pdf=TRUE)
\usepackage{times}
\usepackage{hyperref}

\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\textit{#1}}}
\newcommand{\Rmethod}[1]{{\texttt{#1}}}
\newcommand{\Rfunarg}[1]{{\texttt{#1}}}
\newcommand{\Rclass}[1]{{\textit{#1}}}

\textwidth=6.2in
\textheight=8.5in
%\parskip=.3cm
\oddsidemargin=.1in
\evensidemargin=.1in
\headheight=-.3in

\newcommand{\R}{\textsf{R}}
\newcommand{\pml}{\Robject{pml}}
\newcommand{\phangorn}{\Rpackage{phangorn}}
\newcommand{\ape}{\Rpackage{ape}}
\newcommand{\multicore}{\Rpackage{multicore}}

\newcommand{\term}[1]{\emph{#1}}
\newcommand{\mref}[2]{\htmladdnormallinkfoot{#2}{#1}}



\begin{document}

% Ross Ihakas extenstion for nicer representation 
\DefineVerbatimEnvironment{Sinput}{Verbatim} {xleftmargin=2em}
\DefineVerbatimEnvironment{Soutput}{Verbatim}{xleftmargin=2em}
\DefineVerbatimEnvironment{Scode}{Verbatim}{xleftmargin=2em}
\fvset{listparameters={\setlength{\topsep}{0pt}}}
\renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}}

<<echo=FALSE>>=
options(width=80)
foo <- packageDescription("phangorn")
@


\title{Special features of phangorn (Version \Sexpr{foo$Version})} %$
\author{\mref{mailto:kschliep@snv.jussieu.fr}{Klaus P. Schliep}}
\maketitle

\nocite{Paradis2006}
\section{Introduction}
This document illustrates some of the \phangorn{} some specialised features which are useful but maybe not as wellknown or just not (yet) described elsewhere. This is mainly interesting for someone who wants to explore different models or set up some simulation studies. We show how to construct data objects for different character states other than nucleotides or amino acids or how to set up different models to estimate transition rate. 
\section{User defined data formats}\label{sec:USER}
The vignette \emph{Trees} describes in detail how to estimate phylogenies from nucleotide or amino acids. 

To better understand how to define our own data type it is useful to know a bit about the internal representation of \Robject{phyDat} objects. The internal representation of \Robject{phyDat} object is very similar to \Robject{factor} objects. 
  
As an example we will show here several possibilities to define nucleotide data with gaps defined as a fifth state. When the number of gaps is low and the are missing at random this may be not important. 
 
Let assume we have given a matrix where each row contains a character vector of a taxonomixcal unit:
<<echo=TRUE>>=
library(phangorn)
data = matrix(c("r","a","y","g","g","a","c","-","c","t","c","g", 
                "a","a","t","g","g","a","t","-","c","t","c","a",                                          
                "a","a","t","-","g","a","c","c","c","t","?","g"), 
                dimnames = list(c("t1", "t2", "t3"),NULL), nrow=3, byrow=TRUE)
data
@
Normally we would transform this matrix into an phyDat object and gaps are handled as ambiguious character like "?".  
<<>>=
gapsdata1 = phyDat(data)
gapsdata1
@
Now we will define a "USER" defined object and have to supply a vector levels of the character states for the new data, in our case the for nucleotide states and the gap. Additional we can define ambigious states which can be any of the states. 
<<echo=TRUE>>=
gapsdata2 = phyDat(data, type="USER", levels=c("a", "c", "g", "t", "-"), ambiguity = c("?", "n"))
gapsdata2
@
This is not yet what we wanted as two sites of our alignment, which contain the abigious characters "r" and "y", got deleted.  
%<<echo=TRUE>>=
%gapsdata3 = phyDat(data, type="USER", levels=c("a", "c", "g", "t", "-"), ambiguity = c("?", "n"))
%gapsdata3 
%@
To define ambigious characters like "r" and "y" explicitly we have to supply a contrast matrix similar to contrasts for factors. 
<<echo=TRUE>>=
contrast = matrix(data = c(1,0,0,0,0,
    0,1,0,0,0,
    0,0,1,0,0,
    0,0,0,1,0,   
    1,0,1,0,0,
    0,1,0,1,0,
    0,0,0,0,1,
    1,1,1,1,0,
    1,1,1,1,1),
    ncol = 5, byrow = TRUE)
dimnames(contrast)  = list(c("a", "c", "g", "t", "r", "y", "-", "n", "?"), c("a", "c", "g", "t", "-"))
contrast
gapsdata3 = phyDat(data, type="USER", contrast=contrast)
gapsdata3 
@
Here we defined "n" as a state which can be any nucleotide but not a gap "-" and "?" can be any state including a gap.

These data can be used in all functions available in \phangorn{} to compute distance matrices or perform parsimony and maximum likelihood analysis.  


\section{Estimation of non-standard transition rate matrices}
In the last section~\ref{sec:USER} we described how to set up user defined data formats. Now we describe how to estimate transition matrices with pml. 

Again for nucleotide data the most common models can be called directly in the optim.pml function (e.g. "JC69", "HKY", "GTR" to name a few). Table \ref{models} lists all the available nucleotide models, which can estimated directly in \Robject{optim.pml}. For amino acids several transition matrices are available ("WAG", "JTT", "Dayhoff" and "LG") or can be estimated with optim.pml, e.g. Mathews et al. (2010) \cite{Mathews2010} used this function to estimate a phytochrome amino acid transition matrix. 


We will now show how to estimate a rate matrix with different transition ($\alpha$) and transversion ratio ($\beta$) and a fixed rate to the gap state ($\gamma$) - a kind of Kimura two-parameter model (JC69) for nucleotide data with gaps as fifth state (see table ~\ref{gaps}). 

\begin{table}[htbp]
   \centering
   \begin{tabular}{l|lllll}   
    & a & c & g & t & - \\
   \hline
   a & & & & & \\
   c & $\beta$ & & & & \\
   g & $\alpha$ & $\beta$ & & & \\
   t & $\beta$ & $\alpha$ & $\beta$ & & \\
   - & $\gamma$ & $\gamma$ & $\gamma$ & $\gamma$ & \\   
   \end{tabular}
   \caption{Rate matrix K to optimise. }\label{gaps} 
\end{table}


The parameters subs accepts a vector of consecutive integers and one (or more) element has to be zero (these gets the reference rate of 1).
<<>>=
tree = unroot(rtree(3))
fit = pml(tree, gapsdata3)
fit = optim.pml(fit, optQ=TRUE, subs=c(1,0,1,2,1,0,2,1,2,2), control=pml.control(trace=0))
fit
@

 
Here are some conventions how the models are estimated: \\

If a model is supplied the base frequencies bf and rate matrix Q are optimised according to the model (nucleotides) or the adequate rate matrix and frequencies are chosen (for amino acids). 
If optQ=TRUE and neither a model or subs are supplied than a symmetric (optBf=FALSE) or reversible model (optBf=TRUE, i.e. the GTR for nucleotides) is estimated.  This can be slow if the there are many character states, e.g. for amino acids.

 
\begin{table}[htbp]
   \centering
   \begin{tabular}{|llll|}
   \hline
   model & optQ & optBf & subs \\
   \hline
         JC & FALSE & FALSE & $c(0, 0, 0, 0, 0, 0)$ \\
         F81 & FALSE & TRUE & $c(0, 0, 0, 0, 0, 0)$ \\
         K80 & TRUE & FALSE & $c(0, 1, 0, 0, 1, 0)$ \\
         HKY & TRUE & TRUE  & $c(0, 1, 0, 0, 1, 0)$ \\
         TrNe & TRUE & FALSE & $c(0, 1, 0, 0, 2, 0)$ \\
         TrN & TRUE & TRUE  & $c(0, 1, 0, 0, 2, 0)$ \\
         TPM1 & TRUE & FALSE & $c(0, 1, 2, 2, 1, 0)$ \\
         K81 & TRUE & FALSE & $c(0, 1, 2, 2, 1, 0)$ \\
         TPM1u & TRUE & TRUE & $c(0, 1, 2, 2, 1, 0)$ \\
         TPM2 & TRUE & FALSE & $c(1, 2, 1, 0, 2, 0)$ \\
         TPM2u & TRUE & TRUE & $c(1, 2, 1, 0, 2, 0)$ \\
         TPM3 & TRUE & FALSE & $c(1, 2, 0, 1, 2, 0)$ \\
         TPM3u & TRUE & TRUE & $c(1, 2, 0, 1, 2, 0)$ \\
         TIM1e & TRUE & FALSE & $c(0, 1, 2, 2, 3, 0)$ \\
         TIM1 & TRUE & TRUE & $c(0, 1, 2, 2, 3, 0)$ \\
         TIM2e & TRUE & FALSE & $c(1, 2, 1, 0, 3, 0)$ \\
         TIM2 & TRUE & TRUE & $c(1, 2, 1, 0, 3, 0)$ \\
         TIM3e & TRUE & FALSE & $c(1, 2, 0, 1, 3, 0)$ \\
         TIM3 & TRUE & TRUE & $c(1, 2, 0, 1, 3, 0)$ \\
         TVMe & TRUE & FALSE & $c(1, 2, 3, 4, 2, 0)$ \\
         TVM & TRUE & TRUE & $c(1, 2, 3, 4, 2, 0)$ \\
         SYM & TRUE & FALSE & $c(1, 2, 3, 4, 5, 0)$ \\
         GTR & TRUE & TRUE & $c(1, 2, 3, 4, 5, 0)$ \\
         \hline
   \end{tabular}
   \caption{DNA models available in phangorn and how they are defined. }\label{models} 
\end{table}


  

\section{Generating trees}
\phangorn{} has a few functions to generate trees, which may are interesting for simulation studies. \Rfunction{allTrees} computes all possible bifurcating tree topologies either rooted or unrooted for up to 10 taxa (one has to keep in mind that the number of trees is growing exponentially). 

<<echo=TRUE>>=
trees = allTrees(5)
@

<<label=plotAll,include=FALSE>>=
par(mfrow =c(3,5), mar = c(2,2,2,2)-2)
for(i in 1:15)plot(trees[[i]], cex=1, type="u")
@
\begin{figure}
\begin{center}
<<label=figAll,fig=TRUE,echo=FALSE>>=
<<plotAll>>
@
\end{center}
\caption{all (15) unrooted trees with 5 taxa}
\label{fig:NJ}
\end{figure}



\Rfunction{rNNI} and \Rfunction{rSPR} generate trees which are a defined number of NNI (nearest neighbor interchange) or SPR (subtree pruning and regrafting) away.
  


\bibliographystyle{plain}
\bibliography{phangorn}

\section{Session Information}
The version number of \R{} and packages loaded for generating the vignette were:
<<echo=FALSE,results=tex>>=
toLatex(sessionInfo())
@

\end{document}
