\documentclass{svmult}
% Springer-Verlag style for multiauthors books
\usepackage{graphicx}
\usepackage{float}
\usepackage{url}
\usepackage{amsfonts}

\SweaveOpts{prefix.string=figs/seqinr,eps=F}

\begin{document}

\newcommand{\seqinr}{\texttt{seqin\bf{R}}}
\newcommand{\Seqinr}{\texttt{Seqin\bf{R}}}
\fvset{fontsize= \scriptsize}


%
% R output options
%
<<options, echo=FALSE, fig=FALSE, results=hide, eval=TRUE>>=
options(prompt=" ", continue=" ", width = 77)
CurFileName <- get("file", env = parent.frame(3))
r <- getOption("repos")
r["CRAN"] <- "http://cran.univ-lyon1.fr"
options(repos = r)
#library(seqinr)
library(seqinr, lib = "/Users/lobry/seqinr.Rcheck")#Users/lobry/seqinr.Rcheck")
library(xtable)
library(ade4)
library(ape)
library(MASS)
myplot <- function( res, ... )
{
  plot(res$li[ , 1], res$li[ , 2], ...)
  text(x = res$li[ , 1], y = res$li[ , 2], labels = 1:3, pos = ifelse(res$li[ , 2] < 0, 1, 3))
  perm <- c(3, 1, 2)
  lines( c(res$li[ , 1], res$li[perm, 1]), c(res$li[ , 2], res$li[perm, 2]))
}
@

