%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{EDASeq: Exploratory Data Analysis and Normalization for RNA-Seq data}
%\VignetteDepends{EDASeq}
%\VignetteDepends{yeastRNASeq}
%\VignetteDepends{edgeR}
%\VignetteDepends{DESeq}
%\VignettePackage{EDASeq}
\documentclass{article}

<<style-knitr, eval=TRUE, echo=FALSE, results="asis">>=
BiocStyle::latex()
@

\usepackage{times}
\usepackage[numbers,sort&compress]{natbib}
\usepackage{subfig}
\usepackage{amsmath}


\title{\Biocpkg{EDASeq}: Exploratory Data Analysis and Normalization for RNA-Seq}
\author{Davide Risso}
\date{Modified: March 25, 2014.  Compiled: \today}
\begin{document}

\maketitle

\tableofcontents

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}

In this document, we show how to conduct Exploratory Data Analysis
(EDA) and normalization for a typical RNA-Seq experiment using the
package \Biocpkg{EDASeq}.

One can think of EDA for RNA-Seq as a two-step process: ``read-level''
EDA helps in discovering lanes with low sequencing depths, quality
issues, and unusual nucleotide frequencies, while ``gene-level'' EDA
can capture mislabeled lanes, issues with distributional assumptions
(e.g., over-dispersion), and GC-content bias.

The package also implements both ``within-lane'' and ``between-lane''
normalization procedures, to account, respectively, for within-lane
gene-specific (and possibly lane-specific) effects on read counts
(e.g., related to gene length or GC-content) and for between-lane
distributional differences in read counts (e.g., sequencing depths).

To illustrate the functionality of the \Biocpkg{EDASeq} package, we
make use of the {\it Saccharomyces cerevisiae} RNA-Seq data from
\citet{lee2008novel}. Briefly, a wild-type strain and three mutant
strains were sequenced using the Solexa 1G Genome Analyzer.  For each
strain, there are four technical replicate lanes from the same library
preparation. The reads were aligned using \software{Bowtie}
\citep{langmead2009ultrafast}, with unique mapping and allowing up to
two mismatches.

The \Biocexptpkg{leeBamViews} package provides a subset of the aligned
reads in BAM format. In particular, only the reads mapped between
bases 800,000 and 900,000 of chromosome XIII are considered. We use
these reads to illustrate read-level EDA.

The \Biocexptpkg{yeastRNASeq} package contains gene-level read counts for
four lanes: two replicates of the wild-type strain (``wt'') and
two replicates of one of the mutant strains (``mut'').  We use
these data to illustrate gene-level EDA.

<<setup, echo=FALSE>>=
library(knitr)
opts_chunk$set(dev="pdf", fig.align="center", cache=FALSE, message=FALSE, out.width=".55\\textwidth", echo=TRUE, results="markup", fig.show="hold")
options(width=60)
@

<<data>>=
require(EDASeq)
require(yeastRNASeq)
require(leeBamViews)
@ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Reading in unaligned and aligned read data}
\label{secRead}

\paragraph{Unaligned reads.}

Unaligned (unmapped) reads stored in FASTQ format may be managed via
the class \Rclass{FastqFileList} imported from \Biocpkg{ShortRead}.
Information related to the libraries sequenced in each lane can be
stored in the \Rcode{elementMetadata} slot of the
\Rclass{FastqFileList} object.

<<import-raw>>=
files <- list.files(file.path(system.file(package = "yeastRNASeq"), 
                              "reads"), pattern = "fastq", full.names = TRUE)
names(files) <- gsub("\\.fastq.*", "", basename(files))
met <- DataFrame(conditions=c(rep("mut",2), rep("wt",2)), 
                 row.names=names(files))
fastq <- FastqFileList(files)
elementMetadata(fastq) <- met
fastq
@ 

\paragraph{Aligned reads.}

