Differential expression analysis based on the Negative Binomial (a.k.a. Gamma-Poisson) distribution

Description

This function performs a default analysis through the steps:

• estimation of size factors: estimateSizeFactors

• estimation of dispersion: estimateDispersions

• Negative Binomial GLM fitting and Wald statistics: nbinomWaldTest
  For complete details on each step, see the manual pages of the respective functions. After the DESeq function returns a DESeqDataSet object, results tables (log2 fold changes and p-values) can be generated using the results function. Shrunken LFC can then be generated using the lfcShrink function. All support questions should be posted to the Bioconductor support site: http://support.bioconductor.org .

Usage

DESeq(object, test = c("Wald", "LRT"), fitType = c("parametric",
  "local", "mean"), sfType = c("ratio", "poscounts", "iterate"),
  betaPrior, full = design(object), reduced, quiet = FALSE,
  minReplicatesForReplace = 7, modelMatrixType, useT = FALSE,
  minmu = 0.5, parallel = FALSE, BPPARAM = bpparam())

Arguments

Details

The differential expression analysis uses a generalized linear model of the form:

$$K{ij} im extrm{NB}( mu{ij}, lpha_i)$$

$$mu{ij} = s_j q{ij}$$

$$log2(q{ij}) = x_{j.} eta_i$$

where counts $K{ij}$ for gene i, sample j are modeled using a Negative Binomial distribution with fitted mean $mu{ij}$ and a gene-specific dispersion parameter $lphai$ . The fitted mean is composed of a sample-specific size factor $s_j$ and a parameter $q{ij}$ proportional to the expected true concentration of fragments for sample j. The coefficients $eta_i$ give the log2 fold changes for gene i for each column of the model matrix $X$ . The sample-specific size factors can be replaced by gene-specific normalization factors for each sample using normalizationFactors .

For details on the fitting of the log2 fold changes and calculation of p-values, see nbinomWaldTest if using test="Wald" , or nbinomLRT if using test="LRT" .

Experiments without replicates do not allow for estimation of the dispersion of counts around the expected value for each group, which is critical for differential expression analysis. Analysis without replicates was deprecated in v1.20 and is no longer supported since v1.22.

The argument minReplicatesForReplace is used to decide which samples are eligible for automatic replacement in the case of extreme Cook's distance. By default, DESeq will replace outliers if the Cook's distance is large for a sample which has 7 or more replicates (including itself). This replacement is performed by the replaceOutliers function. This default behavior helps to prevent filtering genes based on Cook's distance when there are many degrees of freedom. See results for more information about filtering using Cook's distance, and the 'Dealing with outliers' section of the vignette. Unlike the behavior of replaceOutliers , here original counts are kept in the matrix returned by counts , original Cook's distances are kept in assays(dds)[["cooks"]] , and the replacement counts used for fitting are kept in assays(dds)[["replaceCounts"]] .

Note that if a log2 fold change prior is used (betaPrior=TRUE) then expanded model matrices will be used in fitting. These are described in nbinomWaldTest and in the vignette. The contrast argument of results should be used for generating results tables.

Value

a DESeqDataSet object with results stored as metadata columns. These results should accessed by calling the results function. By default this will return the log2 fold changes and p-values for the last variable in the design formula. See results for how to access results for other variables.

Seealso

link{results} , lfcShrink , nbinomWaldTest , nbinomLRT

Author

Michael Love

References

Love, M.I., Huber, W., Anders, S. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15:550. https://doi.org/10.1186/s13059-014-0550-8

Examples

# see vignette for suggestions on generating
# count tables from RNA-Seq data
cnts <- matrix(rnbinom(n=1000, mu=100, size=1/0.5), ncol=10)
cond <- factor(rep(1:2, each=5))

# object construction
dds <- DESeqDataSetFromMatrix(cnts, DataFrame(cond), ~ cond)

# standard analysis
dds <- DESeq(dds)
res <- results(dds)

# moderated log2 fold changes
resultsNames(dds)
resLFC <- lfcShrink(dds, coef=2, type="apeglm")

# an alternate analysis: likelihood ratio test
ddsLRT <- DESeq(dds, test="LRT", reduced= ~ 1)
resLRT <- results(ddsLRT)

DESeq2 package for differential analysis of count data

Description

The main functions for differential analysis are DESeq and results . See the examples at DESeq for basic analysis steps. Two transformations offered for count data are the variance stabilizing transformation, vst , and the "regularized logarithm", rlog . For more detailed information on usage, see the package vignette, by typing vignette("DESeq2") , or the workflow linked to on the first page of the vignette. All support questions should be posted to the Bioconductor support site: http://support.bioconductor.org .

Author

Michael Love, Wolfgang Huber, Simon Anders

References

DESeq2 reference:

Love, M.I., Huber, W., Anders, S. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15:550. http://dx.doi.org/10.1186/s13059-014-0550-8

DESeq reference:

Simon Anders, Wolfgang Huber (2010) Differential expression analysis for sequence count data. Genome Biology, 11:106. http://dx.doi.org/10.1186/gb-2010-11-10-r106

DESeqDataSet object and constructors

Description

DESeqDataSet is a subclass of RangedSummarizedExperiment , used to store the input values, intermediate calculations and results of an analysis of differential expression. The DESeqDataSet class enforces non-negative integer values in the "counts" matrix stored as the first element in the assay list. In addition, a formula which specifies the design of the experiment must be provided. The constructor functions create a DESeqDataSet object from various types of input: a RangedSummarizedExperiment, a matrix, count files generated by the python package HTSeq, or a list from the tximport function in the tximport package. See the vignette for examples of construction from different types.

Usage

DESeqDataSet(se, design, ignoreRank = FALSE)
DESeqDataSetFromMatrix(countData, colData, design, tidy = FALSE,
  ignoreRank = FALSE, ...)
DESeqDataSetFromHTSeqCount(sampleTable, directory = ".", design,
  ignoreRank = FALSE, ...)
DESeqDataSetFromTximport(txi, colData, design, ...)

Arguments

Details

Note on the error message "assay colnames() must be NULL or equal colData rownames()": this means that the colnames of countData are different than the rownames of colData. Fix this with: colnames(countData) <- NULL

Value

A DESeqDataSet object.

References

See http://www-huber.embl.de/users/anders/HTSeq for htseq-count

Examples

countData <- matrix(1:100,ncol=4)
condition <- factor(c("A","A","B","B"))
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), ~ condition)

DESeqResults object and constructor

Description

This constructor function would not typically be used by "end users". This simple class extends the DataFrame class of the IRanges package to allow other packages to write methods for results objects from the DESeq2 package. It is used by results to wrap up the results table.

Usage

DESeqResults(DataFrame, priorInfo = list())

Arguments

Value

a DESeqResults object

DESeqTransform object and constructor

Description

This constructor function would not typically be used by "end users". This simple class extends the RangedSummarizedExperiment class of the SummarizedExperiment package. It is used by rlog and varianceStabilizingTransformation to wrap up the results into a class for downstream methods, such as plotPCA .

Usage

DESeqTransform(SummarizedExperiment)

Arguments

Value

a DESeqTransform object

Extract a matrix of model coefficients/standard errors

Description

Note: results tables with log2 fold change, p-values, adjusted p-values, etc. for each gene are best generated using the results function. The coef function is designed for advanced users who wish to inspect all model coefficients at once.

Usage

list(list("coef"), list("DESeqDataSet"))(object, SE = FALSE, ...)

Arguments

Details

Estimated model coefficients or estimated standard errors are provided in a matrix form, number of genes by number of parameters, on the log2 scale. The columns correspond to columns of the model matrix for final GLM fitting, i.e., attr(dds, "modelMatrix") .

Author

Michael Love

Examples

dds <- makeExampleDESeqDataSet(m=4)
dds <- DESeq(dds)
coef(dds)[1,]
coef(dds, SE=TRUE)[1,]

Collapse technical replicates in a RangedSummarizedExperiment or DESeqDataSet

Description

