%\VignetteIndexEntry{04. RNA-seq}
%\VignettePackage{useR2011}

\SweaveOpts{keep.source=TRUE, eps=FALSE}
<<setup, echo=FALSE>>=
options(digits=2)
library(useR2011)
library(edgeR)
library(goseq)
@ 
\section{RNA-seq}
\label{sec:rnaseq}

\subsection{Varieties of RNA-seq}

RNA-seq experiments typically ask about differences in representation
of genes or other features across experimental groups. The analysis of
designed experiments is of course statistical, and hence an ideal task
for \R. The overall structure of the analysis, with tens of thousands
of features and tens of samples, is also reminiscent of microarray
analysis; one might hope that insights from the microarray domain will
apply, at least conceptually, to the analysis of RNA-seq experiments.

The most straight-forward RNA-seq experiments quantify abundance of
known gene models. The known models are derived from reference
databases, reflecting the accumulated wisdom of the community
responsible for the data. The `knownGenes' track of the UCSC genome
browser represents one source of data. It contains, for each gene, the
transcripts and exons that are thought through experimental or
computational approaches to exist. The \Biocpkg{GenomicFeatures}
package allows ready access to this information, as we have seen.  The
data base of known genes is coupled with high throughput sequence data
by counting or otherwise estimating the number of reads associated
with each gene.

A more ambitious approach to RNA-seq attempts to identify novel
transcripts. This requires that sequenced reads be assembled into
contigs that, presumably, correspond to expressed transcripts that are
then located in the genome. Transcripts identified in this way may
correspond to known transcripts, to novel organization of known exons
(e.g., through alternative splicing), or to completely novel
constructs. We will not address the identification of completely novel
transcripts here, but note that having quantified transcript
abundances in several samples one is still interested in the analysis
of designed experiments -- do transcript abundances, novel or
otherwise, differ between experimental groups?

\Bioconductor{} packages play a role in several stages of an RNA-seq
analysis. The \Biocpkg{GenomicRanges} infrastructure we have already
been exposed to can be effectively employed to quantify known exon or
transcript abundances. Quantified abundances are in essence a matrix
of counts, with rows representing features and columns samples. The
\Biocpkg{edgeR}~\cite{pmid19910308} and
\Biocpkg{DESeq}~\cite{AndersHuber2010} packages facilitate analysis of
this data in the context of designed experiments, and are appropriate
when the questions of interest involve between-sample comparisons of
relative abundance. The \Biocpkg{DEXSeq} package extends the approach
in \Biocpkg{edgeR} and \Biocpkg{DESeq} to ask about within-gene,
between group differences in exon use, i.e., for a given gene, do
groups differ in their exon use?

\subsection{Data preparation}

\paragraph{Counting reads aligning to genes}

An essential step is to arrive at some measure of gene representation
amongst the aligned reads. A straight-forward and commonly used
approach is to count the number of times a read overlaps exons. Nuance
arises when a read only partly overlaps an exon; when two exons
overlap (and hence a read appears to be `double counted'); when reads
are aligned with gaps, and the gaps are inconsistent with known exon
boundaries; etc. The \Rfunction{countGenomicOverlaps} function in the
\Biocpkg{GenomicRanges} package provides facilities for implementing
different count strategies, using arguments such as \Rfunarg{type} to
determine the nature of the exon / read overlap, and
\Rfunarg{resolution} to select strategies for counting reads
overlapping multiple subjects. The behavior of the function is also
influenced by whether the queried region is an \Rclass{GRanges} or
\Rclass{GRangesList} instance. 

Software other than \R{} can also be used to summarize count data. An
important point is that the desired input is often raw count data,
rather than normalized (e.g., reads per kilobase of (gene) model per
million mapped reads) values. This is because counts allow information
about uncertainty of estimates to propagate to later stages in the
analysis.

\paragraph{Object creation and filtering}

The following exercise illustrates key steps in data preparation.

