%\VignetteIndexEntry{02. Ranges and strings}
%\VignettePackage{useR2011}

\SweaveOpts{keep.source=TRUE}
<<setup,echo=FALSE>>=
library(useR2011)
plotRanges <- function(x, xlim = x, main = deparse(substitute(x)),
                       col = "black", sep = 0.5, ...)
{
  height <- 1
  if (is(xlim, "Ranges"))
    xlim <- c(min(start(xlim)), max(end(xlim)))
  bins <- disjointBins(IRanges(start(x), end(x) + 1))
  plot.new()
  par(mai=c(0.6, 0.2, 0.6, 0.2))
  plot.window(xlim, c(0, max(bins)*(height + sep)))
  ybottom <- bins * (sep + height) - height
  rect(start(x)-0.5, ybottom, end(x)+0.5, ybottom + height, col = col, ...)
  title(main, cex.main=2.8, font.main=1)
  axis(1)
}
@ 

\section{Ranges and strings}

This section introduces two essential ways in which sequence data are
manipulated. Ranges describe both aligned reads and features of
interest on the genome. Sets of DNA strings represent the reads
themselves and the nucleotide sequence of reference genomes.

\subsection{Reads and genomic features as ranges}

Next-generation sequencing data consists of a large number of short
reads. These are, typically, aligned to a reference genome. Basic
operations are performed on the alignment, e.g., how many reads are
aligned in a genomic range defined by nucleotide coordinates (e.g., in
the exons of a gene), or how many nucleotides from all the aligned
reads cover a set of genomic coordinates. How is this type of data,
the aligned reads and the reference genome, to be represented in \R{}
in a way that allows for effective computation?  

The \Biocpkg{IRanges}, \Biocpkg{GenomicRanges}, and
\Biocpkg{GenomicFeatures} \Bioconductor{} packages provide the
essential infrastructure for these operations; we start with the
\Rclass{GRanges} class, defined in \Biocpkg{GenomicRanges}.  

\paragraph{\Rclass{GRanges}}

Instances of \Rclass{GRanges} are used to specify genomic
coordinates. Suppose we wished to represent two \Dmel{} genes. The
first is located on the positive strand of chromosome 3R, from
position 19967116 to 19973212. The second is on the minus strand of
the X chromosome, with `left-most' base at 18962305, and right-most
base at 18962925. The coordinates are \emph{1-based} (i.e., the first
nucleotide on a chromosome is numbered 1, rather than 0),
\emph{left-most} (i.e., reads on the minus strand are defined to
`start' at the left-most coordinate, rather than the 5' coordinate),
and \emph{closed} (the start and end coordinates are included in the
range; a range with identical start and end coordinates has width 1, a
0-width range is represented by the special construct where the end
coordinate is one less than the start coordinate).

A complete definition of these genes as \Rclass{GRanges} is:
<<GRanges-genes>>=
genes <- GRanges(seqnames=c("3R", "X"), 
                 ranges=IRanges(
                     start=c(19967116, 18962305),
                     end=c(19973212, 18962925)),
                 strand=c("+", "-"),
                 seqlengths=c(`3R`=27905053L, `X`=22422827L))
@ 
The components of a gene are defined as vectors, e.g., of
\Rfunarg{seqnames}, much as one would define a
\Rclass{data.frame}. The start and end coordinates are grouped into an
\Rclass{IRanges} instance. The optional \Rfunarg{seqlengths} argument
specifies the maximum size of each sequence, in this case the
lengths of chromosomes 3R and X in \Dmel. This data is displayed as
<<GRanges-display>>=
genes
@ 
For the curious, the gene coordinates and sequence lengths are derived
from the \Biocdatapkg{org.Dm.eg.db} for genes with Flybase identifiers
FBgn0039155 and FBgn0039155, using the annotation facilities described
in section~\ref{sec:annotation} on annotation.

