%\VignetteIndexEntry{Differential expression for RNA-seq data with dispersion shrinkage}
%\VignettePackage{DSS}


\documentclass{article}
<<style, eval=TRUE, echo=FALSE, results=tex>>=
    BiocStyle::latex(use.unsrturl=FALSE)
@

\renewcommand{\baselinestretch}{1.3}
\usepackage{Sweave}
\SweaveOpts{keep.source=TRUE,eps=FALSE,include=TRUE,width=4,height=4}
%\newcommand{\Robject}[1]{\texttt{#1}}
%\newcommand{\Rpackage}[1]{\textit{#1}}
%\newcommand{\Rclass}[1]{\textit{#1}}
%\newcommand{\Rfunction}[1]{{\small\texttt{#1}}}

\author{Hao Wu \\[1em]Department of Biostatistics and Bioinformatics\\ Emory University\\
  Atlanta, GA 303022 \\ [1em] \texttt{hao.wu@emory.edu}}

\title{\textsf{\textbf{Differential analyses with DSS}}}


\begin{document}
\maketitle
\tableofcontents

%% abstract
\begin{abstract}
This vignette introduces the use of the Bioconductor package
DSS ({\underline D}ispersion {\underline S}hrinkage
for {\underline S}equencing data), which is designed for
differential analysis based on high-throughput sequencing data.
It performs differential expression analyses for RNA-seq,
and differential methylation analyses for bisulfite sequencing (BS-seq) data.
The core of DSS is a procedure based on Bayesian hierarchical model
to estimate and shrink gene- or CpG site-specific dispersions,
then conduct Wald tests for detecting differential expression/methylation.

%Compared with existing methods, DSS provides excellent
%statistical and computational performance.
\end{abstract}


\section{Introduction}
Recent advances in various high-throughput sequencing technologies have
revolutionized genomics research. Among them,
RNA-seq is designed to measure the the abundance of RNA products,
and Bisulfite sequencing (BS-seq) is for measuring DNA methylation.
A fundamental question in functional genomics research
is whether gene expression or DNA methylation
vary under different biological contexts.
Thus, identifying differential expression genes (DEGs)
or differential methylation loci/regions (DML/DMRs)
are key tasks in RNA-seq or BS-seq data analyses.

The differential expression (DE) or differential methylation (DM) analyses
are often based on gene- or CpG-specific statistical test.
A key limitation in RNA- or BS-seq experiments
is that the number of biological replicates is usually limited due to cost constraints.
This can lead to unstable estimation of within group variance,
and subsequently undesirable results from hypothesis testing.
Variance shrinkage methods have been widely applied in DE analyses in
microarray data to improve the estimation of gene-specific within group variances.
These methods are typically based on a Bayesian hierarchical model,
with a prior imposed on the gene-specific variances to provide
a basis for information sharing across all genes.
%In these models, shrinkage is achieved for variance estimation.
%Using shrunk variance in hypothesis tests has been shown to provide better results.

A distinct feature of RNA-seq or BS-seq data is that the measurements are
in the form of counts and have to be modeld by discrete distributions.
%. These data are often assumed to be from the Poisson (for RNA-seq)
%or Binomial (for BS-seq) distributions.
Unlike continuous distributions (such as Gaussian),
the variances depend on means in these discrete distributions.
This implies that the sample variances do not account for biological variation,
and shrinkage cannot be applied on variances directly.
In DSS, we assume that the count data are
from the Gamma-Poisson (for RNA-seq) or Beta-Binomial (for BS-seq) distribution.
We then parameterize the  distributions by a mean and a dispersion parameters.
The dispersion parameters, which represent the biological variation for
replicates within a treatment group, play a central role in the differential analyses.

DSS implements a series of DE/DM detection algorithms based on the
dispersion shrinkage method followed by Wald statistical test
to test each gene/CpG site for differential expression/methylation.
It provides functions for RNA-seq DE analysis for both two group comparision and multi-factor design,
BS-seq DM analysis for two group comparision, multi-factor design, and data without biological replicate.
Simulation and real data results show that the methods provides excellent performance
compared to existing methods, especially when the overall dispersion level is high
or the number of replicates is small.

For more details of the data model, the shrinkage method, and test procedures,
please read \cite{DE} for differential expression from RNA-seq,
\cite{DML} for differential methylation for two-group comparison from BS-seq,
\cite{DMR1} for differential methylation for data without biological replicate,
and \cite{DML-general} for differential methylation for general experimental design.


\section{Using {\tt DSS} for differential expression analysis}

\subsection{Input data preparation}
DSS requires a count table (a matrix of {\bf integers}) for gene expression values
(rows are for genes and columns are for samples).
This is different from the isoform expression based analysis such as in cufflink/cuffdiff,
where the gene expressions are represented as non-integers values.
There are a number of ways to obtain the count table from raw sequencing data (fastq file),
here we provide some example codes using several Bioconductor packages
(the codes require installation of {\tt GenomicFeatures}, {\tt Rsamtools},
and {\tt GenomicRanges} packages).

\begin{enumerate}
\item Sequence alignment. There are several RNA-seq aligner, for example, {\tt tophat}
or {\tt STAR}. Assume the alignment result is saved in a BAM file {\tt data.bam}.
\item Choose a gene annotation. {\tt GenomicFeatures} package provides a convenient way to
access different gene annotations. For example, if one wants to use RefSeq annotation for human
genome build hg19, one can use following codes:
<<eval=FALSE>>=
    library(GenomicFeatures)
    txdb = makeTranscriptDbFromUCSC(genom="hg19",tablename="refGene")
    genes = genes(txdb)
@
\item Obtain count table based on the alignment results and gene annotation.
This can be done in several steps. First read in the BAM file using the {\tt Rsamtools} package:
<<eval=FALSE>>=
    bam=scanBam("data.bam")
@
Next, create {\tt GRanges} object for the aligned sequence reads.
<<eval=FALSE>>=
    IRange.reads=GRanges(seqnames=Rle(bam$rname), ranges=IRanges(bam$pos, width=bam$qwidth))
@
Finally, use the {\tt countOverlaps} function in {\tt GenomicRanges} function to obtain the
read counts overlap each gene.
<<eval=FALSE>>=
    counts = countOverlaps(genes, IRange.reads)
@
\end{enumerate}

There are other ways to obtain the counts, for example, using
{\tt QuasR} or {\tt easyRNASeq} Bioconductor package.
Please refer to the package vignettes for more details.


\subsection{Single factor experiment}
In single factor RNA-seq experiment, DSS also requires
a vector representing experimental designs. The length of the
design vector must match the number of columns of the count table.
Optionally, normalization factors or additional annotation for genes
can be supplied.

The basic data container in the package is {\tt SeqCountSet} class,
which is directly inherited from {\tt ExpressionSet} class
defined in {\tt Biobase}. An object of the class contains all necessary
information for a DE analysis: gene expression values, experimental designs,
and additional annotations.

A typical DE analysis contains the following simple steps.
\begin{enumerate}
\item Create a {\tt SeqCountSet} object using {\tt newSeqCountSet}.
\item Estimate normalization factor using {\tt estNormFactors}.
\item Estimate and shrink gene-wise dispersion using {\tt estDispersion}
\item Two-group comparison using {\tt waldTest}.
\end{enumerate}

The usage of DSS is demonstrated in the simple simulation below.
\begin{enumerate}
\item First load in the library, and make a {\tt SeqCountSet}
object from some counts for 2000 genes and 6 samples.
<<echo=TRUE, result=TRUE>>=
library(DSS)
counts1=matrix(rnbinom(300, mu=10, size=10), ncol=3)
counts2=matrix(rnbinom(300, mu=50, size=10), ncol=3)
X1=cbind(counts1, counts2) ## these are 100 DE genes
X2=matrix(rnbinom(11400, mu=10, size=10), ncol=6)
X=rbind(X1,X2)
designs=c(0,0,0,1,1,1)
seqData=newSeqCountSet(X, designs)
seqData
@
\item Estimate normalization factor.
<<echo=TRUE, result=TRUE>>=
seqData=estNormFactors(seqData)
@
\item Estimate and shrink gene-wise dispersions
<<>>=
seqData=estDispersion(seqData)
@
\item With the normalization factors and dispersions ready, the two-group comparison can be
conducted via a Wald test:
<<>>=
result=waldTest(seqData, 0, 1)
head(result,5)
@
\end{enumerate}

A higher level wrapper function {\tt DSS.DE} is provided
for simple RNA-seq DE analysis in a two-group comparison.
User only needs to provide a count matrix and a vector of 0's and 1's representing the
design, and get DE test results in one line. A simple example is listed below:
<<>>=
counts = matrix(rpois(600, 10), ncol=6)
designs = c(0,0,0,1,1,1)
result = DSS.DE(counts, designs)
head(result)
@


\subsection{Multifactor experiment}
{\tt DSS} provides functionalities for dispersion shrinkage for multifactor experimental designs.
Downstream model fitting (through genearlized linear model)
and hypothesis testing can be performed using other packages such as {\tt edgeR},
with the dispersions estimated from DSS.

Below is an example, based a simple simulation, to illustrate the DE analysis of
a crossed design.
\begin{enumerate}
\item First simulate data for a 2x2 crossed experiments. Note the
counts are randomly generated.
<<echo=TRUE, result=TRUE>>=
library(DSS)
library(edgeR)
counts=matrix(rpois(800, 10), ncol=8)
design=data.frame(gender=c(rep("M",4), rep("F",4)), strain=rep(c("WT", "Mutant"),4))
X=model.matrix(~gender+strain, data=design)
@

\item make SeqCountSet, then estimate size factors and  dispersion
<<echo=TRUE>>=
seqData=newSeqCountSet(counts, as.data.frame(X))
seqData=estNormFactors(seqData)
seqData=estDispersion(seqData)
@


\item Using edgeR's function to do glm model fitting, but plugging in the estimated size factors and  dispersion
from DSS.
<<>>=
fit.edgeR <- glmFit(counts, X, lib.size=normalizationFactor(seqData),
             dispersion=dispersion(seqData))
@

\item Using edgeR's function to do hypothesis testing on the second parameter of the model (gender).
<<echo=TRUE, result=TRUE>>=
lrt.edgeR <- glmLRT(glmfit=fit.edgeR, coef=2)
head(lrt.edgeR$table)
@

\end{enumerate}


%%% DML detection
\section{Using {\tt DSS} for differential methylation analysis}
\subsection{Overview}
To detect differential methylation, statistical tests are conducted at each CpG site,
and then the differential methylation loci (DML) or differential methylation regions (DMR)
are called based on user specified threshold.
A rigorous statistical tests should account for
biological variations among replicates and the sequencing depth.
Most existing methods for DM analysis are based on {\it ad hoc} methods.
For example, using Fisher's exact ignores the biological variations,
using t-test on estimated methylation levels ignores the sequencing depth.
Sometimes arbitrary filtering are implemented: loci with  depth
lower than an arbitrary threshold are filtered out, which results in information loss

The DM detection procedure implemented in DSS is based on
a rigorous Wald test for beta-binomial distributions.
The test statistics depend on the biological variations (characterized by dispersion parameter)
as well as the sequencing depth. An important part of the algorithm
is the estimation of dispersion parameter, which is achieved through a
shrinkage estimator based on a Bayesian hierarchical model \cite{DML}.
An advantage of DSS is that the test can be performed even when
there is no biological replicates. That's because by smoothing,
the neighboring CpG sites can be viewed as ``pseudo-replicates", and the dispersion
can still be estimated with reasonable precision.

DSS also works for general experimental design,  based on a beta-binomial
regression model with ``arcsine'' link function.
Model fitting is performed on transformed data with generalized least
square method, which achieves much improved computational performance compared
with methods based on generalized linear model.

DSS depends on {\tt bsseq} Bioconductor package, which has neat definition of
data structures and many useful utility functions. In order to use the DM detection functionalities,
{\tt bsseq} needs to be pre-installed.

\subsection{Input data preparation}
DSS requires data from each BS-seq experiment to be summarized
into following information for each CG position:
chromosome number, genomic coordinate, total number of reads,
and number of reads showing methylation. For a sample, this information
are saved in a simple text file, with each row representing a CpG site.
Below shows an example of a small part of such a file:
\begin{verbatim}
chr     pos     N       X
chr18   3014904 26      2
chr18   3031032 33      12
chr18   3031044 33      13
chr18   3031065 48      24
\end{verbatim}

One can follow below steps to obtain such data from raw sequence file (fastq file),
using {\tt bismark} (version 0.10.0, commands for newer versions could be different)
for BS-seq alignment and count extraction. These steps require installation
of {\tt bowtie} or {\tt bowtie2}, {\tt bismark}, and the fasta file for reference genome.

\begin{enumerate}

\item Prepare Bisulfite reference genome. This can be done using the
{\tt bismark\_genome\_preparation} function (details in bismark manual). Example command is:\\
{\tt bismark\_genome\_preparation --path\_to\_bowtie /usr/local/bowtie/ --verbose /path/to/refgenomes/}

\item BS-seq alignment. Example command is:\\
{\tt bismark -q -n 1 -l 50  --path\_to\_bowtie /path/bowtie/ BS-refGenome reads.fastq}\\
This step will produce two text files {\tt reads.fastq\_bismark.sam} and
{\tt reads.fastq\_bismark\_SE\_report.txt}.

\item Extract methylation counts using {\tt bismark\_methylation\_extractor} function: \\
{\tt bismark\_methylation\_extractor -s --bedGraph reads.fastq\_bismark.sam}.
This will create multiple txt files to summarize methylation call and cytosine context,
a bedGraph file to display methylation percentage, and a coverage file containing counts information.
The count file contain following columns:{\tt chr, start, end, methylation\%, count methylated, count unmethylated}.
This file can be modified to make the input file for DSS.

\end{enumerate}

A typical DML detection contains two simple steps. First one conduct
DM test at each CpG site, then DML/DMR are called based on the test result
and user specified threshold.

\subsection{DML/DMR detection from two-group comparison}
Below are the steps to call DML or DMR for BS-seq data in two-group comparison setting.

\begin{enumerate}
\item Load in library. Read in text files and create an object of {\tt BSseq} class, which is
defined in {\tt bsseq} Bioconductor package.
This step requires {\tt bsseq} Bioconductor package. {\tt BSseq} class is defined in that package.

<<>>=
library(DSS)
require(bsseq)
path <- file.path(system.file(package="DSS"), "extdata")
dat1.1 <- read.table(file.path(path, "cond1_1.txt"), header=TRUE)
dat1.2 <- read.table(file.path(path, "cond1_2.txt"), header=TRUE)
dat2.1 <- read.table(file.path(path, "cond2_1.txt"), header=TRUE)
dat2.2 <- read.table(file.path(path, "cond2_2.txt"), header=TRUE)
BSobj <- makeBSseqData( list(dat1.1, dat1.2, dat2.1, dat2.2),
     c("C1","C2", "N1", "N2") )[1:10000,]
BSobj
@

\item Perform statistical test for DML by calling {\tt DMLtest} function.
This function basically performs following steps: (1) estimate mean methylation levels
for all CpG site; (2) estimate dispersions at each CpG sites; (3) conduct Wald test.
For the first step, there's an option for smoothing or not.  Because the methylation levels
show strong spatial correlations, smoothing can help obtain better estimates
of mean methylation when the CpG sites are dense in the data
(such as from the whole-genome BS-seq). However for data with sparse CpG,
such as from RRBS or hydroxyl-methylation, smoothing is not recommended.

To perform DML test without smoothing, do:
<<>>=
dmlTest <- DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2"))
head(dmlTest)
@

To perform statistical test for DML with smoothing, do:
<<>>=
dmlTest.sm <- DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2"), smoothing=TRUE)
@

There are two options for smoothing: a simple moving average, or the BSmooth method
implemented in {\tt bsseq} package. The BSmooth method produces much smoother curve,
which is good for visualization purpose. However, it is very computationally intensive,
and the results are not very different from moving average in terms of DMR calling.
So we recommend using moving average.
Smoothing span is an important parameter in smoothing procedure and have  non-trivial
impact on DMR calling. We use 500 bp as default, and think that it performs well in real data tests.


\item With the test results, one can call DML by using {\tt callDML} function.
The results DMLs are sorted by the significance.
<<>>=
  dmls <- callDML(dmlTest, p.threshold=0.001)
  head(dmls)
@

By default, the test is based on the null hypothesis that the difference in methylation levels is 0.
Alternatively, users can specify a threshold for difference. For example, to
detect loci with difference greater than 0.1, do:
<<>>=
  dmls2 <- callDML(dmlTest, delta=0.1, p.threshold=0.001)
  head(dmls2)
@

When delta is specified, the function will compute the posterior probability that the
difference of the means is greater than delta. So technically speaking,
the threshold for p-value here actually refers to the threshold for 1-posterior probability,
or the local FDR. Here we use the same parameter name for the sake of
the consistence of function syntax.

\item DMR detection is also Based on the DML test results, by calling {\tt callDMR} function.
Regions with many statistically significant CpG sites are identified as DMRs.
Some restrictions are provided by users, including the minimum
length, minimum number of CpG sites, percentage of CpG site being significant
in the region, etc. There are some {\it post hoc} procedures to merge nearby DMRs into longer ones.
<<>>=
dmrs <- callDMR(dmlTest, p.threshold=0.01)
head(dmrs)
@
Here the DMRs are sorted by ``areaStat", which is defined in {\tt bsseq}
as the sum of the test statistics of all CpG sites within the DMR.

Similarly, users can specify a threshold for difference. For example, to
detect regions with difference greater than 0.1, do:
<<>>=
  dmrs2 <- callDMR(dmlTest, delta=0.1, p.threshold=0.05)
  head(dmrs2)
@

Note that the distribution of test statistics (and p-values) depends on
the differences in methylation levels and biological variations,
as well as technical factors such as coverage depth. It is very difficulty
to select a natural and rigorous threshold for defining DMRs. We recommend
users try different thresholds in order to obtain satisfactory results.

\item The DMRs can be visualized using {\tt showOneDMR} function,
This function provides more information than the {\tt plotRegion} function in {\tt bsseq}.
It plots the methylation percentages as well as the coverage depths
at each CpG sites, instead of just the smoothed curve.
So the coverage depth information will be available in the figure.

To use the function, do
<<eval=FALSE>>=
  showOneDMR(dmrs[1,], BSobj)
@

The result figure looks like the following.
{\bf Note that the figure below is not generated from the above example.
The example data are from RRBS experiment so the DMRs are much shorter.}

\begin{figure}[h!]
\centerline{\includegraphics[width=5.5in]{aDMR.pdf}}
\end{figure}

\end{enumerate}

\newpage
\subsection{DML/DMR detection from general experimental design}
In DSS, BS-seq data from a general experimental design (such as crossed experiment,
or experiment with covariates) is modeled through a generalized linear model framework.
We use ``arcsine'' link function instead of the typical logit link for it better deals with data
at boundaries (methylation levels close to 0 or 1).  Linear model fitting is done
through ordinary least square on transformed methylation levels. Standard errors
for the estimates are derived with consideration of count data distribution and transformation.
A Wald test is applied to perform hypothesis testing.

\begin{enumerate}

\item Load in data distributed with {\tt DSS}. This is a small portion of a set of
RRBS experiments. There are 5000 CpG sites and 16 samples.
The experiment is a $2\time 2$ design (2 cases and 2 cell types).
There are 4 replicates in each case:cell combination.

<<>>=
data(RRBS)
RRBS
design
@

\item Fit a linear model using {\tt DMLfit.multiFactor} function, include
case, cell, and case:cell interaction.
<<>>=
DMLfit = DMLfit.multiFactor(RRBS, design=design, formula=~case+cell+case:cell)
@

\item Use {\tt DMLtest.multiFactor} function to test the cell effect.
It is important to note that the {\tt coef} parameter is the index
of the coefficient to be tested for being 0. Because the model
(as specified by {\tt formula} in {\tt DMLfit.multiFactor}) include intercept,
the cell effect is the 3rd column in the design matrix, so we use
{\tt coef=3} here.

<<>>=
DMLtest.cell = DMLtest.multiFactor(DMLfit, coef=3)
@

Result from this step is a data frame with chromosome number, CpG site position,
test statistics, p-values (from normal distribution), and FDR.
Rows are sorted by chromosome/position of the CpG sites.
To obtain top ranked CpG sites, one can sort the data frame using following codes:

<<>>=
ix=sort(DMLtest.cell[,"pvals"], index.return=TRUE)$ix
head(DMLtest.cell[ix,])
@

Below is a figure showing the distributions of test statistics and p-values from this example dataset

<<fig=TRUE, height=4, width=8>>=
par(mfrow=c(1,2))
hist(DMLtest.cell$stat, 50, main="test statistics", xlab="")
hist(DMLtest.cell$pvals, 50, main="P values", xlab="")
@

\end{enumerate}

These procedures are computationally very efficient. For a typical RRBS dataset with 4 million CpG sites,
it takes around 20 minutes. In comparison, other methods such as RADMeth or BiSeq
takes at least 10 times longer.


\newpage
\section{Session Info}
<<echo=TRUE, result=TRUE>>=
sessionInfo()
@

\begin{thebibliography}{99}

\bibitem{DML}
\textsc{Hao Feng, Karen Conneely and Hao Wu}. (2014).
\newblock A bayesian hierarchical model to detect differentially methylated
  loci from single nucleotide resolution sequencing data.
\newblock {\em Nucleic Acids Research.\/}~\textbf{42}(8), e69--e69.

\bibitem{DMR1}
\textsc{Hao Wu, Tianlei Xu, Hao Feng, Li Chen, Ben Li, Bing Yao, Zhaohui Qin, Peng Jin and Karen N. Conneely}. (2015).
\newblock {Detection of differentially methylated regions from whole-genome bisulfite sequencing data without replicates.}
\newblock {\em Nucleic Acids Research.\/}~doi: 10.1093/nar/gkv715.

\bibitem{DML-general}
\textsc{Yongseok Park, Hao Wu}. (2016).
\newblock {Differential methylation analysis for BS-seq data under general experimental design.}
\newblock {\em Bioinformatics.\/}~doi:10.1093/bioinformatics/btw026.

\bibitem{DE}
\textsc{Hao Wu, Chi Wang and Zhijing Wu}. (2013).
\newblock {A new shrinkage estimator for dispersion improves differential
  expression detection in RNA-seq data.}
\newblock {\em Biostatistics.\/}~\textbf{14}(2), 232--243.

\end{thebibliography}

\end{document}

