\documentclass{article}
%\VignettePackage{OutbreakTools}
%\VignetteIndexEntry{Handling disease outbreak data}
%\VignetteEngine{knitr}

\usepackage{graphicx}
\usepackage[colorlinks=true,urlcolor=blue]{hyperref}
\usepackage{array}
\usepackage{color}
\usepackage{geometry}
\geometry{verbose,tmargin=3cm,bmargin=3cm,lmargin=3cm,rmargin=3cm}

\usepackage[utf8]{inputenc} % for UTF-8/single quotes from sQuote()
\newcommand{\code}[1]{{{\tt #1}}}
\title{Handling disease outbreak data using \textit{OutbreakTools} \Sexpr{packageDescription("OutbreakTools", fields = "Version")}}
\author{Thibaut Jombart \and Xavier Didelot \and Rolf Ypma \and Lulla Opatowski \and Anne Cori \and Marc A.~Suchard}
\date{\today}




\sloppy
\hyphenpenalty 10000


\begin{document}
%\SweaveOpts{concordance=TRUE}


<<eval=TRUE,echo=FALSE,results='hide'>>=
library(knitr)
opts_chunk$set(fig.path='figs/OutbreakTools-', fig.keep='last', dev='pdf', fig.width=7, fig.height=7,
               tidy=FALSE, warning=FALSE, fig.show='asis', fig.align='center', cache=FALSE,
               out.width=".6\\textwidth")
options(width=80)
@

\color{black}

\maketitle

\begin{abstract}
  This vignette introduces the main functionalities of \textit{OutbreakTools}, a package implementing core
  tools for the analysis of outbreak data.  Disease outbreak data can be diverse and complex, and
  the purpose of \textit{OutbreakTools} is to simplify the handling of this information.  The main feature
  of the package lies in the formal (S4) class \texttt{obkData} (for ``outbreak data''), which
  offers a coherent way of handling data on individuals, samples, contact networks, clinical events,
  as well as phylogenies and genomic sequences. Beyond introducing this data structure, this
  tutorial illustrates how these objects can be handled and visualized in R.
\end{abstract}

\newpage

\tableofcontents


\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Storing outbreak data}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In this section, we first detail the structure of the classes of objects used in
\textit{OutbreakTools}, and then explain how to import data into the package.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Class definitions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Data collected during outbreaks can be hugely diverse and complex.  In \textit{OutbreakTools}, our purpose
is to have a general class of objects which can store virtually any information sampled during
an outbreak, without the user worrying about storage issues and consistency amongst different types of
data.  For most purposes, the core class \texttt{obkData} can be taken as a black box, with which the user can interact using specific
functions called \textit{accessors}. However, a basic understanding of what type of information is stored in these objects will be useful.



% % % % % % % % % % % % % % % % %
\subsubsection{\code{obkData}: storage of outbreak data}
% % % % % % % % % % % % % % % % %
The main class of objects in \textit{OutbreakTools} is \texttt{obkData}.  This formal (S4) class is used
to store various types of information gathered during outbreaks.  The definition of the class in
terms of R objects can be obtained by:
<<results='hide'>>=
library(OutbreakTools)
@
<<>>=
getClassDef("obkData")
@
One can also examine a structure using an empty object:
<<>>=
new("obkData")
@

Each slot of an \texttt{obkData} object is optional. By convention, empty slots are always \texttt{NULL}.
The slots respectively contain:
\begin{itemize}
\item \texttt{@individuals}: a \texttt{data.frame} storing individual data, such as age, sex,
  or onset of symptoms. If not \texttt{NULL}, this \texttt{data.frame} will have exactly one row per
  individual, with row names providing unique identifiers for individuals.
\item \texttt{@records}: a named list of \texttt{data.frame}s storing any time-stamped data gathered at an
  individual level; there is no constraint on the number of \texttt{data.frame}s stored, but each one must contain
  columns named \texttt{individualID} (unique identifiers for individuals) and
  \texttt{date}. Examples: swab data, fever, onset of symptoms, etc.
\item \texttt{@dna}: DNA sequences of one or more genes, stored as an \texttt{obkSequences}
  object. See section below for details on \texttt{obkSequences} objects.
\item \texttt{@contacts}: dynamic contact network between the individuals, stored as an
  \texttt{obkContacts} object. See section below for details on \texttt{obkContacts} objects.
\item \texttt{@context}: a list of \texttt{data.frame}s storing any time-stamped data at a
  non-individual level. Examples: climatic variables, school closures, vaccination campaign, etc.
\item \texttt{@trees}: a list of phylogenetic trees with the class \texttt{multiPhylo} (from the
  \textit{ape} package); can be used to store e.g. a posterior distribution of trees from a Bayesian
  phylogenetic reconstruction using BEAST.
\end{itemize}

The slots of an object \texttt{foo} can be accessed using \texttt{foo@[name-of-the-slot]}.
Let us use the toy outbreak dataset \texttt{ToyOutbreak} and examine its content:
<<>>=
data(ToyOutbreak)
class(ToyOutbreak)
slotNames(ToyOutbreak)
head(ToyOutbreak)
summary(ToyOutbreak)
@

\noindent \texttt{ToyOutbreak} is an \texttt{obkData} object containing information on individuals (\texttt{@individuals}),
samples/records made on individuals (\texttt{@records}), DNA sequences (\texttt{@dna}), a contact
network (\texttt{@contacts}) and one or more phylogenetic trees (\texttt{@trees}).
Accessing a given slot is as easy as:
<<>>=
head(ToyOutbreak@individuals)
head(ToyOutbreak@records$Fever)
ToyOutbreak@trees
@
However, we will see how retrieving information from \texttt{obkData} objects can be made
more powerful using accessors in the following sections.



% % % % % % % % % % % % % % % % %
\subsubsection{\texttt{obkSequences}: storage of DNA sequences for different genes}
% % % % % % % % % % % % % % % % %

Pathogen sequence data can typically be obtained for different genes, making the handling of such
information not entirely trivial: different individuals may have been sequenced for different genes,
at different points in time, etc.
The class \texttt{obkSequences} stores such information.
\texttt{obkSequences} objects contain two slots: \texttt{@dna} and \texttt{@meta}.


The slot \texttt{@dna} is a list of matrices of aligned DNA sequences (in rows), stored using \textit{ape}'s class
\texttt{DNAbin} for efficiency, with each item of the list corresponding to a different gene.
Gene names are the names of the list.
The row names in each matrix contain unique identifiers for the sequences, typically accession numbers.


The slot \texttt{@meta} is a \texttt{data.frame} containing some meta-information about the sequences.
It contains at least two columns for sampled individuals (\texttt{individualID}) and collection
dates (\texttt{date}).
The row names corrspond to sequence labels used in \texttt{@dna}, and respect the same ordering.
\\

Let us examine the DNA information stored in \texttt{ToyOutbreak}:
<<>>=
class(ToyOutbreak@dna)
ToyOutbreak@dna
slotNames(ToyOutbreak@dna)
is.list(ToyOutbreak@dna@dna)
names(ToyOutbreak@dna@dna)
ToyOutbreak@dna@dna$gene1
class(ToyOutbreak@dna@dna$gene1)
class(ToyOutbreak@dna@meta)
head(ToyOutbreak@dna@meta)
@

\noindent \texttt{ToyOutbreak@dna} is an \texttt{obkSequences} object containing DNA sequences for two genes.
The slot \texttt{ToyOutbreak@dna@dna} is a list of \texttt{DNAbin} matrices, each containing
sequences for a given gene.





% % % % % % % % % % % % % % % % %
\subsubsection{\code{obkContacts}: storage of dynamics contact networks}
% % % % % % % % % % % % % % % % %

\textit{obkData} objects can also store contact data between individuals, in the slot \texttt{@contacts}.
These contacts can be fixed or vary in time, in which case data are stored as a dynamic contact network.
The slot \texttt{@contacts} is an instance of the class \texttt{obkContacts}, which currently
contains either a \texttt{network} object (static graph, from the \textit{network} package), or a \texttt{networkDynamic} object, for
contacts varying in time (from the \textit{networkDynamic} package).
These objects are fully documented in their respective vignettes.
Here, we detail a simple toy example from the documentation of \texttt{obkContacts}:
<<>>=
cf <- c("a", "b", "a", "c", "d")
ct <- c("b", "c", "c", "d", "b")
oc.static <- new("obkContacts", cf, ct, directed=FALSE)
slotNames(oc.static)
oc.static
@

\texttt{oc.static} contains a static, non-directed contact network (slot \texttt{@contacts}, class \texttt{network}).
It can be plotted easily using:
<<graphstat,out.width=".7\\textwidth">>=
plot(oc.static, main="Static contact network")
@
<<>>=
onset <- c(1, 2, 3, 4, 5)
terminus <- c(1.2, 4, 3.5, 4.1, 6)
oc.dynamic <- new("obkContacts",cf,ct, directed=FALSE,
                  start=onset, end=terminus)
slotNames(oc.dynamic)
oc.dynamic
@

\texttt{oc.dynamic} is a dynamic graph, i.e. a graph whose vertices and edges can change over time.
By default, plotting the object collapses the graph so that all vertices and edges that exist at
some point are displayed; however, sections of the graph for given time intervals can be obtained
using \texttt{get.contacts} (or alternatively, \texttt{network.extract} on the
\texttt{networkDynamic} object directly).
As a reminder, here is the input of the graph \texttt{oc.dynamic}:
<<>>=
as.data.frame(oc.dynamic)
@
And here are various plots, first of the full (collapsed) contact network, then for different time
intervals (0--2, 2--4, 4--6):
<<dynNet,out.width=".9\\textwidth">>=
par(mfrow=c(2,2))
plot(oc.dynamic@contacts,main="oc.dynamic - collapsed graph",
     displaylabels=TRUE)
plot(get.contacts(oc.dynamic, from=0, to=2),
     main="oc.dynamic - time 0--2", displaylabels=TRUE)
plot(get.contacts(oc.dynamic, from=2, to=4),
     main="oc.dynamic - time 2--4", displaylabels=TRUE)
plot(get.contacts(oc.dynamic, from=4, to=6),
     main="oc.dynamic - time 4--6", displaylabels=TRUE)
@

\textit{networkDynamic} allows for extensive manipulation of dynamic networks.
For more information, refer to the vignette distributed with the package (\texttt{vignette("networkDynamic")}).







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Getting data into \textit{OutbreakTools}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Storing data in \textit{OutbreakTools} requires the following, fairly simple steps:
\begin{enumerate}
\item read data into R
  \begin{enumerate}
  \item read \texttt{data.frame}s storing individuals, samples, and clinical information in R from a
    text file, typically using \texttt{read.table} or \texttt{read.csv} for comma-separated
    files. Every standard spreadsheet software can export data to these formats.
  \item read DNA sequences from separate files containing alignments (one file per gene), typically using \texttt{read.dna} from the ape
    package. While phylogenies can be obtained in R, annotated trees
    produced by Bayesian software such as BEAST can now be imported using
    \texttt{read.annotated.nexus}.
  \end{enumerate}
\item use this information as input to the \texttt{obkData} constructor (\texttt{new("obkData",\ldots)}) to create an
  \texttt{obkData} object.
\end{enumerate}

In the following, we assume that step 1 is sorted and focus on step 2: using the constructor.



% % % % % % % % % % % % % % % % %
\subsubsection{The \code{obkData} constructor}
% % % % % % % % % % % % % % % % %

New objects are created using \texttt{new}, with these slots as arguments.
If no argument is provided, an empty object is created, as seen before:
<<>>=
new("obkData")
@

\noindent This function accepts the following arguments, which mirror to some extent the structure
of the object (see \texttt{?obkData} for more information):
\begin{itemize}
 \item \texttt{individuals}: {a data.frame with a mandatory column named
    \texttt{individualID}, providing unique identifiers for the individuals; if
    missing, row names are used as identifiers.}

  \item \texttt{records}: a list of data.frames, each of which has 2 mandatory
    fields, \texttt{individualID} and \texttt{date}. Dates can be specified as
    \texttt{Date} or \texttt{characters}, in which case they will be
    converted to dates. Most sensible formats will be detected
    automatically and processed. Unusual formats should be provided
    through the argument \texttt{date.format}.  Each item of the list
    should be named according to the type of information recorded,
    e.g. 'swabs', 'temperature', or 'hospitalisation' (admission /
    discharge events).

  \item \texttt{dna}: a list matrices of DNA sequences in \texttt{DNAbin} or
    \texttt{character} format, each component of the list being a
    different gene. A matrix can be provided if there is a single gene.

  \item \texttt{dna.date}: a vector of collection dates for the DNA sequences;
    see \texttt{obkSequences} manpage for more information.

  \item \texttt{dna.individualID}: a vector of individual from which DNA
    sequences where obtained; see \texttt{obkSequences} manpage for
    more information.

  \item \texttt{dna.date.format}: a character string indicating the format of the
    date in \texttt{dna.date} if ambiguous; see \texttt{obkSequences}
    manpage for more information.

  \item \texttt{dna.sep}: the character string used to separate fields
    (e.g. sequenceID/individualID/date) in sequences labels; see
    \texttt{obkSequences} manpage for more information.

  \item \texttt{contacts}: a matrix of characters indicating contacts using two
    columns; if contacts are directed, the first column is 'from', the
    second is 'to'; values should match individual IDs (as returned by
    get.individuals(x)); if numeric values are provided, these are
    converted to integers and assumed to correspond to individuals
    returned by get.individuals(x).

  \item \texttt{context}: a list of data.frames, each of which has 1 mandatory
    field: \texttt{date}. Each item of the list should be named according to
    the type of information recorded, e.g. 'intervention',
    'vaccination', 'climat' (temperature, humidity, etc.), or schools
    (opening/closure).

  \item \texttt{contacts.start}: a vector of dates indicating the beginning of
    each contact.

  \item \texttt{contacts.end}: a vector of dates indicating the end of each
    contact.

  \item \texttt{contacts.duration}: another way to specify contacts.end, as
    duration of contact in days.

  \item \texttt{contacts.directed}: a logical indicating if contacts are
    directed; defaults to FALSE.

  \item \texttt{trees}: a list of phylogenetic trees with the class
    \texttt{multiPhylo} (from the ape package)

  \item \texttt{date.format}: a character string indicating the date format (see
    \texttt{as.Date}); if \texttt{NULL}, date format is detected
    automatically, which is usually a sensible option.

\end{itemize}

\noindent We can now show how to create a new \texttt{obkData} from multiple inputs, using the
dataset \texttt{ToyOutbreakRaw}:
<<>>=
data(ToyOutbreakRaw)
class(ToyOutbreakRaw)
names(ToyOutbreakRaw)
@

Here is an overview of the inputs, including data on individuals:
<<>>=
head(ToyOutbreakRaw$individuals)
@
various time-stamped records:
<<>>=
lapply(ToyOutbreakRaw$records, head)
@
contact information:
<<>>=
head(ToyOutbreakRaw$contacts)
head(ToyOutbreakRaw$contacts.start)
head(ToyOutbreakRaw$contacts.end)
@
DNA sequences:
<<>>=
ToyOutbreakRaw$dna
@
and phylogenetic trees:
<<>>=
ToyOutbreakRaw$trees
@

All this information will be compiled into a single object by:
<<>>=
attach(ToyOutbreakRaw)

x <- new ("obkData", individuals=individuals, records=records,
          contacts=contacts, contacts.start=contacts.start,
          contacts.end=contacts.end, dna=dna,
          dna.individualID=dna.info$individualID,
          dna.date=dna.info$date, sample=dna.info$sample, trees=trees)

detach(ToyOutbreakRaw)

head(x)
summary(x)
@

\texttt{x} is a new, coherent representation of the data.
This representation ensures, amongst other things, that:
\begin{itemize}
  \item individual labels are unique and consistent across records, contacts, DNA sequences and
    patient information
  \item every item is dated using actual dates (\texttt{Date} objects)
  \item every sample/record refers to an individual and a date
  \item every DNA sequence refers to an individual and a date
  \item every DNA sequence belongs to a gene
  \item DNA sequences from the same gene have the same length
  \item every tip of the trees refers to a DNA sequence
  \item every contact refers to documented individuals
\end{itemize}

Having all these items connected allows to simplify data manipulation fairly drastically (see section below on data handling).
For instance, it will be easy to isolate a subset of individuals from the data, which will impact
not only patient information but also the phylogenies, DNA sequences, records and contacts.
It will also be straightforward to add e.g. patient information onto phylogenies.

<<eval=TRUE,echo=FALSE,results='hide'>>=
opts_chunk$set(eval=FALSE, echo=FALSE)
@


% % % % % % % % % % % % % % % %
\subsubsection{Using other constructors: \texttt{obkSequences} and \texttt{obkContacts}}
% % % % % % % % % % % % % % % % %
The classes \texttt{obkSequences} and \texttt{obkContacts}, both used in \texttt{obkData} objects,
also have constructors and can be created independently from \texttt{obkData} objects.
However, the risk is that one would replace e.g. the DNA sequences stored in an
\texttt{obkData} object by a new \texttt{obkSequences}, which would bypass the consistency checks
made by the \texttt{obkData} constructor and possibly lead to an invalid object.
This practice is therefore discouraged for the moment.



%% % % % % % % % % % % % % % % % % %
%% \subsubsection{The \texttt{obkSequences} constructor}
%% % % % % % % % % % % % % % % % % %




%% In most cases, one will not need to construct \texttt{obkSequences} directly, this task being done
%% implicitly will creating \texttt{obkData} objects.
%% However, one might want to modify the DNA sequences stored in an existing \texttt{obkData}, thus
%% needing to build a new \texttt{obkSequences}.
%% As for \texttt{obkData}, \texttt{obkSequences} objects can be created using the constructor
%% \texttt{new("obkSequences", \ldots)}, where ``\ldots'' can be the following arguments:

%% \begin{itemize}
%% \item \texttt{dna}: a list of DNA sequences (not necessarily from the same gene) in \texttt{DNAbin}
%%   or \texttt{character} format; matrices will be accepted too if only one locus is
%%   provided. Sequences must be named using unique identifiers, typically accession
%%   numbers. Typically, this information will be obtained by reading sequence data into R using
%%   \textit{ape}'s \texttt{read.dna} function.
%% \item \texttt{locus}: an optional vector indicating the locus of each sequences; its length must match that of the list of sequences.
%% \end{itemize}

%% \noindent Using these inputs, the \texttt{obkSequences} constructor will sort out sequences per gene and store
%% them as matrices, using one matrix per gene and checking that sequences from the same gene are
%% actually of the same length.

%% Here, we illustrate the creation of \texttt{obkSequences} objects using a dataset of influenza
%% sequences from \textit{adegenet}, first read using \texttt{read.dna}:
%% <<>>=
%% path.file <- system.file("files/usflu.fasta",package="adegenet")
%% path.file
%% flu <- read.dna(path.file, format="fasta")
%% flu
%% @
%% The object is simply created using:
%% <<>>=
%% x <- new("obkSequences", dna=flu)
%% x
%% @
%% As locus information is not provided, the constructor assumed (rightfully so) that all sequences are
%% from the same (unnamed) locus.
%% Here, the sequenced segment is actually hemagglutinin (HA), so we can add this information:
%% <<>>=
%% x <- new("obkSequences", dna=flu, locus=rep("HA",80))
%% x
%% @
%% Now, if we assume that for instance, the first 70 sequences were HA, followed by 8 neuraminidase
%% (NA) and 2 nucleoprotein (NP), then we would use:
%% <<>>=
%% x <- new("obkSequences", dna=flu, locus=rep(c("HA","NA","NP"), c(70,8,2)))
%% x
%% @

%% \noindent Note that sequences do not have to be ordered by locus; the only thing that matters is
%% that the argument \texttt{locus} matches the sequences provided in \texttt{dna}.

%% Replacing the \texttt{@dna} slot of an \texttt{obkData} object is as simple as:
%% <<>>=
%% obj <- new("obkData")
%% obj@dna
%% obj@dna <- x
%% obj
%% @
%% Note however, that this operation does not ensure matching of sequences IDs in \texttt{@dna} with
%% the information provided in \texttt{@sample}.




%% % % % % % % % % % % % % % % % % %
%% \subsubsection{The \texttt{obkContacts} constructor}
%% % % % % % % % % % % % % % % % % %


<<eval=TRUE,echo=FALSE,results='hide'>>=
opts_chunk$set(eval=TRUE, echo=TRUE)
@





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Data handling using \texttt{obkData} objects}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Accessors}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The philosophy underlying formal (S4) classes is that the internal representation of the data can be
complex as long as accessing the information is simple.
This is made possible by decoupling storage and accession: the user is not meant to access the
content of the object directly, but has to use \textit{accessors} to retrieve the information.
In this section, we detail the existing accessors for object classes implemented in \textit{OutbreakTools}.
We use the notation ``\texttt{[\textit{possible-values}]}'' to list or describe possible values of an argument;
the symbols ``\texttt{[]}'' should be omitted from the actual command line. For instance:\\

