This vignette introduces the use of the Bioconductor package DSS (Dispersion Shrinkage for Sequencing data), which is designed for differential analysis on high-throughput sequencing data. It provide functionalities for RNA-seq differential expression, and bisulfite sequencing (BS-seq) differential methylation. 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.
DSS 2.48.0
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.
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. Unlike continuous distributions such as Gaussian, the variances depend on means in these 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 negative-binomial (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.
For more details of the data model, the shrinkage method, and test procedures, please read (Wu, Wang, and Wu 2012) for differential expression from RNA-seq, (Feng, Conneely, and Wu 2014) for differential methylation for two-group comparison from BS-seq, (Wu et al. 2015) for differential methylation for data without biological replicate, and (Park and Wu 2016) for differential methylation for general experimental design.
DSS requires a count table (a matrix of 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). Aligners such as STAR automatically output a count table. Several Bioconductor packages serve this purpose, for example, Rsubread, QuasR, and easyRNASeq. Other tools like Salmon, Kallisto, or RSEM can also be used. Please refer to the package manuals for more details.
In single factor RNA-seq experiment, DSS 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 SeqCountSet class,
which is directly inherited from ExpressionSet class
defined in 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.
- Create a SeqCountSet object using newSeqCountSet.
- Estimate normalization factor using estNormFactors.
- Estimate and shrink gene-wise dispersion using estDispersion.
- Two-group comparison using waldTest.
The usage of DSS is demonstrated in the simple simulation below.
SeqCountSet
object from simulated counts for 2000 genes and 6 samples.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)
X2 = matrix(rnbinom(11400, mu=10, size=10), ncol=6)
X = rbind(X1,X2) ## these are 100 DE genes
designs = c(0,0,0,1,1,1)
seqData = newSeqCountSet(X, designs)
seqData## SeqCountSet (storageMode: lockedEnvironment)
## assayData: 2000 features, 6 samples
## element names: exprs
## protocolData: none
## phenoData
## sampleNames: 1 2 ... 6 (6 total)
## varLabels: designs
## varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation:seqData = estNormFactors(seqData)seqData = estDispersion(seqData)result=waldTest(seqData, 0, 1)
head(result,5)## geneIndex muA muB lfc difExpr stats pval
## 22 22 6.000000 59.66176 -2.225235 -53.66176 -5.493187 3.947447e-08
## 92 92 4.333333 47.51348 -2.295945 -43.18015 -5.471977 4.450425e-08
## 93 93 8.666667 68.35784 -2.016470 -59.69118 -5.441543 5.282107e-08
## 44 44 7.666667 62.31740 -2.040171 -54.65074 -5.420574 5.940810e-08
## 76 76 6.000000 54.24877 -2.130953 -48.24877 -5.415795 6.101681e-08
## fdr
## 22 2.116807e-05
## 92 2.116807e-05
## 93 2.116807e-05
## 44 2.116807e-05
## 76 2.116807e-05A higher level wrapper function 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)## geneIndex muA muB lfc difExpr stats pval
## 55 55 16.857221 7.666667 0.7539478 9.190554 2.771702 0.005576405
## 9 9 6.040833 11.666667 -0.6206354 -5.625834 -2.079509 0.037570573
## 45 45 14.149721 8.333333 0.5058888 5.816388 1.866238 0.062008016
## 100 100 12.829999 7.666667 0.4899562 5.163332 1.752798 0.079636733
## 66 66 13.591942 8.666667 0.4300295 4.925276 1.592430 0.111288096
## 34 34 12.163332 7.666667 0.4386497 4.496665 1.558904 0.119019137
## fdr
## 55 0.5576405
## 9 0.9703365
## 45 0.9703365
## 100 0.9703365
## 66 0.9703365
## 34 0.9703365DSS 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 edgeR,
with the dispersions estimated from DSS.
Below is an example, based a simple simulation, to illustrate the DE analysis of a crossed design.
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)seqData = newSeqCountSet(counts, as.data.frame(X))
seqData = estNormFactors(seqData)
seqData = estDispersion(seqData)fit.edgeR = glmFit(counts, X, dispersion=dispersion(seqData))lrt.edgeR = glmLRT(glmfit=fit.edgeR, coef=2)
head(lrt.edgeR$table)## logFC logCPM LR PValue
## 1 -0.38745557 13.38124 0.9676174224 0.3252751
## 2 0.29009591 13.66414 0.6740642538 0.4116382
## 3 -0.35527698 13.50859 0.9198178075 0.3375228
## 4 0.54582085 13.66542 2.3389074260 0.1261780
## 5 0.01930945 13.44545 0.0025426398 0.9597840
## 6 0.01209868 13.38154 0.0009337767 0.9756222To 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 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 (Feng, Conneely, and Wu 2014). 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 bsseq
Bioconductor package, which has neat definition of
data structures and many useful utility functions. In order to use the DM detection functionalities,
bsseq needs to be pre-installed.
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:
chr pos N X
chr18 3014904 26 2
chr18 3031032 33 12
chr18 3031044 33 13
chr18 3031065 48 24One can follow below steps to obtain such data from raw sequence file (fastq file),
using bismark
(version 0.10.0, commands for newer versions could be different)
for BS-seq alignment and count extraction. These steps require installation
of bowtie
or bowtie2,
bismark, and the fasta file for reference genome.
bismark_genome_preparation function (details in bismark manual). Example command is:bismark_genome_preparation --path_to_bowtie /usr/local/bowtie/ \
--verbose /path/to/refgenomes/bismark -q -n 1 -l 50 --path_to_bowtie \
/path/bowtie/ BS-refGenome reads.fastqThis step will produce two text files reads.fastq_bismark.sam and
reads.fastq_bismark_SE_report.txt.
bismark_methylation_extractor function:bismark_methylation_extractor -s --bedGraph reads.fastq_bismark.samThis 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: chr, start, end, methylation%, count methylated, count unmethylated.
This file can be modified to make the input file for DSS.
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.
Step 1. Load in library. Read in text files and create an object of BSseq class, which is
defined in bsseq Bioconductor 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:1000,]
BSobj## An object of type 'BSseq' with
## 1000 methylation loci
## 4 samples
## has not been smoothed
## All assays are in-memoryStep 2. Perform statistical test for DML by calling 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 often help obtain better estimates
of mean methylation.
To perform DML test without smoothing, do:
dmlTest = DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2"))To perform statistical test for DML with smoothing, do:
dmlTest.sm = DMLtest(BSobj, group1=c("C1", "C2"), group2=c("N1", "N2"),
smoothing=TRUE)User has the option to smooth the methylation levels or not. A simple moving average algorithm is implemented for smoothing. For WGBS data, smoothing is always recommended so that information from nearby CpG sites can be combined to improve the estimation of methylation levels. For data with sparse CpG coverage such as RRBS or hydroxymethylation, smoothing might not alter the results much, but is still recommended. In RRBS, CpG sites are likely to be clustered locally within small genomic regions, so smoothing can potentially help the methylation estimation.
If smoothing is requested, smoothing span is an important parameter which has non-trivial impact on DMR calling. We use 500 bp as default, because it performs well in real data tests according to our experience.
Step 3. With the test results, one can call DML by using callDML function.
The results DMLs are sorted by the significance.
dmls = callDML(dmlTest, p.threshold=0.001)
head(dmls)## chr pos mu1 mu2 diff diff.se stat
## 450 chr18 3976129 0.01027497 0.9390339 -0.9287590 0.06544340 -14.19179
## 451 chr18 3976138 0.01027497 0.9390339 -0.9287590 0.06544340 -14.19179
## 638 chr18 4431501 0.01331553 0.9430566 -0.9297411 0.09273779 -10.02548
## 639 chr18 4431511 0.01327049 0.9430566 -0.9297862 0.09270080 -10.02997
## 710 chr18 4564237 0.91454619 0.0119300 0.9026162 0.05260037 17.15988
## 782 chr18 4657576 0.98257334 0.0678355 0.9147378 0.06815000 13.42242
## phi1 phi2 pval fdr postprob.overThreshold
## 450 0.052591567 0.02428826 1.029974e-45 2.499403e-43 1
## 451 0.052591567 0.02428826 1.029974e-45 2.499403e-43 1
## 638 0.053172411 0.07746835 1.177826e-23 1.429096e-21 1
## 639 0.053121697 0.07746835 1.125518e-23 1.429096e-21 1
## 710 0.009528898 0.04942849 5.302004e-66 3.859859e-63 1
## 782 0.010424723 0.06755651 4.468885e-41 8.133371e-39 1By 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)## chr pos mu1 mu2 diff diff.se stat
## 450 chr18 3976129 0.01027497 0.9390339 -0.9287590 0.06544340 -14.19179
## 451 chr18 3976138 0.01027497 0.9390339 -0.9287590 0.06544340 -14.19179
## 638 chr18 4431501 0.01331553 0.9430566 -0.9297411 0.09273779 -10.02548
## 639 chr18 4431511 0.01327049 0.9430566 -0.9297862 0.09270080 -10.02997
## 710 chr18 4564237 0.91454619 0.0119300 0.9026162 0.05260037 17.15988
## 782 chr18 4657576 0.98257334 0.0678355 0.9147378 0.06815000 13.42242
## phi1 phi2 pval fdr postprob.overThreshold
## 450 0.052591567 0.02428826 1.029974e-45 2.499403e-43 1
## 451 0.052591567 0.02428826 1.029974e-45 2.499403e-43 1
## 638 0.053172411 0.07746835 1.177826e-23 1.429096e-21 1
## 639 0.053121697 0.07746835 1.125518e-23 1.429096e-21 1
## 710 0.009528898 0.04942849 5.302004e-66 3.859859e-63 1
## 782 0.010424723 0.06755651 4.468885e-41 8.133371e-39 1When 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.
Step 4. DMR detection is also Based on the DML test results, by calling 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 post hoc procedures to merge nearby DMRs into longer ones.
dmrs = callDMR(dmlTest, p.threshold=0.01)
head(dmrs)## chr start end length nCG meanMethy1 meanMethy2 diff.Methy areaStat
## 27 chr18 4657576 4657639 64 4 0.506453 0.318348 0.188105 14.34236Here the DMRs are sorted by areaStat, which is defined in 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)## chr start end length nCG meanMethy1 meanMethy2 diff.Methy areaStat
## 31 chr18 4657576 4657639 64 4 0.5064530 0.3183480 0.188105 14.34236
## 19 chr18 4222533 4222608 76 4 0.7880276 0.3614195 0.426608 12.91667Note 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.
The DMRs can be visualized using showOneDMR function,
This function provides more information than the plotRegion function in 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
showOneDMR(dmrs[1,], BSobj)The result figure looks like the following. 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.