The package can deal with aligned (mapped) reads in BAM format, using
the class \Rclass{BamFileList} from \Biocpkg{Rsamtools}.  Again, the
\Rcode{elementMetadata} slot can be used to store lane-level sample
information.

<<import-aligned>>=
files <- list.files(file.path(system.file(package = "leeBamViews"), "bam"),
                    pattern = "bam$", full.names = TRUE)
names(files) <- gsub("\\.bam", "", basename(files))

gt <- gsub(".*/", "", files)
gt <- gsub("_.*", "", gt)
lane <- gsub(".*(.)$", "\\1", gt)
geno <- gsub(".$", "", gt)

pd <- DataFrame(geno=geno, lane=lane,
                row.names=paste(geno,lane,sep="."))

bfs <- BamFileList(files)
elementMetadata(bfs) <- pd
bfs
@ 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Read-level EDA}


\paragraph{Numbers of unaligned and aligned reads.}

One important check for quality control is to look at the total number
of reads produced in each lane, the number and the percentage of reads mapped to a
reference genome. A low total
number of reads might be a symptom of low quality of the input RNA,
while a low mapping percentage might indicate poor quality of the
reads (low complexity), problems with the reference genome, or
mislabeled lanes.

<<plot-total, fig.cap="Per-lane number of mapped reads and quality scores", fig.subcap=c("Number of mapped reads", "Mean per-base quality of mapped reads"), out.width='.49\\linewidth'>>=
colors <- c(rep(rgb(1,0,0,alpha=0.7),2), 
            rep(rgb(0,0,1,alpha=0.7),2),
            rep(rgb(0,1,0,alpha=0.7),2),
            rep(rgb(0,1,1,alpha=0.7),2))
barplot(bfs,las=2,col=colors)

plotQuality(bfs,col=colors,lty=1)
legend("topright",unique(elementMetadata(bfs)[,1]), fill=unique(colors))
@ 

Figure \ref{fig:plot-total}a, produced using the \Rcode{barplot} method for the
\Rclass{BamFileList} class, displays the number of mapped reads for
the subset of the yeast dataset included in the package
\Biocpkg{leeBamViews}.  Unfortunately, \Biocpkg{leeBamViews} does
not provide unaligned reads, but barplots of the total number of reads
can be obtained using the \Rcode{barplot} method for the
\Biocpkg{FastqFileList} class. Analogously, one can plot the percentage
of mapped reads with the \Rcode{plot} method with signature
\texttt{c(x="BamFileList", y="FastqFileList")}. See the manual pages for
details.

\paragraph{Read quality scores.}

As an additional quality check, one can plot the mean per-base (i.e.,
per-cycle) quality of the unmapped or mapped reads in every lane
(Figure \ref{fig:plot-total}b).

\paragraph{Individual lane summaries.}

If one is interested in looking more thoroughly at one lane, it is
possible to display the per-base distribution of quality scores for
each lane (Figure \ref{fig:plot-qual}a) and the number of mapped reads
stratified by chromosome (Figure \ref{fig:plot-qual}b) or strand. As expected,
all the reads are mapped to chromosome XIII.

<<plot-qual, fig.cap="Quality scores and number of mapped reads for lane 'isowt5\\_13e.'", fig.subcap=c("Per-base quality of mapped reads", "Number of mapped reads per-chromosome"), out.width='.49\\linewidth'>>=
plotQuality(bfs[[1]],cex.axis=.8)
barplot(bfs[[1]],las=2)
@ 

\paragraph{Read nucleotide distributions.}

A potential source of bias is related to the sequence composition of
the reads. The function \mbox{\Rfunction{plotNtFrequency}} plots the
per-base nucleotide frequencies for all the reads in a given
lane (Figure \ref{fig:plot-nt}).

<<plot-nt, fig.cap="Per-base nucleotide frequencies of mapped reads for lane 'isowt5\\_13e.'", out.width='.49\\linewidth'>>=
plotNtFrequency(bfs[[1]])
@ 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Gene-level EDA}