\begin{Exercise}
  The \Rpackage{user2011} package contains a data set with pre-computed count
  data. Use the \Rfunction{data} command to load it. Create a variable
  \Robject{grp} to define the groups associated with each column,
  using the column names as a proxy for more authoritative metadata.
  
  Create a \Rfunction{DGEList} object (defined in the \Biocpkg{edgeR}
  package) from the count matrix and group information. Calculate
  relative library sizes using the \Rfunction{calcNormFactors}
  function.
  
  A lesson from the microarray world is to discard genes that cannot
  be informative (e.g., because of lack of variation). Filter reads to
  remove those that are represented at less than 1 per million mapped
  reads, in fewer than 2 samples.
  
  Use \Rfunction{plotMDS} on the filtered reads to perform
  multi-dimensional scaling. Interpret the resulting plot.
\end{Exercise}

\begin{Solution}
  Here we load the data (a matrix of counts) and create treatment
  group names from the column names of the counts matrix.
<<counts>>=
data(counts)
dim(counts)
grp <- factor(sub("[1-4].*", "", colnames(counts)),
              levels=c("untreated", "treated"))
@ 

We use the \Biocpkg{edgeR} package, creating a \Rclass{DGEList} object
from the count and group data. The \Rfunction{calcNormFactors}
function estimates relative library sizes for use as offsets in the
generalized linear model.
<<DGEList>>=
library(edgeR)
dge <- DGEList(counts, group=grp)
dge <- calcNormFactors(dge)
@ 

To filter reads, we scale the counts by the library sizes and express
the results on a per-million read scale. We require that the gene be
represented at a frequency of at least 1 read per million mapped in
two or more of each sample, and use this criterion to subset the
\Rclass{DGEList} instance.
<<DEGList-filter>>=
m <- 1e6 * t(t(dge$counts) / dge$samples$lib.size)
ridx <- rowSums(m > 1) >= 2
table(ridx)   # number filtered / retained
dge <- dge[ridx,]
@

Multi-dimensional scaling takes data in high dimensional space (in our
case, the dimension is equal to the number of genes in the filtered
\Rclass{DGEList} instance) and reduces it to fewer (e.g., 2)
dimensions, allowing easier assessment. The plot is shown in
Figure~\ref{fig:mds}; that the samples separate into distinct groups
provides some reassurance that the data differ according to
treatment. Nonetheless, there appears to be considerable heterogeneity
within groups. Any guess, perhaps from looking at the quality report
generated early, what the within-group differences are due to?
<<mds, fig=TRUE, include=FALSE>>=
plotMDS.dge(dge)
@ 
\end{Solution}
\begin{figure}
  \centering
  \includegraphics[width=.6\textwidth]{RNASeq-mds}
  \caption{MDS plot of lanes from the Pasilla data set.}
  \label{fig:mds}
\end{figure}

\subsection{Differential representation}