Collapses the columns in object by summing within levels of a grouping factor groupby . The purpose of this function is to sum up read counts from technical replicates to create an object with a single column of read counts for each sample. Note: by "technical replicates", we mean multiple sequencing runs of the same library, in constrast to "biological replicates" in which multiple libraries are prepared from separate biological units. Optionally renames the columns of returned object with the levels of the grouping factor. Note: this function is written very simply and can be easily altered to produce other behavior by examining the source code.

Usage

collapseReplicates(object, groupby, run, renameCols = TRUE)

Arguments

Value

the object with as many columns as levels in groupby . This object has assay/count data which is summed from the various columns which are grouped together, and the colData is subset using the first column for each group in groupby .

Examples

dds <- makeExampleDESeqDataSet(m=12)

# make data with two technical replicates for three samples
dds$sample <- factor(sample(paste0("sample",rep(1:9, c(2,1,1,2,1,1,2,1,1)))))
dds$run <- paste0("run",1:12)

ddsColl <- collapseReplicates(dds, dds$sample, dds$run)

# examine the colData and column names of the collapsed data
colData(ddsColl)
colnames(ddsColl)

# check that the sum of the counts for "sample1" is the same
# as the counts in the "sample1" column in ddsColl
matchFirstLevel <- dds$sample == levels(dds$sample)[1]
stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))

Accessors for the 'counts' slot of a DESeqDataSet object.

Description

The counts slot holds the count data as a matrix of non-negative integer count values, one row for each observational unit (gene or the like), and one column for each sample.

Usage

list(list("counts"), list("DESeqDataSet"))(object, normalized = FALSE,
  replaced = FALSE)
list(list("counts"), list("DESeqDataSet,matrix"))(object) <- value

Arguments

Seealso

sizeFactors , normalizationFactors

Author

Simon Anders

Examples

dds <- makeExampleDESeqDataSet(m=4)
head(counts(dds))

dds <- estimateSizeFactors(dds) # run this or DESeq() first
head(counts(dds, normalized=TRUE))

Accessors for the 'design' slot of a DESeqDataSet object.

Description

The design holds the R formula which expresses how the counts depend on the variables in colData . See DESeqDataSet for details.

Usage

list(list("design"), list("DESeqDataSet"))(object)
list(list("design"), list("DESeqDataSet,formula"))(object) <- value
list(list("design"), list("DESeqDataSet,matrix"))(object) <- value

Arguments

Examples

dds <- makeExampleDESeqDataSet(m=4)
design(dds) <- formula(~ 1)

Accessors for the 'dispersionFunction' slot of a DESeqDataSet object.

Description

The dispersion function is calculated by estimateDispersions and used by varianceStabilizingTransformation . Parametric dispersion fits store the coefficients of the fit as attributes in this slot.

Usage

dispersionFunction(object, ...)
dispersionFunction(object, ...) <- value
list(list("dispersionFunction"), list("DESeqDataSet"))(object)
list(list("dispersionFunction"), list("DESeqDataSet,`function`"))(object,
  estimateVar = TRUE) <- value

Arguments

Details

Setting this will also overwrite mcols(object)$dispFit and the estimate the variance of dispersion residuals, see estimateVar below.

Seealso

estimateDispersions

Examples

dds <- makeExampleDESeqDataSet(m=4)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dispersionFunction(dds)

Accessor functions for the dispersion estimates in a DESeqDataSet object.

Description

The dispersions for each row of the DESeqDataSet. Generally, these are set by estimateDispersions .

Usage

dispersions(object, ...)
dispersions(object, ...) <- value
list(list("dispersions"), list("DESeqDataSet"))(object)
list(list("dispersions"), list("DESeqDataSet,numeric"))(object) <- value

Arguments

Seealso

estimateDispersions

Author

Simon Anders

Steps for estimating the beta prior variance

Description

These lower-level functions are called within DESeq or nbinomWaldTest . End users should use those higher-level function instead. NOTE: estimateBetaPriorVar returns a numeric vector, not a DESEqDataSet! For advanced users: to use these functions, first run estimateMLEForBetaPriorVar and then run estimateBetaPriorVar .

Usage

estimateBetaPriorVar(object, betaPriorMethod = c("weighted", "quantile"),
  upperQuantile = 0.05, modelMatrix = NULL)
estimateMLEForBetaPriorVar(object, maxit = 100, useOptim = TRUE,
  useQR = TRUE, modelMatrixType = NULL)

Arguments

Value

for estimateMLEForBetaPriorVar , a DESeqDataSet, with the necessary information stored in order to calculate the prior variance. for estimateBetaPriorVar , the vector of variances for the prior on the betas in the DESeq GLM

Estimate the dispersions for a DESeqDataSet

Description

This function obtains dispersion estimates for Negative Binomial distributed data.

Usage

list(list("estimateDispersions"), list("DESeqDataSet"))(object,
  fitType = c("parametric", "local", "mean"), maxit = 100,
  quiet = FALSE, modelMatrix = NULL, minmu = 0.5)

Arguments

• parametric - fit a dispersion-mean relation of the form:

via a robust gamma-family GLM. The coefficients asymptDisp and extraPois are given in the attribute coefficients of the dispersionFunction of the object.

• local - use the locfit package to fit a local regression of log dispersions over log base mean (normal scale means and dispersions are input and output for dispersionFunction ). The points are weighted by normalized mean count in the local regression.

• mean - use the mean of gene-wise dispersion estimates. |maxit | control parameter: maximum number of iterations to allow for convergence| |quiet | whether to print messages at each step| |modelMatrix | an optional matrix which will be used for fitting the expected counts. by default, the model matrix is constructed from design(object)| |minmu | lower bound on the estimated count for fitting gene-wise dispersion|

Details

Typically the function is called with the idiom:

dds <- estimateDispersions(dds)

The fitting proceeds as follows: for each gene, an estimate of the dispersion is found which maximizes the Cox Reid-adjusted profile likelihood (the methods of Cox Reid-adjusted profile likelihood maximization for estimation of dispersion in RNA-Seq data were developed by McCarthy, et al. (2012), first implemented in the edgeR package in 2010); a trend line capturing the dispersion-mean relationship is fit to the maximum likelihood estimates; a normal prior is determined for the log dispersion estimates centered on the predicted value from the trended fit with variance equal to the difference between the observed variance of the log dispersion estimates and the expected sampling variance; finally maximum a posteriori dispersion estimates are returned. This final dispersion parameter is used in subsequent tests. The final dispersion estimates can be accessed from an object using dispersions . The fitted dispersion-mean relationship is also used in varianceStabilizingTransformation . All of the intermediate values (gene-wise dispersion estimates, fitted dispersion estimates from the trended fit, etc.) are stored in mcols(dds) , with information about these columns in mcols(mcols(dds)) .

The log normal prior on the dispersion parameter has been proposed by Wu, et al. (2012) and is also implemented in the DSS package.

In DESeq2, the dispersion estimation procedure described above replaces the different methods of dispersion from the previous version of the DESeq package.

estimateDispersions checks for the case of an analysis with as many samples as the number of coefficients to fit, and will temporarily substitute a design formula ~ 1 for the purposes of dispersion estimation. Note that analysis of designs without replicates will be removed in the Oct 2018 release: DESeq2 v1.22.0, after which DESeq2 will give an error.

The lower-level functions called by estimateDispersions are: estimateDispersionsGeneEst , estimateDispersionsFit , and estimateDispersionsMAP .

Value

The DESeqDataSet passed as parameters, with the dispersion information filled in as metadata columns, accessible via mcols , or the final dispersions accessible via dispersions .

References

• Simon Anders, Wolfgang Huber: Differential expression analysis for sequence count data. Genome Biology 11 (2010) R106, http://dx.doi.org/10.1186/gb-2010-11-10-r106

• McCarthy, DJ, Chen, Y, Smyth, GK: Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Research 40 (2012), 4288-4297, http://dx.doi.org/10.1093/nar/gks042

