\documentclass{article}

\usepackage[colorlinks=true,linkcolor=blue]{hyperref}
\usepackage{theorem}
\usepackage{underscore}
\usepackage{color}

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

\newcommand{\code}[1]{{\texttt{#1}}}
\newcommand{\file}[1]{{\texttt{#1}}}

\newcommand{\software}[1]{\textsf{#1}}
\newcommand{\R}{\software{R}}
\newcommand{\bioc}{\software{BioConductor}}
\newcommand{\ELAND}{\software{ELAND}}
\newcommand{\MAQ}{\software{MAQ}}
\newcommand{\Bowtie}{\software{Bowtie}}

%% Excercises and Questions
\theoremstyle{break} \newtheorem{Ex}{Exercise}
\theoremstyle{break} \newtheorem{Q}{Question}
%% And solution or answer
\newenvironment{solution}{\color{blue}}{\bigskip}

\begin{document}

\SweaveOpts{keep.source=TRUE,width=9,height=8,eps=FALSE, include=FALSE, prefix.string=docs/IOQAFig}

%% \setkeys{Gin}{width=0.95\textwidth}

\title{I/0 and Quality Assessment using the \Rpackage{ShortRead} package}
\date{May 31, 2009}
\maketitle

\tableofcontents

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

\section{Introduction}

This portion of the course uses the \Rpackage{ShortRead} package to
input aligned and other short read data files, and illustrates some of
the available Solexa-based quality assessment tools.  Activities
during the lab are posed as exercises. So, as a first exercise:
\begin{Ex}
  Start an \R{} session, and load the \Rpackage{ShortRead}
  package. 
<<ShortRead,results=hide,echo=FALSE>>=
library("ShortRead")
<<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
  \Rpackage{ShortRead}.
\end{Ex}
%% 
The course also requires access to sample data. 
\begin{Ex}
  Copy the data from the distribution media to your local hard
  drive. In \R{} change the working directory to point to the data
  location (which should contain the folders \code{extdata},
  \code{scripts}, etc.), along the lines of
<<setwd, eval=FALSE>>=
setwd("c:/Documents and Settings/Desktop/Course/labs")
@ 
<<setwd, echo=FALSE>>=
oldwd = setwd("..")
@ 
%% 
  and confirm that the files have been copied correctly.
\end{Ex}

\section{Aligned read input}

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

\subsection{\Rclass{SolexaPath}: navigating Solexa output}

This section introduces a way to conveniently navigate the hierarchy
of files produced by \ELAND{}; this simplifies subsequent activities,
but the full \ELAND{} output is not required to use
\Rpackage{ShortRead}.

Solexa software processes data in a pipeline. Raw images are extracted
to image intensity files by the Firecrest software component. Image
intensity files are summarized as base calls using the Bustard base
caller. Subsequent analysis is performed by a diversity of software
components called Gerald; one example of a Gerald program is \ELAND{},
a whole-genome aligner.

The pipepline provides (or can be configured to provide -- the people
running the machine have quite a bit of control over this) exquisite
detail about each stage of the process. For instance, Firecrest scans
each lane of a flow cell as 300 \emph{tiles}, arranged in three
serpentine columns. The Firecrest output is summarized in two
different file types for each tile, so there are $2 \times 300 \times
8 = 4800$ files produced by Firecrest alone.

A portion of a file hierarchy is provided as course data.
\begin{Ex}
  Consult the help page for \code{SolexaPath}, and create an instance
  of this object, e.g.,
<<SolexaPath>>=
sp <- SolexaPath("extdata/ELAND/080828_HWI-EAS88_0003")
@ 
\end{Ex}
This command scans the file system rooted at the specified path (the
final directory name given above is a typical top level name for a
Solexa run, encoding the date and machine used, for instance),
identifying likely paths associated with each stage of the pipeline.
\begin{Ex}
  Display this object, and query it for the paths were Gerald analysis
  results are stored.
<<sp-display>>=
sp
analysisPath(sp)
@
% 
  There are two analysis paths, because the data were generated by two
  separate runs of the \ELAND{} software. The analysis paths are
  nested inside \code{baseCallPath(sp)} as a Solexa convention, to
  indicate that the Gerald analyses both use the results of the same
  application of the base calling software.
\end{Ex}
Note how the display of \Robject{sp} is compact, and the names in the
display hint at how to navigate the object. Short read objects can be
very large. Most objects in \Rpackage{ShortRead} and related packages
have been designed to display a summary, rather than the `big'
data. Accessing components of objects, such as with
\Rfunction{analysisPath} in the example above, often returns the data
in all its glory -- frequently, you'll want to adopt a strategy of
assigning such big data to a variable, and inspecting it with
functions like \Rfunction{str} or \Rfunction{head}.

A key functionality provided by \Rpackage{ShortRead} is input of a
diversity of file types, both of files from the Solexa 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, e.g., to
read all files containing image intensity produced by Firecrest. 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 files you are interested in.
\begin{Ex}
  Use \Rfunction{list.files} to display all files in the first Gerald
  analysis directory of the \Robject{sp} object.
<<filePattern>>=
list.files(analysisPath(sp)[[1]])
@ 
\end{Ex}
A full Gerald data set would contain hundreds of files. We'll select
just one file to use in subsequent analysis, based on files matching
the regular expression \texttt{.*_export_head.txt}. These files are
produced using \ELAND{} software run in \texttt{eland-extended}
mode. This mode produces files, one for each lane (and `end' of paired
end reads) that summarize diverse features of \emph{all} reads, and is
a very convenient starting point for analysis.
\begin{Ex}
  We'll use an abbreviated file for most parts of this lab. The
  abbreviated file contains the first 500,000 reads from a single lane
  of a Solexa paired-end run. It is from lane 1, and we will use the
  first end only. The name of the file is
  \file{s_1_1_export_head.txt}.  We will use this as the `pattern'
  to match, and check that we have specified the pattern appropriately
<<pattern>>=
pattern <- "s_1_1_export_head.txt"
list.files(analysisPath(sp), pattern)
@ 
%%
  Our success shouldn't be too surprising in this case, but it often
  pays to check. For instance, the pattern above would also match a file
  with the same name but with \texttt{.tar.gz} appended!
\end{Ex}

\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. By the default for this argument, all
files will be read. An optional third argument allows the user to specify which
type of file is to be read in.
\begin{Ex}
  Use \Rfunction{readAligned} to read in our abbreviated version of
  lane 1. \Rfunction{readAligned} is smart enough to know where
  alignemnt files are located in the Solexa file hierarchy.
<<readAligned>>=
aln <- readAligned(sp, pattern)
@ 
%%
  The \Rfunarg{sp} argument could have been replaced by a directory
  path, e.g., \code{"."} (if the file were in the current working
  directory) or \code{analysisPath(sp)}.

  The third argument (unspecified in the above) allows input of
  diverse types of Solexa alignment files, in addition to input of \MAQ{} and
  \Bowtie{} alignment files.  See the help page for
  \Rfunction{readAligned} for additional details.
\end{Ex}

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{Ex}
  Display the object we've 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 the number of positions recorded as \code{NA}.
<<aln-head>>=
head(sread(aln))
table(strand(aln), useNA="ifany")
sum(is.na(position(aln)))
@ 
%%
  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{Ex}
What are all the \code{NA} values returned by \Rfunction{strand} and
\Rfunction{position}? These correspond to reads that did not align to
the reference genome used by \ELAND{}; that about 1/2 the reads do not
align is below normal, making this an interesting opportunity for
quality assessment. The \Rfunction{strand} function returns a factor
with three levels. The first two describe reads aligned to the plus
and minus strands, the third (\code{*}) is available for successful
alignments where strand information is irrelevant.

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 \ELAND{}
alignment; details of calibration are to be found in Illumina
documentation. The alignments themselves also have qualities
associated with them, with the details of alignment quality differing
between alignment algorithms. 
\begin{Ex}
  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 Solexa. The details of the encoding can be obtained
  by querying \code{quality(aln)} for its \Rfunction{alphabet}; the
  letter \code{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{Ex}
 
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{Ex}
  Retrieve the alignment quality scores, determine how many align
  poorly, and visualize the distribution (Figure~\ref{fig:alignqual})
  of scores.
<<alignQuality,fig=TRUE,height=5>>=
alignQuality(aln)
q <- quality(alignQuality(aln))
sum(q==0)
print(densityplot(q[q>1], plot.points=FALSE, 
                  xlab="Alignment quality"))
@ 
% 
\begin{figure}
  \centering
  \includegraphics[width=.95\textwidth]{IOQAFig-alignQuality}
  \caption{Alignment quality}
  \label{fig:alignqual}
\end{figure}
\end{Ex}

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. For instance, \ELAND{} includes
information about whether the read passed a base-calling filter (based
on strength and consistency of early bases in the read), in addition
to the lane, tile, x and y coordinate of each read.
\begin{Ex}
  Use the \Rfunction{alignData} function to extract the additional
  information in the \ELAND{} alignment file. The underlying data in
  this object can be accessed as though it were a data frame, for
  instance to tally the number of reads passing Solexa base calling
  filter.
<<alignData>>=
alignData(aln)
table(alignData(aln)$filtering)
@ 
\end{Ex}

\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{Ex}
  Select just the aligned reads passing Solexa filtering, and aligning
  to the reference genome.
<<select>>=
filtIdx <- alignData(aln)$filtering=="Y"
alignedIdx <- !is.na(strand(aln))
aln[filtIdx & alignedIdx]
@ 
\end{Ex}
A different approach to subsetting is to use objects of class
\Rclass{SRFilter}. These can be particularly useful as an argument to
\Rfunction{readAligned}, in addition to use in interactive sessions.
\begin{Ex}
  Construct instances of built-in filters to select reads passing the
  Solexa filtering criterion, and uniquely aligning to a fully
  assembled chromosomes. These can be `composed' into a single overall
  filter, and applied to restrict available reads.
<<filter>>=
filt1 <- alignDataFilter(expression(filtering=="Y"))
filt2 <- chromosomeFilter("chr[0-9XYM]+.fa")
filt <- compose(filt1, filt2)
caln <- aln[filt(aln)]
caln
@ 
% 
\end{Ex}
The filters developed above could be used to filter reads while being
read in to R, e.g,. with
<<readAliged-filter,eval=FALSE>>=
readAligned(sp, pattern, filter=filt)
@ 
% 

The \Rfunction{srFilter} function can be used to create custom
filters. The idea is that filter functions accept a single argument
\Rfunarg{x} that is an object to be filtered, and returns a logical
vector that can be used to select elements of the object. 
\begin{Ex}
  As a first example, write and use a filter to select only a single
  read from all that align to a particular chromosome, position, and
  strand.
<<ualignFilter,keep.source=TRUE>>=
ualignFilter <- srFilter(function(x) {
    ## create a numerical index of reads. Divide the index, position,
    ## and strand information between chromosomes. Select the index of
    ## a single read at each unique position and strand. Return the
    ## selected index as a logical vector with the same length as x
    oindex <- seq_len(length(x))
    index <- tapply(oindex, chromosome(x), c)
    pdup <- tapply(position(x), chromosome(x), duplicated)
    sdup <- tapply(strand(x), chromosome(x), duplicated)
    keep <- oindex  %in% unlist(mapply(function(i, p, s) {
        i[!(p & s)]
    }, index, pdup, sdup))
}, name="select only one read per position & strand ")
caln[ualignFilter(caln)]
@ 
\end{Ex}
The filter functions built-in to \Rpackage{ShortRead} use a `factory'
pattern to create instances of each filter that `remember' how the
filters were created. For instance, \code{chromosomeFilter("chr2.fa")}
creates an instance of the chromosome filter to select only
chromosomes matching \code{chr2.fa}.
\begin{Ex}
  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)
caln[sample100(caln)]
@ 
\end{Ex}

\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.
\item Character-encoded base quality scores are interpreted as the
  default for the software package whose output is being parsed, e.g., as `Solexa
  fastq' for \ELAND{}. The object returned by \Rfunction{quality}
  applied to an \Rclass{AlignedRead} object is either
  \Rclass{FastqQuality} or \Rclass{SFastqQuality}.
\end{itemize}
Alignment programs sometimes offer the opportunity to customize
output; such customization needs to be accommodated when reads are
input using \Rpackage{ShortRead}.

\section{Additional input functions}

\Rpackage{ShortRead}, \Rpackage{Biostrings}, and the standard input
functions from \R{} provide additional tools for reading Solexa and
other alignment formats. The \Rfunction{readXStringColumns} function
provides a convenient way to read DNA and quality sequences into
compact data structures. \Rfunction{readFasta} and its counterpart
\Rfunction{readFastq} provide tools for reading FASTA- or FASTQ-
(i.e., including quality annotation) formatted files.
\begin{Ex}
  Files \file{_sequence.txt} contain fastq-formatted sequence and
  quality scores. A sample of this file type is available. Read these
  in to data structured defined in \Rpackage{ShortRead}.  The content
  of \Robject{reads} can be retrieved with the accessor functions \Rfunction{id},
  \Rfunction{sread}, and \Rfunction{quality}.
<<fastq>>=
ap <- analysisPath(sp)[[1]]
reads <- readFastq(ap, "s_1_1_sequence_head.txt$")
head(id(reads))
@ 
\end{Ex}
The \Rfunction{readPrb} functions reads `raw' base call quality scores
from \file{_prb} files in the \code{baseCallPath} directory; the
result is a \Rclass{BStringSet} object that compactly represents the
quality scores in a way analogous to the results of
\Rfunction{quality} applied to \code{aln}.
\begin{Ex}
  The \file{_export.txt} files read by \Rfunction{readAligned} are
  tab-delimited text files. The goal of this exercise is to read the
  DNA sequence and quality score columns in to \R{} as
  \Rclass{DNAStringSet} and \Rclass{BStringSet} objects.  The
  \Rclass{DNAStringSet} and \Rclass{BStringSet} classes represent DNA
  or `biological' strings; they extend the base class
  \Rclass{XString}.

  Start by parsing the first line of a \file{_prb} file to get a sense
  of its content. Specify the \code{colClasses} to be imported, using
  \code{NULL} to indicate that a column should be skipped, and
  \code{DNAString} or \code{BString} to indicate columns that are to
  be read as the corresponding data types (classes). Read the file with
  \Rfunction{readXStringColumns}.
<<readXStringColumns>>=
fl <- list.files(ap, pattern, full=TRUE)
cols <- strsplit(readLines(fl, 1), "\t")[[1]]
length(cols)
cols[9:10]

colClasses <- rep(list(NULL), 22)
colClasses[9:10] <- c("DNAString", "BString")
strings <- readXStringColumns(ap, pattern, colClasses=colClasses)
head(strings[[2]])
@ 
%% 
  For each column in \Robject{colClasses},
  \Rfunction{readXStringColumns} parses the corresponding column in
  all files matching the function argument \Rfunarg{pattern} into a
  single \Rclass{XStringSet}, i.e., concatenating the content of all
  files into a single object.
\end{Ex}
Many of the files produced by Solexa are simple text files. These can
be read in using standard \R{} input commands.
\begin{Ex}
  The \file{_int} files in the \Rfunction{imagaAnalysisPath} are tab-
  and space-delimited records of intensities measured along cycles of
  reads. Each line corresponds to a single cluster (i.e., putative
  read). The first four numbers indicate the lane, tile, x, and
  y. Subsequent numbers come in groups of 4, corresponding the
  intensities of the A, C, G, and T nucleotides over successive
  cycles. Here are coordinates of the first read, followed by the
  first 3 sets of intensities.
<<lowlevel-int-columns>>=
fl <- list.files(imageAnalysisPath(sp),".*_int.*", full=TRUE)
strsplit(readLines(gzfile(fl, open="rb"), 1), "\t")[[1]][1:7]
@ 
  % 
  The files are compressed, so we wrap the file in a \Rfunction{gzfile}
  function call to unzip the file prior to parsing. We now read the
  files into an object, and then convert the intensities into a standard
  array:
<<lowlevel-int>>=
int <- readIntensities(sp)
arr <- as(intensity(int), "array")
@ 
  % 
  Perhaps surprisingly, the intensities for the four nucleotides are
  not independent of one another (Figure~\ref{fig:intensities}). At
  later cycles, the amplitude of the intensities shrink toward zero
  and become less orthogonal. This presumably contributes to decreased
  quality of base calls.
<<intensities-plot-early,fig=TRUE,height=10>>=
print(splom(arr[,,5], pch="."))
@ 
\end{Ex}
\begin{figure}
  \centering
  \includegraphics[width=.95\textwidth]{IOQAFig-intensities-plot-early}
  \caption{Intensities from cycle 5}
  \label{fig:intensities}
\end{figure}

\section{Quality assessment} 

This part of the course addresses \Rpackage{ShortRead} facilities for
assessing quality, primarily of Solexa data. The \Rpackage{ShortRead}
functionality is mean to complement rather than replace QA tools
provided by the \ELAND{} pipeline.

\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 \file{_export.txt} file, and extracts information on reads
and their qualities. Evaluation of the function is straight-forward, e.g.,
\verb|qa <- qa(sp)|,
to visit all files in the \Robject{sp} \Rclass{SolexaPath}. 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{Ex}
  Rather than collating information during the lab, we load the data
  from a previously stored instance.
<<qa-input>>=
load(file.path("data", "qa_080828_081110.rda"))
qa
@ 
\end{Ex}
One feature of \Rpackage{ShortRead} that can speed up this stage of the
operation is the use of clustered computer resources and the
\Rpackage{Rmpi} 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 a subsequent
section. 

The QA information collated from the \file{_export.txt} files is
summarized into a PDF report using the \Rfunction{report}
function. This function currently requires a \LaTeX{} installation,
although the intention is that the report will eventually be generated
in different formats. The command to create the report is
\verb|rpt <- report(qa, dest=tempfile())|.
This creates a PDF file at the location specified by the argument
\Rfunarg{dest}.
\begin{Ex}
  Rather than create a report, we provide a sample, derived from the
  \Robject{qa} object loaded in the previous exercise. The sample is
  at \Sexpr{file.path("docs", "qa_080828_081110.pdf")}.
\end{Ex}
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.

\subsection{Exploring \Robject{qa}}

The \Robject{qa} object is a list-like structure with several
entities.
\begin{Ex}
  The \Robject{readCounts} element of \Robject{qa} is a simple data
  frame summarizing, on a per-lane bases, the total number of reads,
  the number of reads passing Solexa internal filtering, and the
  number of aligned reads.
<<qa-readCount>>=
qa[["readCounts"]]
@ 
% 
  Lanes 1-4 and 6-8 correspond to biologically interesting samples;
  lane 5 is the Solexa $\varphi$X-174 control lane. The number of reads
  (between 2.8 and 4.5 million) is low for a typical experiment
  (official guidelines are provided in Solexa documentation).

  It can be difficult to scan large numbers, so the QA report template
  defines functions that help to display the information in a more
  comprehensible fashion. Source these files into your current \R{}
  session, and read the second and third columns as a proportion of
  the first.
<<qa-template-source>>=
source(file.path("scripts", "qa_solexa.R"))
ppnCount(qa[["readCounts"]])
@ 
% 
  Refer to the QA report for further commentary on this and other
  aspects of the report
\end{Ex}

\subsection{Frequent sequences}

A feature of raw reads, and of many subsequent stages of short read
analysis, is a power law-like relationship between the number of times
a read occurs, and the number of occurrences of a particular sequence
in the sample. This information is contained in the
\Robject{sequenceDistribution} element of \Robject{qa}, and is the
result of running the \Rfunction{tables} command (defined in
\Rpackage{ShortRead}) on a \Rclass{DNAStringSet}
object.
\begin{Ex}
  Retrieve the \Robject{sequenceDistribution} element from
  \Robject{qa}; it is a simple data frame. Look at the contents of the
  data frame using \Rfunction{head}, and select just lane 5 raw reads.
  plot the power-law relationship between the number of reads and
  number of times reads occur. As an
  alternative display of the same information, plot the cumulative
  number of reads as a function of the number of times a read occurs.
<<hoover,fig=TRUE,width=4.5,height=5>>=
df <- qa[["sequenceDistribution"]]
df5raw <- df[df$lane=="s_5_1_export.txt" & df$type=="read",]
head(df5raw)
print(xyplot(log10(nReads)~log10(nOccurrences), df5raw,
             xlab="Copies per read (log 10)",
             ylab="Unique reads (log 10)"))
@ 
<<choover, fig=TRUE,width=4.5,height=5>>=
csum <- with(df5raw, cumsum(nReads * nOccurrences))
csum <- csum / csum[length(csum)]
print(xyplot(csum ~log10(nOccurrences), df5raw,
             xlab="Copies per read (log 10)",
             ylab="Cumulative proportion of reads",
             type="l"))
@
%% 
  Results appear in Figure~\ref{fig:hoover}. The cumulative form of
  this figure appears in the QA report.
\end{Ex}
\begin{figure}
  \includegraphics[width=.5\textwidth]{IOQAFig-hoover}
  \includegraphics[width=.5\textwidth]{IOQAFig-choover}
  \caption{Number of copies of unique reads}
  \label{fig:hoover}
\end{figure}

The power-law relationship between copies per read and number of
unique reads in the control lane consists of three components. At the
left of the graph are $>10^5$ reads that are each represented by only
one copy. These likely correspond to sequencing, base calling, or
other errors associated with the technology. At the right of the
figure are a small number of reads represented many times. These
`frequent' sequences are summarized in the \Robject{frequentSequences}
element of \Robject{qa}; frequent sequences are also reported by the
\Rfunction{tables} function of \Rpackage{ShortRead}.
\begin{Ex}
  Discover the frequent sequences amongst the raw and aligned reads of
  lane 5.
<<freq-seqs>>=
df <- qa[["frequentSequences"]]
head(df[df$lane=="s_5_1_export.txt" & df$type=="read",1:2])
@ 
\end{Ex}
Frequent sequences include poly-A reads, reads where only a few bases
were called, and reads with close similarity to the Solexa primer or
adapter sequence used in sample preparation.
\begin{Ex}
  Many of the primer sequences are filtered out by Solexa criteria,
  but it is worth discovering how many reads are `similar' to this
  sequence. Use the \Rfunction{srdistance} function to identify such
  reads, e.g., amongst those reads that contain no \texttt{N}
  nucleotides.
<<srdistance>>=
seq <- "CGGTTCAGCAGGAATGCCGAGATCGGAAGAGCGGT"
dist <- srdistance(clean(aln), seq)[[1]]
head(table(dist))
@ 
% 
  \Rfunction{srdistance} returns the edit distance between each read
  and the reference sequence, where the edit distance is defined so
  that each base mismatch represents an additional increment of
  1. There are \Sexpr{sum(table(dist)[1:3])} reads that differ at 2 or fewer
  locations, from amongst the \Sexpr{length(clean(aln))} reads in the
  cleaned sample used in this exercise.
\end{Ex}
Reads represented many times may be problematic for downstream
analysis. For instance, sample preparation protocols may involve a PCR
step that results in differential amplification, whereas the analysis
assumes reads are represented in proportion to their occurrence. The
\Rfunction{ualignFilter} function defined above, and the
\Rfunction{srduplicated} function in \Rpackage{ShortRead} represent
two approaches to dealing with this problem by ensuring that reads are
represented exactly once (by some definition of `once'!).

Sequences in the middle portion of the graph in
Figure~\ref{fig:hoover} will often, depending on the nature of the
investigation, represent the sequences of main biological
interest. These are sequences represented an intermediate number of
times, as might be required for reasonable coverage in a SNP discovery
or ChIP-seq experiment. The right-hand graph in
Figure~\ref{fig:hoover} shows a relatively abrupt transition between
reads represented rarely and those represented many times. 

The sample QA report shows that the non-control lanes show much
broader transitions from rarely to frequently represented reads. This
could represent technical shortcomings of this run (e.g., inadequate
enrichment of sample DNA) or features of intrinsic biological interest
(e.g., wide variability in ChIP abundance between binding
sites). Regardless of ultimate source, the broad distribution of read
occurrences implies significant effort may be required to distinguish
`noise' (reads corresponding to those in the left and right portions
of the control lane graph) from signal.

Finally, while the control lane shows a relatively abrupt transition
between reads that occur rarely and those that are common
(Figure~\ref{fig:hoover}), the distribution is in fact 3- or 4-fold
broader than expected under a naive model of random read starts along
the $\varphi$X-174 genome. This is reinforced by alignments to the
reference genome, where clear patterns (e.g., unequal representation
on plus and minus strands; non-uniform coverage) are apparent.
\begin{Ex}
  As an advanced exercise, simulate reads selected uniformly along
  both strands of the 5200bp long $\varphi$X-174 genome. Compare the
  times each read is represented in your sample with those from the
  actual control lane. Hint: use \Rfunction{sample} with
  \Rfunarg{replace=TRUE} to generate the reads, and
  \code{table(table(reads))} to summarize their occurrence; this should
  take less than 5 lines of \R{} code.
\end{Ex}

\subsection{Cycle-specific qualities and base calls}

As a final foray into the details of quality assessment, consider base
calls and quality, and how these change across cycles (see the second
table and Section 4 of the QA report).

The table in the QA report suggests that the control lane (lane 5) is
enriched for A and T; this is confirmed by the figure in section
4. This likely reflects underlying differences in the genomic regions
represented in each lane.
\begin{Ex}
  Use \Rfunction{alphabetFrequency} to summarize nucleotide use in the
  short reads in \Robject{aln}, from the first part of this lab.
<<alphabetFreqeuncy>>=
alphabetFrequency(sread(aln), collapse=TRUE, baseOnly=TRUE, freq=TRUE)
@
%
  The frequency of `other' (i.e., uncalled) nucleotides ($>9\%$) is
  very high; typical runs are $<3\%$; recent runs with GAII
  technologies after 36 cycles are $<1\%$.
\end{Ex}
An unexpected aspect of the figure in Section 4 of the QA report is
apparent trends in nucleotide frequency with cycle. For instance, all
lanes show a marked decrease in A and increase in C across
cycles. This is unexpected in this experiment, where the \emph{a
  priori} expectation is that sequences start at essentially arbitrary
locations in the sequenced DNA: an A is expected as frequently at the
beginning of the sequence as at the end.
\begin{Ex}
  Use \Rfunction{alphabetByCycle} to extract the number of each
  nucleotide sequenced at each cycle. Convert the matrix into a data
  frame containing only the called nucleotides, and plot these counts
  as a function of cycle.
<<abc,fig=TRUE,height=5>>=
abc <- alphabetByCycle(sread(aln))
dim(abc)
abc[1:4,1:5]
abc <- abc[rowSums(abc)!=0,]
df <- as.data.frame(t(abc[1:4,]))
print(xyplot(A+C+G+T~1:nrow(df), df, type="l",
             auto.key=list(x=.75, y=.95, points=FALSE, lines=TRUE),
             xlab="Cycle", ylab="Count"))
@ 
\end{Ex}
\begin{figure}
  \centering
  \includegraphics[width=.95\textwidth]{IOQAFig-abc}
  \caption{Nucleotide frequency per cycle, lane 1}
  \label{fig:abc}
\end{figure}
There are a number of possible contributors to cycle-dependent
nucleotide frequencies, including inadequate reagent volume, and
nucleotide-specific differential accumulation of fluorescent
dyes. Figure~\ref{fig:abc} and the figure in Section 4 of the report
contain additional features that are moderately unexpected, and
unexplained. For instance, the frequency of a nucleotide such as A
changes very systematically across cycles, first increasing and then
decreasing; this seems more regular than expected. Patterns of
nucleotide change also seem to echo one another at similar cycles but
in different lanes, even when the lanes have different biological
material. This occurs for instance in lanes 1-4 of the QA report,
where the last 5 cycles of the C nucleotide seem to change in
(comparative!) unison.

\section{An advanced note}

Several functions are designed to use the \Rpackage{Rmpi} package for
distributed computation, if installed. For instance
<<Rmpi,eval=FALSE>>=
library("Rmpi")
mpi.spawn.Rslaves(nsl=8)
qa <- qa(sp)
mpi.close.Rslaves()
@ 
%% 
distributes the calculation of quality assurance summaries across 8
processors. The \Rfunction{srapply} function can be used to distribute
your own calculations.

\section{Summary}

This portion of the course has provided you with an overview of the input and
quality assessment functionality available in the \Rpackage{ShortRead}
package. A summary of the insights learned might be reflected in
simple, and somewhat naive, work flows.

The first work flow simply performs quality assessment on \ELAND{}
aligned data.
<<workflow-qa, eval=FALSE>>=
sp <- SolexaPath("extdata/ELAND/080828_HWI-EAS88_0003")
rpt <- report(qa(sp), dest="reports/my_report.pdf")
@ 
% 
Performing QA on \ELAND{} data does not commit us to using \ELAND{}
alignments in subsequent steps.

The second work flow reads aligned data in to \R{}. The work flow
might start with aligned reads created by one of many aligners; we
start with \ELAND{} \file{_export.txt} files. There are many possible
issues highlighted in the forgoing discussion. We choose to establish
a series of filters to eliminate some reads at the very start of our
work flow. We eliminate reads with ambiguous base calls and failing
Solexa's internal filtering criteria. Reads aligning to multiple
locations in the genome are not straight-forward to deal with, and are
not essential for ChIP-seq style experiments (this is not the case for
expression or RNA-seq experiments), 
% FIXME (wh 6 June 009) I wonder whether this is really true - why would protein-DNA 
%   interactions in repetitive parts of the genome be more or less interesting
%   than the DNA transcription?
so we remove these. Most close
matches to the Solexa primer sequence are flagged as not passing
Solexa base call filters, but we eliminate reads near to this as
well. Finally, we restrict our attention to those reads that align to
assembled nuclear chromosomes, putting aside for the moment those reads
aligning to organelle genomes (the X and Y chromosomes also require
special consideration, and we might often eliminate these from a 
first-pass work flow, too). 
% FIXME (wh 6 June 009) But the example below does read the reads that mapped to X and Y.
Our second work flow is thus:
<<filter,keep.source=TRUE,eval=FALSE>>=
filt1 <- nFilter()
filt2 <- alignDataFilter(expression(filtering=="Y"))
filt3 <- alignQualityFilter(threshold=1)
filt4 <- srdistanceFilter("CGGTTCAGCAGGAATGCCGAGATCGGAAGAGCGGT", 4)
filt5 <- chromosomeFilter("chr[0-9XY]+.fa")

filt <- compose(filt1, filt2, filt3, filt4, filt5)
aln <- readAligned(sp, "s_1_1_export.txt$", filter=filt)
@ 
% 
This work flow applies equally to MAQ aligned data, with the exception
that Solexa filtering criteria are not available and the chromosome
naming convention in the MAQ-aligned reads in our sample are
different:
<<filter-MAQ,eval=FALSE>>=
maqDir <- file.path("extdata", "MAQ")
filt5 <- chromosomeFilter("chr[0-9XY]+$")
filt <- compose(filt1, filt3, filt4, filt5)
maq <- readAligned(maqDir, "s_8.map", "MAQMap", filter=filt)
@ 
%% 
Each of the filters represents a decision. The decision may be
inappropriate for particular analyses, and may be revisited as
understanding of the data matures.

\section{Session information}

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

\end{document}