RNA-seq differential representation experiments, like classical
microarray experiments, consist of a single statistical design (e.g,
comparing expression of samples assigned to `Treatment' versus
`Control' groups) applied to each feature for which there are aligned
reads. While one could naively perform simple tests (e.g., t-tests) on
all features, it is much more informative to identify important
aspects of RNAseq experiments, and to take a flexible route through
this part of the work flow. Key steps involve formulation of a model
matrix to capture the experimental design, estimation of a test static
to describe differences between groups, and calculation of a $P$ value
or other measure as a statement of statistical significance.

\paragraph{Experimental design} 

In \R, an experimental design is specified with the
\Rfunction{model.matrix} function. The function takes as its first
argument a \Robject{formula} describing the independent variables and
their relationship, and as a second argument a \Robject{data.frame}
containing the (phenotypic) data that the formula describes. A simple
formula might read \verb|~ 1 + grp|, which says that the response is a
linear function involving an intercept (1) plus a term encoded in by
the variable \Robject{grp}. If (as in our case) \Robject{grp} is a
factor, then the first coefficient (column) of the model matrix
corresponds to the first level of \Robject{grp}, and subsequent terms
correspond to \emph{deviations} of each level from the first. If
\Robject{grp} were \Rclass{numeric} rather than \Rclass{factor}, the
formula would represent linear regressions with an intercept. Formulas
are very flexible, allowing representation of models with one, two, or
more factors as main effects, models with or without interaction, and
with nested effects.  

\begin{Exercise}
  To be a more concrete, use the \Rfunction{model.matrix} function and
  a formula involving \Robject{grp} to create the model matrix for our
  experiment.
\end{Exercise}

\begin{Solution}
  Here is the experimental design; it's worth discussing with your
  neighbor the interpretation of the \Robject{design} instance.
<<design>>=
(design <- model.matrix( ~ grp ))
@ 
\end{Solution}

The coefficient (column) labeled `Intercept' corresponds to the first
level of \Robject{grp}, i.e., `untreated'. The coefficient
`grptreated' represents the deviation of the treated gropu from
untreated. Eventually, we will test whether the second coefficient is
significantly different from zero, i.e., whether samples with a `1' in
the second column are, on average different from samples with a `0'.
On the one hand, use of \Rfunction{model.matrix} to specify
experimental design implies that the user is comfortable with
something more than elementary statistical concepts, while on the
other it provides great flexibility in the type of experiment that can
be analyzed.

\paragraph{Negative binomial error}

RNA-seq count data are often described by a negative binomial
model. This model includes a `dispersion' parameter that describes
biological variation beyond the expectation under a Poisson model. The
simplest approach estimates a disperssion parameter from all the
data. The estimate needs to be conducted in the context of the
experimental design, so that variability between experimental factors
is not mistaken for variability in counts. The square root of the
estimated dispersion represents the coefficient of variation between
biological samples.
<<common.dispersion>>=
dge <- estimateCommonDisp(dge, design)
sqrt(dge$common.dispersion)
@ 

This approach assumes that a common dispersion parameter is shared by
all genes. A different approach, appropriate when there are more
samples in the study, is to estimate a dispersion parameter that is
specific to each tag (using \Rfunction{estimateTagwiseDisp} in the
\Biocpkg{edgeR} package). As another alternative, Anders and Huber
\cite{AndersHuber2010} note that dispersion increases as the mean
number of reads per gene decreases. One can estimate the relationship
between dispersion and mean using \Rfunction{estimateGLMTrendedDisp}
in \Rpackage{edgeR}, using a fitted relationship across all genes to
estimate the dispersion of individual genes. Because in our case
sample sizes (biological replicates) are small, gene-wise estimates of
dispersion are likely imprecise. One approach is to moderate these
estimates by calculating a weighted average of the gene-specific and
common dispersion; \Rfunction{estimateGLMTagwiseDisp} performs this
calculation, requiring that the user provide the weight \emph{a
  priori}.

\paragraph{Differential representation}

The final steps in estimating differential representation are to fit
the full model; to perform the likelihood ratio test comparing the
full model to a model in which one of the coefficients has been
removed; and to summarize, from the likelihood ratio calculation,
genes that are most differentially represented. The result is a `top
table' whose row names are the Flybase gene ids used to label the
elements of the \Robject{ex} \Rclass{GRangesList}.

\begin{Exercise}
  Use \Rfunction{glmFit} to fit the general linear model. This
  function requires the input data \Robject{dge}, the experimental
  design \Robject{design}, and as estimate of dispersion.
  
  Use \Rfunction{glmLRT} to form the likelihood ratio test. This
  requires the original data \Robject{dge} and the fitted model from
  the previous part of this question. Which coefficient of the design
  matrix do you wish to test?
  
  Finally, create a `top table' of differentially represented genes
  using \Rfunction{topTags}.
\end{Exercise}

\begin{Solution}
  Here we fit a glm to our data and experimental design, using the
  common dispersion estimate.
<<glmFit>>=
fit <- glmFit(dge, design, dispersion=dge$common.dispersion)
@ 

The fit can be used to calculate a likelihood ratio test, comparing
the full model to a reduced version with the second coefficient
removed. The second coefficient captures the difference between
treated and untreated groups, and the likelihood ratio test asks
whether this term contributes meaningfully to the overall fit.
<<lrt>>=
lrTest <- glmLRT(dge, fit, coef=2)
@ 

Here \Rfunction{topTags} function summarizes results across the experiment.
<<topTags>>=
(tt <- topTags(lrTest))
@ 

As a 'sanity check', summarize the original data for the first several
probes
<<sanity>>=
sapply(rownames(tt$table)[1:4], 
       function(x) tapply(counts[x,], grp, mean))
@ 
\end{Solution}

%% \subsection{Gene set enrichment}

%% \Biocpkg{goseq}

%% issues