• Wu, H., Wang, C. & Wu, Z. A new shrinkage estimator for dispersion improves differential expression detection in RNA-seq data. Biostatistics (2012). http://dx.doi.org/10.1093/biostatistics/kxs033

Examples

dds <- makeExampleDESeqDataSet()
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
head(dispersions(dds))

Low-level functions to fit dispersion estimates

Description

Normal users should instead use estimateDispersions . These low-level functions are called by estimateDispersions , but are exported and documented for non-standard usage. For instance, it is possible to replace fitted values with a custom fit and continue with the maximum a posteriori dispersion estimation, as demonstrated in the examples below.

Usage

estimateDispersionsGeneEst(object, minDisp = 1e-08, kappa_0 = 1,
  dispTol = 1e-06, maxit = 100, quiet = FALSE, modelMatrix = NULL,
  niter = 1, linearMu = NULL, minmu = 0.5)
estimateDispersionsFit(object, fitType = c("parametric", "local",
  "mean"), minDisp = 1e-08, quiet = FALSE)
estimateDispersionsMAP(object, outlierSD = 2, dispPriorVar,
  minDisp = 1e-08, kappa_0 = 1, dispTol = 1e-06, maxit = 100,
  modelMatrix = NULL, quiet = FALSE)
estimateDispersionsPriorVar(object, minDisp = 1e-08,
  modelMatrix = NULL)

Arguments

Value

a DESeqDataSet with gene-wise, fitted, or final MAP dispersion estimates in the metadata columns of the object.

estimateDispersionsPriorVar is called inside of estimateDispersionsMAP and stores the dispersion prior variance as an attribute of dispersionFunction(dds) , which can be manually provided to estimateDispersionsMAP for parallel execution.

Seealso

estimateDispersions

Examples

dds <- makeExampleDESeqDataSet()
dds <- estimateSizeFactors(dds)
dds <- estimateDispersionsGeneEst(dds)
dds <- estimateDispersionsFit(dds)
dds <- estimateDispersionsMAP(dds)
plotDispEsts(dds)

# after having run estimateDispersionsFit()
# the dispersion prior variance over all genes
# can be obtained like so:

dispPriorVar <- estimateDispersionsPriorVar(dds)

Estimate the size factors for a DESeqDataSet

Description

This function estimates the size factors using the "median ratio method" described by Equation 5 in Anders and Huber (2010). The estimated size factors can be accessed using the accessor function sizeFactors . Alternative library size estimators can also be supplied using the assignment function sizeFactors<- .

Usage

list(list("estimateSizeFactors"), list("DESeqDataSet"))(object, type = c("ratio",
  "poscounts", "iterate"), locfunc = stats::median, geoMeans,
  controlGenes, normMatrix, quiet = FALSE)

Arguments

Details

Typically, the function is called with the idiom:

dds <- estimateSizeFactors(dds)

See DESeq for a description of the use of size factors in the GLM. One should call this function after DESeqDataSet unless size factors are manually specified with sizeFactors . Alternatively, gene-specific normalization factors for each sample can be provided using normalizationFactors which will always preempt sizeFactors in calculations.

Internally, the function calls estimateSizeFactorsForMatrix , which provides more details on the calculation.

Value

The DESeqDataSet passed as parameters, with the size factors filled in.

Seealso

estimateSizeFactorsForMatrix

Author

Simon Anders

References

Reference for the median ratio method:

Simon Anders, Wolfgang Huber: Differential expression analysis for sequence count data. Genome Biology 2010, 11:106. http://dx.doi.org/10.1186/gb-2010-11-10-r106

Examples

dds <- makeExampleDESeqDataSet(n=1000, m=4)
dds <- estimateSizeFactors(dds)
sizeFactors(dds)

dds <- estimateSizeFactors(dds, controlGenes=1:200)

m <- matrix(runif(1000 * 4, .5, 1.5), ncol=4)
dds <- estimateSizeFactors(dds, normMatrix=m)
normalizationFactors(dds)[1:3,]

geoMeans <- exp(rowMeans(log(counts(dds))))
dds <- estimateSizeFactors(dds,geoMeans=geoMeans)
sizeFactors(dds)

Low-level function to estimate size factors with robust regression.

Description

Given a matrix or data frame of count data, this function estimates the size factors as follows: Each column is divided by the geometric means of the rows. The median (or, if requested, another location estimator) of these ratios (skipping the genes with a geometric mean of zero) is used as the size factor for this column. Typically, one will not call this function directly, but use estimateSizeFactors .

Usage

estimateSizeFactorsForMatrix(counts, locfunc = stats::median, geoMeans,
  controlGenes)

Arguments

Value

a vector with the estimates size factors, one element per column

Seealso

estimateSizeFactors

Author

Simon Anders

Examples

dds <- makeExampleDESeqDataSet()
estimateSizeFactorsForMatrix(counts(dds))
geoMeans <- exp(rowMeans(log(counts(dds))))
estimateSizeFactorsForMatrix(counts(dds),geoMeans=geoMeans)

FPKM: fragments per kilobase per million mapped fragments

Description

The following function returns fragment counts normalized per kilobase of feature length per million mapped fragments (by default using a robust estimate of the library size, as in estimateSizeFactors ).

Usage

fpkm(object, robust = TRUE)

Arguments

Details

The length of the features (e.g. genes) is calculated one of two ways: (1) If there is a matrix named "avgTxLength" in assays(dds) , this will take precedence in the length normalization. This occurs when using the tximport-DESeq2 pipeline. (2) Otherwise, feature length is calculated from the rowRanges of the dds object, if a column basepairs is not present in mcols(dds) . The calculated length is the number of basepairs in the union of all GRanges assigned to a given row of object , e.g., the union of all basepairs of exons of a given gene. Note that the second approach over-estimates the gene length (average transcript length, weighted by abundance is a more appropriate normalization for gene counts), and so the FPKM will be an underestimate of the true value.

Note that, when the read/fragment counting has inter-feature dependencies, a strict normalization would not incorporate the basepairs of a feature which overlap another feature. This inter-feature dependence is not taken into consideration in the internal union basepair calculation.

Value

a matrix which is normalized per kilobase of the union of basepairs in the GRangesList or GRanges of the mcols(object), and per million of mapped fragments, either using the robust median ratio method (robust=TRUE, default) or using raw counts (robust=FALSE). Defining a column mcols(object)$basepairs takes precedence over internal calculation of the kilobases for each row.

Seealso

fpm

Examples

# create a matrix with 1 million counts for the
# 2nd and 3rd column, the 1st and 4th have
# half and double the counts, respectively.
m <- matrix(1e6 * rep(c(.125, .25, .25, .5), each=4),
ncol=4, dimnames=list(1:4,1:4))
mode(m) <- "integer"
se <- SummarizedExperiment(list(counts=m), colData=DataFrame(sample=1:4))
dds <- DESeqDataSet(se, ~ 1)

# create 4 GRanges with lengths: 1, 1, 2, 2.5 Kb
gr1 <- GRanges("chr1",IRanges(1,1000)) # 1kb
gr2 <- GRanges("chr1",IRanges(c(1,1001),c( 500,1500))) # 1kb
gr3 <- GRanges("chr1",IRanges(c(1,1001),c(1000,2000))) # 2kb
gr4 <- GRanges("chr1",IRanges(c(1,1001),c(200,1300))) # 500bp
rowRanges(dds) <- GRangesList(gr1,gr2,gr3,gr4)

# the raw counts
counts(dds)

# the FPM values
fpm(dds)

# the FPKM values
fpkm(dds)

FPM: fragments per million mapped fragments

Description

Calculates either a robust version (default) or the traditional matrix of fragments/counts per million mapped fragments (FPM/CPM). Note: this function is written very simply and can be easily altered to produce other behavior by examining the source code.

Usage

fpm(object, robust = TRUE)

Arguments

Value

a matrix which is normalized per million of mapped fragments, either using the robust median ratio method (robust=TRUE, default) or using raw counts (robust=FALSE).

Seealso

fpkm

Examples