The \Rclass{GRanges} class has many useful methods defined on
it. Consult the help page
<<GRanges-help, eval=FALSE>>=
?GRanges
@ 
and package vignettes (especially `An Introduction to
\Biocpkg{GenomicRanges}')
<<GRanges-vignettes, eval=FALSE>>=
browseVignettes("GenomicRanges")
@ 
for a comprehensive introduction. A \Rclass{GRanges} instance can be
subset, with accessors for setting, updating, and querying information.
<<GRanges-basics>>=
genes[2]
strand(genes)
width(genes)
length(genes)
names(genes) <- c("FBgn0039155", "FBgn0085359")
genes # now with names
@ 
\Rfunction{strand} returns the strand information in a compact
representation called a \emph{run-length encoding}, this is introduced
in greater detail below.  The `names' could have been specified when
the instance was constructed; once named, the \Rclass{GRanges}
instance can be subset by name like a regular vector.

As the \Rfunction{GRanges} function suggests, the \Rclass{GRanges}
class extends the \Rclass{IRanges} class by adding information about
\Rfunarg{seqname}, \Rfunarg{strand}, and other information
particularly relevant to representing ranges that are on genomes. The
\Rclass{IRanges} class and related data structures (e.g.,
\Rclass{RangedData}) are meant as a more general description of ranges
defined in an arbitrary space. Many methods implemented on the
\Rclass{GRanges} class are `aware' of the consequences of genomic
location, for instance treating ranges on the minus strand differently
(reflecting the 5' orientation imposed by DNA) from ranges on the plus
strand.

\paragraph{Operations on ranges}

The \Rclass{GRanges} class inherits many useful methods from the
\Rclass{IRanges} class; some of these methods are illustrated here. We
use \Rclass{IRanges} to illustrate these operations to avoid
complexities associated with strand and seqname, but the operations
are comparable on \Rclass{GRanges}. We begin with a simple set of
ranges:
<<ranges-ir>>=
ir <- IRanges(start=c(7, 9, 12, 14, 22:24), 
              end=c(15, 11, 12, 18, 26, 27, 28))
@
%% 
<<ranges-ir-plot, results=hide, echo=FALSE, fig=FALSE>>=
png("ranges-ir-plot.png", width=800, height=160)
plotRanges(ir, xlim=c(5, 35), main="Original")
dev.off()
@
These are illustrated in the upper left panel of
Figure~\ref{fig:ranges}.

\begin{figure}
\centering
\includegraphics[width=.5\textwidth,height=!]{ranges-ir-plot}\\
\includegraphics[width=.5\textwidth,height=!]{ranges-shift-ir-plot}\\
\includegraphics[width=.5\textwidth,height=!]{ranges-reduce-ir-plot}
\caption{Ranges}
\label{fig:ranges}
\end{figure}

Methods on ranges can be grouped as follows:
\begin{description}
\item[Intra-range] methods act on each range independently. These
  include \Rfunction{flank}, \Rfunction{narrow}, \Rfunction{reflect},
  \Rfunction{resize}, \Rfunction{restrict}, \Rfunction{shift}. An
  illustration is \Rfunction{shift}, which translates each range by
  the amount specified by the \Rfunarg{shift} argument.  Positive
  values shift to the right, negative to the left; \Rfunarg{shift} can
  be a vector, with each element of the vector shifting the
  corresponding element of the \Rclass{IRanges} instance. Here we
  shift all ranges to the right by 5, with the result illustrated in
  the middle panel of Figure~\ref{fig:ranges}
<<ranges-shift-ir>>=
shift(ir, 5)
@
<<ranges-shift-ir-plot, results=hide, echo=FALSE, fig=FALSE>>=
png("ranges-shift-ir-plot.png", width=800, height=160)
plotRanges(shift(ir, 5), xlim=c(5, 35), main="Shift")
dev.off()
@

\item[Inter-range] methods act on the collection of ranges as a
  whole. These include \Rfunction{disjoin}, \Rfunction{reduce},
  \Rfunction{gaps}, and \Rfunction{range}. An illustration is
  \Rfunction{reduce}, which reduces overlapping ranges into a single
  range, as illustrated in the lower panel of
  Figure~\ref{fig:ranges}.
<<ranges-reduce-ir>>=
reduce(ir)
@
<<ranges-reduce-ir-plot, results=hide, echo=FALSE, fig=FALSE>>=
png("ranges-reduce-ir-plot.png", width=800, height=160)
plotRanges(reduce(ir), xlim=c(5, 35), main="Reduce")
dev.off()
@

  \Rfunction{coverage} is an inter-range operation that calculates how
  many ranges overlap individual positions. Rather than returning
  ranges, \Rfunction{coverage} returns a compressed representation of
  an integer vector known as \Rclass{Rle} (run-length encoding)
<<coverage>>=
coverage(ir)
@    
  The run-length encoding can be interpreted as `a run of length 6 of
  nucleotides covered by 0 ranges, followed by a run of length 2 of
  nucleotides covered by 1 range\ldots'.

\item[Between] methods act on two (or sometimes more) \Rclass{IRanges}
  instances. These include \Rfunction{intersect},
  \Rfunction{setdiff}, \Rfunction{union}, \Rfunction{pintersect},
  \Rfunction{psetdiff}, and \Rfunction{punion}.
  
  \Rfunction{countOverlaps} and \Rfunction{findOverlaps} also operate
  on two sets of ranges. \Rfunction{countOverlaps} takes its first
  argument (the \Rfunarg{query}) and determines how many of the ranges
  in the second argument (the \Rfunarg{subject}) each overlaps. The
  result is an integer vector with one element for each member of
  \Rfunarg{query}. \Rfunction{findOverlaps} performs a similar
  operation but returns a more general matrix-like structure that
  identifies each pair of query / subject overlap. Both arguments
  allow some flexibility in the definition of `overlap'.
\end{description}

\paragraph{\Rfunction{elementMetadata} and \Rfunction{metadata}}

The \Rclass{GRanges} class (actually, most of the data structures
defined or extending those in the \Biocpkg{IRanges} package) has two
additional very useful data components. The
\Rfunction{elementMetadata} function (or its synonym
\Rfunction{values}) allows information on each range to be
stored and manipulated (e.g., subset) along with the \Rclass{GRanges}
instance. The element metadata is represented as a \Rclass{DataFrame},
defined in \Biocpkg{IRanges} and acting like a standard \R{}
\Rclass{data.frame} but with the ability to hold more complicated data
structures as columns (and with element metadata of its own, providing
an enhanced alternative to the \Biocpkg{Biobase} class
\Rclass{AnnotatedDataFrame}).
<<GRanges-elementMetadata>>=
elementMetadata(genes) <- 
    DataFrame(EntrezId=c("42865", "42865"),
              Symbol=c("kal-1", "CG6173"))
@ 
\Rfunction{metadata} allows addition of information to the entire
object. The information is in the form of a list; any data can be
provided.
<<GRanges-metadata>>=
metadata(genes) <- 
    list(CreatedBy="Martin Morgan", Date=date())
@ 

\subsection{Genomic features}

The \Rfunction{GRanges} class is extremely useful for representing
simple ranges. Some next-generation sequence data and genomic features
are more hierarchically structured. A gene may be represented by
several exons within it. An aligned read may be represented by
discontinuous ranges of alignment to a reference. 

The \Rclass{GRangesList} class represents this type of information. It
is a list-like data structure, which each element of the list itself a
\Rclass{GRanges} instance. The gene FBgn0039155 contains several
exons, and can be represented as a length 1 list, where the element of
the list contains a \Rclass{GRanges} object with 7 elements:
<<GRangesList-eg-setup, echo=FALSE>>=
txdbFile <- list.files(bigdata(), "sqlite", full=TRUE)
txdb <- loadFeatures(txdbFile)
fbgn <- exonsBy(txdb, "gene")["FBgn0039155"]
seqlevels(fbgn) <- "chr3R"
@ 
<<GRangesList-eg, echo=FALSE>>=
fbgn
@ 

The \Rclass{GRangesList} object has methods one would expect for lists
(e.g., \Rfunction{length}, subsetting). Many of the methods
introduced for working with \Rclass{IRanges} are also available,
with the method applied element-wise.

\paragraph{The \Biocpkg{GenomicFeatures} package}

Many public resources provide annotations about genomic features. For
instance, the UCSC genome browser maintains the `knownGenes' track of
established exons, transcripts, and coding sequences of many model
organisms. The \Biocpkg{GenomicFeatures} package provides a way to
retrieve, save, and query these resources. The underlying
representation is as sqlite data bases, but the data are available in
\R{} as \Rclass{GRangesList} objects. The following exercise explores
the \Biocpkg{GenomicFeatures} package and some of the functionality
for the \Rclass{IRanges} family introduced above.

\begin{Exercise}
  Use the helper function \Rfunction{bigdata} and
  \Rfunction{list.files} to identify the path to a data base created
  by \Rfunction{makeTranscriptDbFromUCSC}.
  
  Load the saved \Rclass{TranscriptDb} object using
  \Rfunction{loadFeatures}.
  
  Extract all exon coordinates, extracted by gene, using
  \Rfunction{exonsBy}. What is the class of this object? How many
  elements are in the object? What does each element correspond to?
  And the elements of each element? Use \Rfunction{elementLengths} and
  \Rfunction{table} to summarize the number of exons in each gene, for
  instance, how many single-exon genes are there?
  
  Select just those elements corresponding to flybase gene ids
  FBgn0002183, FBgn0003360, FBgn0025111, and FBgn0036449. Use
  \Rfunction{reduce} to simplify gene models, so that exons that
  overlap are considered `the same'.
\end{Exercise}

\begin{Solution}
<<txdb>>=
txdbFile <- list.files(bigdata(), "sqlite", full=TRUE)
txdb <- loadFeatures(txdbFile)
ex0 <- exonsBy(txdb, "gene")
head(table(elementLengths(ex0)))
ids <- c("FBgn0002183", "FBgn0003360", "FBgn0025111", "FBgn0036449")
ex <- reduce(ex0[ids])
@ 
\end{Solution}

\begin{Exercise}
  (Independent) Create a \Rclass{TranscriptDb} instance from UCSC,
  using \Rfunction{makeTranscriptDbFromUCSC}.
\end{Exercise}

\begin{Solution}
<<txdb, eval=FALSE>>=
txdb <- makeTranscriptDbFromUCSC("dm3", "ensGene")
saveFeatures(txdb, "my.dm3.ensGene.txdb.sqlite")
@ 
\end{Solution}

\subsection{Working with strings}

Underlying the ranges of alignments and features are DNA
sequences. The \Biocpkg{Biostrings} package provides tools for working
with this data. The essential data structures are \Rclass{DNAString}
and \Rclass{DNAStringSet}, for working with one or multiple DNA
sequences. The \Biocpkg{Biostrings} package contains additional
classes for representing amino acid and general biological strings.
The \Biocpkg{BSgenome} and related packages (e.g.,
\Biocdatapkg{BSgenome.Dmelanogaster.UCSC.dm3}) are used to represent
whole-genome sequences. The following exercise explores these
packages.

\begin{Exercise}
  The objective of this exercise is to calculate the GC content of the
  exons of a single gene, whose coordinates are specified by the
  \Robject{ex} object of the previous exercise.
  
  Load the \Biocdatapkg{BSgenome.Dmelanogaster.UCSC.dm3} data package,
  containing the USCS representation of \Dmel{} genome assembly dm3.
  
  Extract the sequence name of the first gene of \Robject{ex}. Use
  this to load the appropriate \Dmel{} chromosome.
  
  Use \Rfunction{Views} to create views on to the chromosome that span
  the start and end coordinates of all exons.
  
  The \Rpackage{useR2011} package defines a helper function
  \Rfunction{gcFunction} (developed in a later exercise) to calculate
  GC content. Use this to calculate the GC content in each of the exons.
\end{Exercise}

\begin{Solution}
<<gc-genome>>=
library(BSgenome.Dmelanogaster.UCSC.dm3)
nm <- as.character(unique(seqnames(ex[[1]])))
chr <- Dmelanogaster[[nm]]
v <- Views(chr, start=start(ex[[1]]), end=end(ex[[1]]))
@ 
Here is the helper function, available in the \Rpackage{useR2011}
package, to calculate GC content:
<<gcFunction>>=
gcFunction
@ 
The subject GC content is
<<gcFunction-genome, eval=FALSE>>=
subjectGC <- gcFunction(v)
@ 
\end{Solution}
