%\VignetteDepends{ShortRead, rtracklayer, GenomicFeatures}
%\VignetteIndexEntry{Sequence Analysis : Lab}
%\VignettePackage{SeattleIntro2010}

\documentclass[a4paper]{article}

\usepackage{BioconductorVignette}
\newcommand\ELAND{\software{ELAND}}
\newcommand\MAQ{\software{MAQ}}
\newcommand\bam{\software{BAM}}
\newcommand\Bowtie{\software{Bowtie}}

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

\title{Introduction to Sequences : Reading and Manipulating Short Reads}
\author{Valerie Obenchain and Martin Morgan \\ Fred Hutchinson Cancer Research Center}
\date{9-10 December, 2010}

\begin{document}

<<setup, echo=FALSE>>=
options(width = 60)
olocale=Sys.setlocale(locale="C")
@ 

\maketitle

\tableofcontents

\section{Introduction}

This practical uses the \Biocpkg{ShortRead} package to input aligned
and other short read data files, and illustrates some of the available
sequence manipulation and quality assessment tools.  Activities during
the lab are posed as exercises. So, as a first exercise:
\begin{Exercise}
  Start an \R{} session, and load the \Biocpkg{ShortRead}
  package. 
<<ShortReadVersion>>=
library(ShortRead)
packageDescription("ShortRead")$Version
@ 
%% 
  Confirm that the version of your package is at least as recent as
  the version in this document. Seek assistance from one of the
  course assistants if you need help getting the current version of
  \Biocpkg{ShortRead}.
\end{Exercise}
%%
The course also requires access to sample data. We will be using yeast
data from an RNA-seq experiment. See
\Rcode{browseVignette("SeattleIntro2010")} for experiment details.
\begin{Exercise}
  Load the \Rpackage{SeattleIntro2010} package. Use
<<embl2010, eval=FALSE>>=
library(SeattleIntro2010)
@ 
\end{Exercise}

\section{Aligned read input}

