%\VignetteDepends{DESeq}
%\VignetteIndexEntry{Differential Representation Lab}
%\VignettePackage{SeattleIntro2010}

\documentclass{article}

\usepackage{BioconductorVignette}
\SweaveOpts{keep.source=TRUE,width=9,height=8,eps=FALSE,include=FALSE}

\title{Differential Representation}
\author{Martin Morgan \\ Fred Hutchinson Cancer Research Center}
\date{9-10 December, 2010}


\begin{document}
\maketitle
<<setup, echo=FALSE>>=
options(width = 72)
library(GenomicRanges)
library(lattice)
jnk <- capture.output(library(DESeq))
@ 

This lab illustrates how \Bioconductor{} tools can be used to
investigate differential respresentation. The data used is derived
from Nagalakshmi et al.~\cite{pmid18451266}. This publication defined
the term `RNA-seq' to describe the use of high-througput sequence
methodolgies for assessing differential mRNA abundance. The data were
collected to explore technical, rather than biological, issues.  The
authors developed two different methods, \emph{random hexamer} and
\emph{oligo(dT)}, to generate DNA from poly(A) RNA. Each method had an
`original', `biological', and `technical' replicate. Samples were run
on an early-generation Solexa GAI; the reads are relatively short
(33bp) and relatively limited in number (3.4-6.8 million reads per
sample). Samples were from yeast strain BY4741; an interesting feature
of this strain is that it has deletions in leu2, ura3, met15, and
his3.

\section{Exploration}

\begin{Exercise}
  Review the pre-processing steps used to generate the counts present
  in \Robject{hitspergene}. These steps are summarized in the
  \emph{Appendix} vignette available in the
  \Rpackage{SeattleIntro2010} package.
\end{Exercise}

\begin{Solution}
<<browseVignette, eval=FALSE>>=
browseVignette("SeattleIntro2010")
@ 
\end{Solution}

\begin{Exercise}
  Load the \Rpackage{SeattleIntro2010} package and use the
  \Rfunction{data} function to load the \Robject{hitspergene}
  object. Display it on the screen. Use \Rfunction{class} to determine
  the class of \Robject{hitspergene}. Use the help system to
  understand the class, its accessors, and how to manipulate it. In
  many respects, the object behaves much like a standard \R{}
  \Robject{data.frame}.
\end{Exercise}

\begin{Solution}
<<hitspergene>>=
library(SeattleIntro2010)
data(hitspergene)
hitspergene[1:5,]
class(hitspergene)
@ 
\end{Solution}

\begin{Exercise}
  It is sometimes necessary to convert objects from one representation
  to another. For instance, the \Rpackage{lattice} function
  \Rfunction{splom} prints a `scatter plot matrix', a matrix of
  scatter plots between the columns of a \Rclass{data.frame}. Often,
  it is possible to coerce an object from one type to another using
  functions named like \Rfunction{as.data.frame} or
  \Rcode{as(hitspergene, "data.frame")}. Coerce \Robject{hitspergene}
  to a \Robject{data.frame}, and plot the data frame using
  \Rfunction{splom}. Use the \Rfunarg{pch} argument of
  \Rfunction{splom} to change the character used to plot. Transform
  all the values in \Robject{hitspergene} using the \Rfunction{asinh}
  function, and plot those. 
\end{Exercise}

\begin{Solution}
<<splom1, fig=TRUE>>=
library(lattice)
df1 <- as.data.frame(hitspergene)
splom(df1, pch=".")
@ 
<<splom2, fig=TRUE>>=
df2 <- asinh(df1)
splom(df2, pch=".")
@ 
\end{Solution}