Examining statistics and quality metrics at a read level can help in
discovering problematic libraries or systematic biases in one or more
lanes. Nevertheless, some biases can be difficult to detect at this
scale and gene-level EDA is equally important.

\paragraph{Classes and methods for gene-level counts.}

There are several Bioconductor packages for aggregating reads over
genes (or other genomic regions, such as, transcripts and exons) given
a particular genome annotation, e.g., \Biocpkg{IRanges},
\Biocpkg{ShortRead}, \Biocpkg{Genominator}, \Biocpkg{Rsubread}. See
their respective vignettes for details.

Here, we consider this step done and load the object
\Robject{geneLevelData} from \Biocexptpkg{yeastRNASeq}, which provides
gene-level counts for 2 wild-type and 2 mutant lanes from the yeast
dataset of \citet{lee2008novel} (see the \Biocpkg{Genominator}
vignette for an example on the same dataset).

<<load-gene-level>>=
data(geneLevelData)
head(geneLevelData)
@ 

Since it is useful to explore biases related to length and GC-content,
the \Biocpkg{EDASeq} package provides, for illustration purposes,
length and GC-content for \emph{S. cerevisiae} genes (based on SGD
annotation, version r64 \citep{sgd}).

<<load-lgc>>=
data(yeastGC)
head(yeastGC)
data(yeastLength)
head(yeastLength)
@ 

First, we filter the non-expressed genes, i.e., we consider only the
genes with an average read count greater than 10 across the four lanes
and for which we have length and GC-content information.

<<filter>>=
filter <- apply(geneLevelData,1,function(x) mean(x)>10)
table(filter)
common <- intersect(names(yeastGC),
                    rownames(geneLevelData[filter,])) 
length(common)
@ 

This leaves us with \Sexpr{length(common)} genes.

The \Biocpkg{EDASeq} package provides the \Rclass{SeqExpressionSet}
class to store gene counts, (lane-level) information on the sequenced
libraries, and (gene-level) feature information.  We use the data
frame \Robject{met} created in Section \ref{secRead} for the
lane-level data.  As for the feature data, we use gene length and
GC-content.

<<create-object>>=
feature <- data.frame(gc=yeastGC,length=yeastLength)
data <- newSeqExpressionSet(counts=as.matrix(geneLevelData[common,]),
                            featureData=feature[common,],
                            phenoData=data.frame(
                              conditions=c(rep("mut",2),rep("wt",2)),
                              row.names=colnames(geneLevelData)))
data
@ 

Note that the row names of \Rcode{counts} and \Rcode{featureData}
must be the same; likewise for the row names of \Rcode{phenoData}
and the column names of \Rcode{counts}.  As in the
\Rclass{CountDataSet} class, the expression values can be accessed
with \Rfunction{counts}, the lane information with \Rfunction{pData},
and the feature information with \Rfunction{fData}.

<<show-data>>=
head(counts(data))
pData(data)
head(fData(data))
@ 

The \Rclass{SeqExpressionSet} class has two additional slots:
\Rcode{normalizedCounts} and \Robject{offset} (matrices of the same dimension as \Rcode{counts}),
which may be used to store a matrix of normalized counts and of
normalization offsets, respectively, to be used for subsequent analyses (see Section
\ref{secDE} and the \Biocpkg{edgeR} vignette for details on the
role of offsets). If not specified, the offset is initialized as a matrix
of zeros.

<<show-offset>>=
head(offst(data))
@ 

\paragraph{Between-lane distribution of gene-level counts.}

One of the main considerations when dealing with gene-level counts is
the difference in count distributions between lanes. The
\Rfunction{boxplot} method provides an easy way to produce boxplots of
the logarithms of the gene counts in each lane (Figure \ref{fig:boxplot-genelevel}).

<<boxplot-genelevel, fig.cap="Between-lane distribution of gene-level counts (log).">>=
boxplot(data,col=colors[1:4])
@ 