\texttt{myFunction(x, y=[\textit{"foo" or "bar"}])}\\

\noindent means that the argument \texttt{y} of function \texttt{myFunction} can be either \texttt{"foo"}
or \texttt{"bar"}, and valid calls would be:
<<eval=FALSE>>=
myFunction(x, y="foo")
@
\noindent or:
<<eval=FALSE>>=
myFunction(x, y="bar")
@


% % % % % % % % % % % % % % % % %
\subsubsection{Accessors for \texttt{obkData} objects}
% % % % % % % % % % % % % % % % %

Available accessors are also documented in \texttt{?obkData}.
These functions are meant to retrieve information that is not trivially accessible.
To simply access slots, use the \texttt{@} operator, e.g. \texttt{x@samples},
\texttt{x@individuals}, etc.
\\

All accessors return \texttt{NULL} when information is missing, except for functions returning
number of items, which will return \texttt{0}.
In the following, we illustrate accessors using a random sample of 5 individuals of the toy dataset \texttt{ToyOutbreak}:
<<>>=
data(ToyOutbreak)
set.seed(1)
toKeep <- sample(get.nindividuals(ToyOutbreak),5)
toKeep
x <- subset(ToyOutbreak, individuals=toKeep)
summary(x)
@

\begin{itemize}
\item \texttt{get.individuals(x, data=[\textit{"all" or "individuals" or "records" or "contacts" or "dna" or "context"}])}: returns the
  individual IDs in different components of the object.

