%\VignetteIndexEntry{03. Data exploration}
%\VignettePackage{useR2011}

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

<<setup, echo=FALSE>>=
stopifnot(getRversion() >= "2.13.0" && getRversion() < "2.14.0")
library(useR2011)
@ 

\section{Exploring sequence data: short reads and alignments}

The following sections introduce core tools for working with
high-throughput sequence data. This section focus on the reads and
alignments that are the raw material for
analysis. Section~\ref{sec:rnaseq} addresses statistical approaches to
assessing differential representation in RNA-seq experiments.

\subsection{The \emph{pasilla} data set}

As a running example, we use the \emph{pasilla} data set, derived from
\cite{Brooks2010}. The authors investigate conservation of RNA
regulation between \Dmel{} and mammals. Part of their study used RNAi
and RNA-seq to identify exons regulated by Pasilla (\emph{ps}), the
\Dmel{} ortholog of mammalian NOVA1 and NOVA2. Briefly, their
experiment compared gene expression as measured by RNAseq in S2-DRSC
cells cultured with, or without, a 444bp dsRNA fragment corresponding
to the \emph{ps} mRNA sequence. Their assessment investigated
differential exon use, but our worked example will focus on gene-level
differences.

In this section we look at a subset of the \emph{ps} data,
corresponding to reads obtained from lanes of their RNA seq
experiment, and to the same reads aligned to a \Dmel{} reference
genome. Reads were obtained from GEO and the Short Read Archive (SRA);
reads were aligned to \Dmel{} reference genome dm3 as described in the
\Rpackage{pasilla} experiment data package (this package is
\href{http://bioconductor.org/packages/devel/data/experiment/}{available}
in the `devel' version of \Bioconductor{}; the devel version will be
released in October, 2011).

\subsection{Short reads}

\paragraph{Sequencer technologies}

The Illumina GAII and HiSeq technologies generate sequencers by
measuring incorporation of florescent nucleotides over successive PCR
cycles. These sequencers produce produce output in a variety of
formats, but \emph{FASTQ} is ubiquitous. Each read is represented by a
record of four components:
<<fastq-format, echo=FALSE>>=
fls <- list.files(file.path(bigdata(), "fastq"), full=TRUE)
cat(noquote(readLines(fls[[1]], 4)), sep="\n")
@ 

The first and third lines (beginning with \verb|@| and \verb|+|
respectively) are unique identifiers. In the sample above the
identifier produces by the sequencer typically includes a machine id
followed by colon-separated information on the lane, tile, x, and y
coordinate of the read. The example illustrated here also includes the
SRA accession number, added when the data was submitted to the
archive. The machine identifier could potentially be used to extract
information relative to batch effects. The spacial coordinates (lane,
tile, x, y) are often used to identify optical duplicates (artifacts
introduced when the sequencer falsely interprets reads in close
proximity spatial as distinct); spacial coordinates can also be used
during quality assessment to identify spatial artifacts of sequencing,
though these spatial effects are rarely pursued.

The second and fourth lines of the FASTQ record are the nucleotides
and qualities of each cycle in the read. It is given in 5' to 3'
orientation as seen by the sequencer. A letter \texttt{N} is used to
signify bases that the sequencer was not able to call.  The fourth
line of the FASTQ record encodes the quality (confidence) of the
corresponding base call. The quality score is encoded following one of
several conventions, with the general notion being that letters later
in the visible ASCII alphabet
<<ascii, echo=FALSE>>=
cat(rawToChar(as.raw(32+1:47)),
    rawToChar(as.raw(32+48:94)), sep="\n")
@ 
are of lower quality; this is developed further below. Both the
sequence and quality scores may span multiple lines.

Technologies other than Illumina use different formats to represent
sequences. Roche 454 sequence data is generated by `flowing' labeled
nucleotides over samples, with greater intensity corresponding to
longer runs of A, C, G, or T. This data is represented as a series of
`flow grams' (a kind of run-length encoding of the read) in Standard
Flowgram Format (SFF). The \Bioconductor{} package
\Biocpkg{R453Plus1Toolbox} has facilities for parsing SFF files, but
after quality control steps the data are frequently represented (with
some loss of information) as FASTQ. SOLiD technologies produce
sequence data using a `color space' model. This data is not easily
read in to \R{}, and much of the error-correcting benefit of the color
space model is lost when converted to FASTQ; SOLiD sequences are not
well-handled by \Bioconductor{} packages.

\paragraph{Short reads in \R{}}

FASTQ files can be read in to \R{} using the \Rfunction{readFastq}
function from the \Biocpkg{ShortRead} package. Use this function by
providing the path to a FASTQ file. There are sample data files
available in the \Rpackage{useR2011} package, each consisting of 1
million reads from a lane of the Pasilla data set.
<<readFastq>>=
fastqDir <- file.path(bigdata(), "fastq")
fastqFiles <- list.files(fastqDir, full=TRUE)
fq <- readFastq(fastqFiles[1], withIds=TRUE)
fq
@ 
The data are represented as an object of class \Rclass{ShortReadQ}
(`short read and \emph{q}uality').
<<sread>>=
head(sread(fq), 3)
head(quality(fq), 3)
head(id(fq), 3)
@ 
The \Rclass{ShortReadQ} class illustrates \emph{class
  inheritance}. It extends the \Rclass{ShortRead} class
<<getClasse>>=
getClass("ShortReadQ")
@ 
Methods defined on \Rclass{ShortRead} are available for
\Rclass{ShortReadQ}.
<<showMethods-ShortRead, eval=FALSE>>=
showMethods(class="ShortRead", where=getNamespace("ShortRead"))
@ 
For instance, the \Rfunction{width} can be used to demonstrate that
all reads consist of 37 nucleotides.
<<width-ShortRead>>=
table(width(fq))
@ 
The \Rfunction{alphabetByCycle} function summarizes use
of nucleotides at each cycle in a (equal width) \Rclass{ShortReadQ} or
\Rclass{DNAStringSet} instance.
<<width-ShortReadQ>>=
abc <- alphabetByCycle(sread(fq))
abc[1:4, 1:8]
@ 

FASTQ files are getting larger. A very common reason for looking at
data at this early stage in the processing pipeline is to explore
sequence quality. In these circumstances it is often not necessary to
parse the entire FASTQ file. Instead create a representative sample
<<FastqSampler>>=
sampler <- FastqSampler(fastqFiles[1], 1000000)
yield(sampler) # sample of 1000000 reads
@ 

\Biocpkg{ShortRead} contains facilities for quality assessment of
FASTQ files. Here we generate a report from a sample of 1 million
reads from each of our files and display it in a web browser
<<qa, eval=FALSE>>=
qas0 <- Map(function(fl, nm) {
    fq <- FastqSampler(fl)
    qa(yield(fq), nm)
}, fastqFiles, 
   sub("_subset.fastq", "", basename(fastqFiles)))
qas <- do.call(rbind, qas0)
rpt <- report(qas, dest=tempfile())
browseURL(rpt)
@ 
A report from a larger subset of the experiment is available
<<report,eval=FALSE >>=
rpt <- system.file("GSM461176_81_qa_report", package="useR2011")
browseURL(rpt)
@ 

\begin{Exercise}
  Use the helper function \Rfunction{bigdata} (defined in the
  \Rpackage{useR2011} package) and the \Rfunction{file.path} and
  \Rfunction{list.files} functions to locate two fastq files from
  \cite{Brooks2010} (the files were obtained as described in the
  appendix and \Rpackage{pasilla} experiment data package, available
  with the development versions of \R{} and \Bioconductor{}).
  
  Input one of the fastq files using \Rfunction{readFastq} from the
  \Biocpkg{ShortRead} package.
  
  Use \Rfunction{alphabetFrequency} to summarize the GC content of all
  reads (hint: use the \Rfunction{sread} accessor to extract the
  reads, and the \Rfunarg{collapse=TRUE} argument to the
  \Rfunction{alphabetFrequency} function). Using the helper function
  \Rfunction{gcFunction} from the \Rpackage{useR2011} package, draw a
  histogram of the distribution of GC frequencies across reads.
  
  Use \Rfunction{alphabetByCycle} to summarize the frequency of each
  nucleotide, at each cycle. Plot the results using
  \Rfunction{matplot}, from the \Rpackage{graphics} package.
  
  As an advanced exercise, and if on Mac or Linux, use the
  \Rpackage{multicore} package and \Rfunction{mclapply} to read and
  summarize the GC content of reads in two files in parallel.
\end{Exercise}

\begin{Solution}
Discovery:
<<fastq-discovery>>=
list.files(bigdata())
fls <- list.files(file.path(bigdata(), "fastq"), full=TRUE)
@ 
Input:
<<fastq-input>>=
fq <- readFastq(fls[1])
@ 
GC content:
<<gcC>>=
alf0 <- alphabetFrequency(sread(fq), as.prob=TRUE, collapse=TRUE)
sum(alf0[c("G", "C")])
@ 
A histogram of the GC content of individual reads is obtained with
<<gc-reads>>=
gc <- gcFunction(sread(fq))
hist(gc)
@ 
Alphabet by cycle:
<<abc>>=
abc <- alphabetByCycle(sread(fq))
matplot(t(abc[c("A", "C", "G", "T"),]), type="l")
@ 
Advanced (Mac, Linux only): processing on multiple cores.
<<abc-mclapply, eval=FALSE>>=
library(multicore)
gc0 <- mclapply(fls, function(fl) {
  fq <- readFastq(fl)
  gc <- gcFunction(sread(fq))
  table(cut(gc, seq(0, 1, .05)))
})
## simplify list of length 2 to 2-D array
gc <- simplify2array(gc0)
matplot(gc, type="s")
@ 
\end{Solution}

\begin{Exercise}
  Use \Rfunction{quality} to extract the quality scores of the short
  reads. Interpret the encoding qualitatively.
  
  Convert the quality scores to a numeric matrix, using
  \Rfunction{as}. Inspect the numeric matrix (e.g., using
  \Rfunction{dim}) and understand what it represents.
  
  Use \Rfunction{colMeans} to summarize the average quality score by
  cycle. Use \Rfunction{plot} to visualize this
\end{Exercise}

\begin{Solution}
<<quality>>=
head(quality(fq))
qual <- as(quality(fq), "matrix")
dim(qual)
plot(colMeans(qual), type="b")
@ 
\end{Solution}

\subsection{Alignments}

Most down-stream analysis of short read sequences is based on reads
aligned to reference genomes. There are many aligners available,
including
\href{http://bio-bwa.sourceforge.net/}{BWA}~\cite{pmid20080505,pmid19451168},
\href{http://bowtie-bio.sourceforge.net/index.shtml}{Bowtie}~\cite{pmid19261174},
\href{http://research-pub.gene.com/gmap/}{GSNAP}, and Illumina's
ELAND; merits of these are discussed in the literature. There are also
alignment algorithms implemented in \Bioconductor{} (e.g.,
\Rfunction{matchPDict} in the \Biocpkg{Biostrings} package, and the
\Biocpkg{Rsubread} package); \Rfunction{matchPDict} is particularly
useful for flexible alignment of moderately sized subsets of data.

\paragraph{Alignment formats}

Most main-stream aligners produce output in SAM (text-based) or BAM
format. A SAM file is a text file, with one line per aligned read, and
fields separated by tabs. Here is an example, split into fields.
<<SAM>>=
fl <- system.file("extdata", "ex1.sam", package="Rsamtools")
strsplit(readLines(fl, 1), "\t")[[1]]
@

Fields in a SAM file are summarized in Table~\ref{tbl:sam}.
\begin{table}
  \centering
  \caption{Fields in a SAM record. From
    \url{http://samtools.sourceforge.net/samtools.shtml}}
  \medskip
  \label{tbl:sam}
  \begin{tabular}{lll}
    Field & Name & Value \\\hline
    1 & QNAME & Query (read) NAME \\
    2 & FLAG & Bitwise FLAG, e.g., strand of alignment \\
    3 & RNAME & Reference sequence NAME \\
    4 & POS & 1-based leftmost POSition of sequence \\
    5 & MAPQ & MAPping Quality (Phred-scaled) \\
    6 & CIAGR & Extended CIGAR string \\
    7 & MRNM & Mate Reference sequence NaMe \\
    8 & MPOS & 1-based Mate POSistion \\
    9 & ISIZE & Inferred insert SIZE \\
    10 & SEQ & Query SEQuence on the reference strand \\
    11 & QUAL & Query QUALity \\
    12$+$ & OPT & OPTional fields, format TAG:VTYPE:VALUE \\
  \end{tabular}
\end{table}
We recognize from the FASTQ file the identifier string, read sequence
and quality. The alignment is to a chromosome `seq1' starting at
position 1. The strand of alignment is encoded in the `flag', field.
The alignment record also includes a measuring of mapping quality and
a CIGAR string describing the nature of the alignment. In this case,
the CIGAR is 35M, indicating that the alignment consisted of 35
\texttt{M}atches or mismatches, with no indels or gaps; indels are
represented by \texttt{I} and \texttt{D}; gaps (e.g., from alignments
spanning introns) by \texttt{N}.

BAM files encode the same information as SAM files, but in a format
that is more efficiently parsed by software; BAM files are the primary
way in which aligned reads are imported in to \R{}. 

\paragraph{Aligned reads in \R{}}

The \Rfunction{readGappedAlignments} function from the
\Biocpkg{GenomicRanges} package reads essential information from a BAM
file in to \R. The result is an instance of the
\Rclass{GappedAlignments} class. The \Rclass{GappedAlignments} class
has been designed to allow useful manipulation of many reads (e.g., 20
million) under moderate memory requirements (e.g., 4 GB).
<<readGappedAlignments>>=
alnFile <- system.file("extdata", "ex1.bam", package="Rsamtools")
aln <- readGappedAlignments(alnFile)
head(aln, 3)
@
The \Rfunction{readGappedAlignments}
function takes an additional parameter, \Rfunarg{which}, allowing the
user to specify regions of the BAM file (e.g., known gene coordinates)
from which to extract alignments.

A \Rclass{GappedAlignments} instance is like a data frame, but with
accessors as suggested by the column names. It is easy to query for,
e.g., the distribution of reads aligning to each strand, the width of
reads, or the cigar strings
<<GappedAlignments-accessors>>=
table(strand(aln))
table(width(aln))
head(sort(table(cigar(aln)), decreasing=TRUE))
@ 

\begin{Exercise}
  Use \Rfunction{bigdata}, \Rfunction{file.path} and
  \Rfunction{list.files} to obtain file paths to the BAM files. These
  are a subset of the aligned reads, overlapping just four genes.
  
  Input the aligned reads from one file using
  \Rfunction{readGappedAlignments}. Explore the reads, e.g., using
  \Rfunction{table} or \Rfunction{xtabs} to summarize which chromosome
  and strand the subset of reads is from.
  
  The object \Robject{ex} created earlier contains coordinates of four
  genes. Use \Rfunction{countOverlaps} to first determine the number
  of genes an individual read aligns to, and then the number of
  uniquely aligning reads overlapping each gene. Since the RNAseq
  protocol was not strand-sensitive, set the strand of \Robject{aln}
  to \code{*}.
  
  Write the sequence of steps required to calculate counts as a simple
  function, and calculate counts on each file. On Mac or Linux, can
  you easily parallelize this operation?
\end{Exercise}

\begin{Solution}
  We discover the location of files using standard \R{} commands:
<<bam-ex-fls>>=
fls <- list.files(file.path(bigdata(), "bam"), "bam$", full=TRUE)
names(fls) <- sub("_.*", "", basename(fls))
@ 
Use \Rfunction{readGappedAlignments} to input data from one of the
files, and standard \R{} commands to explore the data.
<<bam-ex-input>>=
## input
aln <- readGappedAlignments(fls[1])
xtabs(~rname + strand, as.data.frame(aln))
@
To count overlaps in regions defined in a previous exercise, load
the regions.
<<bam-ex-roi>>=
data("ex")           # from an earlier exercise
@ 
Many RNA-seq protocols are not
strand aware, i.e., reads align to the plus or minus strand regardless
of the strand on which the corresponding gene is encoded. Adjust the
strand of the aligned reads to indicate that strand is not known.
<<bam-ex-strand>>=
strand(aln) <- "*"   # protocol not strand-aware
@ 
For simplicity, we are interested in reads that align to only a single
gene. Count the number of genes a read aligns to\ldots
<<bam-ex-hits>>=
hits <- countOverlaps(aln, ex)
table(hits)
@ 
and reverse the operation to count the number of times each region of
interest aligns to a uniquely overlapping alignment.
<<bam-ex-cnt>>=
cnt <- countOverlaps(ex, aln[hits==1])
@ 
A simple function for counting reads is
<<bam-count-fun>>=
counter <- 
    function(filePath, range)
{
    aln <- readGappedAlignments(filePath)
    strand(aln) <- "*"
    hits <- countOverlaps(aln, range)
    cnt <- countOverlaps(range, aln[hits==1])
    names(cnt) <- names(range)
    cnt
}
@ 
This can be applied to all files using \Rfunction{sapply}
<<bam-count-all>>=
counts <- sapply(fls, counter, ex)
@ 
The counts in one BAM file are independent of counts in another BAM
file. This encourages us to count reads in each BAM file in parallel,
decreasing the length of time required to execute our program. On
Linux and Mac OS, a straight-forward way to parallelize this operation
is:
<<bam-count-mclapply, eval=FALSE>>=
if (require(multicore))
    simplify2array(mclapply(fls, counter, ex))
@ 
\end{Solution}

The \Rclass{GappedAlignments} class inputs only some of the fields of
a BAM file, and may not be appropriate for all uses. In these cases
the \Rfunction{scanBam} function in \Biocpkg{Rsamtools} provides
greater flexibility. The idea is to view BAM files as a kind of data
base. Particular regions of interest can be selected, and the
information in the selection restricted to particular fields. These
operations are determined by the values of a \Rclass{ScanBamParam}
object, passed as the named \Rfunarg{param} argument to
\Rfunction{scanBam}.

\begin{Exercise}
  Consult the help page for \Rfunction{ScanBamParam}, and construct an
  object that restricts the information returned by a
  \Rfunction{scanBam} query to the aligned read DNA sequence. Your
  solution will use the \Rfunarg{what} parameter to the
  \Rfunction{ScanBamParam} function.

  Use the \Rfunction{ScanBamParam} object to query a BAM file, and
  calculate the GC content of all aligned reads. Summarize the
  GC content as a histogram (Figure~\ref{fig:gc-read}).
\end{Exercise}
\begin{figure}
  \centering
  \includegraphics[width=.6\textwidth]{DataExploration-gc-read}
  \caption{GC content in aligned reads}
  \label{fig:gc-read}
\end{figure}
\begin{Solution}
<<gc-read, fig=TRUE, include=FALSE>>=
param <- ScanBamParam(what="seq")
seqs <- scanBam(fls[[1]], param=param)
readGC <- gcFunction(seqs[[1]][["seq"]])
hist(readGC)
@ 
\end{Solution}