The \Rfunction{MDPlot} method produces a mean-difference plot (MD-plot)
of read counts for two lanes (Figure \ref{fig:md-plot}).

<<md-plot, fig.cap="Mean-difference plot of the gene-level counts (log) of lanes 'mut\\_1' and 'wt\\_1.'">>=
MDPlot(data,c(1,3))
@ 

\paragraph{Over-dispersion.}

Although the Poisson distribution is a natural and simple way to model
count data, it has the limitation of assuming equality of the mean and
variance. For this reason, the negative binomial distribution has been
proposed as an alternative when the data show over-dispersion. The
function \Rfunction{meanVarPlot} can be used to check whether the
count data are over-dispersed (for the Poisson distribution, one would
expect the points in Figures \ref{fig:plot-mean-var} to be evenly
scattered around the black line).

<<plot-mean-var, fig.cap="Mean-variance relationship for the two mutant lanes and all four lanes: the black line corresponds to the Poisson distribution (variance equal to the mean), while the red curve is a lowess fit.", fig.subcap=c("Mutant lanes", "All four lanes"), out.width='.49\\linewidth'>>=
meanVarPlot(data[,1:2], log=TRUE, ylim=c(0,16))
meanVarPlot(data, log=TRUE, ylim=c(0,16))
@ 

Note that the mean-variance relationship should be examined within
replicate lanes only (i.e., conditional on variables expected to
contribute to differential expression). For the yeast dataset, it is
not surprising to see no evidence of over-dispersion for the two
mutant technical replicate lanes (Figure \ref{fig:plot-mean-var}a); likewise for
the two wild-type lanes. However, one expects over-dispersion in the
presence of biological variability, as seen in Figure
\ref{fig:plot-mean-var}b when considering at once all four mutant and wild-type lanes
\citep{anders2010differential,bullard2010evaluation,robinson2010edger}.


\paragraph{Gene-specific effects on read counts.}

Several authors have reported selection biases related to sequence
features such as gene length, GC-content, and mappability
\citep{bullard2010evaluation,hansen2011removing,oshlack2009transcript,risso2011gc}.

In Figure \ref{fig:plot-gc}, obtained using \Rfunction{biasPlot}, one can
see the dependence of gene-level counts on GC-content. The same plot
could be created for gene length or mappability instead of GC-content.

<<plot-gc, fig.cap="Lowess regression of the gene-level counts (log) on GC-content for each lane, color-coded by experimental condition.">>=
biasPlot(data, "gc", log=TRUE, ylim=c(1,5))
@ 

To show that GC-content dependence can bias differential expression
analysis, one can produce stratified boxplots of the log-fold-change
of read counts from two lanes using the \Rfunction{biasBoxplot} method
(Figure \ref{fig:boxplot-gc}). Again, the same type of plots can be created
for gene length or mappability.

<<boxplot-gc, fig.cap="Boxplots of the log-fold-change between 'mut\\_1' and 'wt\\_1' lanes stratified by GC-content.">>=
lfc <- log(counts(data)[,3]+0.1) - log(counts(data)[,1]+0.1)
biasBoxplot(lfc, fData(data)$gc)
@ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Normalization}

Following \citet{risso2011gc}, we consider two main types of effects
on gene-level counts: (1) within-lane gene-specific (and possibly
lane-specific) effects, e.g., related to gene length or GC-content,
and (2) effects related to between-lane distributional differences,
e.g., sequencing depth. Accordingly,
\Rfunction{withinLaneNormalization} and
\Rfunction{betweenLaneNormalization} adjust for the first and second
type of effects, respectively.  We recommend to normalize for
within-lane effects prior to between-lane normalization.

We implemented four within-lane normalization methods, namely: loess
robust local regression of read counts (log) on a gene feature such as
GC-content (\Rcode{loess}), global-scaling between feature strata
using the median (\Rcode{median}), global-scaling between feature
strata using the upper-quartile (\Rcode{upper}), and full-quantile
normalization between feature strata (\Rcode{full}). For a discussion
of these methods in context of GC-content normalization see
\citet{risso2011gc}.

