%\VignetteDepends{IRanges, GenomicRanges, Biostrings, BSgenome, GenomicFeatures, BSgenome.Scerevisiae.UCSC.sacCer2, SeattleIntro2010}
%\VignetteIndexEntry{Working with DNA strings and ranges: Exercises}
%\VignettePackage{SeattleIntro2010}

\documentclass{article}

\usepackage{BioconductorVignette}
\usepackage{underscore}

\SweaveOpts{keep.source=TRUE, eps=FALSE, include=FALSE}

\title{Working with DNA strings and ranges: Exercises}
\author{Herv\'e Pag\`es\footnote{Fred Hutchinson Cancer Research Center, Seattle, WA 98008}}
\date{9-10 December, 2010}

\begin{document}

<<setup, echo=FALSE>>=
options(width = 60)
library(GenomicFeatures)
library(BSgenome.Scerevisiae.UCSC.sacCer2)
library(SeattleIntro2010)
@

\maketitle

\tableofcontents


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Preliminaries}

First let's make sure you have a working \R{}-2.12 installation with all
the packages needed for the course installed. To install any missing
package (or to update the installed packages), start \R{} and use the
following command:

<<preliminaries,eval=FALSE>>=
source("http://bioconductor.org/course-packages/install-SeattleIntro2010.R")
@


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Use case I: Extracting sequences from a reference genome}

In this introductory use case, we learn how to extract DNA sequences
from a \textit{BSgenome data package} for a set of given locations.
In particular we use the transcript, exon and CDS locations stored
in a \Rclass{TranscriptDb} object to extract the sequences of those
features.

The \Bioconductor{} data repositories provide a \textit{BSgenome data
package} for the sacCer2 genome (Yeast): the
\Biocpkg{BSgenome.\-Scerevisiae.\-UCSC.sacCer2} package.
You should normally have it installed.
It contains the full DNA sequences of the sacCer2 genome.

Also the \Biocpkg{SeattleIntro2010} package contains a \Rclass{TranscriptDb}
object corresponding to this genome. Note that a \Rclass{TranscriptDb} object
contains positional (and relational) information about features but it
does not contain sequences. If we need to extract the sequence of a given
feature, an easy way is to query the
\Biocpkg{BSgenome.\-Scerevisiae.\-UCSC.sacCer2} package with the
\Rfunction{getSeq} function.

Note that the result of this query is meaningful only if the
\Rclass{TranscriptDb} object contains positional annotations relative
to the genome stored in the \textit{BSgenome data package}.
For example using a \Rclass{TranscriptDb} object based on
\Biocpkg{BSgenome.\-Hsapiens.\-UCSC.hg18} to extract sequences from
\Biocpkg{BSgenome.\-Hsapiens.\-UCSC.hg19} would not make sense.

\begin{Exercise}
  \begin{itemize}

  \item Start \R{} and load the \Biocpkg{SeattleIntro2010} package.

  \item Use \Rcode{system.file(package="SeattleIntro2010")} to get the
  full path to the top-level folder of the installed package.
  (The top-level folder of an installed package should always be
  treated as read-only).

  \item Use \Rfunction{list.files} on the previous result.

  \item The \Rclass{TranscriptDb} object that we are looking for is in the
  \code{extdata} subfolder. Use \Rcode{list.files(system.file("extdata",
  package="SeattleIntro2010"))} to see it. It's the
  \code{sacCer2_sgdGene.sqlite} file.

  \item The \code{sacCer2_sgdGene.sqlite} file is an SQLite database
  that stores the transcript, exon and CDS locations relative to the
  sacCer2 genome as well as the relations between those features and
  their corresponding genes. This information was extracted from the
  ``SGD Genes'' track for sacCer2 at the UCSC Genome Browser, and
  formatted into an SQLite database that can be loaded in \R{} with the
  \Rfunction{loadFeatures} function from the \Biocpkg{GenomicFeatures}
  package.

  \item Load the \Biocpkg{GenomicFeatures} package and use the
  \Rfunction{loadFeatures} function to load the \code{sacCer2_sgdGene.sqlite}
  db into your session.
  Let's call \Robject{txdb} the returned object.

  \end{itemize}
\end{Exercise}