\begin{Exercise}
  One feature of the \Rclass{DataFrame} class is that the elements
  (columns) can retain `metadata' describing their content. Extract
  the element metadata to discover which smaples correspond to each
  methodology and replicate. Subset \Robject{hitspergene} so that it
  contains only columns corresponding to the `Original' and
  `Biological' replicate. 
  
  Now remove rows for which none of the samples had any hits. Do this
  by converting \Robject{hitspergene} to a \Rclass{data.frame} using
  \Rfunction{as.data.frame}, and using \Rfunction{rowSums} to sum up
  the number of hits per row. Take a quick look at the distribution of
  the number of hits per row using the \Rfunction{histogram} function,
  perhaps transforming the data (e.g., using the \Rfunction{asinh}
  function). Finally, remove rows with zero counts using logical
  subsetting. How many rows have been removed
\end{Exercise}

\begin{Solution}
<<hitspergene-reduce, fig=TRUE>>=
elementMetadata(hitspergene)
biolReps <- elementMetadata(hitspergene)$Replicate != "Technical"
count <- hitspergene[, biolReps]

df1 <- as.data.frame(count)
hitsPerRow <- rowSums(df1)
histogram(asinh(hitsPerRow))

count <- count[hitsPerRow != 0,]
dim(count)
sum(hitsPerRow==0)  ## number of rows removed
@ 
\end{Solution}

\section{Differential Representation}

\begin{Exercise}
We'll look at differential representation using the \Biocpkg{DESeq}
package~\cite{pmid20979621}. Load the package and take a few minutes
to browse its vignette.
\end{Exercise}
\begin{Solution}
<<DESeq, eval=FALSE>>=
library(DESeq)
browseVignettes("DESeq")
@ 
\end{Solution}

We'll take a `fast track' through the analysis. There are four steps:
\begin{enumerate}
\item Convert our data into a format understood by \Rpackage{DESeq}
\item Calculate scale factors so that libraries can be treated as
  though they have comparable numbers of reads.
\item Estimate functions to describe the relationship between mean and
  variance
\item Evaluate significance of differences between the random hexamer
  and oligo(dT) methods.
\end{enumerate}

\begin{Exercise}
  The \Rpackage{DESeq} package requires two pieces of information: a
  \Rclass{data.frame} of counts, and a vector describing the
  treatments to which each column of the count data frame
  belongs. Extract this information from the \Robject{count} object
  created above, and create a \Rclass{CountDataSet} using
  \Rfunction{newCountDataSet}. A slight transformation is to convert
  the \Rclass{numeric} data in \Robject{counts} to \Rclass{integer},
  rounding up using the \Rfunction{ceiling} function.
\end{Exercise}
\begin{Solution}
<<newCountDataSet>>=
df1 <- ceiling(as.data.frame(count))
sample <- elementMetadata(count)$Sample
cds <- newCountDataSet(df1, sample)
@ 
\end{Solution}

\begin{Exercise}
  Estimate the relative contributions of each sample with
  \Rfunction{estimateSizeFactors}. Read the help page for this
  function (with \Rcode{?estimateSizeFactors}) to understand what it
  is doing. Extract the estimated size factors from the result with
  the \Rfunction{sizeFactors} function; how do the size factors
  compare with, say, a naive weight based on reads per sample? What
  factors are likely to contribute to this difference?
\end{Exercise}
\begin{Solution}
<<estimateSizeFactors>>=
cds <- estimateSizeFactors(cds)
sizeFactors(cds)
colSums(df1) / mean(colSums(df1))
@ 
\end{Solution}

\begin{Exercise}
  Estimate the variance functions using
  \Rfunction{estimateVarianceFunctions}; see its help page for details
\end{Exercise}
\begin{Solution}
<<estimateVarianceFunctions>>=
cds <- estimateVarianceFunctions(cds)
@ 
\end{Solution}

\begin{Exercise}
  Finally, perform a test of signficance between groups using
  \Rfunction{nbinomTest}. Look at the first few rows of the result
  using \Rfunction{head}, and compare what you see with the
  description of return value on the \Rfunction{nbinomTest} help
  page. Use \Rfunction{order} to determine the row indexes that will
  place the rows of the result into increasing order, based on the
  adjusted p-value, \Robject{padj}.
\end{Exercise}
\begin{Solution}
<<topTable>>=
nbTopTable <- nbinomTest(cds, "dT", "RH")
head(nbTopTable, 3)
o <- order(nbTopTable$padj)
head(nbTopTable[o,], 3)
@ 
\end{Solution}

\section{Evaluation}

This section is self-directed, asking you to use the facilities of
\Rpackage{DESeq} to assess whether the steps we performed lead to
reasonable results. In particular:

\begin{Exercise}
  Follow the vignette \emph{Analysing RNA-Seq data with the "DESeq"
    package} in the \Rpackage{DESeq} package to assess whether
  estimated variance functions are reasonable. Do this using the
  \Rfunction{scvPlot} function. Note that `shot noise' dominates at
  low count number; what are the consequences of this for the results
  that can be inferred? Further explore your fit as described in the
  \Biocpkg{DESeq} vignette.
\end{Exercise}

\begin{Exercise}
  Follow the vignette to explore your results. Plot the log2 fold
  change against the base mean, analogous to the \Rfunction{plotDE}
  function defined in the vignette. Identify the genes with largest
  and smallest fold change. 
\end{Exercise}

\section{Resources}

\bibliographystyle{plain}
\bibliography{SeattleIntro2010}

\end{document}