\item  \texttt{get.nindividuals(x, data=[\textit{"all" or "individuals" or "records" or "contacts" or "dna" or "context"}])}:
  returns the number of individuals in different components of the object.
<<>>=
get.nindividuals(x)
get.nindividuals(x, "records")
get.nindividuals(x, "dna")
get.nindividuals(x, "contacts")
@
There are 5 individuals in the data, except for contact information; this is because contacts were
only recorded between the first 20 individuals of \texttt{ToyOutbreak}:
<<>>=
get.individuals(ToyOutbreak, "contacts")
@

\item \texttt{get.nlocus(x)}: returns the number of loci.
\item \texttt{get.locus(x)}: returns the names of the loci in the data.
<<>>=
get.nlocus(x)
get.locus(x)
@

\item \texttt{get.nsequences(x, what=[\textit{"total" or "bylocus"}])}: returns the number of sequences in \texttt{@dna}.
\item \texttt{get.sequences(x)}: returns the IDs of the sequences in \texttt{@dna}.
<<>>=
get.nsequences(x)
get.nsequences(x, "bylocus")
get.sequences(x)
@


\item \texttt{get.trees(x)}: returns the content of \texttt{x@trees}.
<<>>=
get.trees(x)
@

\item \texttt{get.dna(x, locus=[\textit{locus IDs}], id=[\textit{sequence IDs}])}: returns a list of
  matrices of DNA sequences; the arguments \texttt{locus} and \texttt{id} are optional; if provided,
  they should be character strings corresponding to the name of the loci and/or sequences to be
  retained. Integers or logical will be treated as indicators based on the results of
  \texttt{get.locus} or \texttt{get.sequences}.