\Robject{txdb} is a \Rclass{TranscriptDb} object. You will learn
more about those objects in the next session (in particular, how to
make your own). For now, we're just going to extract all the exon
locations from it and then query the
\Biocpkg{BSgenome.\-Scerevisiae.\-UCSC.sacCer2} with \Rfunction{getSeq}
to extract their sequences.

\begin{Exercise}
  \begin{itemize}

  \item Extract all the exon locations from \Robject{txdb} with the
  \Rfunction{exons} function. The result is a \Rclass{GRanges} object.

  \item How many exons? What \Rclass{GRanges} accessor would you use
  to extract the lengths of those exons?

  \item What's the longest exon? Where is it located?

  \item Are there exons on the 2micron plasmid?

  \item Load the \Biocpkg{BSgenome.\-Scerevisiae.\-UCSC.sacCer2} package.

  \item There is only one symbol defined in this package, the
  \Robject{Scerevisiae} object (you can check this with
  \Rcode{ls("package:BSgenome.Scerevisiae.UCSC.sacCer2")}). Display it.

  \item Get the lengths of the sequences with \Rcode{seqlengths(Scerevisiae)}.

  \item Try to load any sequence with
  \Rcode{Scerevisiae[["some sequence name"]]}.

  \item Have a quick look at the man page for the \Rfunction{getSeq} function
  (with \Rcode{?getSeq}) and use it on \Robject{Scerevisiae} and the
  \Rclass{GRanges} object created previously to extract the exons sequences.
  An important caveat here is to use \Rcode{as.character=FALSE}.

  \end{itemize}
\end{Exercise}

We end up with a \Rclass{DNAStringSet} object containing our exon sequences.
Later we will learn more about \Rclass{DNAString} and \Rclass{DNAStringSet}
objects.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Use case II: Importing and manipulating a \Rclass{GappedAlignments}
         object}

An high-throughput sequencing experiment produces reads that need to be
aligned against a reference genome. Although the \Bioconductor{} software
provides some fast pattern matching tools that can be used for aligning
the reads, this step is typically done with a third-party software like
Bowtie, BWA, Eland, etc...

The \Biocpkg{SeattleIntro2010} package contains reads from the Nagalakshmi
et al.~\cite{pmid18451266} experiment (Yeast RNA-seq). They have been aligned
against the sacCer2 reference genome (using the BWA software) and stored
in a BAM file. In order to keep the package to a reasonable size, only
the reads from a single lane (oligo(dT)-primed, original) with hits on
chromosome I to V have been kept. They are in the
\code{SRR002051.chrI-V.bam} file.

\Bioconductor{} provides several tools to load a BAM file:
\begin{itemize}
  \item The low-level \Rfunction{scanBam} function from the
  \Biocpkg{Rsamtools} package. Returns a list of lists.
  \item The middle-level \Rfunction{readGappedAlignments} function
  from the \Biocpkg{GenomicRanges} package. Returns a
  \Rclass{GappedAlignments} object.
  \item The high-level \Rfunction{readAligned} function from
  the \Biocpkg{ShortRead} package. Returns an \Rclass{AlignedRead}
  object.
\end{itemize}

Here we will focus on the middle-level solution. As we will see,
\Rclass{GappedAlignments} objects don't store as much information
as \Rclass{AlignedRead} objects (e.g. the read sequences, read qualities
and alignment scores are not stored), but, unlike \Rclass{AlignedRead}
objects, they can store alignments with indels and gaps.

\begin{Exercise}
  \begin{itemize}

  \item Load the \code{SRR002051.chrI-V.bam} file (located in the
  \code{extdata} subfolder of the \Biocpkg{SeattleIntro2010} package)
  with the \Rfunction{readGappedAlignments} function from the
  \Biocpkg{GenomicRanges} package. Let's call the result \Robject{galn}.

  \item Display \Robject{galn}. What's its class? What's its length?

  \end{itemize}
\end{Exercise}

\Robject{galn} is a \Rclass{GappedAlignments} object. According
to the man page for those objects (\Rcode{?GappedAlignments}):

\vspace{\baselineskip}
{\addtolength{\leftskip}{0.5in}\addtolength{\rightskip}{0.5in}
  \noindent A GappedAlignments object is a vector-like object where each
  element describes an \textit{alignment}, that is, how a given sequence
  (called ``query'' or ``read'', typically short) aligns to a reference
  sequence (typically long).
\par}
\vspace{\baselineskip}