# generate a dataset with size factors: .5, 1, 1, 2
dds <- makeExampleDESeqDataSet(m = 4, n = 1000,
interceptMean=log2(1e3),
interceptSD=0,
sizeFactors=c(.5,1,1,2),
dispMeanRel=function(x) .01)

# add a few rows with very high count
counts(dds)[4:10,] <- 2e5L

# in this robust version, the counts are comparable across samples
round(head(fpm(dds), 3))

# in this column sum version, the counts are still skewed:
# sample1 < sample2 & 3 < sample 4
round(head(fpm(dds, robust=FALSE), 3))

# the column sums of the robust version
# are not equal to 1e6, but the
# column sums of the non-robust version
# are equal to 1e6 by definition

colSums(fpm(dds))/1e6
colSums(fpm(dds, robust=FALSE))/1e6

Shrink log2 fold changes

Description

Adds shrunken log2 fold changes (LFC) and SE to a results table from DESeq run without LFC shrinkage. For consistency with results , the column name lfcSE is used here although what is returned is a posterior SD. Three shrinkage estimators for LFC are available via type (see the vignette for more details on the estimators). The apeglm publication demonstrates that 'apeglm' and 'ashr' outperform the original 'normal' shrinkage estimator.

Usage

lfcShrink(dds, coef, contrast, res, type = c("normal", "apeglm", "ashr"),
  lfcThreshold = 0, svalue = FALSE, returnList = FALSE,
  format = c("DataFrame", "GRanges", "GRangesList"), apeAdapt = TRUE,
  apeMethod = "nbinomCR", parallel = FALSE, BPPARAM = bpparam(),
  quiet = FALSE, ...)

Arguments

Details

Specifying apeglm passes along DESeq2 MLE log2 fold changes and standard errors to the apeglm function in the apeglm package, and re-estimates posterior LFCs for the coefficient specified by coef . Specifying ashr passes along DESeq2 MLE log2 fold changes and standard errors to the ash function in the ashr package, with arguments mixcompdist="normal" and method="shrink" . See vignette for a comparison of shrinkage estimators on an example dataset. For all shrinkage methods, details on the prior is included in priorInfo(res) , including the fitted_g mixture for ashr.

For normal , and design as a formula, shrinkage cannot be applied to coefficients in a model with interaction terms. For normal and user-supplied model matrices, shrinkage is only supported via coef . For normal with numeric- or list-style contrasts, it is possible to use lfcShrink , but likely easier to use DESeq with betaPrior=TRUE followed by results , because the numeric or list should reference the coefficients from the expanded model matrix. These coefficients will be printed to console if 'contrast' is used with normal .

Value

a DESeqResults object with the log2FoldChange and lfcSE columns replaced with shrunken LFC and SE. For consistency with results (and similar to the output of bayesglm ) the column name lfcSE is used here, although what is returned is a posterior SD. For normal and for apeglm the estimate is the posterior mode, for ashr it is the posterior mean. priorInfo(res) contains information about the shrinkage procedure, relevant to the various methods specified by type .

References

Publications for the following shrinkage estimators:

type="normal" :

Love, M.I., Huber, W., Anders, S. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15:550. https://doi.org/10.1186/s13059-014-0550-8

type="apeglm" :

Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895

type="ashr" :

Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2. https://doi.org/10.1093/biostatistics/kxw041

Related work, the bayesglm function in the arm package:

Gelman, A., Jakulin, A., Pittau, M.G. and Su, Y.-S. (2009). A Weakly Informative Default Prior Distribution For Logistic And Other Regression Models. The Annals of Applied Statistics 2:4. http://www.stat.columbia.edu/~gelman/research/published/ priors11.pdf

Examples

set.seed(1)
dds <- makeExampleDESeqDataSet(n=500,betaSD=1)
dds <- DESeq(dds)
res <- results(dds)

# these are the coefficients from the model
# we can specify them using 'coef' by name or number below
resultsNames(dds)

res.shr <- lfcShrink(dds=dds, coef=2)
res.shr <- lfcShrink(dds=dds, contrast=c("condition","B","A"))
res.ape <- lfcShrink(dds=dds, coef=2, type="apeglm")
res.ash <- lfcShrink(dds=dds, coef=2, type="ashr")

Make a simulated DESeqDataSet

Description

Constructs a simulated dataset of Negative Binomial data from two conditions. By default, there are no fold changes between the two conditions, but this can be adjusted with the betaSD argument.

Usage

makeExampleDESeqDataSet(n = 1000, m = 12, betaSD = 0,
  interceptMean = 4, interceptSD = 2, dispMeanRel = function(x) 4/x +
  0.1, sizeFactors = rep(1, m))

Arguments

Value

a DESeqDataSet with true dispersion, intercept and beta values in the metadata columns. Note that the true betas are provided on the log2 scale.

Examples

dds <- makeExampleDESeqDataSet()
dds

Likelihood ratio test (chi-squared test) for GLMs

Description

This function tests for significance of change in deviance between a full and reduced model which are provided as formula . Fitting uses previously calculated sizeFactors (or normalizationFactors ) and dispersion estimates.

Usage

nbinomLRT(object, full = design(object), reduced, betaTol = 1e-08,
  maxit = 100, useOptim = TRUE, quiet = FALSE, useQR = TRUE,
  minmu = 0.5)

Arguments

Details

The difference in deviance is compared to a chi-squared distribution with df = (reduced residual degrees of freedom - full residual degrees of freedom). This function is comparable to the nbinomGLMTest of the previous version of DESeq and an alternative to the default nbinomWaldTest .

Value

a DESeqDataSet with new results columns accessible with the results function. The coefficients and standard errors are reported on a log2 scale.

Seealso

DESeq , nbinomWaldTest

Examples

dds <- makeExampleDESeqDataSet()
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomLRT(dds, reduced = ~ 1)
res <- results(dds)

Wald test for the GLM coefficients

Description

This function tests for significance of coefficients in a Negative Binomial GLM, using previously calculated sizeFactors (or normalizationFactors ) and dispersion estimates. See DESeq for the GLM formula.

Usage

nbinomWaldTest(object, betaPrior = FALSE, betaPriorVar,
  modelMatrix = NULL, modelMatrixType, betaTol = 1e-08, maxit = 100,
  useOptim = TRUE, quiet = FALSE, useT = FALSE, df, useQR = TRUE,
  minmu = 0.5)

Arguments

Details

The fitting proceeds as follows: standard maximum likelihood estimates for GLM coefficients (synonymous with "beta", "log2 fold change", "effect size") are calculated. Then, optionally, a zero-centered Normal prior distribution ( betaPrior ) is assumed for the coefficients other than the intercept.

Note that this posterior log2 fold change estimation is now not the default setting for nbinomWaldTest , as the standard workflow for coefficient shrinkage has moved to an additional function link{lfcShrink} .

For calculating Wald test p-values, the coefficients are scaled by their standard errors and then compared to a standard Normal distribution. The results function without any arguments will automatically perform a contrast of the last level of the last variable in the design formula over the first level. The contrast argument of the results function can be used to generate other comparisons.

The Wald test can be replaced with the nbinomLRT for an alternative test of significance.

Notes on the log2 fold change prior:

The variance of the prior distribution for each non-intercept coefficient is calculated using the observed distribution of the maximum likelihood coefficients. The final coefficients are then maximum a posteriori estimates using this prior (Tikhonov/ridge regularization). See below for details on the prior variance and the Methods section of the DESeq2 manuscript for more detail. The use of a prior has little effect on genes with high counts and helps to moderate the large spread in coefficients for genes with low counts.

The prior variance is calculated by matching the 0.05 upper quantile of the observed MLE coefficients to a zero-centered Normal distribution. In a change of methods since the 2014 paper, the weighted upper quantile is calculated using the wtd.quantile function from the Hmisc package. The weights are the inverse of the expected variance of log counts, so the inverse of $1/ar{mu} + lpha_{tr}$ using the mean of normalized counts and the trended dispersion fit. The weighting ensures that noisy estimates of log fold changes from small count genes do not overly influence the calculation of the prior variance. See estimateBetaPriorVar . The final prior variance for a factor level is the average of the estimated prior variance over all contrasts of all levels of the factor.