<<normalization>>=
dataWithin <- withinLaneNormalization(data,"gc", which="full")
dataNorm <- betweenLaneNormalization(dataWithin, which="full")
@ 

Regarding between-lane normalization, the package implements three of
the methods introduced in \citet{bullard2010evaluation}:
global-scaling using the median (\Rcode{median}), global-scaling using
the upper-quartile (\Rcode{upper}), and full-quantile normalization
(\Rcode{full}).

Figure \ref{fig:plot-gc-norm} shows how after full-quantile within- and
between-lane normalization, the GC-content bias is reduced and the
distribution of the counts is the same in each lane.

<<plot-gc-norm, fig.cap="Full-quantile within- and between-lane normalization. (a) Lowess regression of normalized gene-level counts (log) on GC-content for each lane. (b) Between-lane distribution of normalized gene-level counts (log).", fig.subcap=c("GC-content", "Count distribution"), out.width='.49\\linewidth'>>=
biasPlot(dataNorm, "gc", log=TRUE, ylim=c(1,5))
boxplot(dataNorm, col=colors)
@ 

\paragraph{Offset.}
\label{secDE}

Some authors have argued that it is better to leave the count data
unchanged to preserve their sampling properties and instead use an
offset for normalization purposes in the statistical model for read
counts
\citep{anders2010differential,hansen2011removing,robinson2010edger}. This
can be achieved easily using the argument \Rcode{offset} in both
normalization functions.

<<norm-offset>>=
dataOffset <- withinLaneNormalization(data,"gc",
                                      which="full",offset=TRUE)
dataOffset <- betweenLaneNormalization(dataOffset,
                                       which="full",offset=TRUE)
@ 

Note that the \Robject{dataOffset} object will have both normalized
counts and offset stored in their respective slots.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Differential expression analysis}

One of the main applications of RNA-Seq is differential expression
analysis.  The normalized counts (or the original counts and the
offset) obtained using the \Biocpkg{EDASeq} package can be supplied
to packages such as \Biocpkg{edgeR} \citep{robinson2010edger} or
\Biocpkg{DESeq} \citep{anders2010differential} to find differentially
expressed genes. This section should be considered only as an
illustration of the compatibility of the results of \Biocpkg{EDASeq}
with two of the most widely used packages for differential expression;
our aim is not to compare differential expression strategies (e.g.,
normalized counts vs. offset).

\subsection{edgeR}\label{edger}

We can perform a differential expression analysis with
\Biocpkg{edgeR} based on the original counts by passing an offset to
the generalized linear model. For simplicity, we estimate a common dispersion
parameter for all genes. See the \Biocpkg{edgeR} vignette for details
about how to perform a differential expression analysis using a
gene-specific dispersion or more complex designs.

<<edger>>=
library(edgeR)
design <- model.matrix(~conditions, data=pData(dataOffset))
disp <- estimateGLMCommonDisp(counts(dataOffset), 
                              design, offset=-offst(dataOffset))

fit <- glmFit(counts(dataOffset), design, disp, offset=-offst(dataOffset))

lrt <- glmLRT(fit, coef=2)
topTags(lrt)
@ 

\subsection{DESeq}

We can perform a differential expression analysis with
\Biocpkg{DESeq} based on the normalized counts by using the coerce
method from the \Rclass{SeqExpressionSet} class to the
\Rclass{CountDataSet} class of \Biocpkg{DESeq}.  When working with
data that have been normalized for both within- and between-lane
effects, we force the size factors to be one, since differences in
lane sequencing depths have already been accounted for in our
between-lane normalization. One could also consider only within-lane
normalization and account for differences in sequencing depth by
estimating the size factors using \Rpackage{DESeq}.