<<>>=
get.dna(x)
@
returns all the DNA sequences, in two matrices corresponding to the different genes.
We can request e.g. only the second gene:
<<>>=
get.dna(x, locus=2)
@
or even just specific sequences, say ("311" and "222"):
<<>>=
get.dna(x, id=c("311","222"))
@
Note that we could also refer to sequences by their index in \texttt{get.sequences}:
<<>>=
get.sequences(x)
identical(get.dna(x, id=c("311","222")), get.dna(x, id=c(2,1)))
@

\item \texttt{get.ncontacts(x, from=NULL, to=NULL)}: returns the number of contacts in \texttt{x@contacts}; the
  optional arguments \texttt{from} and \texttt{to} can be used, in the case of dynamic networks, to
  specify the range of dates for which contacts should be kept.

\item \texttt{get.contacts(x, from=NULL, to=NULL)}: returns the contacts in \texttt{x@contacts}; the
  optional arguments \texttt{from} and \texttt{to} can be used, in the case of dynamic networks, to
  specify the range of dates for which contacts should be kept.
  Here, the object \texttt{x} contains no contact information, as the individuals of the samples
  retained were had no documented contacts:
<<>>=
get.ncontacts(ToyOutbreak)
get.individuals(ToyOutbreak@contacts)
get.individuals(x)
get.ncontacts(x)
@