When a log2 fold change prior is used (betaPrior=TRUE), then nbinomWaldTest will by default use expanded model matrices, as described in the modelMatrixType argument, unless this argument is used to override the default behavior. This ensures that log2 fold changes will be independent of the choice of reference level. In this case, the beta prior variance for each factor is calculated as the average of the mean squared maximum likelihood estimates for each level and every possible contrast.

Value

a DESeqDataSet with results columns accessible with the results function. The coefficients and standard errors are reported on a log2 scale.

Seealso

DESeq , nbinomLRT

Examples

dds <- makeExampleDESeqDataSet()
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)
res <- results(dds)

Normalized counts transformation

Description

A simple function for creating a DESeqTransform object after applying: f(count(dds,normalized=TRUE) + pc) .

Usage

normTransform(object, f = log2, pc = 1)

Arguments

Seealso

varianceStabilizingTransformation , rlog

Accessor functions for the normalization factors in a DESeqDataSet object.

Description

Gene-specific normalization factors for each sample can be provided as a matrix, which will preempt sizeFactors . In some experiments, counts for each sample have varying dependence on covariates, e.g. on GC-content for sequencing data run on different days, and in this case it makes sense to provide gene-specific factors for each sample rather than a single size factor.

Usage

normalizationFactors(object, ...)
normalizationFactors(object, ...) <- value
list(list("normalizationFactors"), list("DESeqDataSet"))(object)
list(list("normalizationFactors"), list("DESeqDataSet,matrix"))(object) <- value

Arguments

Details

Normalization factors alter the model of DESeq in the following way, for counts $K{ij}$ and normalization factors $NF{ij}$ for gene i and sample j:

$$K{ij} im extrm{NB}( mu{ij}, lpha_i)$$

$$mu{ij} = NF{ij} q_{ij}$$

Note

Normalization factors are on the scale of the counts (similar to sizeFactors ) and unlike offsets, which are typically on the scale of the predictors (in this case, log counts). Normalization factors should include library size normalization. They should have row-wise geometric mean near 1, as is the case with size factors, such that the mean of normalized counts is close to the mean of unnormalized counts. See example code below.

Examples

dds <- makeExampleDESeqDataSet(n=100, m=4)

normFactors <- matrix(runif(nrow(dds)*ncol(dds),0.5,1.5),
ncol=ncol(dds),nrow=nrow(dds),
dimnames=list(1:nrow(dds),1:ncol(dds)))

# the normalization factors matrix should not have 0's in it
# it should have geometric mean near 1 for each row
normFactors <- normFactors / exp(rowMeans(log(normFactors)))
normalizationFactors(dds) <- normFactors

dds <- DESeq(dds)

Normalize for gene length

Description

Normalize for gene length using the output of transcript abundance estimators

Usage

normalizeGeneLength(...)

Arguments

Details

This function is deprecated and moved to a new general purpose package, tximport, which will be added to Bioconductor.

Plot of normalized counts for a single gene

Description

Normalized counts plus a pseudocount of 0.5 are shown by default.

Usage

plotCounts(dds, gene, intgroup = "condition", normalized = TRUE,
  transform = TRUE, main, xlab = "group", returnData = FALSE,
  replaced = FALSE, pc, ...)

Arguments

Examples

dds <- makeExampleDESeqDataSet()
plotCounts(dds, "gene1")

Plot dispersion estimates

Description

A simple helper function that plots the per-gene dispersion estimates together with the fitted mean-dispersion relationship.

Usage

list(list("plotDispEsts"), list("DESeqDataSet"))(object, ymin, CV = FALSE,
  genecol = "black", fitcol = "red", finalcol = "dodgerblue",
  legend = TRUE, xlab, ylab, log = "xy", cex = 0.45, ...)

Arguments

Author

Simon Anders

Examples

dds <- makeExampleDESeqDataSet()
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
plotDispEsts(dds)

MA-plot from base means and log fold changes

Description

A simple helper function that makes a so-called "MA-plot", i.e. a scatter plot of log2 fold changes (on the y-axis) versus the mean of normalized counts (on the x-axis).

Usage

list(list("plotMA"), list("DESeqDataSet"))(object, alpha = 0.1, main = "",
  xlab = "mean of normalized counts", ylim, MLE = FALSE, ...)
list(list("plotMA"), list("DESeqResults"))(object, alpha, main = "",
  xlab = "mean of normalized counts", ylim, MLE = FALSE, ...)

Arguments

Details

This function is essentially two lines of code: building a data.frame and passing this to the plotMA method for data.frame from the geneplotter package. The code of this function can be seen with: getMethod("plotMA","DESeqDataSet") If users wish to modify the graphical parameters of the plot, it is recommended to build the data.frame in the same manner and call plotMA . If the object contains a column svalue then these will be used for coloring the points (with a default alpha=0.005 ).

Author

Michael Love

Examples

dds <- makeExampleDESeqDataSet()
dds <- DESeq(dds)
plotMA(dds)
res <- results(dds)
plotMA(res)

Sample PCA plot for transformed data

Description

This plot helps to check for batch effects and the like.

Usage

list(list("plotPCA"), list("DESeqTransform"))(object, intgroup = "condition",
  ntop = 500, returnData = FALSE)

Arguments

Value

An object created by ggplot , which can be assigned and further customized.

Note

See the vignette for an example of variance stabilization and PCA plots. Note that the source code of plotPCA is very simple. The source can be found by typing DESeq2:::plotPCA.DESeqTransform or getMethod("plotPCA","DESeqTransform") , or browsed on github at https://github.com/mikelove/DESeq2/blob/master/R/plots.R Users should find it easy to customize this function.

Author

Wolfgang Huber

Examples

# using rlog transformed data:
dds <- makeExampleDESeqDataSet(betaSD=1)
rld <- rlog(dds)
plotPCA(rld)

# also possible to perform custom transformation:
dds <- estimateSizeFactors(dds)
# shifted log of normalized counts
se <- SummarizedExperiment(log2(counts(dds, normalized=TRUE) + 1),
colData=colData(dds))
# the call to DESeqTransform() is needed to
# trigger our plotPCA method.
plotPCA( DESeqTransform( se ) )

Sparsity plot

Description

A simple plot of the concentration of counts in a single sample over the sum of counts per gene. Not technically the same as "sparsity", but this plot is useful diagnostic for datasets which might not fit a negative binomial assumption: genes with many zeros and individual very large counts are difficult to model with the negative binomial distribution.

Usage

plotSparsity(x, normalized = TRUE, ...)

Arguments

Examples

dds <- makeExampleDESeqDataSet(n=1000,m=4,dispMeanRel=function(x) .5)
dds <- estimateSizeFactors(dds)
plotSparsity(dds)

Accessors for the 'priorInfo' slot of a DESeqResults object.

Description

The priorInfo slot contains details about the prior on log fold changes

Usage

priorInfo(object, ...)
priorInfo(object, ...) <- value
list(list("priorInfo"), list("DESeqResults"))(object)
list(list("priorInfo"), list("DESeqResults,list"))(object) <- value

Arguments

Replace outliers with trimmed mean

Description

Note that this function is called within DESeq , so is not necessary to call on top of a DESeq call. See the minReplicatesForReplace argument documented in link{DESeq} .

Usage

replaceOutliers(object, trim = 0.2, cooksCutoff, minReplicates = 7,
  whichSamples)
replaceOutliersWithTrimmedMean(object, trim = 0.2, cooksCutoff,
  minReplicates = 7, whichSamples)

Arguments

Details

This function replaces outlier counts flagged by extreme Cook's distances, as calculated by DESeq , nbinomWaldTest or nbinomLRT , with values predicted by the trimmed mean over all samples (and adjusted by size factor or normalization factor). This function replaces the counts in the matrix returned by counts(dds) and the Cook's distances in assays(dds)[["cooks"]] . Original counts are preserved in assays(dds)[["originalCounts"]] .