<<deseq>>=
library(DESeq)
counts <- as(dataNorm, "CountDataSet")
sizeFactors(counts) <- rep(1,4)
counts <- estimateDispersions(counts)
res <- nbinomTest(counts, "wt", "mut")
head(res)
@ 

\section{Definitions and conventions}

\subsection{Rounding}
After either within-lane or between-lane normalization, the expression
values are not counts
anymore. However, their distribution still shows some typical features
of counts distribution (e.g., the variance depends on the mean).
Hence, for most applications, it is useful to round the normalized
values to recover count-like values, which we refer to as ``pseudo-counts''. 

By default, both
\Rfunction{withinLaneNormalization} and
\Rfunction{betweenLaneNormalization} round the normalized values to
the closest integer. This behavior can be changed by specifying
\texttt{round=FALSE}. This gives the user more flexibility and assures
that rounding approximations do not affect subsequent computations
(e.g., recovering the offset from the normalized counts).

\subsection{Zero counts}
To avoid problems in the computation of logarithms (e.g. in log-fold-changes), we add a small positive constant (namely $0.1$) to the
counts. For instance, the log-fold-change between $y_1$ and $y_2$ is
defined as 
\begin{equation*}
  \frac{\log(y_1 + 0.1)}{\log(y_2 + 0.1)}.
\end{equation*}

\subsection{Offset}
We define an offset in the normalization as
\begin{equation*}
  o = \log(y_{norm} + 0.1) - \log(y_{raw} + 0.1),
\end{equation*}
where $y_{norm}$ and $y_{raw}$ are the normalized and raw counts, respectively.

One can easily recover the normalized data from the raw counts and
offset, as shown here:

<<unrounded>>=
dataNorm <- betweenLaneNormalization(data, round=FALSE, offset=TRUE)

norm1 <- normCounts(dataNorm)
norm2 <- exp(log(counts(dataNorm) + 0.1 ) + offst(dataNorm)) - 0.1

head(norm1 - norm2)
@ 

Note that the small constant added in the definition of offset does
not matter when pseudo-counts are considered, i.e.,

<<rounded>>= 
head(round(normCounts(dataNorm)) - round(counts(dataNorm) * exp(offst(dataNorm))))
@ 

We defined the offset as the log-ratio between normalized and raw
counts. However, the \Biocpkg{edgeR} functions expect as offset
argument the log-ratio between raw and normalized counts. One must use
\texttt{-offst(offsetData)} as the offset argument of \Biocpkg{edgeR}
(see Section \ref{edger}).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{SessionInfo}