\item \texttt{get.data(x, data=[\textit{name of data seeked}], where=NULL,
  drop=[\textit{TRUE/FALSE}], showSource=[\textit{TRUE/FALSE}])}: multi-purpose accessor
  seeking a data field with a given name in the entire dataset; \texttt{data} can be the name of a
  slot, or the name of a column in \texttt{x@individuals}, in the data.frames in \texttt{x@records}, or
  \texttt{x@context}, or in the \texttt{@dna@meta}. The optional argument \texttt{where} allows one to specify in which slot the
  information should be looked for. The argument \texttt{drop} states whether to return a vector
  (\texttt{TRUE}), or a one-column \texttt{data.frame} (\texttt{FALSE}), while \texttt{showSource}
  allows to put information in context (i.e., adding \texttt{individualID}, \texttt{date} and
  \texttt{source} where applicable.

\noindent For instance, we can retrieve temperature measurements using:
<<getData1>>=
get.data(x,"temperature")
get.data(x,"temperature", showSource=TRUE)
@
or the sex of the different individuals:
<<>>=
get.data(x, "Sex")
@
Several fields can be requested, so long as they are stored in the same slot; for instance:
<<>>=
get.data(x, c("Sex","Age","infector"))
@
The source (where matching fields were found) will be indicated if \texttt{showSource} is \texttt{TRUE}:
<<>>=
get.data(x, c("Sex","Age","infector"), showSource=TRUE)
@
This is especially useful when the same field appears in different slots, such as \texttt{date}:
<<>>=
get.data(x, "date")
@
actually corresponds to:
<<>>=
get.data(x, "date", showSource=TRUE)
@
as there are dates in both \texttt{@records} and \texttt{@dna}.
To retain only the latter, we use the argument \texttt{where}:
<<>>=
get.data(x, "date", where="records", showSource=TRUE)
@

A failed search will return NULL with a warning; for instance, we can try searching for ``sugarman'':
<<sugarman, warning=TRUE>>=
get.data(x, "sugarman")
@

\end{itemize}



% % % % % % % % % % % % % % % % %
\subsubsection{Accessors for \texttt{obkSequences} objects}
% % % % % % % % % % % % % % % % %

Accessors of \texttt{obkSequences} objects are basically a subset of what is available for \texttt{obkData}.
They work in the same way, and use the same arguments; they include:
\begin{itemize}
  \item \texttt{get.locus}
  \item \texttt{get.nlocus}
  \item \texttt{get.sequences}
  \item \texttt{get.nsequences}
  \item \texttt{get.dna}
  \item \texttt{get.individuals}
  \item \texttt{get.nindividuals}
  \item \texttt{get.dates}
  \item \texttt{get.ndates}
\end{itemize}




% % % % % % % % % % % % % % % % %
\subsubsection{Accessors for \texttt{obkContacts} objects}
% % % % % % % % % % % % % % % % %
Accessors of \texttt{obkContacts} objects are basically a subset of what is available for \texttt{obkData}.
They work in the same way, and use the same arguments; they include:
\begin{itemize}
  \item \texttt{get.nindividuals}
  \item \texttt{get.individuals}
  \item \texttt{get.ncontacts}
  \item \texttt{get.contacts}
  \item \texttt{get.dates}
  \item \texttt{get.ndates}
\end{itemize}

Another useful function is \texttt{as.matrix}, which converts the object into an adjacency matrix
(by default)), a matrix of incidence, or a matrix listing edges.
For instance, using a graph derived from the first 10 individuals in \texttt{ToyOutbreak}:
<<out.width==".8\\textwidth">>=
x <- subset(ToyOutbreak, individuals=1:10)
get.ncontacts(x)
plot(x@contacts, main="Contacts in x", label.cex=1.25, vertex.cex=2)
@