This section illustrates input of aligned reads. It focuses on aligned
reads in \bam{} format; 
reading data produced by software such as \MAQ{} or \Bowtie{}
is described in the \Biocpkg{ShortRead} `Overview'
vignette and on the \Rfunction{readAligned} help page.

\subsection{Preliminaries}

Vendor and third-party software is likely to process raw images, base
calling (\Biocpkg{Rolexa} is a \Bioconductor{} package providing
alternative base calling), and alignment to a reference genome (see
\Biocpkg{BSgenome} and the \Rfunction{matchPDict} function of
\Biocpkg{Biostrings}). The following work flow starts with reads aligned
with \software{BWA} software in \bam{} format.

We'll start by creating a variable, \Rcode{extdataDir}, containing the
directory holding the sample data, and check to make sure that we
have defined the right path.
<<extdataDir>>=
extdataDir <- 
    system.file("extdata", package="SeattleIntro2010")
file.exists(extdataDir)
stopifnot(file.exists(extdataDir))
@ 

A key functionality provided by \Biocpkg{ShortRead} is input of a
diversity of file types, both of files from the Illumina pipeline and
from other software and in formats appropriate for other
technologies. The interface to these input functions is meant to
facilitate reading one or more files into a single object. The
interface is like that for \Rfunction{list.files}: provide a directory
path where relevant files are to be found, plus a regular expression
to select the files you are interested in.
\begin{Exercise}
  Use \Rfunction{list.files} to display all files in the
  \Rcode{extdataDir} directory.
<<filePattern>>=
list.files(extdataDir)
@ 
\end{Exercise}
We'll select just one file to use in this analysis, based on files
matching the regular expression \texttt{SRR002051.chrI-V.bam\$}. 
The use of the special character `\$' ensures that we don't also match
the \bam{} index file (extention \texttt{.bai}) in the same directory. 
Other inputs (e.g., Solexa, \MAQ, \Bowtie, or \software{BWA}
files) are also supported.
\begin{Exercise}
  We'll use an abbreviated file for most parts of this lab. The
  abbreviated file contains chromosomes 1 through 5 from a single lane
  of an Illumina run. The name of the file is
  \texttt{SRR002051.chrI-V.bam}.  We will match this `pattern'
  and check that we have specified the pattern appropriately.
<<pattern>>=
pattern <- "SRR002051.chrI-V.bam$"
list.files(extdataDir, pattern)
@ 
%%
\end{Exercise}

\subsection{\Rfunction{readAligned} and the \Rclass{AlignedRead} class}

The \Rfunction{readAligned} function can be used to input aligned
reads. The first argument is a directory path where alignment files
are to be found. The second argument is the regular expression to
select files to be read. An optional third argument allows the user to
specify which type of file is to be read in.
\begin{Exercise}
  Use \Rfunction{readAligned} to read in the data, 
  specifying \texttt{BAM} as the type of file.
<<readAligned>>=
aln <- readAligned(extdataDir, pattern, type="BAM")
@ 
%%
  See the help page for \Rfunction{readAligned} for additional
  details about supported file types.
\end{Exercise}

What does \Rfunction{readAligned} input? It inputs the short read
sequences and base call qualities, and the chromosome, position, and
strand information associated with short read alignments. This
information is expected to be provided by all short read alignemnt
software.
\begin{Exercise}
  Display the object we have read in.
<<aln>>=
aln
@ 
%
  There are \Sexpr{length(aln)} reads in the object, each read
  consisting of \Sexpr{width(aln)} nucleotides. View the first several
  reads and query information about, e.g., the number of reads that
  align to each strand, or to each chromosome.
<<aln-head>>=
head(sread(aln))
table(strand(aln), useNA="ifany")
table(chromosome(aln), useNA="ifany")
@ 
%%
  Notice how elements of the \Robject{aln} object are extracted using
  \emph{accessors} such as \Rfunction{sread} and \Rfunction{strand};
  these are described on the help page for the class of the
  \Robject{aln} object (indicated in the display of \Robject{aln},
  above, as class \Rclass{AlignedRead}); note that the help page
  refers to the help page for \Rfunction{accessors} to enumerate
  additional ways of accessing the data.
\end{Exercise}

The \Rfunction{strand} function returns a factor
with three levels. The first two describe reads aligned to the plus
and minus strands, the third (\Rcode{*}) is available for successful
alignments where strand information is irrelevant.
This \bam{} file was created by filtering out reads that did not
align to chromosomes I-V so there are no \Rcode{NA} values in the 
data. The approximately equal distribution of reads between 
the positive and negative strands indicate there is no bias. 

Aligned reads contain several different kinds of information about
`quality'. Individual bases are assessed for quality during base
calling. These `raw' base qualities are `calibrated' during 
alignment. The alignments themselves also have qualities
associated with them, with the details of alignment quality differing
between alignment algorithms.
 
\subsection{Base calls}
\begin{Exercise}
  Use \Rfunction{alphabetFrequency} to summarize nucleotide use in the
  short reads in \Robject{aln}.
<<alphabetFreqeuncy>>=
alphabetFrequency(sread(aln), collapse=TRUE, baseOnly=TRUE, 
                  as.prob=TRUE)