The DESeq function calculates a diagnostic measure called Cook's distance for every gene and every sample. The results function then sets the p-values to NA for genes which contain an outlying count as defined by a Cook's distance above a threshold. With many degrees of freedom, i.e. many more samples than number of parameters to be estimated-- it might be undesirable to remove entire genes from the analysis just because their data include a single count outlier. An alternate strategy is to replace the outlier counts with the trimmed mean over all samples, adjusted by the size factor or normalization factor for that sample. The following simple function performs this replacement for the user, for samples which have at least minReplicates number of replicates (including that sample). For more information on Cook's distance, please see the two sections of the vignette: 'Dealing with count outliers' and 'Count outlier detection'.

Value

a DESeqDataSet with replaced counts in the slot returned by counts and the original counts preserved in assays(dds)[["originalCounts"]]

Seealso

DESeq

Extract results from a DESeq analysis

Description

results extracts a result table from a DESeq analysis giving base means across samples, log2 fold changes, standard errors, test statistics, p-values and adjusted p-values; resultsNames returns the names of the estimated effects (coefficents) of the model; removeResults returns a DESeqDataSet object with results columns removed.

Usage

results(object, contrast, name, lfcThreshold = 0,
  altHypothesis = c("greaterAbs", "lessAbs", "greater", "less"),
  listValues = c(1, -1), cooksCutoff, independentFiltering = TRUE,
  alpha = 0.1, filter, theta, pAdjustMethod = "BH", filterFun,
  format = c("DataFrame", "GRanges", "GRangesList"), test,
  addMLE = FALSE, tidy = FALSE, parallel = FALSE,
  BPPARAM = bpparam(), minmu = 0.5)
resultsNames(object)
removeResults(object)

Arguments

• a character vector with exactly three elements: the name of a factor in the design formula, the name of the numerator level for the fold change, and the name of the denominator level for the fold change (simplest case)

• a list of 2 character vectors: the names of the fold changes for the numerator, and the names of the fold changes for the denominator. these names should be elements of resultsNames(object) . if the list is length 1, a second element is added which is the empty character vector, character() . (more general case, can be to combine interaction terms and main effects)

• a numeric contrast vector with one element for each element in resultsNames(object) (most general case) If specified, the name argument is ignored. |name | the name of the individual effect (coefficient) for building a results table. Use this argument rather than contrast for continuous variables, individual effects or for individual interaction terms. The value provided to name must be an element of resultsNames(object) .| |lfcThreshold | a non-negative value which specifies a log2 fold change threshold. The default value is 0, corresponding to a test that the log2 fold changes are equal to zero. The user can specify the alternative hypothesis using the altHypothesis argument, which defaults to testing for log2 fold changes greater in absolute value than a given threshold. If lfcThreshold is specified, the results are for Wald tests, and LRT p-values will be overwritten.| |altHypothesis | character which specifies the alternative hypothesis, i.e. those values of log2 fold change which the user is interested in finding. The complement of this set of values is the null hypothesis which will be tested. If the log2 fold change specified by name or by contrast is written as $ eta $ , then the possible values for altHypothesis represent the following alternate hypotheses: |

|* greaterAbs: $|eta| > extrm{lfcThreshold} $ , and p-values are two-tailed |

|* lessAbs: $ |eta| < extrm{lfcThreshold} $ , p-values are the maximum of the upper and lower tests. The Wald statistic given is positive, an SE-scaled distance from the closest boundary |

• greater: $ eta > extrm{lfcThreshold} $

• less: $ eta < - extrm{lfcThreshold} $ |listValues | only used if a list is provided to contrast : a numeric of length two: the log2 fold changes in the list are multiplied by these values. the first number should be positive and the second negative. by default this is c(1,-1)| |cooksCutoff | theshold on Cook's distance, such that if one or more samples for a row have a distance higher, the p-value for the row is set to NA. The default cutoff is the .99 quantile of the F(p, m-p) distribution, where p is the number of coefficients being fitted and m is the number of samples. Set to Inf or FALSE to disable the resetting of p-values to NA. Note: this test excludes the Cook's distance of samples belonging to experimental groups with only 2 samples.| |independentFiltering | logical, whether independent filtering should be applied automatically| |alpha | the significance cutoff used for optimizing the independent filtering (by default 0.1). If the adjusted p-value cutoff (FDR) will be a value other than 0.1, alpha should be set to that value.| |filter | the vector of filter statistics over which the independent filtering will be optimized. By default the mean of normalized counts is used.| |theta | the quantiles at which to assess the number of rejections from independent filtering| |pAdjustMethod | the method to use for adjusting p-values, see ?p.adjust| |filterFun | an optional custom function for performing independent filtering and p-value adjustment, with arguments res (a DESeqResults object), filter (the quantitity for filtering tests), alpha (the target FDR), pAdjustMethod . This function should return a DESeqResults object with a padj column.| |format | character, either "DataFrame" , "GRanges" , or "GRangesList" , whether the results should be printed as a DESeqResults DataFrame, or if the results DataFrame should be attached as metadata columns to the GRanges or GRangesList rowRanges of the DESeqDataSet . If the rowRanges is a GRangesList , and GRanges is requested, the range of each gene will be returned| |test | this is automatically detected internally if not provided. the one exception is after nbinomLRT has been run, test="Wald" will generate Wald statistics and Wald test p-values.| |addMLE | if betaPrior=TRUE was used, whether the "unshrunken" maximum likelihood estimates (MLE) of log2 fold change should be added as a column to the results table (default is FALSE). This argument is preserved for backward compatability, as now the recommended pipeline is to generate shrunken MAP estimates using lfcShrink . This argument functionality is only implemented for contrast specified as three element character vectors.| |tidy | whether to output the results table with rownames as a first column 'row'. the table will also be coerced to data.frame| |parallel | if FALSE, no parallelization. if TRUE, parallel execution using BiocParallel , see next argument BPPARAM| |BPPARAM | an optional parameter object passed internally to bplapply when parallel=TRUE . If not specified, the parameters last registered with register will be used.| |minmu | lower bound on the estimated count (used when calculating contrasts)|

Details

The results table when printed will provide the information about the comparison, e.g. "log2 fold change (MAP): condition treated vs untreated", meaning that the estimates are of log2(treated / untreated), as would be returned by contrast=c("condition","treated","untreated") . Multiple results can be returned for analyses beyond a simple two group comparison, so results takes arguments contrast and name to help the user pick out the comparisons of interest for printing a results table. The use of the contrast argument is recommended for exact specification of the levels which should be compared and their order.

If results is run without specifying contrast or name , it will return the comparison of the last level of the last variable in the design formula over the first level of this variable. For example, for a simple two-group comparison, this would return the log2 fold changes of the second group over the first group (the reference level). Please see examples below and in the vignette.

The argument contrast can be used to generate results tables for any comparison of interest, for example, the log2 fold change between two levels of a factor, and its usage is described below. It can also accomodate more complicated numeric comparisons. The test statistic used for a contrast is:

$$c^t eta / qrt{c^t Sigma c }$$

The argument name can be used to generate results tables for individual effects, which must be individual elements of resultsNames(object) . These individual effects could represent continuous covariates, effects for individual levels, or individual interaction effects.

Information on the comparison which was used to build the results table, and the statistical test which was used for p-values (Wald test or likelihood ratio test) is stored within the object returned by results . This information is in the metadata columns of the results table, which is accessible by calling mcols on the DESeqResults object returned by results .

On p-values:

By default, independent filtering is performed to select a set of genes for multiple test correction which maximizes the number of adjusted p-values less than a given critical value alpha (by default 0.1). See the reference in this man page for details on independent filtering. The filter used for maximizing the number of rejections is the mean of normalized counts for all samples in the dataset. Several arguments from the filtered_p function of the genefilter package (used within the results function) are provided here to control the independent filtering behavior. In DESeq2 version >= 1.10, the threshold that is chosen is the lowest quantile of the filter for which the number of rejections is close to the peak of a curve fit to the number of rejections over the filter quantiles. 'Close to' is defined as within 1 residual standard deviation. The adjusted p-values for the genes which do not pass the filter threshold are set to NA .