\noindent (note: see \texttt{?plot.network} to customize such graphics).
<<>>=
as.matrix(x@contacts)
as.matrix(x@contacts, "edgelist")
@

Lastly, for dynamic graphs, the function \texttt{as.data.frame} returns all the relevant information:
<<>>=
as.data.frame(x@contacts)
@




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Subsetting the data}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A lot of data handling lies in creating subsets of the data based on some given criteria.
The method \texttt{subset} for \texttt{obkData} objects allows for a range of manipulations.
The syntax is as follows:
<<eval=FALSE, tidy=FALSE>>=
subset(x, individuals=NULL, locus=NULL, sequences=NULL,
       date.from=NULL, date.to=NULL, date.format=NULL, ...)
@
See \texttt{?subset.obkData} for the details of these arguments.
The function works in a fairly intuitive way.
The arguments \texttt{individuals}, \texttt{locus} and \texttt{sequences} are
vectors of characters indicating items to be kept.
If integers or logicals are provided, these are assumed to match the output of
\texttt{get.[\textit{...}]}.
For instance, these two formulations are equivalent:
<<subset1>>=
data(ToyOutbreak)
x1 <- subset(ToyOutbreak, individuals=1:10)
x2 <- subset(ToyOutbreak, get.individuals(ToyOutbreak)[1:10])
identical(x1,x2)
@

Another, non-exclusive way of subsetting the data is using dates.
The arguments \texttt{date.from} and \texttt{date.to} are used for indicating the range of dates of
samples to be retained.
For instance, the range of data in the influenza H1N1 pandemic dataset \texttt{FluH1N1pdm2009} is:
<<>>=
data(FluH1N1pdm2009)
attach(FluH1N1pdm2009)

x <- new("obkData", individuals = individuals, dna = FluH1N1pdm2009$dna,
      dna.individualID = samples$individualID, dna.date = samples$date,
      trees = FluH1N1pdm2009$trees)

detach(FluH1N1pdm2009)

range(get.data(x, "date"))
@
We can retain data collected during the first month using:
<<subsetdate>>=
min.date <- min(get.dates(x))
min.date
min.date+31
x1 <- subset(x, date.to=min.date+31)
summary(x)
summary(x1)
@
Note that dates can also be provided as character strings in any sensible format, in which case
\texttt{subset} detects it automatically.
\\

Finally, note that several filters can be specified at the same time.
For instance, in the following we extract European data collected between the 1st June and the 31st August:
<<lastsubset>>=
temp <- get.data(x, "location", showSource=TRUE)
head(temp)
toKeep <- temp$individualID[temp$location=="Europe"]
x.summerEur <- subset(x, date.from="01/06/2009", date.to="31/08/2009",
                      indiv=toKeep)
summary(x.summerEur)
head(x.summerEur)
@





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Obtaining phylogenies from genetic sequences}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The package \textit{ape} implements a wide range of genetic distances (see \texttt{?dist.dna}) and
most usual algorithms for distance-based phylogenetic reconstruction.
In \textit{OutbreakTools}, the function \texttt{make.phylo} is a wrapper for these methods, allowing to
derive trees for a selection or all the genes present in an \texttt{obkData} object.
Trees can be stored in the \texttt{obkData} (\texttt{result='obkData'}) or returned as a
\texttt{multiPhylo} object (\texttt{result='multiPhylo'}).
We illustrate this procedure using \texttt{x.summerEur}, the data of pandemic H1N1 influenza collected in Europe during
the summer 2009 (see previous section):
<<>>=
x.summerEur@trees <- NULL
get.nsequences(x.summerEur)
@
\texttt{make.phylo} admits a range of arguments allowing to select which genes (\texttt{locus}),
model of evolution (\texttt{model}), and tree reconstruction method (\texttt{method}) should be used.
By default, a Neighbour-Joining tree based on Hamming distances (number of differing nucleotides) is
derived for every gene, and the resulting trees are plotted:
<<makephylo, fig.keep="all">>=
x2 <- make.phylo(x.summerEur)
summary(x2)
@