% SeqinR adds
\mainmatter              % start of the contributions
%
\title{\Seqinr{} \Sexpr{packageDescription("seqinr")$Version}: a contributed package to the R project 
for statistical computing devoted to biological sequences retrieval and analysis}
%
\titlerunning{\Seqinr{} \Sexpr{packageDescription("seqinr")$Version}}  % abbreviated title (for running head)
%                                     also used for the TOC unless
%                                     \toctitle is used
%
\author{Delphine Charif\inst{1}, Jean R. Lobry\inst{1}, Leonor Palmeira\inst{1} }
%
\authorrunning{D. Charif \& J.R. Lobry \& L. Palmeira}   % abbreviated author list (for running head)
%
%%%% modified list of authors for the TOC (add the affiliations)
\tocauthor{Delphine Charif (Universit\'{e} Claude Bernard - Lyon I),
  Jean R. Lobry (Universit\'{e} Claude Bernard - Lyon I), Leonor
  Palmeira (Universit\'{e} Claude Bernard - Lyon I)}
%
\institute{Universit\'{e} Claude Bernard - Lyon I\\
Laboratoire de Biom\'{e}trie, Biologie \'{E}volutive\\
CNRS UMR 5558 - INRIA Helix project\\
43 Bd 11/11/1918\\F-69622 VILLEURBANNE CEDEX, FRANCE\\
\texttt{http://pbil.univ-lyon1.fr/members/lobry/}
}

\maketitle              % typeset the title of the contribution

\begin{abstract}
The \seqinr{} package for the R environment is a library of
utilities to retrieve and analyse biological sequences. It provides an
interface between i) the R language and environment for statistical
computing and graphics and ii) the ACNUC sequence retrieval system for
nucleotide and protein sequence databases such as GenBank, EMBL,
SWISS-PROT. ACNUC is very efficient in providing direct access to
subsequences of biological interest (\textit{e.g.} protein coding regions, tRNA
or rRNA coding regions) present in GenBank and in EMBL. Thanks to a
simple query language, it is then easy under R to select sequences of
interest and then use all the power of the R environment to analyze
them. The ACNUC databases can be locally installed but they are more
conveniently accessed through a web server to take advantage of
centralized daily updates. The aim of this paper is to provide a
handout on basic sequence analyses under \seqinr{} with a special focus
on multivariate methods.
\end{abstract}

\section{Introduction}

\subsection{About R and CRAN}

R \cite{R, RfromR} is a \emph{libre} language and environment for statistical computing and graphics 
which provides a wide variety of statistical and graphical techniques: linear and 
nonlinear modelling, statistical tests, time series analysis, classification, clustering, etc. 
Please consult the R project homepage at \texttt{http://www.R-project.org/} for 
further information. 

<<nmirrors, echo = FALSE, eval=T>>=
nmirrors <- function( url = paste("http://cran.r-project.org", "mirrors.html", sep="/") )
{
  readLines(url) -> tmp
  dtab <- grep("<table", tmp)
  ftab <- grep("</table", tmp)
  if( length(dtab) != length(ftab) )
  {
    stop("Unbalanced table while reading file")
  }
  if( ! all(ftab - dtab > 0) )
  {
    stop("overlapping table marks")
  }
  nm <- 0
  for( i in 1:length(dtab) )
  {
    nm <- nm + length(grep("http:", tmp[dtab[i]:ftab[i]]))
  }
  nm <- nm / 2
  return( list( nm = nm, nc = length(dtab) ) )
}
@

The Comprehensive R Archive Network, CRAN, is a network of servers 
around the world that store identical, up-to-date, versions of code and documentation 
for R. At compilation time of this document, there were
\Sexpr{ifelse(exists("nmirrors"), nmirrors()$nm, "???")} 
mirrors available 
from \Sexpr{ifelse(exists("nmirrors"), nmirrors()$nc, "???")} countries.
Please use the CRAN mirror nearest to you to minimize network load, they are
listed at \texttt{http://cran.r-project.org/mirrors.html}, and can be directly
selected with the function \texttt{chooseCRANmirror()}.

\subsection{About this document}

In the terminology of the R project \cite{R, RfromR}, this document is a package \emph{vignette}.
The examples given thereafter were run under \texttt{\Sexpr{R.version.string}}
on \Sexpr{date()} with Sweave \cite{Sweave}.
The last compiled version of this document is distributed along with the \seqinr{}
package in the \texttt{/doc} folder. Once \seqinr{} has been installed, the
full path to the package is given by the following R code :

\begin{Schunk}
\begin{Sinput}
 .find.package("seqinr")
\end{Sinput}
\begin{Soutput}
[1] "/Users/lobry/Library/R/library/seqinr"
\end{Soutput}
\end{Schunk}


\subsection{About sequin and \seqinr{}}

Sequin is the well known sofware used to submit sequences to GenBank, \seqinr{}
has definitively no connection with sequin. \seqinr{} is just a shortcut, with
no google hit, for "Sequences in R".

However, as a mnemotechnic tip, you may think about the \seqinr{} package
as the {\bf{R}}eciprocal function of sequin: with sequin you can submit sequences
to Genbank, with \seqinr{} you can {\bf{R}}etrieve sequences from Genbank. This is
a very good summary of a major functionality of the \seqinr{} package: to
provide an efficient access to sequence databases under R.

\subsection{About getting started}

You need a computer connected to the Internet. First, install R on your computer.
There are distributions for Linux, Mac and Windows users
on the CRAN (\texttt{http://cran.r-project.org}). Then, install the \texttt{ape}, 
\texttt{ade4} and \texttt{seqinr} packages. This can be done directly in an R console
with for instance the command \texttt{install.packages("seqinr")}. 
Last, load the \seqinr{} package with:

\begin{Schunk}
\begin{Sinput}
library(seqinr)
\end{Sinput}
\end{Schunk}

The command \texttt{lseqinr()} lists all what is defined in the package \seqinr{}:

<<lseqinr,eval=T>>=
lseqinr()[1:9]
@

We have printed here only the first 9 entries because they are too numerous.
To get help on a specific function, say \texttt{aaa()}, just prefix its name
with a question mark, as in \texttt{?aaa} and press enter.

\subsection{About running R in batch mode}

Although R is usually run in an interactive mode, some data pre-processing 
and analyses could be too long. You can run your R code in batch mode
in a shell with a command that typically looks like :

\begin{verbatim}
unix$ R CMD BATCH input.R results.out &
\end{verbatim}

where \texttt{input.R} is a text file with the R code you want to run and
\texttt{results.out} a text file to store the outputs. Note that in batch mode,
the graphical user interface is not active so that some graphical devices 
(\textit{e.g.} \texttt{x11}, \texttt{jpeg}, \texttt{png}) are not
available (see the R FAQ \cite{RFAQ} for further details).

It's worth noting that R uses the XDR representation of binary objects in binary saved files, 
and these are portable across all R platforms. The \texttt{save()} and \texttt{load()}
functions are very efficient (because of their binary nature) for saving and restoring any 
kind of R objects, in a platform independent way. To give a striking real example, at a given time
on a given platform, it was about $4$ minutes long to import a numeric table with 70000 lines and 64 columns
with the defaults settings of the \texttt{read.table()} function. Turning it into binary format,
it was then about $8$ \emph{seconds} to restore it with the \texttt{load()} function.
It is therefore advisable in the \texttt{input.R} batch file to save important data or
results (with something like \texttt{save(mybigdata, file = "mybigdata.RData")})
so as to be able to restore them later efficiently in the interactive mode (with something
like \texttt{load("mybigdata.RData")}).


\subsection{About the learning curve}

If you are used to work with a purely graphical user interface, you may feel frustrated in the
beginning of the learning process because apparently simple things are not so easily
obtained (\textit{ce n'est que le premier pas qui co{\^{u}te} !}).
In the long term, however, you are a winner for the following reasons.

<<numberofpackages, fig=F, echo=F, results = hide, eval=T>>=
numberofpackages <- nrow(available.packages())
@

\begin{description}
\item{\bf Wheel (the):} do not re-invent (there's a patent \cite{wheel} on it anyway).
At the compilation time of this document there were 
\Sexpr{ifelse(exists("numberofpackages"), numberofpackages, "???")}
contributed packages available. Even if you don't want to be spoon-feed 
\textit{{\`a} bouche ouverte}, 
it's not a bad
idea to look around there just to check what's going on in your own application field.
Specialists all around the world are there.

\item{\bf Hotline:} there is a very reactive discussion list to help you, just make sure to
read the posting guide there: \url{http://www.R-project.org/posting-guide.html}
before posting. Because of the high traffic on this list, we strongly suggest to answer \emph{yes} at the
question \emph{Would you like to receive list mail batched in a daily  digest?} when
subscribing at \url{https://stat.ethz.ch/mailman/listinfo/r-help}. Some \textit{bons mots}
from the list are archived in the R \texttt{fortunes} package.
\item{\bf Automation:} consider the 178 pages of figures in the additional data file 1
(\url{http://genomebiology.com/2002/3/10/research/0058/suppl/S1}) from \cite{lobrysueoka}. 
They were produced in part automatically (with a proprietary
software that is no more maintained) and manually, involving a lot of
tedious and repetitive manipulations (such as italicising species names by hand in subtitles).
In few words, a waste of time. The advantage of the R environment is that once you are
happy with the outputs (including graphical outputs) of an analysis for species x, it's very
easy to run the same analysis on n species. 

\item{\bf Reproducibility:} if you do not consider the reproducibility of scientific results
to be a serious problem in practice, then the paper by Jonathan Buckheit and David Donoho
\cite{repro} is a must read. Molecular data are available in public databases, this is
a necessary but not sufficient condition to allow for the reproducibility of results.
Publishing the R source code that was used in your analyses is a simple way
to greatly facilitate the reproduction of your results at the expense of no extra cost. 
At the expense of a little extra cost, you may consider to set up a RWeb server
so that even the laziest reviewer may reproduce your results just by clicking on
the "do it again" button in his web browser (\textit{i.e.} without installing any
software on his computer). For an example involving the \seqinr{} pacakage, 
follow this link \url{http://pbil.univ-lyon1.fr/members/lobry/repro/bioinfo04/}
to reproduce on-line the results from \cite{fifine}.

\item{\bf Fine tuning:} you have full control on everything, even the source code
for all functions is available. 
The following graph was specifically designed to illustrate
the first experimental evidence \cite{chargaff} that, on average, we have also [A]=[T] and [C]=[G] 
in single-stranded DNA. These data from Chargaff's lab give the base composition of the L (Ligth) 
strand for 7 bacterial chromosomes.

\setkeys{Gin}{width=0.5\textwidth}
<<chargaff,fig=T, results=hide,eval=T>>=
example(chargaff)
@
\setkeys{Gin}{width=0.8\textwidth}

This is a very specialised graph. The filled areas correspond to non-allowed values beause the sum 
of the four bases frequencies cannot exceed 100 \%. The white areas correspond to possible values 
(more exactly to the projection from $\mathbb{R}^4$ to the corresponding $\mathbb{R}^2$ 
planes of the region of allowed values). The lines correspond to the very small subset of allowed 
values for which we have in addition [A]=[T] and [C]=[G]. Points represent observed values in
the $7$ bacterial chromosomes. The whole graph is entirely defined by the code given in the
example of the \texttt{chargaff} dataset (\texttt{?chargaff} to see it).

Another example of highly specialised graph is given by the function \texttt{tablecode()} to
display a genetic code as in textbooks :

\setkeys{Gin}{width=0.6\textwidth}
<<tablecode1,fig=T, results=hide,eval=T>>=
tablecode(dia=F)
@
\setkeys{Gin}{width=0.8\textwidth}

It's very convenient in practice to have a genetic code at hand, and moreover here,
all genetic code variants are available :

\setkeys{Gin}{width=0.6\textwidth}
<<tablecode2,fig=T, results=hide,eval=T>>=
tablecode(numcode = 2, dia=F)
@
\setkeys{Gin}{width=0.8\textwidth}

As from \seqinr{} 1.0-4, it is possible to export the table of a genetic code into a \LaTeX~document,
for instance table \ref{code3.tex} and table \ref{code4.tex} were automatically generated with the following R code:

<<latextablecode,fig=F,eval=T>>=
tablecode(numcode = 3, urn.rna = s2c("TCAG"), latexfile = "code3.tex")
tablecode(numcode = 4, urn.rna = s2c("TCAG"), latexfile = "code4.tex")
@

\input{code3.tex}
\input{code4.tex}

The tables were then inserted in the \LaTeX~file with:
\begin{verbatim}
\input{code3.tex}
\input{code4.tex}
\end{verbatim}

\item{\bf Data as fast moving targets:} in research area, data are not always stable. compare the following
graph :
\setkeys{Gin}{width=0.5\textwidth}
<<dbg, fig = TRUE, eval=T>>=
dbg <- get.db.growth()
plot(x = dbg$date,
  y = log10(dbg$Nucl),
  las = 1,
  main = "The growth of DNA databases",
  xlab = "Year",
  ylab = "Log10 number of nucleotides")
@
\setkeys{Gin}{width=0.8\textwidth}

with figure 1 in \cite{lobrylncs}, data have been updated since then but the same R code
was used to produce the figure, ensuring an automatic update. For \LaTeX~users, it's
worth mentioning the fantastic tool contributed by Friedrich Leish \cite{Sweave}
called \texttt{Sweave()} that allows for the automatic insertion
of R outputs (including graphics) in a \LaTeX~document. In the same spirit, there
is a package called \texttt{xtable} to coerce R data into \LaTeX~tables, for instance
table \ref{df} here was produced this way, enforcing a complete coherence between
the R code example and the table.

\end{description}


%
% How to get sequences
%
\section{How to get sequence data}

\subsection{Importing raw sequence data from fasta files}

The fasta format is very simple and widely used for simple import of
biological sequences. It begins with a single-line description starting
with a character \texttt{>}, followed by lines of sequence data
of maximum 80 character each. Examples of files in fasta format
are distributed with the \seqinr{} package in the \texttt{sequences}
directory:

<<fastafiles,eval=T>>=
list.files(path = system.file("sequences", package = "seqinr"), pattern = ".fasta")
@

The function \texttt{read.fasta()} imports sequences from fasta files
into your workspace, for example:

<<readfasta, eval=T>>=
seqaa <- read.fasta(File = system.file("sequences/seqAA.fasta", package = "seqinr"), seqtype="AA")
seqaa
@

A more consequent example is given in the fasta file \texttt{ct.fasta} which
contains the complete genome of \textit{Chlamydia trachomatis} that was
used in \cite{oriloc}. You should be able to reproduce figure 1b from this
paper with the following code:

<<oriloc, fig=TRUE, results = hide, eval=T>>=
out <- oriloc(seq.fasta = system.file("sequences/ct.fasta", package ="seqinr"),
      g2.coord = system.file("sequences/ct.coord", package = "seqinr"),
     oldoriloc = TRUE)
plot(out$st, out$sk/1000, type="l", xlab = "Map position in Kb",
         ylab = "Cumulated composite skew in Kb", 
         main = expression(italic(Chlamydia~~trachomatis)~~complete~~genome), las = 1)
abline(h = 0, lty = 2)
text(400, -4, "Terminus")
text(850, 9, "Origin")
@

Note that the algorithm has been improved since then and that it's
more advisable to use the default option \texttt{oldoriloc = FALSE}
if you are interested in the prediction of origins and terminus of
replication from base composition biases (more on this at
\url{http://pbil.univ-lyon1.fr/software/oriloc.html}). See also \cite{smorfland}
for a recent review on this topic.

<<oriloc2, fig=TRUE, results = hide, eval=T>>=
out <- oriloc(seq.fasta = system.file("sequences/ct.fasta", package ="seqinr"),
      g2.coord = system.file("sequences/ct.coord", package = "seqinr"))
plot(out$st, out$sk/1000, type="l", xlab = "Map position in Kb",
         ylab = "Cumulated composite skew in Kb", 
         main = expression(italic(Chlamydia~~trachomatis)~~complete~~genome), las = 1)
mtext("New version")
abline(h = 0, lty = 2)
text(400, -4, "Terminus")
text(850, 9, "Origin")
@

\subsection{Importing aligned sequence data}

Aligned sequence data are very important in evolutionary studies,
in this representation all vertically aligned positions are supposed
to be homologous, that is sharing a common ancestor. This is a
mandatory starting point for comparative studies. 
There is a function in \seqinr{} called \texttt{read.alignment()} to 
read aligned sequences data from various formats  (\texttt{mase}, 
\texttt{clustal}, \texttt{phylip}, \texttt{fasta} or \texttt{msf})
produced by common external programs for multiple sequence alignment.

%The data returned by \texttt{read.alignment()} are of class
%alignment. Whereas sequences are stored as vector of character for the
%class \texttt{"SeqFastadna"}, \texttt{"SeqFastaAA"} and  \texttt{"SeqAcnucWeb"},
%they are stored as vector of strings for the class \texttt{"alignment"}.

Let's give an example. The gene coding for the mitochondrial cytochrome oxidase I 
is essential and therefore often used in phylogenetic studies because of its
ubiquitous nature. The following two sample tests
of aligned sequences of this gene (extracted from ParaFit \cite{parafit}), 
are distributed along with the \seqinr{} package:

<<readaln,eval=T>>=
louse <- read.alignment(system.file("sequences/louse.fasta", package = "seqinr"), format = "fasta")
louse$nam
gopher <- read.alignment(system.file("sequences/gopher.fasta", package = "seqinr"), format = "fasta")
gopher$nam
@

\begin{figure}[htbp]
   \begin{center}
   \begin{tabular}{cc}
      \includegraphics[width=0.5\textwidth]{figs/louse}&
      \includegraphics[width=0.5\textwidth]{figs/gopher}\\
   \end{tabular}
   \end{center}
   \caption{Louse (left) and gopher (right). 
   Images are from the wikipedia (\protect\url{http://www.wikipedia.org/}).
   The picture of the chewing louse \textit{Damalinia limbata} found on Angora goats
   was taken by Fiorella Carnevali (ENEA, Italy). The gopher drawing is from
   Gustav Mützel, Brehms Tierleben, Small Edition 1927.}
   \label{lousegopher}
\end{figure}

The aligned sequences are now imported in your R environment.
The $8$ genes of the first sample are from various species of louse (insects parasitics
on warm-blooded animals) and the $8$ genes of the second sample are from their corresponding
gopher hosts (a subset of rodents), see figure \ref{lousegopher} :

<<names,eval=T>>=
l.names <- readLines(system.file("sequences/louse.names", package = "seqinr"))
l.names
g.names <- readLines(system.file("sequences/gopher.names", package = "seqinr"))
g.names
@

\Seqinr{} has very few methods devoted to phylogenetic analyses but many are
available in the \texttt{ape} package. This allows for a very fine tuning
of the graphical outputs of the analyses thanks to the power of the R
facilities. For instance, a natural question
here would be to compare the topology of the tree of the hosts and their
parasites to see if we have congruence between host and parasite evolution.
In other words, we want to display two phylogenetic trees face to face. This
would be tedious with a program devoted to the display of a single phylogenetic
tree at time, involving a lot of manual copy/paste operations, hard to reproduce,
and then boring to maintain with data updates.

How does it looks under R? First, we need to \emph{infer} the tree topologies
from data. Let's try as an \emph{illustration} the famous neighbor-joining tree estimation 
of Saitou and Nei \cite{nj} with Jukes and Cantor's correction \cite{JC}
for multiple substitutions.

<<calculnjsurJC,eval=T>>=
library(ape)
louse.JC <- dist.dna(x =  lapply(louse$seq, s2c) , model = "JC69")
gopher.JC <- dist.dna(x =  lapply(gopher$seq, s2c) , model = "JC69")

l <- nj(louse.JC)
g <- nj(gopher.JC)
@ 

Now we have an estimation for \emph{illustrative} purposes of the tree topology for the parasite 
and their hosts. We want to plot the two trees face to face, and for this we must change 
R graphical parameters. The first thing to do is to save the current graphical parameter
settings so as to be able to restore them later:

<<savegraphicalparameters,eval=T>>=
op <- par(no.readonly = TRUE)
@

The meaning of the \texttt{no.readonly = TRUE} option here is that graphical
parameters are not all settable, we just want to save those we can change at will. Now,
we can play with graphics :

\setkeys{Gin}{width=\textwidth}
<<face2face, fig=TRUE, width=16, height=8,eval=T>>=
g$tip.label <- paste(1:8, g.names)
l$tip.label <- paste(1:8, l.names)

layout(matrix(data = 1:2, nrow = 1, ncol = 2), width=c(1.4, 1))
par(mar=c(2,1,2,1))
plot(g, adj = 0.8, cex = 1.4, use.edge.length=FALSE, 
  main = "gopher (host)", cex.main = 2)
plot(l,direction="l", use.edge.length=FALSE, cex = 1.4,
  main = "louse (parasite)", cex.main = 2)                                         
@
\setkeys{Gin}{width=0.8\textwidth}

We now restore the old graphical settings that were previously saved:

<<restoregraphicalparameters,eval=T>>=
par(op)
@

OK, this may look a little bit obscure if you are not fluent in programming, but please
try the following experiment. In your current working directory, that is in the
directory given by the \texttt{getwd()} command, create a text file called
\texttt{essai.r} with your favourite text editor, and copy/paste the previous R
commands, that is :

\tiny
\begin{verbatim}
louse <- read.alignment(system.file("sequences/louse.fasta", package = "seqinr"), format = "fasta")
gopher <- read.alignment(system.file("sequences/gopher.fasta", package = "seqinr"), format = "fasta")
l.names <- readLines("http://pbil.univ-lyon1.fr/software/SeqinR/Datasets/louse.names")
g.names <- readLines("http://pbil.univ-lyon1.fr/software/SeqinR/Datasets/gopher.names")
louse.JC <- dist.dna(x =  lapply(louse$seq, s2c), model = "JC69" )
gopher.JC <- dist.dna(x =  lapply(gopher$seq, s2c), model = "JC69" )
l <- nj(louse.JC)
g <- nj(gopher.JC)
g$tip.label <- paste(1:8, g.names)
l$tip.label <- paste(1:8, l.names)
layout(matrix(data = 1:2, nrow = 1, ncol = 2), width=c(1.4, 1))
par(mar=c(2,1,2,1))
plot(g, adj = 0.8, cex = 1.4, use.edge.length=FALSE, 
  main = "gopher (host)", cex.main = 2)
plot(l,direction="l", use.edge.length=FALSE, cex = 1.4,
  main = "louse (parasite)", cex.main = 2)        
\end{verbatim} 
\normalsize

Make sure that your text has been saved and then go back to R console to enter
the command :

\scriptsize
\begin{verbatim}
source("essai.r")
\end{verbatim}
\normalsize

This should reproduce the previous face-to-face phylogenetic trees in your R graphical device. 
Now, your boss is unhappy with working with the Jukes and Cantor's model \cite{JC}
and wants you to use the Kimura's 2-parameters distance \cite{K80} instead.
Go back to the text editor to change \texttt{model = "JC69"} by \texttt{model = "K80"},
save the file, and in the R console \texttt{source("essai.r")} again, you should
obtain the following graph :

\setkeys{Gin}{width=\textwidth}
<<K80, fig=TRUE, echo=F, results=hide, width=16, height=8,eval=T>>=
louse.JC <- dist.dna(x =  lapply(louse$seq, s2c))
gopher.JC <- dist.dna(x =  lapply(gopher$seq, s2c))
l <- nj(louse.JC)
g <- nj(gopher.JC)
g$tip.label <- paste(1:8, g.names)
l$tip.label <- paste(1:8, l.names)
layout(matrix(data = 1:2, nrow = 1, ncol = 2), width=c(1.4, 1))
par(mar=c(2,1,2,1))
plot(g, adj = 0.8, cex = 1.4, use.edge.length=FALSE, 
  main = "gopher (host)", cex.main = 2)
plot(l,direction="l", use.edge.length=FALSE, cex = 1.4,
  main = "louse (parasite)", cex.main = 2)        
@
\setkeys{Gin}{width=0.8\textwidth}

Nice congruence, isn't it? Now, something even worst, there was a error in the
aligned sequence set : the first base in the first sequence in the file
\texttt{louse.fasta} is not a C but a T. To locate the file on your system,
enter the following command:

<<whereislouse,fig=F,eval=T>>=
system.file("sequences/louse.fasta", package = "seqinr")
@

Open the \texttt{louse.fasta} file
in your text editor, fix the error, go back to the R console to
\texttt{source("essai.r")} again. That's all, your graph is now consistent with
the updated dataset.


\subsection{Complex queries in ACNUC databases}


As a rule of thumb, after compression one nucleotide needs one octet
of disk space storage (because you need also the annotations corresponding
to the sequences), so that most likely you won't have enough space on
your computer to work with a local copy of a complete DNA database.
The idea is to import under R only the subset of sequences you are
interested in. This is done in three steps:

\subsubsection{Choose a bank}

Select the database from which you want to extract sequences with the \texttt{choosebank()} function.
This function initiates a remote access to an ACNUC database. Called without arguments,
\texttt{choosebank()} returns the list of available databases:

<<choixbanque1, eval=T>>=
choosebank()
@

Biological sequence databases are fast moving targets, and for publication purposes it is
recommended to specify on which release you were working on when you made the job.
To get more informations about available databases on the server, just set 
the \texttt{infobank} parameter to \texttt{TRUE}. For
instance, here is the result for the three first databases on the default server 
at the compilation time (\today) of this document:

<<choixbanquemoreinfo, eval=T>>=
choosebank(infobank = TRUE)[1:3, ]
@

Note that there is a \texttt{status} column because a database could be unavailable
for a while during updates. If you try call \texttt{choosebank(bank = "bankname")} 
when the bank called \texttt{bankname} is off from server, you will get an explicit 
error message stating that this bank is temporarily unavailable, for instance:

%
% Ca on est obligé de le mettre en dur parce que ca génère une erreur !
%
\begin{Schunk}
\begin{Sinput}
 choosebank("off")
\end{Sinput}
\begin{Soutput}
Error in choosebank("off") : Database with name -->off<-- is currently off for maintenance, please try again later.
\end{Soutput}
\end{Schunk}

Some special purpose databases are not listed by default. These are \textit{tagged} databases
that are only listed if you provide an explicit \texttt{tagbank} argument to the \texttt{choosebank()}
function. Of special interest for teaching purposes is the \texttt{TP} tag, an acronym for
\textit{Travaux Pratiques} which means "practicals", and corresponds to \emph{frozen}
databases so that you can set up a practical whose results are stable from year to year. Currently
available frozen databases at the default server are:

<<frozen, eval=T>>=
choosebank(tagbank = "TP", infobank = TRUE)
@

Now, if you want to work with a given database, say GenBank, just call \texttt{choosebank()}
with \texttt{"genbank"} as its first argument and store the result in a variable
in the workspace, called for instance \texttt{mybank} in the example thereafter:

<<choixbanque2, eval=T>>=
mybank <- choosebank("genbank")
str(mybank)
@

The list returned by \texttt{choosebank()} here means that in the database
called \texttt{\Sexpr{ifelse(exists("mybank"), mybank$bankname, "???")}} at the compilation time
of this document there were 
\texttt{\Sexpr{ifelse(exists("mybank"), formatC(as.integer(mybank$totseqs), big.mark=","), "???")}}
sequences from
\texttt{\Sexpr{ifelse(exists("mybank"), formatC(as.integer(mybank$totspecs), big.mark=","), "???")}}
species and a total of
\texttt{\Sexpr{ifelse(exists("mybank"), formatC(as.integer(mybank$totkeys), big.mark=","), "???")}}
keywords. The status of the bank was
\texttt{\Sexpr{ifelse(exists("mybank"), mybank$status, "???")}}, 
and the release information was
\texttt{\Sexpr{ifelse(exists("mybank"), mybank$release, "???")}}.
For specialized databases, some relevant informations are also given in the
\texttt{details} component, for instance:

<<exdetails,fig=F, eval=T>>=
choosebank("taxobacgen")$details
@

The previous command has a side-effect that is worth mentioning. 
As from \seqinr~1.0-3, the result of the \texttt{choosebank()} function is automatically
stored in a global variable named \texttt{banknameSocket}, so that if no socket argument
is given to the \texttt{query()} function, the last opened database will be used by default
for your requests.
This is just a matter of convenience so that you don't have to explicitly specify the details of the
socket connection when working with the last opened database. You have, however,
full control of the process since \texttt{choosebank()} returns (invisibly) all the
required details. There is no trouble to open \emph{simultaneously} many databases.
You are just limited by the number of simultaneous connections your build of R is
allowed\footnote{
There is a very convenient function called \texttt{closeAllConnections()} in the R base package if
you want to close all open connections at once.}.

For advanced users who may wish to access to more than one database at time, a good advice
is to close them with the function \texttt{closebank()} as soon as possible so that the maximum
number of simultaneous connections is never reached. In the example below, we want to
display the number of taxa (\textit{i.e.} the number of nodes) in the species taxonomy associated
with each available database (including frozen databases). For this, we loop over available databases and 
close them as soon as the information has been retrieved.

<<taxaperbank,fig=T,eval=F, eval=T>>=
banks <- c(choosebank(), choosebank(tagbank="TP"))
ntaxa <- numeric(0)
for(i in banks){
  ntaxa[i] <- as.numeric(choosebank(i)$totspecs)
  closebank()
}
dotchart(log10(ntaxa[order(ntaxa)]), pch = 19,
main = "Number of taxa in available databases",
xlab = "Log10(number of taxa)")
@

\subsubsection{Make your query}

For this section, set up the default bank to GenBank, so that you don't have to provide the sockets details
for the \texttt{query()} function:

<<settogenbankbeforequery, eval=T>>=
choosebank("genbank")
@

Then, you have to say what you want, that is to compose a query
to select the subset of sequences you are interested in. The way to do this is
documented under \texttt{?query}, we just give here a simple example. 
In the query below, we want to select all the coding sequences 
(\texttt{t=cds}) from cat (\texttt{sp=felis catus}) that are not 
(\texttt{et no}) partial sequences (\texttt{k=partial}). 
We want the result to be stored in an object called \texttt{completeCatsCDS}.
<<query1,eval=T>>=
query("completeCatsCDS", "sp=felis catus et t=cds et no k=partial")
@

Now, there is in the workspace an object called \texttt{completeCatsCDS}, which
does not contain the sequences themselves but the \emph{sequence names} (and various relevant informations
such as the genetic code and the frame) that fit 
the query. They are stored in the \texttt{req} component of the object,
let's see the name of the first ten of them:

<<getNames,eval=T>>=
sapply(completeCatsCDS$req[1:10], getName)
@

The first sequence that fit our request is \texttt{\Sexpr{ifelse(exists("completeCatsCDS"), getName(completeCatsCDS$req[[1]]), "???")}},
the second one is \texttt{\Sexpr{ifelse(exists("completeCatsCDS"), getName(completeCatsCDS$req[[2]]), "???")}}, and so on. Note that
the sequence name may have an extension, this corresponds to \emph{subsequences},
a specificity of the ACNUC system that allows to handle easily a
subsequence with a biological meaning, typically a gene. The list of available subsequences
in a given database is given by the function \texttt{getType()}, for example the list
of available subsequences in GenBank is given in table \ref{genbank}.

%
% Besoin d'edition manuelle du fichier genbank.tex pour virer les caracteres spéciaux Latex, ici "_"
%
<<xtablegenbank, fig = FALSE, echo = FALSE,eval=FALSE>>=
choosebank("genbank") -> bank
tmp <- getType(bank$s)
tmp <- t(data.frame(tmp))
row.names(tmp)<-1:nrow(tmp)
names(tmp)<-NULL
colnames(tmp) <- c("Type","Description")
print(xtable(tmp, digits = rep(0,3), caption = paste("Available subsequences in", bank$bankname), label = "genbank"), 
file = "genbank.tex")
@
\input{genbank.tex}


The component \texttt{call} of \texttt{completeCatsCDS} keeps automatically a 
trace of the way you have selected the sequences: 

<<list1call,fig=F,eval=T>>=
completeCatsCDS$call
@

At this stage you can quit your R 
session saving the workspace image. The next time an R session is opened with the 
workspace image restored, there will be an object called \texttt{completeCatsCDS}, and 
looking into its \texttt{call} component will tell you that it contains the names 
of complete coding sequences from \textit{Felis catus}.

In practice, queries for sequences are rarely done in one step and are more likely
to be the result of an iterative, progressively refining, process. An important point
is that a list of sequences can be re-used. For instance, we can re-use \texttt{completeCatsCDS}
to get only the list of sequences that were published in 2004:

<<query2,eval=T>>=
query("ccc2004", "completeCatsCDS et y=2004")
length(ccc2004$req)
@

Hence, there were \Sexpr{ifelse(exists("ccc2004"), length(ccc2004$req), "???")} complete coding sequences published in 2004 for
\textit{Felis catus} in GenBank.

As from release 1.0-3 of the \seqinr{} package, there is new parameter \texttt{virtual}
which allows to disable the automatic retrieval of information for all list elements. This is interesting for list
with many elements, for instance :

<<queryvirtual,eval=T>>=
query("allcds", "t=cds", virtual = TRUE)
allcds$nelem
@

There are therefore \texttt{\Sexpr{ifelse(exists("allcds"), formatC(as.integer(allcds$nelem), big.mark=","), "???")}} coding
sequences in this version of GenBank\footnote{
which is stored in the \texttt{release} component of the object \texttt{banknameSocket}
and current value is today (\today): \texttt{banknameSocket\$release = 
\Sexpr{ifelse(exists("banknameSocket"), banknameSocket$release, "???")}}.
}. 
It would be long to get all the informations for the elements
of this list, so we have set the parameter \texttt{virtual} to \texttt{TRUE} and the \texttt{req}
component of the list has not been documented:

<<nodoc,eval=T>>=
allcds$req
@

However, the list can still be re-used\footnote{
of course, as long as the socket connection with the server has not been lost: virtual lists details are only
known by the server.}, 
for instance we may extract from this list all the sequences
from, say, \textit{Mycoplasma genitalium}:

<<chtouille,eval=T>>=
query("small", "allcds et sp=mycoplasma genitalium", virtual = TRUE)
small$nelem
@

There are then \texttt{\Sexpr{ifelse(exists("small"), formatC(as.integer(small$nelem), big.mark=","), "???")}} elements in
the list \texttt{small}, so that we can safely repeat the previous query without asking for a
virtual list:

<<chtouille2,eval=T>>=
query("small", "allcds et sp=mycoplasma genitalium")
sapply(small$req, getName)[1:10]
@

Here are some illustrations of using virtual list to answer simple questions about the
current GenBank release.

\begin{description}
\item[\textbf{Man.}] How many sequences are available for our species?
<<man, eval=T>>=
query("man","sp=homo sapiens",virtual=T)
man$nelem
@
There are \texttt{\Sexpr{ifelse(exists("man"), formatC(man$nelem, big.mark=","), "???")}} sequences from \textit{Homo sapiens}.

\item[\textbf{Sex.}] How many sequences are annotated with a keyword starting by sex?
<<sex, eval=T>>=
query("sex","k=sex@",virtual=T)
sex$nelem
@
There are \texttt{\Sexpr{ifelse(exists("sex"), formatC(sex$nelem, big.mark=","), "???")}} such sequences.

\item[\textbf{tRNA.}] How many complete tRNA sequences are available?
<<trnacplt, eval=T>>=
query("trna","t=trna et no k=partial",virtual=T)
trna$nelem
@
There are \texttt{\Sexpr{ifelse(exists("trna"), formatC(trna$nelem, big.mark=","), "???")}} complete tRNA sequences.

\item[\textbf{Nature vs. Science.}] In which journal were the more sequences published?
<<natvsscience, eval=T>>= 
query("nature","j=nature",virtual=T)
nature$nelem
query("science","j=science",virtual=T)
science$nelem
@
There are \texttt{\Sexpr{ifelse(exists("nature"), formatC(nature$nelem, big.mark=","), "???")}} sequences published
in \textit{Nature} and
\texttt{\Sexpr{ifelse(exists("science"), formatC(science$nelem, big.mark=","), "???")}} sequences published in
\textit{Science}, so that the winner is 
\textit{\Sexpr{ifelse(exists("nature"), ifelse(nature$nelem < science$nelem, "Science", "Nature"), "???")}}.

%
% \item[TriTryp] quand la ref de Science/309/404 sera dans genbank
%

\item[\textbf{Smith.}] How many sequences have Smith (last name) as author?
<<smith, eval=T>>=
query("smith","au=smith",virtual=T)
smith$nelem
@
There are \texttt{\Sexpr{ifelse(exists("smith"), formatC(smith$nelem, big.mark=","), "???")}} such sequences.

\item[\textbf{YK2.}] How many sequences were published after year 2000 (included)?
<<yk2, eval=T>>=
query("yk2","y>2000",virtual=T)
yk2$nelem
@
There are \texttt{\Sexpr{ifelse(exists("yk2"), formatC(yk2$nelem, big.mark=","), "???")}} sequences published after year 2000.

\item[\textbf{Organelle contest.}] Do we have more sequences from chloroplast genomes or from mitochondion genomes?
<<organelles, eval=T>>=
query("chloro","o=chloroplast",virtual=T)
chloro$nelem
query("mito","o=mitochondrion",virtual=T)
mito$nelem
@
There are \texttt{\Sexpr{ifelse(exists("chloro"), formatC(chloro$nelem, big.mark=","), "???")}} sequences from
chloroplast genomes and
\texttt{\Sexpr{ifelse(exists("mito"), formatC(mito$nelem, big.mark=","), "???")}} sequences from mitochondrion
genomes, so that the winner is 
\Sexpr{ifelse(exists("mito"), ifelse(chloro$nelem < mito$nelem, "mitochondrion", "chloroplast"), "???")}.


\end{description}

\subsubsection{Extract sequences of interest}

The sequence itself is obtained with the function \texttt{getSequence()}.
For example, the first 50 nucleotides of the first sequence of our request are:

<<getSequence, eval=T>>=
myseq <- getSequence(completeCatsCDS$req[[1]])
myseq[1:50]
@
They can also be coerced as string of character with the function \texttt{c2s()}:
<<SequenceAsString, eval=T>>=
c2s(myseq[1:50])
@
Note that what is done by \texttt{getSequence()} is much more complex
than a substring extraction because subsequences of biological interest are
not necessarily contiguous or even on the same DNA strand. Consider for
instance the following coding sequence from sequence \texttt{AE003734}:

\scriptsize
\begin{verbatim}
AE003734.PE35        Location/Qualifiers    (length=1833 bp)
     CDS             join(complement(162997..163210),
                     complement(162780..162919),complement(161238..162090),
                     146568..146732,146806..147266)
                     /gene="mod(mdg4)"
                     /locus_tag="CG32491"
                     /note="CG32491 gene product from transcript CG32491-RT;
                     trans-splicing"
\end{verbatim}
\normalsize

To get the coding sequence manually you would have join 5 different pieces 
from \texttt{AE003734} and some of them are in the complementary strand. 
With \texttt{getSequence()} you don't have to think about this. Just make a
query with the sequence name:

<<transplicing1, eval=T>>=
query("transspliced", "N=AE003734.PE35")
length(transspliced$req)
getName(transspliced$req[[1]])
@

Ok, now there is in your workspace an object called \texttt{transspliced} which \texttt{req}
component is of length one (because you have asked for just one sequence) and the name of the
single element of the req component is \Sexpr{ifelse(exists("transspliced"), getName(transspliced$req[[1]]), "???")} (because this
is the name of the sequence you wanted). Let see the first 50 base of this sequence:

<<transsplicing2, eval=T>>=
getSequence(transspliced$req[[1]])[1:50]
@

All the complex transsplicing operations have been done here. You can check that there is no
in-frame stop codons\footnote{
Stop codons are represented by the character \texttt{*} when translated into protein.} 
with the \texttt{getTrans()} function to translate this coding sequence into protein:

<<transsplicing3, eval=T>>=
getTrans(transspliced$req[[1]])[1:50]
table(getTrans(transspliced$req[[1]]))
@

In a more graphical way:

<<transp4,fig=T, eval=T>>=
aacount <- table(getTrans(transspliced$req[[1]]))
aacount <- aacount[order(aacount)]
names(aacount) <- aaa(names(aacount))
dotchart(aacount, pch = 19, xlab = "Stop and amino-acid counts",
main = "There is only one stop codon in AE003734.PE35")
abline(v=1, lty = 2)
@

Note that the relevant variant of the genetic code was automatically set up during the translation
of the sequence into protein. This is because the \texttt{transspliced\$req[[1]]} object belongs to the 
\texttt{SeqAcnucWeb} class:

<<transsplicing4, eval=T>>=
class(transspliced$req[[1]])
@

Therefore, when you are using the \texttt{getTrans()} function, you are automatically redirected
to the \texttt{getTrans.SeqAcnucWeb()} function which knows how to take into account the relevant frame
and genetic code for your coding sequence.


\section{How to deal with sequences}

\subsection{Sequence classes}


There are at present three classes of sequences, depending on the way they were obtained:

\begin{itemize}
      \item {\bfseries seqFasta} is the class for the sequences that were imported from a fasta file
      \item {\bfseries seqAcnucWeb} is the class for the sequences coming from an ACNUC database server
      \item {\bfseries seqFrag} is the class for the sequences that are fragments of other sequences
\end{itemize}

\subsection{Generic methods for sequences}

All sequence classes are sharing a common interface, so that there are very few method names we have to remember. 
In addition, all classes have their specific as.ClassName method that return an instance of the class,
and is.ClassName method to check whether an object belongs or not to the class. 
Available methods are: 
\\
\\
\begin{tabular}{|@{} c @{}|@{} c @{}|@{} c @{}|}
\hline
{\bfseries Methods} & {\bfseries Result} & {\bfseries Type of result} \\
\hline \hline
{\bfseries getFrag} & a sequence fragment & a sequence fragment \\
\hline
{\bfseries getSequence} & the sequence & vector of characters \\
\hline
{\bfseries getName} & the name of a sequence & string \\
\hline
{\bfseries getLength} & the length of a sequence & numeric vector \\
\hline
{\bfseries getTrans} & translation into amino-acids & vector of characters \\
\hline
{\bfseries getAnnot} & sequence annotations & vector of string \\
\hline
{\bfseries getLocation} & position of a Sequence on its parent sequence & list of numeric vector \\
\hline
\end{tabular}

\subsection{Internal representation of sequences}

The current mode of sequence storage is done with vectors of characters instead of strings.
This is very convenient for the user because all R tools to manipulate vectors are immediatly available. 
The price to pay is that this storage mode is extremly expensive in terms of memory.
They are two utilities called \texttt{s2c()} and \texttt{c2s()} that allows to convert strings into 
vector of characters, and \textit{vice versa}, respectively.

\subsubsection{Sequences as vectors of characters}

In the vectorial representation mode, all the very convenient R tools for indexing vectors
are at hand.
\begin{enumerate}
\item Vectors can be indexed by a vector of \emph{positive} integers saying which
elements are to be selected. As we have already seen, the first 50 elements of a sequence
are easily extracted thanks to the binary operator \texttt{from:to}, as in:

<<fromto,eval=T>>=
1:50
myseq[1:50]
@

The \texttt{seq()} function allows to build more complexe integer vectors. For instance
in coding sequences it is very common to focus on third codon positions where
selection is weak. Let's extract bases from third codon positions:

<<seqtcp,eval=T>>=
tcp <- seq(from = 3, to = length(myseq), by = 3)
tcp[1:10]
myseqtcp <- myseq[tcp]
myseqtcp[1:10]
@
 
 \item Vectors can also be indexed by a vector of \emph{negative} integers saying which
elements have to be removed. For instance, if we want to keep first and second codon positions,
the easiest way is to remove third codon positions:

<<seqfscp, eval=T>>=
-tcp[1:10]
myseqfscp <- myseq[-tcp]
myseqfscp[1:10]
@

\item Vectors are also indexable by a vector of \emph{logicals} whose \texttt{TRUE}
values say which elements to keep. Here is a different way to extract all third coding positions
from our sequence. First, we define a vector of three logicals with only the last one true:

<<ind, eval=T>>=
ind <- c(F, F, T)
ind
@

This vector seems too short for our purpose because our sequence is much more longer
with its \Sexpr{ifelse(exists("myseq"), length(myseq), "???")} bases. But under R vectors are automatically \emph{recycled}
when they are not long enough:

<<myseqtcp2, eval=T>>=
(1:30)[ind]
myseqtcp2 <- myseq[ind]
@

The result should be the same as previously:

<<identical, eval=T>>=
 identical(myseqtcp, myseqtcp2)
@

This recycling rule is extremely convenient in practice but may have surprising
effects if you assume (incorrectly) that there is a stringent dimension control for R vectors
as in linear algebra.

\end{enumerate}

Another advantage of working with vector of characters is that most R functions
are vectorized so that many things can be done without explicit looping. Let's
give some very simple examples:

<<vectorized1, eval=T>>=
tota <- sum(myseq == "a")
@

The total number of \texttt{a} in our sequence is \Sexpr{ifelse(exists("tota"), tota, "???")}. Let's compare
graphically the different base counts in our sequence :

<<vecto2, height = 3, fig=TRUE, eval=T>>=
basecount <- table(myseq)
myseqname <- getName(completeCatsCDS$req[[1]])
dotchart(basecount, xlim = c(0, max(basecount)), pch = 19,
  main = paste("Base count in",  myseqname))
@

<<vecto3, fig=TRUE, eval=T>>=
dinuclcount <- count(myseq, 2)
dotchart(dinuclcount[order(dinuclcount)], xlim = c(0, max(dinuclcount)), pch = 19,
  main = paste("Dinucleotide count in",  myseqname))
@

<<vecto4, height = 9, fig=TRUE, eval=T>>=
codonusage <- uco(myseq)
dotchart.uco(codonusage, main = paste("Codon usage in",  myseqname))
@


\subsubsection{Sequences as strings}

If you are interested in (fuzzy) pattern matching, then it is advisable to work with
sequence as strings to take advantage of \emph{regular expression} implemented
in R. The function \texttt{words.pos()} returns the positions of all occurrences
of a given regular expression. Let's suppose we want to know where are the trinucleotides
"cgt" in a sequence, that is the fragment CpGpT in the direct strand:

<<cgt, eval=T>>=
mystring <- c2s(myseq)
words.pos("cgt", mystring)
@

We can also look for the fragment CpGpTpY to illustrate fuzzy matching because
Y (IUPAC code for pyrimidine) stands C or T:

<<fuzzy, eval=T>>=
words.pos("cgt[ct]", mystring)
@

To look for all CpC dinucleotides separated by 3 or 4 bases:
<<fuzzy2, eval=T>>=
words.pos("cc.{3,4}cc", mystring)
@

Virtually any pattern is easily encoded with a regular expression. This is
especially useful at the protein level because many functions can be attributed 
to short linear motifs.

%
% Multivariate analyses
%
\section{Multivariate analyses}

\subsection{Correspondence analysis}

This is the most popular multivariate data analysis technique for amino-acid
and codon count tables, its application, however, is not without pitfalls \cite{misuse}.
Its primary goal is to transform a table of counts
into a graphical display, in which each gene (or protein) and each codon (or amino-acid)
is depicted as a point. Correspondence analysis (CA) may be defined as a special 
case of principal components analysis (PCA) with a different underlying metrics.
The interest of the metrics in CA, that is the way we measure the distance between
two individuals, is illustrated bellow with a very simple example (Table \ref{toyaa} inspired from \cite{CG}) with only
three proteins having only three amino-acids, so that we can represent exactly
on a map the consequences of the metric choice.

<<toyaa,eval=T>>=
data(toyaa)
toyaa
@

<<toyaa, fig = FALSE, echo = FALSE,eval=T>>=
print(xtable(toyaa, digits = rep(0,4), caption = "A very simple example of amino-acid counts in three proteins to be loaded with \\texttt{data(toyaa).}", label = "toyaa"), 
file = "toyaa.tex")
@
\input{toyaa.tex}

Let's first use the regular Euclidian metrics between two proteins $i$ and $i'$,
\begin{equation}
d^2(i,i') = \sum_{j=1}^{J}(n_{ij} - n_{i'j})^2
\label{euclidian}
\end{equation}
to visualize this small data set:

<<euclidian, fig = TRUE,eval=T>>=
library(ade4)
pco <- dudi.pco(dist(toyaa), scann = F, nf = 2)
myplot <- function( res, ... )
{
  plot(res$li[ , 1], res$li[ , 2], ...)
  text(x = res$li[ , 1], y = res$li[ , 2], labels = 1:3, pos = ifelse(res$li[ , 2] < 0, 1, 3))
  perm <- c(3, 1, 2)
  lines( c(res$li[ , 1], res$li[perm, 1]), c(res$li[ , 2], res$li[perm, 2]))
}
myplot(pco, main = "Euclidian distance", asp = 1, pch = 19, xlab = "", ylab = "", las = 1)
@

From this point of view, the first individual is far away from the two others. But
thinking about it, this is a rather trivial effect of protein size:

<<protsize,eval=T>>=
rowSums(toyaa)
@

With \Sexpr{ifelse(exists("toyaa"), rowSums(toyaa)[1], "???")} amino-acids, the first protein is two times bigger 
than the others so that when computing the Euclidian distance (\ref{euclidian}) its $n_{ij}$ entries
are on average bigger, sending it away from the others.
To get rid of this trivial effect, the first
obvious idea is to divide counts by protein lengths so as to work with 
\emph{protein profiles}. The corresponding distance is,

\begin{equation}
d^2(i,i') = \sum_{j=1}^{J}(\frac{n_{ij}}{n_{i\bullet}} - \frac{n_{i'j}}{n_{i'\bullet}})^2
\label{euclprof}
\end{equation}

where $n_{i\bullet}$ and $n_{i'\bullet}$ are the total number of amino-acids
in protein $i$ and $i'$, respectively.

<<profile, fig = TRUE,eval=T>>=
profile <- toyaa/rowSums(toyaa)
profile
dudi.pco(dist(profile), scann = F, nf = 2) -> pco1
myplot(pco1, main = "Euclidian distance on protein profiles", asp = 1, pch = 19, xlab = "", ylab = "",
  ylim = range(pco1$li[ , 2])*1.2)
@

The pattern is now completely different with the three protein equally spaced.
This is normal because in terms of relative amino-acid composition they are
all differing two-by-two by $5\%$ at the level of two amino-acids only. We
have clearly removed the trivial protein size effect, but this is still not completely
satisfactory. The proteins are differing by $5\%$ for all amino-acids but the situation 
is somewhat different for \texttt{Cys} because this amino-acid is very rare.
A difference of $5\%$ for a rare amino-acid has not the same significance than
a difference of $5\%$ for a common amino-acid such as \texttt{Ala} in our
example. To cope with this, CA make use of a variance-standardizing
technique to compensate for the larger variance in high frequencies and the 
smaller variance in low frequencies. This is achieved with the use of the 
\emph{chi-square distance} ($\chi^2$) which differs from the previous Euclidean distance 
on profiles (\ref{euclprof}) in that each square is weighted by the inverse of 
the frequency corresponding to each term,

\begin{equation}
d^2(i,i') = n_{\bullet\bullet}\sum_{j=1}^{J}\frac{1}{n_{{\bullet}j}}(\frac{n_{ij}}{n_{i\bullet}} - \frac{n_{i'j}}{n_{i'\bullet}})^2
\label{chi}
\end{equation}

where $n_{{\bullet}j}$ is the total number of amino-acid of kind $j$ and
$n_{\bullet\bullet}$ the total number of amino-acids. With this point
of view, the map is now like this:

<<afc, fig = TRUE,eval=T>>=
coa <- dudi.coa(toyaa, scann = FALSE, nf = 2)
myplot(coa, main = expression(paste(chi^2," distance")), 
  asp = 1, pch = 19, xlab = "", ylab = "")
@

The pattern is completely different with now protein number 3 which is far away from
the others because it is enriched in the rare amino-acid \texttt{Cys} as compared to
others.

The purpose of this small example was to demonstrates that the metric choice
is not without dramatic effects on the visualisation of data. Depending on your
objectives, you may agree or disagree with the $\chi^2$ metric choice, that's
not a problem, the important point is that you should be aware that there is
an underlying model there, \textit{chacun a son go{\^u}t} ou 
\textit{chacun {\`a} son go{\^u}t}, it's up to you.

Now, if you agree with the  $\chi^2$ metric choice, there's a nice
representation that may help you for the interpretation of results. 
This is a kind of "biplot" representation in which the lines and
columns of the dataset are simultaneously represented, in the
right way, that is as a graphical \textit{translation} of a 
mathematical theorem, but let's see how does it look like in practice: 

<<scatter, fig = TRUE,eval=T>>=
scatter(coa, clab.col = 0.8, clab.row = 0.8, posi = "none")
@

What is obvious is that the Cys content has a major effect on protein
variability here, no scoop. Please note how the information is well
summarised here: protein number 3 differs because it's enriched in
in Cys ; protein number 1 and 2 are almost the same but there is a
small trend protein number 1 to be enriched in Ala. As compared to
to table \ref{toyaa} this graph is of poor information here, so let's
try a more big-rooom-sized example (with $20$ columns so as to
illustrate the dimension reduction technique).

Data are from \cite{lobrygautier}, a sample of the proteome of
\textit{Escherichia coli}. According to the title of this paper,
the most important factor for the between-protein variability is
hydrophilic - hydrophobic gradient. Let's try to reproduce this
assertion :

<<lobrygautier,fig=T,eval=T>>=
download.file(url="ftp://pbil.univ-lyon1.fr/pub/datasets/NAR94/data.txt", destfile = "data.txt")

ec <- read.table(file = "data.txt", header = TRUE, 
    row.names = 1)
    
ec.coa <- dudi.coa(ec, scann = FALSE, nf = 1)
F1 <- ec.coa$li[,1]
hist(F1, proba = TRUE, xlab = "First factor for amino-acid variability",
col = grey(0.8), border = grey(0.5), las = 1, ylim = c(0,6),
        main="Protein distribution on first factor")
lines(density(F1, adjust = 0.5), lwd = 2)
@

There is clearly a bimodal distribution of proteins on the first factor. What are the
the amino-acid coordinates on this factor?

<<lobrygautier2,fig=T,height=5, eval=T>>=
aacoo <- ec.coa$co[ , 1]
names(aacoo) <- rownames(ec.coa$co)
aacoo <- sort( aacoo)
dotchart(aacoo, pch = 19, xlab = "Coordinate on first factor",
main = "Amino acid coordinates on first factor")
@

Aliphatic and aromatic amino-acids have positive values while charged amino-acids
have negative values\footnote{The physico-chemical classes for amino acids are
given in the component \texttt{AA.PROPERTY} of the \texttt{SEQINR.UTIL}
object.}. Let's try to compute the GRAVY score (\textit{i.e.} the Kyte and Doolittle 
hydropathic index\cite{KD}) of our proteins to compare this with their coordinates 
on the first factor. We need first the amino-acid \emph{relatives} frequencies in the
proteins, for this we divide the all the amino-acid counts by the total by row:

<<lobrygautier3,fig=F, eval=T>>=
ecfr <- ec/rowSums(ec)
ecfr[1:5, 1:5]
@

We need also the coefficients corresponding to the GRAVY score:

<<lobrygautier4,fig=F, eval=T>>=
gravy <- read.table(file ="ftp://pbil.univ-lyon1.fr/pub/datasets/NAR94/gravy.txt")
gravy[1:5, ]
coef <- gravy$V2
@

The coefficient are given in the alphabetical order of the three letter code for
the amino acids, that is in a different order than in the object \texttt{ecfr}:

<<lobrygautier5,fig=F, eval=T>>=
names(ecfr)
@

We then re-order the columns of the data set and check that everthing is OK:

<<lobrygautier6,fig=F, eval=T>>=
ecfr <- ecfr[ , order(names(ecfr))]
ecfr[1:5,1:5]
all(names(ecfr) == tolower(as.character(gravy$V1)))
@

Now, thanks to R build-in matrix multiplication, it's only one line to compute
the GRAVY score:

<<lobrygautier5,fig=T,eval=T>>=
as.matrix(ecfr) %*% coef -> gscores
plot(gscores,F1,xlab="GRAVY Score", ylab="F1 Score",las=1,main="The first factor is protein hydrophaty")
@

The proteins with high GRAVY scores are integral membrane proteins, and those
with low scores are cytoplasmic proteins. Now, suppose that we want to adjust
a mixture of two normal distributions to get an estimate of the proportion of
cytoplasmic and integral membrane proteins. We first have a look on the predefined
distributions (Table \ref{dpqr}), but there is apparently not an out of the box
solution.
<<dpqrtable, fig = FALSE, results = hide, echo = FALSE, eval=T>>=
noms <- ls("package:stats")
startwith <- function(nom, char)
{
  substr(nom, 1, 1) == char
}
(dnoms <- noms[sapply(noms, startwith, "d")])
(pnoms <- noms[sapply(noms, startwith, "p")])
(qnoms <- noms[sapply(noms, startwith, "q")])
(rnoms <- noms[sapply(noms, startwith, "r")])
suffix <- function(nom)
{
  substr(nom, 2, nchar(nom))
}
(dnoms <- sapply(dnoms, suffix, USE.NAMES = FALSE))
(pnoms <- sapply(pnoms, suffix, USE.NAMES = FALSE))
(qnoms <- sapply(qnoms, suffix, USE.NAMES = FALSE))
(rnoms <- sapply(rnoms, suffix, USE.NAMES = FALSE))
(dpqr <- dnoms[dnoms %in% pnoms & dnoms %in% qnoms & dnoms %in% rnoms])
(ddpqr <-  paste("d",dpqr, sep = ""))
(pdpqr <-  paste("p",dpqr, sep = ""))
(qdpqr <-  paste("q",dpqr, sep = ""))
(rdpqr <-  paste("r",dpqr, sep = ""))
data.frame(cbind(ddpqr,pdpqr, qdpqr, rdpqr), row.names = dpqr) -> data
names(data) <- c("d","p","q","r")
data
print(xtable(data, 
caption = "Density, distribution function, quantile function and random generation for the predefined distributions under R",
label = "dpqr"), 
file = "dpqrtable.tex")
@
\input{dpqrtable.tex}
We then define our own probability density function and then use \texttt{fitdistr} from package
\texttt{MASS} to get a maximum likelihood estimate of the parameters:

<<logfitdis,fig=T,eval=T>>=
dmixnor <- function(x, p, m1, sd1, m2, sd2){
  p*dnorm(x, m1, sd1) + (1 - p)*dnorm(x, m2, sd2)
}
library(MASS)
fitdistr(F1, dmixnor, list(p=0.88, m1=-0.04, sd1=0.076, m2=0.34, sd2=0.07))$estimate -> e
e
hist(F1, proba = TRUE, col = grey(0.8), 
main = "Ajustement with a mixture of two normal distributions",
xlab = "First factor for amino-acid variability", las = 1)
xx <- seq(from = min(F1), to = max(F1), length = 200)
lines(xx, dmixnor(xx,e[1],e[2],e[3],e[4],e[5]), lwd = 2)
@

\subsection{Synonymous and non-synonymous analyses}

Genetic codes are surjective applications from the set codons ($n=64$)
into the set of amino-acids ($n=20$) :

\setkeys{Gin}{width=\textwidth}
<<surjective, echo=F, fig=T, eval=T>>=
#
# Insert 2 (surjective genetic codes)
#
numcode <- 1 # To choose the genetic code

#
# General layout
#
symbols(x = rep(0,3), y = rep(0,3), 
        circles = c(1, 0.75, 0.45), 
        inches = FALSE,
        bg = c("pink", "white", "lightblue"), 
        xlim = c(-1, 1),
        ylim = c(-1, 1),
        bty = "n",
        asp = 1,
        main = paste("The surjective nature of genetic codes\nGenetic code number", numcode),
        xlab = "", ylab= "",
        xaxt ="n", yaxt = "n")
title( sub = "Adapted from insert 2 in Lobry & Chessel (2003) JAG 44:235")

words() -> codons
unlist(lapply(lapply(codons,s2c),translate, numcode = numcode)) -> aa
aaa(aa) -> aa3
#
# sort by alphabetical order of three-letter code of amino-acids
#

neworder <- order(aa3)
aa3 <- aa3[neworder]
aa <- aa[neworder]
codons <- codons[neworder]
#
# Text for codons
#
cangles <- seq(0, 2*pi, le = 65)[1:64]
text(x = sin(cangles)*0.9, y = cos(cangles)*0.9, labels = codons, cex = 0.65)
#
# Text for aa3
#
aangles <- seq(0, 2*pi, le = 22)[1:21]
text(x = sin(aangles)*0.35, y = cos(aangles)*0.35, 
     labels = unique(aa3), cex = 0.8)
#
# Text for aa
#
text(x = sin(aangles)*0.25, y = cos(aangles)*0.25, 
     labels = unique(aa), cex = 0.8)
#
# Draw lines
#
for( i in 1:64 )
{
  target <- aaa(translate(s2c(codons[i]), numcode = numcode))
  n <- which( unique(aa3) == target)
  lines(x = c(sin(cangles[i])*0.85, sin(aangles[n])*0.4), 
        y = c(cos(cangles[i])*0.85, cos(aangles[n])*0.4) )
}
@
\setkeys{Gin}{width=0.8\textwidth}

Two codons encoding the same amino-acid are said synonymous while
two codons encoding a different amino-acid are said non-synonymous.
The distinction between the synonymous and non-synonymous level are
very important in evolutionary studies because most of the selective
pressure is expected to work at the non-synonymous level, because the
amino-acids are the components of the proteins, and therefore more likely
to be subject to selection.

$K_s$ and $K_a$ are an estimation of the number of substitutions per synonymous site 
and per non-synonymous site, respectively, between two protein-coding genes \cite{kaks}.
The $\frac{K_{a}}{K_{s}}$ ratio is used as tool to evaluate selective pressure (see \cite{hurst}
for a nice back to basics). Let's give a simple illustration with three orthologous genes of the 
thioredoxin familiy from \textit{Homo sapiens}, \textit{Mus musculus},
and \textit{Rattus norvegicus} species: 

<<ortho, eval=T>>=
ortho <- read.alignment(system.file("sequences/ortho.fasta", package = "seqinr"), format="fasta")
kaks.ortho <- kaks(ortho)
kaks.ortho$ka/kaks.ortho$ks
@

The  $\frac{K_{a}}{K_{s}}$ ratios are less than $1$, suggesting a selective 
pressure on those proteins during evolution.

For transversal studies (\textit{i.e.} codon usage studies in a genome at the time it was sequenced)
there is little doubt that the strong requirement to distinguish between synonymous and an non-synonymous
variability was the source of many mistakes \cite{misuse}. We have just shown here with a scholarship
example that the metric choice is not neutral. If you consider that the $\chi^{2}$ metric is not too bad,
with respect to your objectives, and that you want to quantify the synonymous and an non-synonymous
variability, please consider reading this paper \cite{lobrychessel}, and follow this link
\url{http://pbil.univ-lyon1.fr/members/lobry/repro/jag03/} for on-line reproducibility.

Let's now use the toy example given in table \ref{toycodon} to illustrate how to study synonymous
and non-synonymous codon usage.

<<toycodon,eval=T>>=
data(toycodon)
toycodon
@

<<xtabletoycodon, fig = FALSE, echo = FALSE,eval=T>>=
print(xtable(toycodon, digits = rep(0,11), caption = "A very simple example of codon counts in three coding sequences to be loaded with \\texttt{data(toycodon).}", label = "toycodon"), 
file = "toycodon.tex")
@
\input{toycodon.tex}

Let's first have a look to global codon usage, we do not take into account the structure
of the genetic code:

<<globaltoycodon, fig=T,eval=T>>=
global <- dudi.coa(toycodon, scann = FALSE, nf= 2)
myplot(global, asp = 1, pch = 19, xlab = "", ylab = "", main = "Global codon usage")
@

From a global codon usage point of view, coding sequence number 3 is away.
To take into account the genetic code structure, we need to know for which amino-acid the codons are coding.
The codons are given by the names of the columns of the object \texttt{toycodon}:

<<toy1,fig=F,eval=T>>=
names(toycodon)
@


Put all codon names into a single string:

<<toy2,fig=F,eval=T>>=
c2s(names(toycodon))
@

Transform this string as a vector of characters:

<<toy3,fig=F,eval=T>>=
s2c(c2s(names(toycodon)))
@

Translate this into amino-acids using the default genetic code:

<<toy4,fig=F,eval=T>>=
translate(s2c(c2s(names(toycodon))))
@

Use the three letter code for amino-acid instead:

<<toy5,fig=F,eval=T>>=
aaa(translate(s2c(c2s(names(toycodon)))))
@

Make this a factor:

<<toy5,fig=F,eval=T>>=
facaa <- factor(aaa(translate(s2c(c2s(names(toycodon))))))
facaa
@

The non synonymous codon usage analysis is the between amino-acid analysis:

<<toy6, fig=T,eval=T>>=
nonsynonymous <- t(between(dudi = t(global), fac = facaa, scann = FALSE, nf = 2))
myplot(nonsynonymous, asp = 1, pch = 19, xlab = "", ylab = "", main = "Non synonymous codon usage")
@

This is reminiscent of something, let's have a look at amino-acid counts:

<<toy7,eval=T>>=
by(t(toycodon), facaa, colSums)
@

This is exactly the same data set that we used previously (table \ref{toyaa}) at the amino-acid
level. The non synonymous codon usage analysis is exactly the same as the amino-acid analysis.
Coding sequence number 3 is far away because it codes for many Cys, a rare amino-acid. Note
that at the global codon usage level, this is also the major visible structure. To get rid of this
amino-acid effect, we use the synonymous codon usage analysis, that is the within amino-acid
analysis:

<<toy8, fig=T,eval=T>>=
synonymous <- t(within(dudi = t(global), fac = facaa, scann = FALSE, nf = 2))
myplot(synonymous, asp = 1, pch = 19, xlab = "", ylab = "", main = "Synonymous codon usage")
@

Now, coding sequence number 2 is away. When the amino-acid effect is removed, the pattern is then
completely different. To interpret the result we look at the codon coordinates on the first factor of
synonymous codon usage:

<<toy9,fig=T,eval=T>>=
synonymous$co[ , 1, drop = FALSE] -> tmp
tmp <- tmp[order(tmp$Axis1), , drop = FALSE]
colcod <- sapply(rownames(tmp), function(x) ifelse(substr(x,3,3) == "c" || substr(x,3,3) == "g", "blue", "red"))
pchcod <- ifelse(colcod=="red",1,19)
dotchart(tmp$Axis1, labels = toupper(rownames(tmp)),
color = colcod, pch = pchcod,
main = "Codon coordinates on first factor\nfor synonymous codon usage")
legend("topleft", inset = 0.02, legend = c("GC ending codons", "AT ending codons"),
text.col = c("blue", "red"), pch = c(19,1), col = c("blue","red"), bg = "white")
@

At the synonymous level, coding sequence number 2 is different because it is enriched in GC-ending codons
as compared to the two others. Note that this is hard to see at the global codon usage level because of the
strong amino-acid effect.

\begin{figure}[htbp]
   \begin{center}
      \includegraphics{figs/lobgau5}
   \end{center}
   \caption{Screenshot of figure 5 from \cite{lobrygautier}. Each point represents
   a protein. This was to show the correlation between the codon adaptation index (CAI Score)
   with the second factor of correspondence analysis at the amino-acid level (F2 Score). Highly
   expressed genes have a high CAI value.
   }
   \label{lobgau5}
\end{figure}


To illustrate the interest of synonymous codon usage analyses, let's use now a more realistic example.
In \cite{lobrygautier} there was an assertion stating that selection for
translation optimisation in \textit{Escherichia coli} was also visible at the amino-acid level.
The argument was in figure 5 of the paper (\textit{cf} fig \ref{lobgau5}), that can be reproduced\footnote{
  the code to reproduce all figures from \cite{lobrygautier} is available at
  \url{http://pbil.univ-lyon1.fr/members/lobry/repro/nar94/}.
} with the following R code:

%
% C'est long celui-la, vaudrait mieux refaire les calculs a partir de ec999, ca illustrerait
% comment on calcule le CAI
%
<<reprolobgau,fig=T, eval=T>>=
ec <- read.table(
     file = "ftp://pbil.univ-lyon1.fr/pub/datasets/NAR94/data.txt",
     header = TRUE,
     row.names = 1)
ec.coa <- dudi.coa(ec, scann = FALSE, nf=3)
tmp <- read.table(
     file = "ftp://pbil.univ-lyon1.fr/pub/datasets/NAR94/ecoli999.cai")
cai <- exp(tmp$V2)
plot(cai, ec.coa$li[,2], pch=20, xlab="CAI Score", ylab="F2 Score",
     main="Fig 5 from Lobry & Gautier (1994) NAR 22:3174")
@


So, there was a correlation between the CAI (Codon Adaptation Index \cite{CAI}) and
the second factor for amino-acid composition variability. However, this
is not completely convincing because the CAI is not completely independent
of the amino-acid composition of the protein. Let's use within amino-acid
correspondence analysis to remove the amino-acid effect. Here is a commented
step-by-step analysis:

<<vlg1, fig = F, eval=T>>=
data(ec999)
class(ec999)
names(ec999)[1:10]
ec999[[1]][1:50]
@

This is to load the data from \cite{lobrygautier} which is available as \texttt{ec999}
in the \seqinr{} package. The letters \texttt{ec} are for the bacterium 
\textit{Escherichia coli} and the number \texttt{999} means that there were
$999$ coding sequences available from this species at that time. The class of the
object \texttt{ec999} is a list, which names are the coding sequence names, for instance
the first coding sequence name is 
\texttt{\Sexpr{ifelse(exists("ec999"), names(ec999)[1], "???")}}.
Each element of the list is a vector of character, we have listed just above the 50 first
character of the first coding sequence of the list with \texttt{ec999[[1]][1:50]}, we
can see that there is a start codon (ATG) at the beginning of the first coding sequence.

<<vlg2, fig = F, eval=T>>=
ec999.uco <- lapply(ec999, uco) # compute codon usage for all CDS
class(ec999.uco)
class(ec999.uco[[1]])
ec999.uco[[1]]
@

This is to compute the codon usage, that is how many times each codon is used
in each coding sequence. Because \texttt{ec999} is a list, we use the function
\texttt{lapply()} to apply the same function, \texttt{uco()}, to all the
elements of the list and we store the result in the object \texttt{ec999.uco}.
The object \texttt{ec999.uco} is a list too, and all its elements belong to
the class table.

<<vlg3, fig = F, eval=T>>=
df <- as.data.frame(lapply(ec999.uco, as.vector)) # put it in a dataframe
dim(df)
df[1:5,1:5]
@

This is to put the codon usage into a data.frame. Note that the codons are in row
and the coding sequences are in columns. This is more convenient for the following
because groups for within and between analyses are usually handled by row.

<<vlg4, fig = F, eval=T>>=
row.names(df) <- names(ec999.uco[[1]]) # add codon names
df[1:5,1:5]
@

This is to keep a trace of codon names, just in case we would like to re-order
the dataframe \texttt{df}. This is important because we can now play with
the data at will without loosing any critical information. 

<<vlg5, fig = F, eval=T>>=
ec999.coa <- dudi.coa(df = df, scannf = FALSE) # run global correspondence analysis
ec999.coa
@

This is to run global correspondence analysis of codon usage.
We have set the \texttt{scannf} parameter to \texttt{FALSE}
because otherwise the eigenvalue bar plot is displayed for the
user to select manually the number of axes to be kept.


<<vlg6, fig = F, eval=T>>=
facaa <- as.factor(aaa(translate(s2c(c2s(rownames(df)))))) # define a factor for amino-acids
facaa
@

This is to define a factor for amino-acids. The function \texttt{translate()} use by
default the standard genetic code and this is OK for \textit{E. coli}.

<<vlg7, fig = F, eval=T>>=
ec999.syn <- within(dudi = ec999.coa, fac = facaa, scannf = FALSE) # run synonymous codon usage analysis
ec999.syn
@

This is to run the synonymous codon usage analysis. The value of the \texttt{ratio} component of
the object \texttt{ec999.syn} shows that most of the variability is at the synonymous
level, a common situation in codon usage studies.

<<vlg8, fig = F, eval=T>>=
ec999.btw <- between(dudi = ec999.coa, fac  = facaa, scannf = FALSE) # run non-sysnonymous codon usage analysis <=> amino-acid usage
ec999.btw
@

This is to run the non-sysnonymous codon usage analysis, or amino-acid usage analysis.

<< vlg9,fig=T,eval=T>>=
x <- ec999.syn$co[,1] 
y <- ec999.btw$co[,2]
kxy <- kde2d(x,y, n = 100)
nlevels <- 25
breaks <- seq(from = min(kxy$z), to = max(kxy$z), length = nlevels + 1)
col <- cm.colors(nlevels)
image(kxy, breaks = breaks, col = col, xlab = "First synonymous factor",
ylab = "Second non-synonymous factor", xlim = c(-0.5, 0.5),
ylim  = c(-0.3, 0.3), las = 1,
main = "The second factor for amino-acid variability is\ncorrelated with gene expressivity")
contour(kxy, add = TRUE, nlevels = nlevels, drawlabels = FALSE)
box()
abline(c(0,1), lty=2)
abline(lm(y~x))
legend("topleft", lty = c(2,1), legend = c("y = x", "y = lm(y~x)"), inset = 0.01, bg = "white")
@

This is to plot the whole thing. We have extracted the coding sequences coordinates
on the first synonymous factor and the second non-synonymous factor within
\texttt{x} and \texttt{y}, respectively. Because we have many points, we
use the two-dimensional kernel density estimation provided by the function
\texttt{kde2d()} from package \texttt{MASS}.


%
% Akashi Gojobori
%

<<xaacost, fig = FALSE, echo = FALSE,eval=T>>=
data(aacost)
print(xtable(aacost, digits = rep(0,8), caption = "Aerobic cost of amino-acids in \\textit{Escherichia coli} and G+C classes to be loaded with \\texttt{data(aacost).}", 
label = "aacost"), 
file = "aacost.tex")
@
\input{aacost.tex}

\clearpage
\SweaveInput{prochlo.rnw}

\clearpage
\section{FAQ: Frequently Asked Question}

\subsection{How do I compute a score on my sequences?}

In the example below we want to compute the G+C content in third codon
positions for complete ribosomal CDS from \textit{Escherichia coli}:

<<gc3, fig=F,eval=T>>=
choosebank("emblTP")
query("ecribo","sp=escherichia coli ET t=cds ET k=ribosom@ ET NO k=partial")
sapply(sapply(ecribo$req, getSequence), GC3)
@

At the amino-acid level, we may get an estimate of the isoelectric point of
the proteins this way:

<<computePI,fig=F,eval=T>>=
sapply( sapply(sapply(ecribo$req, getSequence), getTrans), computePI)
@

Note that some pre-defined vectors to compute linear forms on sequences are
available in the \texttt{EXP} data.

As a matter of convenience, you may encapsulate the computation of your favorite score 
within a function this way :

<<encapsulate,fig=F,eval=T>>=
GC3m <- function(list, ind = 1:list$nelem) sapply(sapply(list$req[ind], getSequence), GC3)
GC3m(ecribo)
GC3m(ecribo, 1:10)
@

\section{Releases notes}

\subsection{release 1.0-4}

\begin{itemize}
\item The scaling factor $n_{\bullet\bullet}$ was missing in equation \ref{chi}.
\item The files \texttt{louse.fasta}, \texttt{louse.names}, \texttt{gopher.fasta}, \texttt{gopher.names}
and \texttt{ortho.fasta} that were used for examples in the previous version of this document are
no more downloaded from the internet since they are now distributed in the \texttt{sequences/} folder
of the package.
\item An example of synonymous and non synonymous codon usage analysis was added to the
vignette along with two toy data sets (\texttt{toyaa} and \texttt{toycodon}).
\item A FAQ section was added to the vignette.
\item A bug in \texttt{getAnnot()} when the number of lines was zero is now fixed.
\item There is now a new argument, \texttt{latexfile}, in \texttt{tablecode()} to export genetic codes
tables in a \LaTeX~document, for instance table \ref{code3.tex} and table \ref{code4.tex} here.
\item Function \texttt{splitseq()} has been entirely rewritten to improve speed.
\end{itemize}

\subsection{release 1.0-3}

\begin{itemize}
\item The new package maintainer is Dr. Simon Penel, PhD, who has now a fixed position in the
laboratory that issued \seqinr~(\texttt{penel@biomserv.univ-lyon1.fr}). Delphine Charif was
successful too to get a fixed position in the same lab, with now a different research task (but who knows?).
Thanks to the close vicinity of our pioneering maintainers the transition was sweet. The DESCRIPTION
file of the \seqinr{} package has been updated to take this into account.

\item The reference paper for the package is now \textit{in press}. We do not have the full
reference for now, you may use \texttt{citation("seqinr")} to check if it is complete now:
<<cite, fig=F, eval=T>>=
citation("seqinr")
@

\item There was a bug when sending a \texttt{gfrag} request to the server for long (Mb range) 
sequences. The length argument was converted to scientific notations that are not understand by the
server. This is now corrected and should work up the the Gb scale.

\item The \texttt{query()} function has been improved by de-looping list element info request,
there are now download at once which is much more efficient. For example, a query from a
researcher-home ADSL connection with a list with about 1000 elements was 60 seconds and
is now only 4 seconds (\textit{i.e.} 15 times faster now).

\item A new parameter \texttt{virtual} has been added to \texttt{query()} 
so that long lists can stay on the server without trying to download
them automatically. A query like \texttt{query(s\$socket,"allcds","t=cds", virtual = TRUE)} is 
now possible.

\item Relevant genetic codes and frames are now automatically propagated.

\item \Seqinr{}~sends now its name and version number to the server.

% Pas arrvié à le repoduire
%
%\item A bug as been reported for intensive \texttt{kaks()} calls.

\item Strict control on ambiguous DNA base alphabet has been relaxed.

\item Default value for parameter \texttt{invisible} of function \texttt{query()} is now \texttt{TRUE}.

\end{itemize}

\section{Acknowledgments}

Please enter \texttt{contibutors()} in your R console.

%
% ---- Bibliography ----
%
\begin{thebibliography}{99}
%

\bibitem{fifine}
Charif, D., Thioulouse, J., Lobry, J.R., Perri{\`e}re, G.:
Online synonymous sodon usage analyses with the ade4 and seqinR packages.
Bioinformatics {\bf 21} (2005) 545--547.
\url{http://pbil.univ-lyon1.fr/members/lobry/repro/bioinfo04/}.

\bibitem{repro}
Buckheit, J., Donoho, D.L.:
Wavelab and reproducible research.
(1995) \textit{In} A. Antoniadis (ed.), \textit{Wavelets  and Statistics}, Springer-Verlag, Berlin, New  York.

\bibitem{CG}
Gautier, C:
Analyses statistiques et {\'e}volution des s{\'e}quences d'acides nucl{\'e}iques.
PhD thesis (1987), Universit{\'e} Claude Bernard - Lyon I.

\bibitem{RFAQ}
Hornik, K.: 
The R FAQ.
ISBN 3-900051-08-9 (2005) \url{http://CRAN.R-project.org/doc/FAQ/}.

\bibitem{KD}
Kyte, J., Doolittle, R.F.: 
A simple method for displaying the hydropathic character of a protein. 
J. Mol. Biol. {\bf 157} (1982) 105--132.

\bibitem{Sweave}
Leisch, F.:
Sweave: Dynamic generation of statistical reports using literate data analysis.
Compstat 2002 --- Proceedings in Computational Statistics (2002) 575--580
ISBN 3-7908-1517-9.

\bibitem{oriloc}
Frank, A.C., Lobry, J.R.:
Oriloc: prediction of replication boundaries in unannotated bacterial chromosomes.
Bioinformatics {\bf 16} (2000) 560--561.

\bibitem{hurst}
Hurst, L.D.: 
The Ka/Ks ratio: diagnosing the form of sequence evolution. 
Trends Genet. {\bf 18} (2002) 486--487.

\bibitem{R}
Ihaka, R., Gentleman, R.:
R: A Language for Data Analysis and Graphics.
J. Comp. Graph. Stat. {\bf 3} (1996) 299--314

\bibitem{JC}
Jukes, T.H., Cantor, C.R.:
Evolution of protein molecules.
(1969) pp. 21--132. \textit{In} H.N. Munro (ed.), \textit{Mammalian Protein Metabolism},
Academic Press, New York.

\bibitem{wheel}
Keogh, J.:
Circular transportation facilitation device.
(2001) Australian Patent Office application number \textit{AU 2001100012 A4}.
\url{www.ipmenu.com/archive/AUI_2001100012.pdf}. 

\bibitem{K80}
Kimura, M.:
A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. 
J. Mol. Evol. {\bf 16} (1980) 111--120.

\bibitem{parafit}
Legendre, P., Desdevises, Y., Bazin, E.:
A statistical test for host-parasite coevolution.
Syst. Biol. {\bf 51} (2002) 217--234.

\bibitem{kaks}
Li, W.-H.:
Unbiased estimation of the rates of synonymous and nonsynonymous substitution.
J. Mol. Evol. {\bf 36} (1993) 96--9.

\bibitem{lobrylncs}
Lobry, J.R.:
Life history traits and genome structure: aerobiosis and G+C content in bacteria. 
Lecture Notes in Computer Sciences {\bf 3039} (2004) 679--686.
\url{http://pbil.univ-lyon1.fr/members/lobry/repro/lncs04/}.

\bibitem{lobrychessel}
Lobry, J.R., Chessel, D.:
Internal correspondence analysis of codon and amino-acid usage in thermophilic bacteria. 
J. Appl. Genet. {\bf 44} (2003) 235--261.
\url{http://jay.au.poznan.pl/html1/JAG/pdfy/lobry.pdf}

\bibitem{lobrygautier}
Lobry, J.R., Gautier, C.:
Hydrophobicity, expressivity and aromaticity are the major trends of amino-acid usage in 
999 \textit{Escherichia coli} chromosome-encoded genes. 
Nucleic Acids Res {\bf 22} (1994) 3174--3180.
\url{http://pbil.univ-lyon1.fr/members/lobry/repro/nar94/}

\bibitem{lobrysueoka}
Lobry, J.R., Sueoka, N.:
Asymmetric directional mutation pressures in bacteria.
Genome Biology {\bf 3} (2002) research0058.1--research0058.14.
\url{http://genomebiology.com/2002/3/10/research/0058}.

\bibitem{smorfland}
Mackiewicz, P., Zakrzewska-Czerwi{\'n}ska, J., Zawilak, A., Dudek, M.R., Cebrat, S.:
Where does bacterial replication start? Rules for predicting the \textit{oriC} region.
Nucleic Acids Res. {\bf 32} (2004) 3781--3791.

\bibitem{misuse}
Perri{\`e}re, G., Thioulouse, J:
Use and misuse of correspondence analysis in codon usage studies.
Nucleic Acids Res. {\bf 30} (2002) 4548--4555.

\bibitem{RfromR}
R Development Core Team:
R: A language and environment for statistical computing
(2004) ISBN 3-900051-00-3, http://www.R-project.org

\bibitem{chargaff}
Rudner, R., Karkas, J.D., Chargaff, E.:
Separation of microbial deoxyribonucleic acids into complementary strands. 
Proc. Natl. Acad. Sci. USA, {\bf 63} (1969) 152--159.

\bibitem{nj}
Saitou, N., Nei, M.:
The neighbor-joining method: a new method for reconstructing phylogenetic trees.
Mol. Biol. Evol. {\bf 4} (1984) 406--425.

\bibitem{CAI}
Sharp, P.M., Li, W.-H.:
The codon adaptation index - a measure of directional synonymous codon usage bias, and its potential applications.
Nucleic Acids Research {\bf 15} (1987) 1281---1295.

\bibitem{Karlin}
Karlin, S., Brendel, V.:
Chance and Statistical Significance in Protein and DNA Sequence Analysis.
Science {\bf 257} (1992) 39--49.

\bibitem{UV}
Palmeira, L., Guéguen L., Lobry, J.R.:
UV-targeted dinucleotides are not depleted in light-exposed Prokaryotic genomes.
\textit{in prep.}

\bibitem{Singer}
Singer, C.E., Ames, B.N.:
Sunlight Ultraviolet and Bacterial DNA Base Ratios.
Science {\bf 170} (1970) 822--826.

\bibitem{Bak}
Bak, A.L., Atkins, J.F., Singer, C.E., Ames, B.N.:
Evolution of DNA Base Compositions in Microorganisms.
Science {\bf 175} (1972) 1391--1393.

\bibitem{Setlow}
Setlow, R. B.:
Cyclobutane-Type Pyrimidine Dimers in Polynucleotides.
Science {\bf 153} (1966) 379--386.

\end{thebibliography}

\end{document}