By default, results assigns a p-value of NA to genes containing count outliers, as identified using Cook's distance. See the cooksCutoff argument for control of this behavior. Cook's distances for each sample are accessible as a matrix "cooks" stored in the assays() list. This measure is useful for identifying rows where the observed counts might not fit to a Negative Binomial distribution.

For analyses using the likelihood ratio test (using nbinomLRT ), the p-values are determined solely by the difference in deviance between the full and reduced model formula. A single log2 fold change is printed in the results table for consistency with other results table outputs, however the test statistic and p-values may nevertheless involve the testing of one or more log2 fold changes. Which log2 fold change is printed in the results table can be controlled using the name argument, or by default this will be the estimated coefficient for the last element of resultsNames(object) .

If useT=TRUE was specified when running DESeq or nbinomWaldTest , then the p-value generated by results will also make use of the t distribution for the Wald statistic, using the degrees of freedom in mcols(object)$tDegreesFreedom .

Value

For results : a DESeqResults object, which is a simple subclass of DataFrame. This object contains the results columns: baseMean , log2FoldChange , lfcSE , stat , pvalue and padj , and also includes metadata columns of variable information. The lfcSE gives the standard error of the log2FoldChange . For the Wald test, stat is the Wald statistic: the log2FoldChange divided by lfcSE , which is compared to a standard Normal distribution to generate a two-tailed pvalue . For the likelihood ratio test (LRT), stat is the difference in deviance between the reduced model and the full model, which is compared to a chi-squared distribution to generate a pvalue .

For resultsNames : the names of the columns available as results, usually a combination of the variable name and a level

For removeResults : the original DESeqDataSet with results metadata columns removed

Seealso

DESeq , lfcShrink

References

Richard Bourgon, Robert Gentleman, Wolfgang Huber: Independent filtering increases detection power for high-throughput experiments. PNAS (2010), http://dx.doi.org/10.1073/pnas.0914005107

Examples

## Example 1: two-group comparison

dds <- makeExampleDESeqDataSet(m=4)

dds <- DESeq(dds)
res <- results(dds, contrast=c("condition","B","A"))

# with more than two groups, the call would look similar, e.g.:
# results(dds, contrast=c("condition","C","A"))
# etc.

## Example 2: two conditions, two genotypes, with an interaction term

dds <- makeExampleDESeqDataSet(n=100,m=12)
dds$genotype <- factor(rep(rep(c("I","II"),each=3),2))

design(dds) <- ~ genotype + condition + genotype:condition
dds <- DESeq(dds)
resultsNames(dds)

# Note: design with interactions terms by default have betaPrior=FALSE

# the condition effect for genotype I (the main effect)
results(dds, contrast=c("condition","B","A"))

# the condition effect for genotype II
# this is, by definition, the main effect *plus* the interaction term
# (the extra condition effect in genotype II compared to genotype I).
results(dds, list( c("condition_B_vs_A","genotypeII.conditionB") ))

# the interaction term, answering: is the condition effect *different* across genotypes?
results(dds, name="genotypeII.conditionB")

## Example 3: two conditions, three genotypes

# ~~~ Using interaction terms ~~~

dds <- makeExampleDESeqDataSet(n=100,m=18)
dds$genotype <- factor(rep(rep(c("I","II","III"),each=3),2))
design(dds) <- ~ genotype + condition + genotype:condition
dds <- DESeq(dds)
resultsNames(dds)

# the condition effect for genotype I (the main effect)
results(dds, contrast=c("condition","B","A"))

# the condition effect for genotype III.
# this is the main effect *plus* the interaction term
# (the extra condition effect in genotype III compared to genotype I).
results(dds, contrast=list( c("condition_B_vs_A","genotypeIII.conditionB") ))

# the interaction term for condition effect in genotype III vs genotype I.
# this tests if the condition effect is different in III compared to I
results(dds, name="genotypeIII.conditionB")

# the interaction term for condition effect in genotype III vs genotype II.
# this tests if the condition effect is different in III compared to II
results(dds, contrast=list("genotypeIII.conditionB", "genotypeII.conditionB"))

# Note that a likelihood ratio could be used to test if there are any
# differences in the condition effect between the three genotypes.

# ~~~ Using a grouping variable ~~~

# This is a useful construction when users just want to compare
# specific groups which are combinations of variables.

dds$group <- factor(paste0(dds$genotype, dds$condition))
design(dds) <- ~ group
dds <- DESeq(dds)
resultsNames(dds)

# the condition effect for genotypeIII
results(dds, contrast=c("group", "IIIB", "IIIA"))

Apply a 'regularized log' transformation

Description

This function transforms the count data to the log2 scale in a way which minimizes differences between samples for rows with small counts, and which normalizes with respect to library size. The rlog transformation produces a similar variance stabilizing effect as varianceStabilizingTransformation , though rlog is more robust in the case when the size factors vary widely. The transformation is useful when checking for outliers or as input for machine learning techniques such as clustering or linear discriminant analysis. rlog takes as input a DESeqDataSet and returns a RangedSummarizedExperiment object.

Usage

rlog(object, blind = TRUE, intercept, betaPriorVar,
  fitType = "parametric")
rlogTransformation(object, blind = TRUE, intercept, betaPriorVar,
  fitType = "parametric")

Arguments

Details

Note that neither rlog transformation nor the VST are used by the differential expression estimation in DESeq , which always occurs on the raw count data, through generalized linear modeling which incorporates knowledge of the variance-mean dependence. The rlog transformation and VST are offered as separate functionality which can be used for visualization, clustering or other machine learning tasks. See the transformation section of the vignette for more details, including a statement on timing. If rlog is run on data with number of samples in [30-49] it will print a message that it may take a few minutes, if the number of samples is 50 or larger, it will print a message that it may take a "long time", and in both cases, it will mention that the vst is a much faster transformation.

The transformation does not require that one has already estimated size factors and dispersions.

The regularization is on the log fold changes of the count for each sample over an intercept, for each gene. As nearby count values for low counts genes are almost as likely as the observed count, the rlog shrinkage is greater for low counts. For high counts, the rlog shrinkage has a much weaker effect. The fitted dispersions are used rather than the MAP dispersions (so similar to the varianceStabilizingTransformation ).

The prior variance for the shrinkag of log fold changes is calculated as follows: a matrix is constructed of the logarithm of the counts plus a pseudocount of 0.5, the log of the row means is then subtracted, leaving an estimate of the log fold changes per sample over the fitted value using only an intercept. The prior variance is then calculated by matching the upper quantiles of the observed log fold change estimates with an upper quantile of the normal distribution. A GLM fit is then calculated using this prior. It is also possible to supply the variance of the prior. See the vignette for an example of the use and a comparison with varianceStabilizingTransformation .

The transformed values, rlog(K), are equal to $rlog(K{ij}) = log_2(q{ij}) = eta{i0} + eta{ij}$ , with formula terms defined in DESeq .

The parameters of the rlog transformation from a previous dataset can be frozen and reapplied to new samples. See the 'Data quality assessment' section of the vignette for strategies to see if new samples are sufficiently similar to previous datasets. The frozen rlog is accomplished by saving the dispersion function, beta prior variance and the intercept from a previous dataset, and running rlog with 'blind' set to FALSE (see example below).

Value

a DESeqTransform if a DESeqDataSet was provided, or a matrix if a count matrix was provided as input. Note that for DESeqTransform output, the matrix of transformed values is stored in assay(rld) . To avoid returning matrices with NA values, in the case of a row of all zeros, the rlog transformation returns zeros (essentially adding a pseudocount of 1 only to these rows).

Seealso

plotPCA , varianceStabilizingTransformation , normTransform

References

Reference for regularized logarithm (rlog):

Michael I Love, Wolfgang Huber, Simon Anders: Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 2014, 15:550. http://dx.doi.org/10.1186/s13059-014-0550-8

Examples