\texttt{x2} now contains a phylogenetic tree derived from the sequences in \texttt{x.summerEur}.
This one can be plotted simply, using:
<<out.width=".75\\textwidth">>=
library(ape)
plot(get.trees(x2)[[1]])
axisPhylo()
@
or alternatively:
<<>>=
plot(x2, "phylo")
@
\noindent Note that we could ask for a different model of evolution, for instance Kimura's 2
parameters distance:
<<tree1,fig.keep="last",out.width=".75\\textwidth">>=
x3 <- make.phylo(x.summerEur, locus=1, ask=FALSE, model="K80")
plot(get.trees(x3)[[1]])
axisPhylo()
@

Finally, note that \textit{OutbreakTools} also integrates functions to read annotated trees with Newick
(\texttt{read.annotated.tree}) or NEXUS (\texttt{read.annotated.nexus}) formats. This will be
particularly useful to process the outputs of Bayesian phylogenetic reconstruction software such as BEAST.
See \texttt{?read.annotated.nexus} for more information.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Simulating outbreak data}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\textit{OutbreakTools} provides some basic functionality for the simulation of outbreak data through the \texttt{simuEpi} function.
A basic SIR (susceptible-infectious-removed) model is assumed, and the result is returned as a list containing the SIR dynamics (\texttt{x\$dynamics}), an
\texttt{obkData} object (\texttt{x\$x}) and an optional \texttt{ggplot} graphic of the SIR dynamics (\texttt{x\$plot}).
\\

The arguments are as follows:
\begin{itemize}
\item \texttt{N}: the size of the population, which remains constant throughout. The simulation will start with one
infectious individual, \texttt{N-1} susceptibles and zero removed. Default is \texttt{N=1000}.
\item \texttt{D}: duration of the simulation, in days. Default is \texttt{D=10}.
\item \texttt{beta}: probability that a susceptible individual becomes infected by a given infectious individual on a given day.
Default is \texttt{beta=0.001}.
\item \texttt{nu}: rate of recovery, ie the probability that an infectious individual becomes removed on a given day.
Default is \texttt{nu=10}.
\item \texttt{L}: length of genetic sequences to be generated. Default is \texttt{L=1000}.
\item \texttt{mu}: rate of mutation per site per transmission event. Default is \texttt{mu=0.001}.
\item \texttt{plot}: logical indicating whether to create a plot of the SIR trajectory over
  time. Default is \texttt{plot=TRUE}. Plot will be a \texttt{ggplot} object stored as the
  \texttt{\$plot} slot of the returned list.
\item \texttt{makePhyloTree}: logical indicating whether to create a neighbor-joining tree from the simulated sequences.
Default is \texttt{makePhyloTree=FALSE}.
\end{itemize}

Let us look at an example in a very small population of size \texttt{N=50} and with the infectious rate \texttt{beta}
raised accordingly to generate a few transmission events:
<<fig.width=10, out.width="\\textwidth">>=
set.seed(1)
x <- simuEpi(N=50, D=20, beta=0.01,plot=TRUE,makePhylo=TRUE)
summary(x)
x$dynamics
summary(x$x)
@

