%\VignetteIndexEntry{Efficient genome searching with Biostrings and the BSgenome data packages}
%\VignetteKeywords{BSgenome, genome, DNA, RNA, Sequence, Biostrings, Sequence alignment, SNPs} 
%\VignettePackage{BSgenome}

%
% NOTE -- ONLY EDIT THE .Rnw FILE!!!  The .tex file is
% likely to be overwritten.
%
\SweaveOpts{keep.source=TRUE}
\documentclass[10pt]{article}

%\usepackage{amsmath}
%\usepackage[authoryear,round]{natbib}
\usepackage{hyperref}
\usepackage{underscore}

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

\newcommand{\scscst}{\scriptscriptstyle}
\newcommand{\scst}{\scriptstyle}

\newcommand{\R}{\textsf{R}}
\newcommand{\code}[1]{\texttt{#1}}
\newcommand{\term}[1]{\emph{#1}}
\newcommand{\Rpackage}[1]{\textsf{#1}}
\newcommand{\Rfunction}[1]{\texttt{#1}}
\newcommand{\Robject}[1]{\texttt{#1}}
\newcommand{\Rclass}[1]{\textit{#1}}
\newcommand{\Rmethod}[1]{\textit{#1}}
\newcommand{\Rfunarg}[1]{\textit{#1}}

\bibliographystyle{plainnat}
 
\begin{document}
%\setkeys{Gin}{width=0.55\textwidth}

\title{Efficient genome searching with Biostrings and the BSgenome data packages}
\author{Herv\'e Pag\`es}
\maketitle

\tableofcontents


% ---------------------------------------------------------------------------

\section{The Biostrings-based genome data packages}

The Bioconductor project provides data packages that contain the full genome
sequences of a given organism. These packages are called
{\it Biostrings-based genome data packages} because the sequences they
contain are stored in some of the basic containers defined in the
\Rpackage{Biostrings} package, like the \Rclass{DNAString},
the \Rclass{DNAStringSet} or the \Rclass{MaskedDNAString} containers.
Regardless of the particular sequence data that they contain, all the
Biostrings-based genome data packages are very similar and can be
manipulated in a consistent and easy way.
They all require the \Rpackage{BSgenome} package in order to work properly.
This package, unlike the Biostrings-based genome data packages,
is a software package that provides the infrastructure needed to
support them (this is why the Biostrings-based genome data packages
are also called {\it BSgenome data packages}).
The \Rpackage{BSgenome} package itself requires the \Rpackage{Biostrings}
package.

See the man page for the \Rfunction{available.genomes} function
(\Rfunction{?available.genomes}) for more information about how to get
the list of all the BSgenome data packages currently available in your
version of Bioconductor (you need an internet connection so that
\Rfunction{available.genomes} can query the Bioconductor package
repositories).

More genomes can be added if necessary. Note that the process of making a
BSgenome data package is not yet documented but you are welcome to ask for
help on the bioc-devel mailing list
(\url{http://bioconductor.org/docs/mailList.html}) if you need a genome
that is not yet available.



% ---------------------------------------------------------------------------

\section{Finding an arbitrary nucleotide pattern in a chromosome}

In this section we show how to find (or just count) the occurences of some
arbitrary nucleotide pattern in a chromosome. The basic tool for this is
the \Rfunction{matchPattern} (or \Rfunction{countPattern}) function
from the \Rpackage{Biostrings} package.

First we need to install and load the BSgenome data package for the
organism that we want to look at. In our case, we want to search
chromosome I of {\it Caenorhabditis elegans}.

UCSC provides several versions of the C. elegans genome: ce1, ce2 and ce4.
These versions correspond to different {\it releases} from WormBase,
which are the WS100, WS120 and WS170 releases, respectively.
See \url{http://genome.ucsc.edu/FAQ/FAQreleases#release1} for the list
of all UCSC genome releases and for the correspondance between UCSC
versions and release names.

The BSgenome data package for the ce2 genome is
\Rpackage{BSgenome.Celegans.UCSC.ce2}. Note that ce1 and ce4 are not
available in Bioconductor but they could be added if there is demand for
them.

% DELETED


% ---------------------------------------------------------------------------

\section{Finding an arbitrary nucleotide pattern in an entire genome}

Now we want to extend our analysis to the {\it forward} and {\it reverse}
strands of all the C. elegans chromosomes.
More precisely, here is the analysis we want to perform:

\begin{itemize}

  \item{The input dictionary: }
  Our input is a dictionary of 50 patterns. Each pattern is a short nucleotide
  sequence of 15 to 25 bases (As, Cs, Gs and Ts only, no Ns).
  It is stored in a FASTA file called \Robject{"ce2dict0.fa"}.
% ce2dict0.fa was generated with
%   > set.seed(23)
%   > ce2dict0 <- sapply(1:50, function(i) { x <- genome[[sample(seqnames(genome), 1)]]; start <- as.integer(runif(1, min=1, length(x)-20)); end <- start + as.integer(runif(1, min=14, max=25)); x <- subXString(x, start, end); if (sample(1:2, 1) == 2) x <- reverseComplement(x); as.character(x) })
%   > names(ce2dict0) <- paste("pattern", 1:50, sep="")
%   > write.XStringViews(XStringViews(ce2dict0, "DNAString"), file="ce2dict0.fa", format="fasta")
  See the {\it Finding all the patterns of a constant width dictionary
  in an entire genome} section of this document for a very efficient way
  to deal with the special case where all the patterns in the input
  dictionary have the same length.

  \item{The target: }
  Our target (or subject) is the {\it forward} and {\it reverse} strands of
  the seven C. elegans chromosomes (14 sequences in total).
  We want to find and report all occurences (or hits) of every pattern
  in the target. Note that a given pattern can have 0, 1 or several hits
  in 0, 1 or 2 strands of 0, 1 or several chromosomes.

  \item{Exact or inexact matching? }
  We are interested in exact matches only (for now).

  \item{The output: }
  We want to put the results of this analysis in a file so we can send
  it to our collaborators for some post analysis work.
  Our collaborators are not necessarily familiar with R or Bioconductor
  so dumping a high-level R object (like a list or a data frame) into an
  .rda file is not an option. For maximum portability (one of our
  collaborators wants to use Microsoft Excel for the post analysis) we
  choose to put our results in a tabulated file where one line describes
  one hit. The columns (or fields) of this file will be (in this order):
  \begin{itemize}
    \item{seqname: }
    the name of the chromosome where the hit occurs.
    \item{start: }
    an integer giving the starting position of the hit.
    \item{end: }
    an integer giving the ending position of the hit.
    \item{strand: }
    a plus (\Robject{+}) for a hit in the positive strand
    or a minus (\Robject{-}) for a hit in the negative strand.
    \item{patternID: }
    we use the unique ID provided for every pattern in the
    \Robject{"ce2dict0.fa"} file.
  \end{itemize}

\end{itemize}

% DELETED


% ---------------------------------------------------------------------------

\section{Some precautions when using \Rmethod{matchPattern}}

Improper use of \Rmethod{matchPattern} (or \Rmethod{countPattern}) can affect
performance.

If needed, the \Rmethod{matchPattern} and \Rmethod{countPattern} methods
convert their first argument (the pattern) to an object of the same class
than their second argument (the subject) before they pass it to the subroutine
that actually implements the fast search algorithm.

So if you need to reuse the same pattern a high number of times,
it's a good idea to convert it {\it before} to pass it to the
\Rmethod{matchPattern} or \Rmethod{countPattern} method.

% DELETED

% The above example is Raphael's use case discussed on BioC on Feb 2006.
% In Biostrings 1, the equivalent would be:
% src <- sapply(1:100,
%               function(i) {
%                 paste(sample(c("A","C","G","T"), 25, replace=TRUE),
%                 collapse="")
%               }
%        )
% tmpseq <- DNAString(src)
% someStats <- function(v)
% {
%     GC <- DNAString("GC")
%     CG <- DNAString("CG")
%     sapply(1:length(v),
%            function(i) {
%                y <- v[i]
%                c(alphabetFrequency(y)[2:5],
%                  GC=length(matchDNAPattern(GC, y)),
%                  CG=length(matchDNAPattern(CG, y)))
%            }
%     )
% }
% someStats(tmpseq[1:10])



% ---------------------------------------------------------------------------

\section{Masking the chromosome sequences}

Starting with Bioconductor 2.2, the chromosome sequences in a {\it BSgenome
data package} can have built-in masks.
Starting with Bioconductor 2.3, there can be up to 4 built-in masks per
sequence. These will always be (in this order): (1) the mask of assembly gaps,
(2) the mask of intra-contig ambiguities, (3) the mask of repeat regions that
were determined by the RepeatMasker software, and (4) the mask of repeat
regions that were determined by the Tandem Repeats Finder software (where only
repeats with period less than or equal to 12 were kept).

For a given package, all the sequences will always have the same number of
masks.

% DELETED

The built-in masks are named consistenly across all the BSgenome data packages
available in Bioconductor:

\begin{table}[ht]
\begin{center}
\begin{tabular}{l|l|l|l}
\hline
Name  & Active by default & Short description & Long description \\
\hline
AGAPS & yes & assembly gaps
      & Masks the big N-blocks that have been placed
        between the contigs during the assembly.
        This mask is consistent with the Gap track from UCSC
        Genome Browser.\\
AMB   & yes & intra-contig ambiguities
      & Masks any IUPAC ambiguity letter that was found in the contig
        regions of the original sequence.
        Note that only As, Cs, Gs and Ts remain unmasked when the
        AGAPS and AMB masks are both active (before SNPs are eventually
        injected, see below). \\
RM    & no & RepeatMasker
      & Masks the repeat regions determined by the RepeatMasker software.
        This mask is consistent with the RepeatMasker track from
        UCSC Genome Browser. \\
TRF   & no & Tandem Repeats Finder
      & Masks the tandem repeat regions that were determined by the
        Tandem Repeats Finder software (with period of 12 or less). \\
\hline
\end{tabular}
\end{center}
\caption{The built-in masks provided by the BSgenome data packages.}
\end{table}

When displaying a masked sequence (here a \Rclass{MaskedDNAString}
object), the {\it masked width} and {\it masked ratio} are reported for
each individual mask, as well as for all the masks together, and for all
the active masks together.
The {\it masked width} is the total number of nucleotide positions
that are masked and the {\it masked ratio} is the {\it masked width}
divided by the length of the sequence.

% DELETED


% ---------------------------------------------------------------------------

\section{Hard masking}

coming soon...



% ---------------------------------------------------------------------------
\section{Injecting known SNPs in the chromosome sequences}

coming soon...



% ---------------------------------------------------------------------------

\section{Finding all the patterns of a constant width dictionary
         in an entire genome}

The \Rfunction{matchPDict} function can be used instead of \Rfunction{matchPattern}
for the kind of analysis described in the {\it Finding an arbitrary nucleotide
pattern in an entire genome} section but it will be much faster (between 100x and
10000x faster depending on the size of the input dictionary).
Note that a current limitation of \Rfunction{matchPDict} is that it only works with
a dictionary of DNA patterns where all the patterns have the same number of
nucleotides (constant width dictionary).
See \Rfunction{?matchPDict} for more information.

Here is how our \Rfunction{runAnalysis1} function can be modified in order to use
\Rfunction{matchPDict} instead of \Rfunction{matchPattern}:
<<e1>>=
runOneStrandAnalysis <- function(dict0, bsgenome, seqnames, strand,
                                 outfile="", append=FALSE)
{
    cat("\nTarget: strand", strand, "of", metadata(bsgenome)$genome,
        "chromosomes", paste(seqnames, collapse=", "), "\n")
    if (strand == "-")
        dict0 <- reverseComplement(dict0)
    pdict <- PDict(dict0)
    for (seqname in seqnames) {
        subject <- bsgenome[[seqname]]
        cat(">>> Finding all hits in strand", strand, "of chromosome", seqname, "...\n")
        mindex <- matchPDict(pdict, subject)
        matches <- extractAllMatches(subject, mindex)
        writeHits(seqname, matches, strand, file=outfile, append=append)
        append <- TRUE
        cat(">>> DONE\n")
    }
}

runAnalysis2 <- function(dict0, outfile="")
{
    library(BSgenome.Celegans.UCSC.ce2)
    genome <- BSgenome.Celegans.UCSC.ce2
    seqnames <- seqnames(genome)
    runOneStrandAnalysis(dict0, genome, seqnames, "+", outfile=outfile, append=FALSE)
    runOneStrandAnalysis(dict0, genome, seqnames, "-", outfile=outfile, append=TRUE)
}
@

Remember that \Rfunction{matchPDict} only works if all the patterns in the
input dictionary have the same length so for this 2nd analysis, we will truncate
the patterns in \Robject{ce2dict0} to 15 nucleotides:

% STRANGE THAT ALSO THIS DOES NOT WORK ... ???
%<<e2>>=
%ce2dict0cw15 <- DNAStringSet(ce2dict0, end=15)
%@

% DELETED


% ---------------------------------------------------------------------------

\section{Session info}

<<g1>>=
sessionInfo()
@

\end{document}