<<sessionInfo, results="asis">>=
toLatex(sessionInfo())
@ 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{thebibliography}{8}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{doi: #1}\else
  \providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi

\bibitem[Lee et~al.(2008)Lee, Hansen, Bullard, Dudoit, and
  Sherlock]{lee2008novel}
A.~Lee, K.D. Hansen, J.~Bullard, S.~Dudoit, and G.~Sherlock.
\newblock Novel low abundance and transient {RNAs} in yeast revealed by tiling
  microarrays and ultra high-throughput sequencing are not conserved across
  closely related yeast species.
\newblock \emph{PLoS Genet}, 4\penalty0 (12):\penalty0 e1000299, 2008.

\bibitem[Langmead et~al.(2009)Langmead, Trapnell, Pop, and
  Salzberg]{langmead2009ultrafast}
B.~Langmead, C.~Trapnell, M.~Pop, and S.L. Salzberg.
\newblock Ultrafast and memory-efficient alignment of short {DNA} sequences to
  the human genome.
\newblock \emph{Genome Biol}, 10\penalty0 (3):\penalty0 R25, 2009.

\bibitem[{Saccharomyces Genome Database}(r64)]{sgd}
{Saccharomyces Genome Database}.
\newblock {\url{http://www.yeastgenome.org}}, r64.

\bibitem[Bullard et~al.(2010)Bullard, Purdom, Hansen, and
  Dudoit]{bullard2010evaluation}
J.H. Bullard, E.~Purdom, K.D. Hansen, and S.~Dudoit.
\newblock Evaluation of statistical methods for normalization and differential
  expression in {mRNA-Seq} experiments.
\newblock \emph{BMC bioinformatics}, 11\penalty0 (1):\penalty0 94, 2010.

\bibitem[Robinson et~al.(2010)Robinson, McCarthy, and Smyth]{robinson2010edger}
M.D. Robinson, D.J. McCarthy, and G.K. Smyth.
\newblock {edgeR}: a {Bioconductor} package for differential expression
  analysis of digital gene expression data.
\newblock \emph{Bioinformatics}, 26\penalty0 (1):\penalty0 139, 2010.

\bibitem[Anders and Huber(2010)]{anders2010differential}
S.~Anders and W.~Huber.
\newblock Differential expression analysis for sequence count data.
\newblock \emph{Genome Biology}, 11\penalty0 (10):\penalty0 R106, 2010.

\bibitem[Oshlack and Wakefield(2009)]{oshlack2009transcript}
A.~Oshlack and M.J. Wakefield.
\newblock Transcript length bias in RNA-seq data confounds systems biology.
\newblock \emph{Biology Direct}, 4\penalty0 (1):\penalty0 14, 2009.

\bibitem[Risso et~al.(2011)Risso, Schwartz, Sherlock, and
  Dudoit]{risso2011gc} D.~Risso, K.~Schwartz, G.~Sherlock, and
  S.~Dudoit.  \newblock {GC-Content Normalization for RNA-Seq Data}.
  \newblock Technical report \#291, University of California,
  Berkeley, Division of Biostatistics, 2011.  \newblock
  {\url{http://www.bepress.com/ucbbiostat/paper291/}}.

\bibitem[Hansen et~al.(2011)Hansen, Irizarry, and
  Wu]{hansen2011removing} K.D. Hansen, R.A. Irizarry, and Z.~Wu.
  \newblock Removing technical variability in {RNA-Seq} data using
  conditional quantile normalization.  \newblock Technical report
  \#227, Johns Hopkins University, Dept. of Biostatistics Working
  Papers, 2011.  \newblock
  {\url{http://www.bepress.com/jhubiostat/paper227/}}.

\end{thebibliography}


\end{document}



% LocalWords:  EDASeq GC Saccharomyces cerevisiae Solexa langmead ultrafast BAM
% LocalWords:  leeBamViews yeastRNASeq mut FASTQ FastqFileList Rsamtools fastq
% LocalWords:  elementMetadata BamFileList gsub basename DataFrame bam geno sep
% LocalWords:  bfs barplot las dataset barplots qual plotQuality lty topright
% LocalWords:  isowt plotNtFrequency nt Bioconductor exons IRanges ShortRead gc
% LocalWords:  Genominator Rsubread geneLevelData SGD sgd lgc yeastGC rownames
% LocalWords:  yeastLength SeqExpressionSet newSeqExpressionSet exprs phenoData
% LocalWords:  featureData colnames ExpressionSet pData fData edgeR offst md
% LocalWords:  boxplot boxplots genelevel MDPlot meanVarPlot  anders bullard
% LocalWords:  robinson lowess mappability hansen oshlack risso biasPlot ylim
% LocalWords:  biasBoxplot lfc withinLaneNormalization betweenLaneNormalization
% LocalWords:  loess quantile dataWithin dataNorm dataOffset DESeq disp glmFit
% LocalWords:  estimateGLMCommonDisp lrt glmLRT coef CountDataSet deseq eval
% LocalWords:  sizeFactors estimateDispersions nbinomTest SessionInfo tex doi
% LocalWords:  sessionInfo toLatex urlstyle Dudoit RNAs microarrays PLoS 
% LocalWords:  Trapnell Salzberg Purdom mRNA BMC bioinformatics Smyth Irizarry
% LocalWords:  Biostatistics