\textit{Alignments} are also called ``hits'' or ``matches''.
It's important to note that the relationship between the original set
of queries to align (i.e. the input of the aligner software) and the
elements of the \Rclass{GappedAlignments} object is generally not one-to-one.
Queries with no hits don't show up in the \Rclass{GappedAlignments} object.
Queries with multiple hits might have been part of the output of the
aligner software and stored in the BAM file (by assigning more than
one BAM record to the same query ID). In that case the resulting
\Rclass{GappedAlignments} object would have more than one element
corresponding to the same query ID.

Like for \Rclass{TranscriptDb}, the information stored in a
\Rclass{GappedAlignments} object is positional only (no sequence is
stored), and positions are always 1-based with respect to the 5' end of
the plus strand of the reference sequence.

The basic accessors for a \Rclass{GappedAlignments} object are
\Rfunction{rname}, \Rfunction{strand}, \Rfunction{cigar}, \Rfunction{qwidth},
\Rfunction{start}, \Rfunction{end}, \Rfunction{width} and \Rfunction{ngap}.
They all return a vector of the same length as the \Rclass{GappedAlignments}
object itself.

Note that the ``r'' in \Rfunction{rname} and the ``q'' in \Rfunction{qwidth}
stand for ``reference'' and ``query'', respectively, that is, the
\Rfunction{rname} accessor gets (or sets) the name of the \textit{reference}
sequence for each alignment and the \Rfunction{qwidth} accessor gets
the width of the \textit{query} (i.e. its number of bases).
Although they are not prefixed with ``r'', \Rfunction{start}, \Rfunction{end}
and \Rfunction{width} are relative to the reference genome.

In the next exercise we look more closely at the difference between
\Rfunction{qwidth} and \Rfunction{width}.