We can see that \Sexpr{get.nindividuals(x$x)} individuals got infected over the time period of \texttt{D=20} days during which the outbreak
was simulated.
The actual transmission tree is stored as contact information:
<<>>=
plot(x$x, "contacts", main="Transmission tree")
@
The object also possesses a Neighbor-Joining tree based on the simulated sequence data:
<<>>=
plot(x$x, "phylo")
@






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Graphics for \texttt{obkData} objects}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Several plotting options are available for \texttt{obkData}, corresponding to different
sub-functions (see \texttt{?plot.obkData}).
The syntax to use is
\texttt{plot(x, y=[\textit{"timeline" or "geo" or "mst" or "phylo" or "contacts"}], ...)} where
\texttt{x} is an \textit{obkData} object, and \texttt{y} indicates the type of graphic to generate.
Further arguments can be passed via \texttt{...}.
The different types of graphics are:
\begin{itemize}
\item `\texttt{timeline}': plots the timeline of the outbreak; the timeline of every case
is plotted in a single window; uses \texttt{plotIndividualTimeline}.
\item `\texttt{geo}' plots the cases on a map. Needs geographical information. Uses \texttt{plotGeo}.
\item `\texttt{mst}': plots a minimal spanning tree of the genetic data. Uses \texttt{plotggMST}.
\item `\texttt{phylo}': plots a phylogenetic tree of the genetic data. Uses \texttt{plotggphy}.
\item `\texttt{contacts}': plots a phylogenetic tree of the genetic data. Uses the plot method for
  \textit{obkContacts}.

\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Plotting a timeline of samples}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This plotting option relies on the function \texttt{plotIndividualTimeline}; see
\texttt{?plotIndividualTimeline} for more information.
Let's plot the outbreak of equine influenza provided in \texttt{HorseFlu}:
<<>>=
data(HorseFlu)
summary(HorseFlu)
@
\noindent The default plot is a timeline showing all time-stamped data
<<plottime1>>=
plot(HorseFlu,'timeline')
@
A problem appears here: in this dataset, there is simply too much information to display on a single
graphic.
This can be improved by passing further arguments to \texttt{plotIndividualTimeline}:
<<>>=
args(plotIndividualTimeline)
@
For instance, we can choose to visualize only vaccination dates:
<<>>=
plot(HorseFlu,'timeline', what="Vac")
@
note that the argument \texttt{what} actually uses regular expressions to find matching fields in
the data, so that here \texttt{'Vac'} allows us to keep \texttt{FirstVac} and \texttt{LastVac}.
\\

Individuals can also be ordered and colored according to individual meta information.
For instance, to visualize collection dates of DNA and sort individuals per yards:
<<>>=
plotIndividualTimeline(HorseFlu, what="dna", colorBy="yardID", orderBy="yardID",plotNames=TRUE)
@
Note that only individuals for which requested information is present (here, DNA sequences) are plotted.
It is also possible to specify a subset of individuals using \texttt{selection}:
<<plotfirst20>>=
plot(HorseFlu,selection=1:20, colorBy="yardID", orderBy="yardID", size=5)
@



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Visualizing samples on a map}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This plotting option relies on the function \texttt{plotGeo}; see
%% \texttt{?plotGeo} for more information.
If geographical information is available, the function \texttt{plotGeo} can be used to visualize the
cases on a map (which is by default downloaded from googlemaps).
\texttt{plotGeo} is the function used by the generic \texttt{plot} of \texttt{obkData} when the
second argument is \texttt{'geo'}.
Geographical information can be provided as longitude/latitudes, or as strings specifying locations (which are converted to lon/lat using googlemaps).
Let us plot the toy outbreak already used before, and which already contains longitudes and latitudes.
<<>>=
data(ToyOutbreak)
head(ToyOutbreak@individuals)
@
We specify the columns holding these data with \texttt{location}, and we have to tell the function
that these are valid lon/lat with 'isLonLat' (which defaults to FALSE):
<<plotgeo1,dev='png',results='hide',message=FALSE>>=
plot(ToyOutbreak,'geo', location=c('lon','lat'), zoom=14)
@

\noindent We can also colour individuals by a certain characteristic using \texttt{colorBy} (here,
by sex), and even centre the map on a given individual using \texttt{center}:
<<dev='png',results='hide',message=FALSE>>=
plot(ToyOutbreak,'geo', location=c('lon','lat'), zoom=15,
     colorBy='Sex', center='11')
@



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Building minimum spanning trees from genetic sequences}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This plotting option relies on the function \texttt{plotggMST}; see
\texttt{?plotggMST} for more information.

It can be useful to plot a minimal spanning tree of the sequences, to quickly visualize the genetic diversity and the relation between sequences.
This can be achieved using \texttt{plotggMST}, or simply \texttt{plot} using \texttt{mst} for the
second argument:
<<plotmst1, cache=TRUE, fig.keep="last">>=
data(HorseFlu)
plot(HorseFlu,'mst')
@

\noindent this is a large tree, we can also look at the diversity within one individual,
e.g. individual 42:
<<>>=
plot(HorseFlu,'mst',individualID=42)
@





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Plotting phylogenetic trees}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This plotting option relies on the function \texttt{plotggphy}; see
%% \texttt{?plotggphy} for more information.

Phylogenies stored in \texttt{obkData} (slot \texttt{@trees}) can be plotted using \texttt{plotggphy}.
This function can be particularly useful as it allows for taking the collection dates into account
and for plotting a time tree (where branch length represent time, rather than quantity of
evolution).
We illustrate this function using data on pandemic influenza stored in \texttt{FluH1N1pdm2009}.
We first create an \texttt{obkData}:
<<>>=
data(FluH1N1pdm2009)
attach(FluH1N1pdm2009)

x <- new("obkData", individuals = individuals, dna = FluH1N1pdm2009$dna,
      dna.individualID = samples$individualID, dna.date = samples$date,
      trees = FluH1N1pdm2009$trees)

detach(FluH1N1pdm2009)

summary(x)
@
The phylogenie(s) contained in \texttt{x} can be extracted by:
<<>>=
get.trees(x)
tre <- get.trees(x)[[1]]
tre
@
and plotted using \textit{ape}'s standard \texttt{plot} function:
<<pdh1n1tree1, out.width="0.8\\textwidth">>=
plot(get.trees(x)[[1]], show.tip=FALSE)
@
However, we are loosing the temporal information about the samples:
<<pdh1n1,dev='png'>>=
plot(x, colorBy="location", orderBy="location")
@

The basic plot of \texttt{plotggphy} gives a tree quite similar to \textit{ape}'s:
<<fig.keep="last", out.width="0.8\\textwidth">>=
plotggphy(x)
@
However, \texttt{plotggphy} is also more flexible and powerful.
In particular, the argument \texttt{build.tip.attribute} allows to derive attributes for the tips
based on information on samples and individuals.
Here, for instance, we can use it to retrieve dates for each tip:
<<pdh1n1tree2, out.width="0.8\\textwidth">>=
p <- plotggphy(x, ladderize = TRUE,  branch.unit = "year")
@

\noindent Note that \texttt{p} is a graphical (\texttt{ggplot}) object, which can be re-used later to generate and
modify the plot.
Importantly, other attributes can also be used and represented by colors on the tips.
For instance, \texttt{x} contains information about the location of different individuals:
<<>>=
head(x@individuals)
@
Which can be exploited by:
<<pdh1n1tree3, out.width="\\textwidth">>=
p <- plotggphy(x, ladderize = TRUE, branch.unit = "year",
               tip.color = "location", tip.size = 3, tip.alpha = 0.75)
@




\end{document}