@

  The frequency of `other' (i.e., uncalled) nucleotides is
  typically $<1\%$.
\end{Exercise}
  The expectation is that sequences start at essentially arbitrary
locations in the sequenced DNA so nucloetide frequencies should be
approximately constant over the cycles (i.e., an A is expected as 
frequently at the beginning of the sequence as at the end).
\begin{Exercise}
  Use \Rfunction{alphabetByCycle} to extract the number of each
  nucleotide A, C, T, G, sequenced at each cycle. Plot the transpose
  of the matrix using \Rfunction{matplot}, to see how nucleotide
  counts change across cycles.
<<abc-lab-plot,fig=TRUE,prefix=FALSE,height=4>>=
abc <- alphabetByCycle(sread(aln))
dim(abc)
abc[1:4,1:5]
matplot(t(abc[1:4,]), type="l", lty=c(1:4), 
col=c("green","red","black","blue"), xlab="Cycle", ylab="Count")
legend("top", dimnames(abc)$alphabet[1:4], lty=c(1:4), 
col=c("green","red","black","blue"), bty="n", horiz=T)
@ 
\end{Exercise}
\begin{figure}[ht]
  \centering
  \includegraphics[width=.95\textwidth]{abc-lab-plot.pdf}
  \caption{Nucleotide frequency per cycle}  
  \label{fig:abc}
\end{figure}
There are a number of possible contributors to cycle-dependent
nucleotide frequencies. Early experiments such as this frequently
suffered from inadequate reagent volume and nucleotide-specific
differential accumulation of fluorescent dyes.

\begin{Exercise}
  Retrieve calibrated base quality from \Robject{aln}. 
<<quality>>=
head(quality(aln))
@ 
% 
  These qualities are string-encoded $-10 \log_{10}$
  probabilities. The encoding in this case follows a convention
  established by Illumina. The details of the encoding can be obtained
  by querying \Rcode{quality(aln)} for its \Rfunction{alphabet}; the
  letter \Rcode{A} corresponds to a $-10 log_{10}$ score of 1. 

  Numeric values are readily retrieved as a matrix, with rows
  corresponding to reads and columns to cycles.  Computations can then
  be performed on them, e.g., to determine average calibrated quality
  scores as a function of cycle.
<<quality-scores>>=
alf <- alphabet(quality(aln))
m <- as(quality(aln), "matrix")
colMeans(m)
@ 
\end{Exercise}
 


\begin{Exercise}
  It may be of interest to identify the GC-rich reads in our sample.
<<gcRegions, fig=TRUE, prefix=FALSE, height=4>>=
res <- alphabetFrequency(sread(aln), baseOnly=TRUE)
head(res)
gcContent <- rowSums(res[,c("G", "C")]) / rowSums(res)
print(histogram(gcContent, xlab="Percent GC", 
                ylab="Percent of total reads"))
gcIdx <- gcContent > 0.6
aln[gcIdx]
@
\begin{figure}[ht]
  \centering
  \includegraphics[width=.95\textwidth]{gcRegions}
  \caption{GC-rich reads}
  \label{fig:gcrich}
\end{figure}
\end{Exercise}

\subsection{Alignment}

Alignment qualities are accessible with \Rfunction{alignQuality}. This
returns an object that can contain quality scores in different
formats; to extract the actual quality scores, use
\Rfunction{quality}. Reads failing to align or to align in multiple
locations have an alignment quality of 0.
\begin{Exercise}
  Retrieve the alignment quality scores, determine how many align
  poorly, and visualize the distribution (Figure~\ref{fig:alignqual})
  of scores.
<<alignQuality,fig=TRUE,include=FALSE,prefix=FALSE,height=4>>=
alignQuality(aln)
q <- quality(alignQuality(aln))
sum(q==0)
print(densityplot(q[q>1], plot.points=FALSE, 
                  xlab="Alignment quality"))
@ 
% 
\begin{figure}[ht]
  \centering
  \includegraphics[width=.95\textwidth]{alignQuality}
  \caption{Alignment quality}
  \label{fig:alignqual}
\end{figure}
\end{Exercise}

Alignment algorithms produce information in addition to basic data
about chromosome, position, and strand alignment. The exact content
varies between algorithms, and is available with
\Rfunction{alignData}. \Rfunction{alignData} returns an
\Rclass{AlignedDataFrame} object that contains these data and a
metadata description of them. 

\subsection{Subsets and filters}

A very common operation is to reduce the number of reads used for
subsequent analysis. This can be done in a coordinated fashion by
creating a subset of \Robject{aln}.
\begin{Exercise}
  Select just the aligned reads from chromosome III and V. 
<<select>>=
chrfilt <- chromosomeFilter("chrIII|chrV") 
aln[chrfilt(aln)]
@ 
\end{Exercise}


The filter functions built-in to \Biocpkg{ShortRead} use a `factory'
pattern to create instances of each filter that `remember' how the
filters were created. For instance, \Rcode{chromosomeFilter("chrV")}
creates an instance of the chromosome filter to select only
chromosomes matching \Rcode{chrV}.
\begin{Exercise}
  As an advanced example, the following filter subsamples a
  (user-specified) number of reads. The \Rfunction{samplingFilter}
  function uses the factory pattern, so filters created with it
  remember how many reads to sample.
<<sampleFunction>>=
samplingFilter <- function(sampleSize) {
    srFilter(function(x) {
        idx <- seq_len(length(x))
        idx %in% sample(idx, sampleSize)
    }, name="Demo sampling filter")
}
sample100 <- samplingFilter(100)
aln[sample100(aln)]
@
\end{Exercise}


\subsection{Cautions}

There are several confusing areas associated with reading data aligned
with various software packages. (1) Some alignment programs and genome
resources start numbering nucleotides of the subject sequence at 0,
whereas others start at 1. (2) Some alignment programs report matches
on the minus strand in terms of the `left-most' position of the read
(i.e., the location of the 3' end of the aligned read), whereas other
report `five-prime' matches (i.e., in terms of the 5' end of the
read), regardless of whether the alignment is on the plus or minus
strand. (3) Some alignment programs reverse complement the sequence of
reads aligned to the minus strand. (4) Base qualities are sometimes
encoded as character strings, but the encoding differs between `fastq'
and `solexa fastq'. It seems that all combinations of these choices
are common `in the wild'.

The help page for \Rfunction{readAligned} attempts to be explicit
about how reads are formatted. Briefly:
\begin{itemize}
\item Subject sequence nucleotides are numbered starting at 1, rather
  than zero. \Rfunction{readAligned} adjusts the coordinate system of
  input reads if necessary (e.g., when reading MAQ alignments).
\item ELAND and Bowtie alignments on the minus strand are reported in
  `left-most' coordinates systems.
\item ELAND and Bowtie alignments on the minus strand are not reverse
  complemented.
\end{itemize}
Alignment programs sometimes offer the opportunity to customize
output; such customization needs to be accommodated when reads are
input using \Biocpkg{ShortRead}.


\section{Quality assessment} 

This part of the course addresses \Biocpkg{ShortRead} facilities for
assessing quality, primarily of \bam{} data. The \Biocpkg{ShortRead}
functionality is mean to complement rather than replace QA tools
provided by manufacturer pipelines.

The data we are using is from an RNA-seq experiment, using either 
random hexamer or oligo(dT) primers to isolate RNA. The data consist 
of six lanes of an early Solexa run, three from random hexamer, three from 
oligo(dT). In each group of three, there is an original, biological replicate 
and technical replicate (of original). 

\subsection{Generating a QA report}

Creating a QA report is a two-step process. The first step is to visit
necessary files to collate information in a compact
representation. The second step is to present the information in a
useful format.

The \Rfunction{qa} function collates information for the QA report. It
visits each \_fastq.bam file, and extracts information on reads
and their qualities. Evaluation of the function is straight-forward, e.g.,
\verb|qa <- qa(extdataDir, pattern, type="BAM")|. The
process of collating files can be time consuming (each export file
must be parsed, taking 3-4 minutes per file) and memory intensive
(lanes are processed independently of one another, but processing a full lane
consumes 2-3 GB of memory). The return value of the \Rfunction{qa}
function is actually quite compact, and easy to work with.
\begin{Exercise}
  Rather than collating information during the lab, we load the data
  from a previously stored instance. The components of the qa object
  displayed with names(qa) give a sense of what has been collated
  during the call to qa but are primarily for internal use.
<<qa-input>>=
data("qa", package="SeattleIntro2010")
names(qa)
@ 
\end{Exercise}
One feature of \Biocpkg{ShortRead} that can speed up this stage of the
operation is the use of clustered computer resources and the
\Rpkg{Rmpi} or \Rpkg{mullticore} package; \Rfunction{qa} uses the
\Rfunction{srapply} function to automatically detect and distribute
collation tasks across pre-established nodes. This is outlined in more
detail in the \Biocpkg{ShortRead} Vignette. 

Primer and replicate type information can be seen in the hitspergene data file included
with the course package.
<<qa-laneinfo>>=
data("hitspergene", package="SeattleIntro2010")
elementMetadata(hitspergene)
@


The QA information collated from all six lanes is
summarized into an html report using the \Rfunction{report}
function. The command to create the report is
\verb|rpt <- report(qa, dest=tempfile())|.  This creates an HTML file
at the location specified by the argument \Rfunarg{dest}.
\begin{Exercise}
  Create a report from the \Robject{qa} object loaded in the previous
  exercise. Do this with the command
<<qa-report>>=
rpt <- report(qa, dest=tempfile())
@ 
%% 
View the report in your browser with
<<qa-browse, eval=FALSE>>=
browseURL(rpt) 
@ 
\end{Exercise}
The QA report provides summary statistics about the numbers of reads
and alignments, base calls and qualities, characteristics of per-lane
and per-tile read quality, and other information. The QA report is
self-documenting, providing a narrative description of each section.
Take some time to look over the report and observe the behavior in
the `Cycle-specific base calls and read quality' section.


\section{Session information}

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

\end{document}