dds <- makeExampleDESeqDataSet(m=6,betaSD=1)
rld <- rlog(dds)
dists <- dist(t(assay(rld)))
# plot(hclust(dists))

Show method for DESeqResults objects

Description

Prints out the information from the metadata columns of the results object regarding the log2 fold changes and p-values, then shows the DataFrame using the standard method.

Usage

list(list("show"), list("DESeqResults"))(object)

Arguments

Author

Michael Love

Accessor functions for the 'sizeFactors' information in a DESeqDataSet object.

Description

The sizeFactors vector assigns to each column of the count matrix a value, the size factor, such that count values in the columns can be brought to a common scale by dividing by the corresponding size factor (as performed by counts(dds, normalized=TRUE) ). See DESeq for a description of the use of size factors. If gene-specific normalization is desired for each sample, use normalizationFactors .

Usage

list(list("sizeFactors"), list("DESeqDataSet"))(object)
list(list("sizeFactors"), list("DESeqDataSet,numeric"))(object) <- value

Arguments

Seealso

estimateSizeFactors

Author

Simon Anders

Summarize DESeq results

Description

Print a summary of the results from a DESeq analysis.

Usage

list(list("summary"), list("DESeqResults"))(object, alpha, list())

Arguments

Author

Michael Love

Examples

dds <- makeExampleDESeqDataSet(m=4)
dds <- DESeq(dds)
res <- results(dds)
summary(res)

Unmix samples using loss in a variance stabilized space

Description

Unmixes samples in x according to pure components, using numerical optimization. The components in pure are added on the scale of gene expression (either normalized counts, or TPMs). The loss function when comparing fitted expression to the samples in x occurs in a variance stabilized space. This task is sometimes referred to as "deconvolution", and can be used, for example, to identify contributions from various tissues. Note: some groups have found that the mixing contributions may be more accurate if very lowly expressed genes across x and pure are first removed. We have not explored this fully. Note: if the pbapply package is installed a progress bar will be displayed while mixing components are fit.

Usage

unmix(x, pure, alpha, shift, power = 1, format = "matrix",
  quiet = FALSE)

Arguments

Value

a matrix, the mixture components for each sample in x (rows). The "pure" samples make up the columns, and so each row sums to 1. If colnames existed on the input matrices they will be propagated to the output matrix. If format="list" , then a list, containing as elements: (1) the matrix of mixture components, (2) the correlations in the variance stabilized space of the fitted samples to the samples in x , and (3) the fitted samples as a matrix with the same dimension as x .

Examples

# some artificial data
cts <- matrix(c(80,50,1,100,
1,1,60,100,
0,50,60,100), ncol=4, byrow=TRUE)
# make a DESeqDataSet
dds <- DESeqDataSetFromMatrix(cts,
data.frame(row.names=seq_len(ncol(cts))), ~1)
colnames(dds) <- paste0("sample",1:4)

# note! here you would instead use
# estimateSizeFactors() to do actual normalization
sizeFactors(dds) <- rep(1, ncol(dds))

norm.cts <- counts(dds, normalized=TRUE)

# 'pure' should also have normalized counts...
pure <- matrix(c(10,0,0,
0,0,10,
0,10,0), ncol=3, byrow=TRUE)
colnames(pure) <- letters[1:3]

# for real data, you need to find alpha after fitting estimateDispersions()
mix <- unmix(norm.cts, pure, alpha=0.01)

Apply a variance stabilizing transformation (VST) to the count data

Description

This function calculates a variance stabilizing transformation (VST) from the fitted dispersion-mean relation(s) and then transforms the count data (normalized by division by the size factors or normalization factors), yielding a matrix of values which are now approximately homoskedastic (having constant variance along the range of mean values). The transformation also normalizes with respect to library size. The rlog is less sensitive to size factors, which can be an issue when size factors vary widely. These transformations are useful when checking for outliers or as input for machine learning techniques such as clustering or linear discriminant analysis.

Usage

varianceStabilizingTransformation(object, blind = TRUE,
  fitType = "parametric")
getVarianceStabilizedData(object)

Arguments

Details

For each sample (i.e., column of counts(dds) ), the full variance function is calculated from the raw variance (by scaling according to the size factor and adding the shot noise). We recommend a blind estimation of the variance function, i.e., one ignoring conditions. This is performed by default, and can be modified using the 'blind' argument.

Note that neither rlog transformation nor the VST are used by the differential expression estimation in DESeq , which always occurs on the raw count data, through generalized linear modeling which incorporates knowledge of the variance-mean dependence. The rlog transformation and VST are offered as separate functionality which can be used for visualization, clustering or other machine learning tasks. See the transformation section of the vignette for more details.

The transformation does not require that one has already estimated size factors and dispersions.

A typical workflow is shown in Section Variance stabilizing transformation in the package vignette.

If estimateDispersions was called with:

fitType="parametric" , a closed-form expression for the variance stabilizing transformation is used on the normalized count data. The expression can be found in the file vst.pdf which is distributed with the vignette.

fitType="local" , the reciprocal of the square root of the variance of the normalized counts, as derived from the dispersion fit, is then numerically integrated, and the integral (approximated by a spline function) is evaluated for each count value in the column, yielding a transformed value.

fitType="mean" , a VST is applied for Negative Binomial distributed counts, 'k', with a fixed dispersion, 'a': ( 2 asinh(sqrt(a k)) - log(a) - log(4) )/log(2).

In all cases, the transformation is scaled such that for large counts, it becomes asymptotically (for large values) equal to the logarithm to base 2 of normalized counts.

The variance stabilizing transformation from a previous dataset can be frozen and reapplied to new samples. See the 'Data quality assessment' section of the vignette for strategies to see if new samples are sufficiently similar to previous datasets. The frozen VST is accomplished by saving the dispersion function accessible with dispersionFunction , assigning this to the DESeqDataSet with the new samples, and running varianceStabilizingTransformation with 'blind' set to FALSE (see example below). Then the dispersion function from the previous dataset will be used to transform the new sample(s).

Limitations: In order to preserve normalization, the same transformation has to be used for all samples. This results in the variance stabilizition to be only approximate. The more the size factors differ, the more residual dependence of the variance on the mean will be found in the transformed data. rlog is a transformation which can perform better in these cases. As shown in the vignette, the function meanSdPlot from the package vsn can be used to see whether this is a problem.

Value

varianceStabilizingTransformation returns a DESeqTransform if a DESeqDataSet was provided, or returns a a matrix if a count matrix was provided. Note that for DESeqTransform output, the matrix of transformed values is stored in assay(vsd) . getVarianceStabilizedData also returns a matrix.

Seealso

plotPCA , rlog , normTransform

Author

Simon Anders

References

Reference for the variance stabilizing transformation for counts with a dispersion trend:

Simon Anders, Wolfgang Huber: Differential expression analysis for sequence count data. Genome Biology 2010, 11:106. http://dx.doi.org/10.1186/gb-2010-11-10-r106

Examples

dds <- makeExampleDESeqDataSet(m=6)
vsd <- varianceStabilizingTransformation(dds)
dists <- dist(t(assay(vsd)))
# plot(hclust(dists))

Quickly estimate dispersion trend and apply a variance stabilizing transformation

Description

This is a wrapper for the varianceStabilizingTransformation (VST) that provides much faster estimation of the dispersion trend used to determine the formula for the VST. The speed-up is accomplished by subsetting to a smaller number of genes in order to estimate this dispersion trend. The subset of genes is chosen deterministically, to span the range of genes' mean normalized count. This wrapper for the VST is not blind to the experimental design: the sample covariate information is used to estimate the global trend of genes' dispersion values over the genes' mean normalized count. It can be made strictly blind to experimental design by first assigning a design of ~1 before running this function, or by avoiding subsetting and using varianceStabilizingTransformation .

Usage

vst(object, blind = TRUE, nsub = 1000, fitType = "parametric")

Arguments

Value

a DESeqTranform object or a matrix of transformed, normalized counts

Examples

dds <- makeExampleDESeqDataSet(n=2000, m=20)
vsd <- vst(dds)