\begin{Exercise}
  \begin{itemize}

  \item Compare the output of \Rcode{table(qwidth(galn))} and
  \Rcode{table(width(galn))}. What we observe is that, even if all the
  queries are 33-base long, the portion of the reference sequence that
  they are aligned to can have a slightly different length.
  This is because those alignments contain \textit{indels}
  (represented with I's or a D's in the CIGAR) or \textit{gaps}
  (represented with N's in the CIGAR).

  \item Do the alignments have gaps?

  \item The \Biocpkg{GenomicRanges} package provides some CIGAR utility
  functions documented in \Rcode{?\`{}cigar-utils\`} (note the use of the
  backticks).
  In order to find the alignments that don't have a \textit{simple cigar}
  (\textit{simple cigar} means ``only M's in it''), we're going to use
  the \Rfunction{cigarOpTable} function.
  Extract the CIGAR strings from the \Robject{galn} object and pass them
  to \Rfunction{cigarOpTable}.

  \item What we have now is a numeric matrix with one row per CIGAR string
  and one column per valid CIGAR letter \footnote{This set of CIGAR letters
  corresponds to the extended CIGAR specification as described in the
  SAMtools spec -- see \url{http://samtools.sourceforge.net/} for more
  information}.
  Drop the M column and use \Rfunction{rowSums} to collapse the remaining
  columns into a single column.

  \item From here it should be easy to extract the index of alignments
  that don't have a \textit{simple cigar}. Use this index to subset
  \Robject{galn}.

  \end{itemize}
\end{Exercise}

One last thing before we move on to the next section.

\begin{Exercise}
  \begin{itemize}

  \item Turn \Robject{galn} into a \Rclass{GRanges} object (let's call
  it \Robject{hits}). Tip: Consult the man page for \Rclass{GappedAlignments}
  objects for how to do this.

  \end{itemize}
\end{Exercise}

We will use this \Rclass{GRanges} object in the next section.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Use case III: Looking at the reads that don't hit a (known) gene}

In this section, we want to identify the reads that don't hit a (known)
gene as well as the regions in the genomes covered by those reads.
Because those reads are coming from an RNA-seq experiment, those regions
would be good candidates for de-novo gene/transcript discovery.

\begin{Exercise}
  \begin{itemize}

  \item Load the \code{sacCer2_sgdGene.sqlite} file (located in the
  \code{extdata} subfolder of the \Biocpkg{SeattleIntro2010} package)
  with \Rfunction{loadFeatures}. Let's call \Robject{txdb} this
  \Rclass{TranscriptDb} object.

  \item Use the \Rfunction{countOverlaps} function to find the elements
  in the \Robject{hits} object (from the previous section) that don't
  hit any gene.
  Note that there are basically 2 approaches to this: one based on the
  start/end of the transcripts and one based on the start/end of the exons.
  Let's call \Robject{hits0} the subset of \Robject{hits} made of
  those elements that don't hit any gene.

  \item Use the \Rfunction{reduce} function to extract the regions covered
  by \Robject{hits0}. Additionally you can perform a sanity check by
  making sure those regions don't intersect with the exons.

  \end{itemize}
\end{Exercise}

We want to refine the way we've determined the above regions by using a
method based on the \textit{depth} of the coverage of the hits. More
precisely we want to find the regions of the genome where the coverage of
\Robject{hits0} is greater than (or equal to) some threshold (e.g. 10).
An alternative way to formulate this is: we want to find the regions of
the genome corresponding to all the bases that receive at least 10 hits.

\begin{Exercise}
  \begin{itemize}

  \item Compute the genome wide coverage of \Robject{hits0}. The result
  is a \Rclass{SimpleRleList} object. This sounds complicated but
  it helps to think of it as a named list of \Rclass{Rle} objects.
  The names of the list are the chromosomes and each top-level element in
  the list is an \Rclass{Rle} object representing the coverage for the
  corresponding chromosome. Let's call this \Rclass{SimpleRleList} object
  \Robject{cov0}. Like you would do with a standard \Rcode{list} object,
  you can use the double-bracket operator (\Rcode{[[}) to extract an
  element from \Robject{cov0}.

  \item Use the \Rfunction{slice} function to \textit{slice} \Robject{cov0}
  horizontally. The result of this is a \Rclass{SimpleRleViewsList}
  object. This sounds even more complicated but here again it helps
  to think of it as a named list of \Rclass{RleViews} objects.
  The names of the list are still the chromosomes and each top-level element
  in the list now is an \Rclass{RleViews} object i.e. a set of views on
  the coverage vector for the corresponding chromosome.

  \item Use the generic \textit{coercion} function \Rfunction{as} to turn
  the result of the previous slicing into a \Robject{GRanges} object.

  \end{itemize}
\end{Exercise}

There is an issue with the above method: we've lost the strand information.
We might want to retain it if our final goal were to identify new
transcripts. In the next exercise we improve our ``pile up and slice''
method to propagate the strand information stored in \Robject{hits0}.

\begin{Exercise}
  \begin{itemize}

  \item Split \Robject{hits0} in two \Rclass{GRanges} object: one
  containing the hits located on the plus strand and one containing
  the hits located on the minus strand.

  \item Apply the ``pile up and slice'' method used in the previous
  exercise to each \Rclass{GRanges} object. This transforms each
  \Rclass{GRanges} object into another \Rclass{GRanges} object.

  \item Note that the two \Rclass{GRanges} objects obtained previously
  are unstranded. Add the strand information to them.

  \item Combine the two \Rclass{GRanges} objects obtained previously
  with the generic function \Rfunction{c}.

  \item One last thing we might want to do is use \Rfunction{reduce}
  on the final result. The only effect of reducing here is to reorder
  the regions first by chromosome and then by strand (this is how
  \Rfunction{reduce} orders the elements of a GRanges object).

  \end{itemize}
\end{Exercise}

Finally, to make it easier to repeat the transformation we've done
to \Robject{hits0} in the previous exercise but now with different
values of the threshold used in the slicing step, we want to wrap
our code in a function.

\begin{Exercise}
  \begin{itemize}

  \item Write a function (\Rfunction{coveredRegions}) that takes a
  set of hits (in a \Rclass{GRanges} object) and a threshold (\Rcode{lower})
  and returns a \Rclass{GRanges} object containing the regions where the
  coverage of the hits is greater than (or equal to) the threshold.
  Make sure the returned \Rclass{GRanges} object is stranded and reduced.

  \item Sanity check: compare \Rcode{coveredRegions(hits0, 1)} with
  \Rcode{reduce(hits0)}.

  \item Extract the DNA sequences corresponding to the regions returned
  by \Rcode{coveredRegions(hits0, 10)}. Try with other threshold
  values.

  \end{itemize}
\end{Exercise}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Use case IV: Measuring the complexity of the reads}

A common task when dealing with HTS data is to filter out reads with a low
complexity like poly-As or reads made of the repetition of the same 2-mer
(dinucleotide) etc... In this section, we implement a simple function that
takes a \Rclass{DNAStringSet} object (the reads) and returns a score for each
read based on its complexity. The approach we use is inspired and adapted
from the DUST algorithm.

The basic idea behind DUST is that a DNA sequence with a ``poor trinucleotide
content'' (i.e. with a small number of \textit{distinct} trinucleotides) is
considered to have a low complexity. For example, the following sequences
have a low complexity:

\begin{itemize}
  \item AAAAAAAAAA: contains only 1 tri-nculeotide: AAA
  \item ATATATATATATATA: contains 2 tri-nucleotide: ATA and TAT
\end{itemize}

On the other hand, a DNA sequence with a ``rich trinucleotide
content'' (i.e. many \textit{distinct} trinucleotides) is considered
to have a high complexity. For example, the following 36-mer:

\begin{itemize}
  \item GGGCTACATGACGGTCCTGTATTTAGCCAGAGGATC
\end{itemize}

has the highest complexity a 36-mer can have because all the
trinucleotides contained in it (34 in total) are distinct.

Here is how we will compute the score of a given DNA sequence:
\begin{itemize}
  \item Count the number of occurences (frequency) of each possible
  trinucleotide:
  $F_{AAA}$, $F_{AAC}$, $F_{AAG}$, $F_{AAT}$, $F_{ACA}$, ..., $F_{TTT}$
  (64 in total).
  \item Substract $1$ to the frequencies that are not zero. That is:
  for each $F_{xxx}$, if $F_{xxx} \geq 1$ then $F_{xxx} = F_{xxx} - 1$.
  \item $score = 1 - \frac{\sqrt{\sum{F_{xxx}^2}}}{L - 3}$ where $L$
  is the length of the sequence.
\end{itemize}

A score of 0 indicates a poly-A, poly-C, poly-G or poly-T.
A score of 1 can only be obtained by a sequence where all trinucleotides
are distinct which implies that the sequence is no more than 66 bases long.
Sequences longer than this cannot obtain a score of 1 and this scoring
algorithm is not expected to give meaningful results on long sequences.
Also, the score of very short sequences (i.e. 4 <= length <= 6) is not
very meaningful either.

The \Biocpkg{SeattleIntro2010} package contains reads from the Nagalakshmi
et al. experiment stored in a FASTQ file (unaligned reads).
Let's start by loading and ``cleaning'' them.

\begin{Exercise}
  \begin{itemize}

  \item Use the \Rfunction{read.DNAStringSet} function from the
  \Biocpkg{Biostrings} package to load the \code{SRR002051.reads1-50k.fastq}
  file located in the \code{extdata} subfolder of the
  \Biocpkg{SeattleIntro2010} package.
  Note that, by default, \Rfunction{read.DNAStringSet} expects a FASTA file.
  Consult the man page for \Rfunction{read.DNAStringSet} to see how to
  read a FASTQ file.

  \item Use the \Rfunction{alphabetFrequency} function (with
  \Rcode{collapse=TRUE}) on those reads.

  \item The scoring algorithm described above doesn't apply to reads
  that contain IUPAC ambiguity letters. Remove those reads.
  Tip: Try to use \Rfunction{alphabetFrequency} with \Rcode{baseOnly=TRUE}.
%Sol: cleanreads <- reads[alphabetFrequency(reads, baseOnly=TRUE)[ , "other"] == 0]

  \end{itemize}
\end{Exercise}

Let's call \Robject{cleanreads} the \Rclass{DNAStringSet} object
containing the result of the previous exercise.
Now we are going to compute the score of those reads.

\begin{Exercise}
  \begin{itemize}

  \item Use the \Rfunction{trinucleotideFrequency} function from the
  \Biocpkg{Biostrings} package to compute the trinucleotide frequencies
  of the clean reads. The result (let's call it \Robject{tnf}) is a matrix
  containing the trinucleotide counts for each input read. The matrix
  has 1 row per read and 64 columns (i.e. one column per each possible
  trinucleotide).

  \item Use the \Robject{tnf} matrix to compute the scores of the clean
  reads by following the steps described above.
  Tip 1: Note that, by performing vectorized arithmetic (a key feature of
  arithmetic operations in \R{}), we can avoid the use of loops and be
  very fast.
  Tip 2: You can use \Rfunction{rowSums} on a numeric matrix to sum all the
  coefficients that belong to the same row, and this for all the rows.

  \item Wrap the previous code in a function i.e. implement the
  \Rfunction{complexityOfReads} function that computes the complexity
  scores of each element of a \Rclass{DNAStringSet} object.

  \item Use this function to display the clean reads with a score of 0,
  then the reads with a score of 1.

  \item Plot the histogram of the complexities of the reads in
  \Robject{cleanreads}.

  \end{itemize}
\end{Exercise}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Use case V: Pattern matching}

In this last use case, we're going to try to align our reads with the
pattern matching tools available in the \Biocpkg{Biostrings} package.

We can divide those tools in 2 categories:

\begin{itemize}

  \item General purpose pattern matching/aligning tools:
  \begin{itemize}

    \item \Rfunction{matchPattern} (and convenience wrapper
    \Rfunction{countPattern});

    \item \Rfunction{vmatchPattern} (and convenience wrapper
    \Rfunction{vcountPattern});

    \item \Rfunction{pairwiseAlignment}.

  \end{itemize}

  \item Specialized pattern matching/finding tools:
  \begin{itemize}

    \item \Rfunction{matchPDict} (and convenience wrappers
    \Rfunction{countPDict}, \Rfunction{whichPDict});

    \item \Rfunction{vmatchPDict} (and convenience wrappers
    \Rfunction{vcountPDict}, \Rfunction{vwhichPDict});

    \item \Rfunction{matchLRPatterns};

    \item \Rfunction{matchProbePair};

    \item \Rfunction{trimLRPatterns};

    \item \Rfunction{matchPWM};

    \item \Rfunction{findPalindromes}.

  \end{itemize}

\end{itemize}

The \Rfunction{matchPDict} function (and family) is the tool of choice
when it comes to matching a high number of short DNA sequences (the
patterns) against a long reference DNA sequence (the subject).
The \Rfunction{countPDict} and \Rfunction{whichPDict} functions are
convenience wrappers for situations where counting the hits only is
desired (\Rfunction{countPDict}), or when we just want to know which
patterns have at least 1 hit (\Rfunction{whichPDict}).
A common feature of all the \Rfunction{*PDict} matching functions is
that, in order to be efficient, they all require to preprocess the
set of patterns (aka ``the pattern dictionary''). This is achieved with
the \Rfunction{PDict} constructor:

<<importance-of-preprocessing-setup, echo=FALSE>>=
library(SeattleIntro2010)
library(Biostrings)
filepath <- system.file("extdata", "SRR002051.reads1-50k.fastq", package="SeattleIntro2010")
reads <- read.DNAStringSet(filepath, format="fastq")
cleanreads <- reads[alphabetFrequency(reads, baseOnly=TRUE)[ , "other"] == 0L]
library(BSgenome.Scerevisiae.UCSC.sacCer2)
@

<<importance-of-preprocessing>>=
system.time(nhits <- countPDict(cleanreads, Scerevisiae$chrI))
system.time(nhits2 <- countPDict(PDict(cleanreads), Scerevisiae$chrI))
identical(nhits, nhits2)
table(nhits)
@

Note that, by default, \Rfunction{matchPDict} and family look for
exact matches only. See \Rcode{?PDict} for how to allow a small numbers
of mismatches (at the price of a significant slowdown though).

Let's start by illustrating the basic usage of
\Rfunction{PDict}/\Rfunction{countPDict}.

\begin{Exercise}
  \begin{itemize}

  \item Preprocess the \Robject{cleanreads} dictionary with \Rfunction{PDict}.

  \item Use \Rfunction{countPDict} to count the number of hits per pattern
  against sacCer2 chromosome IV (the longest chromosome in Yeast).

  \end{itemize}
\end{Exercise}

Note that, without taking any special precaution, the hit counts we
got are for hits on the plus strand of chrIV only. However, a simple
and efficient way to count the hits on the minus strand is to repeat
the 2 steps above on the reverse complement of \Robject{cleanreads}.

\begin{Exercise}
  \begin{itemize}

  \item Use the \Rfunction{reverseComplement} function on
  \Robject{cleanreads}.

  \item Count the hits on the minus strand of chrIV.

  \item Finally add the 2 vectors of hit counts you got.

  \end{itemize}
\end{Exercise}

Note that one way to make this faster (at the cost of more memory usage)
is to preprocess the patterns and their reverse complements in the same
\Rclass{PDict} object, but this goes beyond the scope of this simple
introduction to fast pattern matching with \Biocpkg{Biostrings}.

In the next exercise, we want to match our clean reads against the
\textit{full} sacCer2 genome, by searching the plus and minus strands
of each reference sequence. Again, we will do exact matching only.

\begin{Exercise}
  \begin{itemize}

  \item Write a function \Rfunction{countSingleStrandGenomeHits}, that
  takes 2 arguments: (1) a pattern dictionary (\Rclass{DNAStringSet})
  and, (2) a \Rclass{BSgenome object} (e.g. \Robject{Scerevisiae}).
  The function will return the number of hits on the entire genome
  for each pattern in the input dictionary. The function should only
  count the hits on the plus strand.
  \begin{itemize}
     \item Tip 1: Use \Rfunction{seqnames} on a \Rclass{BSgenome} object
     to get its sequence names.
     \item Tip 2: Loop over the sequence names with a \Rcode{for} loop.
     \item Tip 3: Before you start the loop, create an integer vector of
     the length of the input dictionary with
     \Rcode{nhits <- integer(length(dict))}. It will be initialized with
     zeroes. Then, inside the loop, add the results for each sequence
     to \Robject{nhits}.
  \end{itemize}

  \item Use \Rfunction{countSingleStrandGenomeHits} to count the hits
  of our clean reads against the full sacCer2 genome. The final result
  (let's call it \Rcode{gwnhits} for ``genome wide nhits'') must take
  both strands into account.

  \item What's the percentage of reads that don't hit the genome?
  How could this be improved?

  \end{itemize}
\end{Exercise}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Solutions and session information}

\subsection{Exercise 1}
\begin{verbatim}
library(SeattleIntro2010)
filepath <- system.file("extdata", "sacCer2_sgdGene.sqlite",
                        package="SeattleIntro2010")
library(GenomicFeatures)
txdb <- loadFeatures(filepath)
\end{verbatim}

\subsection{Exercise 2}
\begin{verbatim}
ex <- exons(txdb)
length(ex)  # number of exons
width(ex)  # lengths of exons
max(width(ex))  # length of longest exon
ex[which.max(width(ex))]  # longest exon
ex[seqnames(ex) == "2micron"]  # exons on the 2micron plasmid
library(BSgenome.Scerevisiae.UCSC.sacCer2)
Scerevisiae  # displays Scerevisiae
seqlengths(Scerevisiae)
Scerevisiae[["chrI"]]
getSeq(Scerevisiae, ex, as.character=FALSE)
\end{verbatim}

\subsection{Exercise 3}
\begin{verbatim}
filepath <- system.file("extdata", "SRR002051.chrI-V.bam",
                        package="SeattleIntro2010")
galn <- readGappedAlignments(filepath)
class(galn)
length(galn)
\end{verbatim}

\subsection{Exercise 4}
\begin{verbatim}
table(qwidth(galn))  # one single value (33)
table(width(galn))  # most are 33 but not all!
any(ngap(galn) != 0)  # no gaps
cigtable <- cigarOpTable(cigar(galn))
dim(cigtable)  # big matrix! (many rows)
colnames(cigtable)  # M is the first column
idx <- which(rowSums(cigtable[ , -1]) != 0)  # final index
galn[idx]
\end{verbatim}

\subsection{Exercise 5}
\begin{verbatim}
hits <- grg(galn)
\end{verbatim}

\subsection{Exercise 6}
\begin{verbatim}
hits0 <- hits[countOverlaps(hits, exons(txdb)) == 0]
reduce(hits0)
intersect(reduce(hits0), exons(txdb))  # sanity check
\end{verbatim}

\subsection{Exercise 7}
\begin{verbatim}
cov0 <- coverage(hits0)  # SimpleRleList object
cov0[["chrV"]]  # Rle object
slice(cov0, lower=200)  # SimpleRleViewsList object
as(slice(cov0, lower=200), "GRanges")
\end{verbatim}

\subsection{Exercise 8}
\begin{verbatim}
plushits0 <- hits0[strand(hits0) == "+"]
minushits0 <- hits0[strand(hits0) == "-"]
plusregions0 <- as(slice(coverage(plushits0), lower=200), "GRanges")
minusregions0 <- as(slice(coverage(minushits0), lower=200), "GRanges")
strand(plusregions0) <- "+"
strand(minusregions0) <- "-"
c(plusregions0, minusregions0)
reduce(c(plusregions0, minusregions0))
\end{verbatim}

\subsection{Exercise 9}
\begin{verbatim}
coveredRegions <- function(hits, lower)
{
  plushits <- hits[strand(hits) == "+"]
  minushits <- hits0[strand(hits) == "-"]
  plusregions <- as(slice(coverage(plushits), lower=lower), "GRanges")
  minusregions <- as(slice(coverage(minushits), lower=lower), "GRanges")
  strand(plusregions) <- "+"
  strand(minusregions) <- "-"
  reduce(c(plusregions, minusregions))
}
identical(coveredRegions(hits0, 1), reduce(hits0))  # TRUE
getSeq(Scerevisiae, coveredRegions(hits0, 10), as.character=FALSE)
getSeq(Scerevisiae, coveredRegions(hits0, 1000), as.character=FALSE)
\end{verbatim}

\subsection{Exercise 10}
\begin{verbatim}
filepath <- system.file("extdata", "SRR002051.reads1-50k.fastq",
                        package="SeattleIntro2010")
reads <- read.DNAStringSet(filepath, format="fastq")
alphabetFrequency(reads, collapse=TRUE)
cleanreads <- reads[alphabetFrequency(reads, baseOnly=TRUE)[ , "other"] == 0]
\end{verbatim}

\subsection{Exercise 11}
\begin{verbatim}
tnf <- trinucleotideFrequency(cleanreads)

tnf <- tnf - 1
tnf[tnf < 0] <- 0
scores <- 1 - sqrt(rowSums(tnf^2)) / (width(cleanreads) - 3)

complexityOfReads <- function(reads)
{
    tnf <- trinucleotideFrequency(reads)
    tnf <- tnf - 1
    tnf[tnf < 0] <- 0
    1 - sqrt(rowSums(tnf^2)) / (width(reads) - 3)
}

cleanreads[complexityOfReads(cleanreads) == 0]
cleanreads[complexityOfReads(cleanreads) == 1]

plot(hist(complexityOfReads(cleanreads), breaks=50))
\end{verbatim}

\subsection{Exercise 12}
\begin{verbatim}
pdict <- PDict(cleanreads)
library(BSgenome.Scerevisiae.UCSC.sacCer2)
nhits <- countPDict(pdict, Scerevisiae[["chrIV"]])
\end{verbatim}

\subsection{Exercise 13}
\begin{verbatim}
rc_cleanreads <- reverseComplement(cleanreads)
rc_pdict <- PDict(rc_cleanreads)
rc_nhits <- countPDict(rc_pdict, Scerevisiae[["chrIV"]])
nhits + rc_nhits
\end{verbatim}

\subsection{Exercise 14}
\begin{verbatim}
countSingleStrandGenomeHits <- function(dict, genome)
{
  pdict <- PDict(dict)
  nhits <- integer(length(dict))
  for (seqname in seqnames(genome)) {
    nhits <- nhits + countPDict(pdict, genome[[seqname]])
  }
  nhits
}
plus_nhits <- countSingleStrandGenomeHits(cleanreads, Scerevisiae)
minus_nhits <- countSingleStrandGenomeHits(rc_cleanreads, Scerevisiae)
gwnhits <- plus_nhits + minus_nhits
sum(gwnhits == 0) / length(gwnhits)  # 80% of the reads have no hit!
\end{verbatim}

\subsection{Session information}

<<sessionInfo,results=tex,echo=FALSE>>=
toLatex(sessionInfo())
@ 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibliographystyle{plain}
\bibliography{SeattleIntro2010}

\end{document}
