deseq2

Why un-normalized counts?

As input, the DESeq2 package expects count data as obtained, e.g. from RNA-seq or another high-throughput sequencing experiment, in the form of a matrix of integer values. The value in the i-th row and the j-th column of the matrix tells how many reads can be assigned to gene i in sample j. Analogously, for other types of assays, the rows of the matrix might correspond e.g. to binding regions (with ChIP-Seq) or peptide sequences (with quantitative mass spectrometry). We will list method for obtaining count matrices in sections below.

The values in the matrix should be un-normalized counts or estimated counts of sequencing reads (for single-end RNA-seq) or fragments (for paired-end RNA-seq). The RNA-seq workflow describes multiple techniques for preparing such count matrices. It is important to provide count matrices as input for DESeq2's statistical model [Love2014] to hold, as only the count values allow assessing the measurement precision correctly. The DESeq2 model internally corrects for library size, so transformed or normalized values such as counts scaled by library size should not be used as input.

The DESeqDataSet

The object class used by the DESeq2 package to store the read counts and the intermediate estimated quantities during statistical analysis is the :class:`~DESeqDataSet.DESeqDataSet`, which will usually be represented in the code here as an object dds.

A technical detail is that the :class:`~DESeqDataSet.DESeqDataSet` class extends the :class:`AnnData` class of the anndata package.

A :class:`~DESeqDataSet.DESeqDataSet` object must have an associated design matrix. The design matrix expresses the variables which will be used in modeling. The matrix is built from a formula starting with a tilde (~) followed by the variables with plus signs between them (it will be coerced into a formula if it is not already). The design can be changed later, however then all differential analysis steps should be repeated, as the design formula is used to estimate the dispersions and to estimate the log2 fold changes of the model.

We will now show 2 ways of constructing a :class:`~DESeqDataSet.DESeqDataSet`, depending on what pipeline was used upstream of DESeq2 to generated counts or estimated counts:

1. From a count matrix :ref:`countmat`
2. From an :class:`AnnData` object :ref:`ad`

Count matrix input

The constructor :meth:`.DESeqDataSet.__init__` can be used if you already have a matrix of read counts prepared from another source. Another method for quickly producing count matrices from alignment files is the featureCounts function [Liao2013] in the Rsubread package. To use the constructor of :meth:`.DESeqDataSet.__init__`, the user should provide the counts matrix, the information about the samples (the rows of the count matrix) as a data frame, and the design formula.

To demonstrate the construction of :class:`~DESeqDataSet.DESeqDataSet` from a count matrix, we will read in count data from the pasilla package, available in inmoose. We read in a count matrix, which we will name cts, and the sample information table, which we will name sample_data. Further below we describe how to extract these objects from, e.g. featureCounts output:

We examine the count matrix and sample data to see if they are consistent in terms of sample order:

Note that these are not in the same order with respect to samples!

It is absolutely critical that the rows of the count matrix and the rows of the sample data (information about samples) are in the same order. DESeq2 will not make guesses as to which row of the count matrix belongs to which row of the sample data, these must be provided to DESeq2 already in consistent order.

As they are not in the correct order as given, we need to re-arrange one or the other so that they are consistent in terms of sample order (if we do not, later functions would produce an error). We additionally need to chop off the "fb" of the row names of sample_data, so the naming is consistent:

If you have used the featureCounts function [Liao2013] in the Rsubread package, the matrix of read counts can be directly provided from the "counts" element in the list output. The count matrix and sample data can typically be read into Python from flat files using import functions from pandas or numpy.

With the count matrix, cts, and the sample information, sample_data, we can construct a :class:`~DESeqDataSet.DESeqDataSet`:

If you have additional feature data, it can be added to the :class:`~DESeqDataSet.DESeqDataSet` by adding to the metadata columns of a newly constructed object. (Here we add redundant data just for demonstration, as the gene names are already the rownames of the dds.):

:class:`AnnData` input

If one has already created or obtained an :class:`AnnData`, it can be easily input into DESeq2 as follows. First we load the module containing the airway dataset:

The constructor function below shows the generation of a :class:`~DESeqDataSet.DESeqDataSet` from an :class:`AnnData` ad:

Pre-filtering

While it is not necessary to pre-filter low count genes before running the DESeq2 functions, there are two reasons which make pre-filtering useful: by removing rows in which there are very few reads, we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2. It can also improve visualizations, as features with no information for differential expression are not plotted.

Here we perform a minimal pre-filtering to keep only rows that have at least 10 reads total. Note that more strict filtering to increase power is automatically applied via :ref:`independent filtering<indfilt>` on the mean of normalized counts within the :meth:`.DESeqDataSet.results` function:

Alternatively, a popular filter is to ensure at least X samples with a count of 10 or more, where X can be chosen as the sample size of the smallest group of samples:

>>> keep = (dds.counts() >= 10).sum(axis=1) >= X
>>> dds = dds[keep,:]

Note on factor levels

By default, Python will choose a reference level for factors based on alphabetical order, it chooses the first value as the reference. Then, if you never tell the DESeq2 functions which level you want to compare against (e.g. which level represents the control group), the comparisons will be based on the alphabetical order of the levels. There are two solutions: you can either explicitly tell :meth:`.DESeqDataSet.results` which comparison to make using the contrast argument (this will be shown later), or you can explicitly set the factors levels by specifying the desired reference value first. In order to see the change of reference levels reflected in the results names, you need to either run :func:`DESeq` or :func:`nbinomWaldTest` / :func:`nbinomLRT` after the re-leveling operation. Setting the factor levels can be done with the:code:reorder_categories function:

If you need to subset the columns of a :class:`~DESeqDataSet.DESeqDataSet`, i.e. when removing certain samples from the analysis, it is possible that all the samples for one or more levels of a variable in the design formula would be removed. In this case, the remove_unused_categories function can be used to remove those levels which do not have samples in the current :class:`~DESeqDataSet.DESeqDataSet`:

Collapsing technical replicates

DESeq2 provides a function :func:`collapseReplicates` which can assist in combining the counts from technical replicates into single columns of the count matrix. The term technical replicate implies multiple sequencing runs of the same library. You should not collapse biological replicates using this function. See the manual page for an example of the use of :func:`collapseReplicates`.

About the pasilla dataset

We continue with the :doc:`/pasilla` data constructed from the count matrix method above. This data set is from an experiment on Drosophila melanogaster cell cultures and investigated the effect of RNAi knock-down of the splicing factor pasilla [Brooks2011]_. The detailed transcript of the production of the pasilla data is provided in the vignette of the data package pasilla.

Log fold change shrinkage for visualization and ranking

Shrinkage of effect size (LFC estimates) is useful for visualization and ranking of genes. To shrink the LFC, we pass the dds object to the function :func:`lfcShrink`. Below we specify to use the apeglm method for effect size shrinkage [Zhu2018], which improves on the previous estimator.

We provide the dds object and the name or number of the coefficient we want to shrink, where the number refers to the order of the coefficient as it appears in dds.resultsNames():

>>> dds.resultsNames()
>>> resLFC = dds.lfcShrink(coef="condition_treated_vs_untreated", type="apeglm")
>>> resLFC

Shrinkage estimation is discussed more in a :ref:`later section<shrink>`.

p-values and adjusted p-values

We can order our results table by the smallest p-value:

We can summarize some basic tallies using the :meth:`.DESeqResults.summary` function:

How many adjusted p-values were less than 0.1:

The :meth:`.DESeqDataSet.results` function contains a number of arguments to customize the results table which is generated. You can read about these arguments by looking up the documentation of :meth:`.DESeqDataSet.results`. Note that the :meth:`.DESeqDataSet.results` function automatically performs independent filtering based on the mean of normalized counts for each gene, optimizing the number of genes which will have an adjusted p-value below a given FDR cutoff, alpha. Independent filtering is further discussed :ref:`below<indfilt>`. By default the argument alpha is set to 0.1. If the adjusted p-value cutoff will be a value other than 0.1, alpha should be set to that value:

MA-plot

In DESeq2, the function :meth:`.DESeqResults.plotMA` shows the log2 fold changes attributable to a given variable over the mean of normalized counts for all the samples in the :class:`~DESeqDataSet.DESeqDataSet`. Points will be colored red if the adjusted p-value is less than 0.1. Points which fall out of the window are plotted as open triangles pointing either up or down:

It is more useful visualize the MA-plot for the shrunken log2 fold changes, which remove the noise associated with log2 fold changes from low count genes without requiring arbitrary filtering thresholds:

Alternative shrinkage estimators

The moderated log fold changes proposed by [Love2014] use a normal prior distribution, centered on zero and with a scale that is fit to the data. The shrunken log fold changes are useful for ranking and visualization, without the need for arbitrary filters on low count genes. The normal prior can sometimes produce too strong of shrinkage for certain datasets. In DESeq2 version 1.18, we include two additional adaptive shrinkage estimators, available via the type argument of :func:`lfcShrink`.

The options for type are:

• apeglm is the adaptive t prior shrinkage estimator from the apeglm package [Zhu2018]. As of version 1.28.0, it is the default estimator.
• ashr is the adaptive shrinkage estimator from the ashr package [Stephens2016]. Here DESeq2 uses the ashr option to fit a mixture of Normal distributions to form the prior, with method="shrinkage".
• normal: is the the original DESeq2 shrinkage estimator, an adaptive Normal distribution as prior.

If the shrinkage estimator apeglm is used in published research, please cite [Zhu2018].

If the shrinkage estimator ashr is used in published research, please cite [Stephens2016].

In the LFC shrinkage code above, we specified coef="condition_treated_vs_untreated". We can also just specify the coefficient by the order that it appears in dds.resultsNames(), in this case coef=2. For more details explaining how the shrinkage estimators differ, and what kinds of designs, contrasts and output is provided by each, see the :ref:`extended section on shrinkage estimators<moreshrink>`:

>>> dds.resultsNames()
>>> # because we are interested in treated vs untreated, we set 'coef=2'
>>> resNorm = dds.lfcShrink(coef=2, type="normal")
>>> resAsh = dds.lfcShrink(coef=2, type="ashr")

More information on results columns

Information about which variables and tests were used can be found by inspecting the attribute description on the :class:`~results.DESeqResults` object:

For a particular gene, a log2 fold change of -1 for condition treated vs untreated means that the treatment induces a multiplicative change in observed gene expression level of 2^− 1 = 0.5 compared to the untreated condition. If the variable of interest is continuous-valued, then the reported log2 fold change is per unit of change of that variable.

Exporting results to CSV files

A plain-text file of the results can be exported using the method :meth:`.DESeqResults.to_csv`. We suggest using a descriptive file name indicating the variable and level which were tested:

Exporting only the results which pass an adjusted p-value threshold can be accomplished by subsetting the results:

Wald test individual steps

The function :func:`DESeq` runs the following functions in order:

>>> dds = dds.estimateSizeFactors()
>>> dds = dds.estimateDispersions(dds)
>>> dds = nbinomWaldTest(dds)

Control features for estimating size factors

In some experiments, it may not be appropriate to assume that a minority of features (genes) are affected greatly by the condition, such that the standard median-ratio method for estimating the size factors will not provide correct inference (the log fold changes for features that were truly un-changing will not centered on zero). This is a difficult inference problem for any method, but there is an important feature that can be used: the controlGenes argument of :meth:`.DESeqDataSet.estimateSizeFactors`. If there is any prior information about features (genes) that should not be changing with respect to the condition, providing this set of features to controlGenes will ensure that the log fold changes for these features will be centered around 0. The paradigm then becomes:

>>> dds = dds.estimateSizeFactors(controlGenes=ctrlGenes)
>>> dds = DESeq(dds)

Contrasts

A contrast is a linear combination of estimated log2 fold changes, which can be used to test if differences between groups are equal to zero. The simplest use case for contrasts is an experimental design containing a factor with three levels, say A, B and C. Contrasts enable the user to generate results for all 3 possible differences: log2 fold change of B vs A, of C vs A, and of C vs B. The contrast argument of :meth:`.DESeqDataSet.results` function is used to extract test results of log2 fold changes of interest, for example:

>>> dds.results(contrast=["condition","C","B"])

Log2 fold changes can also be added and subtracted by providing a list to the contrast argument which has two elements: the names of the log2 fold changes to add, and the names of the log2 fold changes to subtract. The names used in the list should come from dds.resultsNames(). Alternatively, a numeric vector of the length of dds.resultsNames() can be provided, for manually specifying the linear combination of terms. A tutorial describing the use of numeric contrasts for DESeq2 explains a general approach to comparing across groups of samples. Demonstrations of the use of contrasts for various designs can be found in the examples section of the help page :meth:`.DESeqDataSet.results`. The mathematical formula that is used to generate the contrasts can be found :ref:`below<theory>`.

Interactions

Interaction terms can be added to the design formula, in order to test, for example, if the log2 fold change attributable to a given condition is different based on another factor, for example if the condition effect differs across genotype.

>>> dds.obs["group"] = Factor([f"{g}{c}" for (g,c) in zip(dds.obs["genotype"], dds.obs["condition"])])
>>> dds.design = "~ group"
>>> dds = DESeq(dds)
>>> dds.resultsNames()
>>> dds.results(contrast=["group", "IB", "IA"])

Time-series experiments

There are a number of ways to analyze time-series experiments, depending on the biological question of interest. In order to test for any differences over multiple time points, once can use a design including the time factor, and then test using the likelihood ratio test as described in the following section, where the time factor is removed in the reduced formula. For a control and treatment time series, one can use a design formula containing the condition factor, the time factor, and the interaction of the two. In this case, using the likelihood ratio test with a reduced model which does not contain the interaction terms will test whether the condition induces a change in gene expression at any time point after the reference level time point (time 0). An example of the later analysis is provided in the RNA-seq workflow.

Likelihood ratio test

DESeq2 offers two kinds of hypothesis tests: the Wald test, where we use the estimated standard error of a log2 fold change to test if it is equal to zero, and the likelihood ratio test (LRT). The LRT examines two models for the counts, a full model with a certain number of terms and a reduced model, in which some of the terms of the full model are removed. The test determines if the increased likelihood of the data using the extra terms in the full model is more than expected if those extra terms are truly zero.

The LRT is therefore useful for testing multiple terms at once, for example testing 3 or more levels of a factor at once, or all interactions between two variables. The LRT for count data is conceptually similar to an analysis of variance (ANOVA) calculation in linear regression, except that in the case of the Negative Binomial GLM, we use an analysis of deviance (ANODEV), where the deviance captures the difference in likelihood between a full and a reduced model.

The likelihood ratio test can be performed by specifying test="LRT" when using the :func:`DESeq` function, and providing a reduced design formula, e.g. one in which a number of terms from dds.design are removed. The degrees of freedom for the test is obtained from the difference between the number of parameters in the two models. A simple likelihood ratio test, if the full design was ~condition would look like:

>>> dds = DESeq(dds, test="LRT", reduced="~1")
>>> res = dds.results()

If the full design contained other variables, such as a batch variable, e.g. ~batch + condition then the likelihood ratio test would look like:

>>> dds = DESeq(dds, test="LRT", reduced="~batch")
>>> res = dds.results()

Extended section on shrinkage estimators

Here we extend the :ref:`discussion of shrinkage estimators<shrink>`. Below is a summary table of differences between methods available in :func:`lfcShrink` via the type argument (and for further technical reference on use of arguments please the documentation of :func:`lfcShrink`):

Beginning with the first row, all shrinkage methods provided by DESeq2 are good for ranking genes by "effect size", that is the log2 fold change (LFC) across groups, or associated with an interaction term. It is useful to contrast ranking by effect size with ranking by a p-value or adjusted p-value associated with a null hypothesis: while increasing the number of samples will tend to decrease the associated p-value for a gene that is differentially expressed, the estimated effect size or LFC becomes more precise. Also, a gene can have a small p-value although the change in expression is not great, as long as the standard error associated with the estimated LFC is small.

The next two rows point out that apeglm and ashr shrinkage methods help to preserve the size of large LFC, and can be used to compute s-values. These properties are related. As noted in the :ref:`previous section<shrink>`, the original DESeq2 shrinkage estimator used a Normal distribution, with a scale that adapts to the spread of the observed LFCs. Because the tails of the Normal distribution become thin relatively quickly, it was important when we designed the method that the prior scaling is sensitive to the very largest observed LFCs. As you can read in the DESeq2 paper, under the section, "Empirical prior estimate", we used the top 5% of the LFCs by absolute value to set the scale of the Normal prior (we later added weighting the quantile by precision). ashr, published in 2016, and apeglm use wide-tailed priors to avoid shrinking large LFCs. While a typical RNA-seq experiment may have many LFCs between -1 and 1, we might consider a LFC of >4 to be very large, as they represent 16-fold increases or decreases in expression. ashr and apeglm can adapt to the scale of the entirety of LFCs, while not over-shrinking the few largest LFCs. The potential for over-shrinking LFC is also why DESeq2's shrinkage estimator is not recommended for designs with interaction terms.

What are s-values? This quantity proposed by [Stephens2016] gives the estimated rate of false sign among genes with equal or smaller s-value. [Stephens2016] points out they are analogous to the q-value of [Storey2003]. The s-value has a desirable property relative to the adjusted p-value or q-value, in that it does not require supposing there be a set of null genes with LFC = 0 (the most commonly used null hypothesis). Therefore, it can be benchmarked by comparing estimated LFC and s-value to the "true LFC" in a setting where this can be reasonably defined. For these estimated probabilities to be accurate, the scale of the prior needs to match the scale of the distribution of effect sizes, and so the original DESeq2 shrinkage method is not really compatible with computing s-values.

The last four rows explain differences in whether coefficients or contrasts can have shrinkage applied by the various methods. All three methods can use coef with either the name or numeric index from dds.resultsNames() to specify which coefficient to shrink. All three methods allow for a positive lfcThreshold to be specified, in which case, they will return p-values and adjusted p-values or s-values for the LFC being greater in absolute value than the threshold (see :ref:`this section<thresh>` for normal). For apeglm and ashr, setting a threshold means that the s-values will give the "false sign or small" rate (FSOS) among genes with equal or small s-value. We found FSOS to be a useful description for when the LFC is either the wrong sign or less than the threshold distance from 0.

TODO

Finally, normal and ashr can be used with arbitrary specified contrast because normal shrinks multiple coefficients simultaneously (apeglm does not), and because ashr does not estimate a vector of coefficients but models estimated coefficients and their standard errors from upstream methods (here, DESeq2's MLE). Although apeglm cannot be used with contrast, we note that many designs can be easily rearranged such that what was a contrast becomes its own coefficient. In this case, the dispersion does not have to be estimated again, as the designs are equivalent, up to the meaning of the coefficients. Instead, one need only run nbinomWaldTest to re-estimate MLE coefficients -- these are necessary for apeglm -- and then run lfcShrink specifying the coefficient of interest in dds.resultsNames().

We give some examples below of producing equivalent designs for use with coef. We show how the coefficients change with patsy.dmatrix, but the user would, for example, either change the levels of dds.obs.condition or replace the :attr:`.DESeqDataSet.design`, then run :func:`nbinomWaldTest` followed by :func:`lfcShrink`.

Three groups:

Three groups, compare condition effects:

Two groups, two individuals per group, compare within-individual condition effects:

Recommendations for single-cell analysis

The DESeq2 developers and collaborating groups have published recommendations for the best use of DESeq2 for single-cell datasets, which have been described first in [Berge2018]. Default values for DESeq2 were designed for bulk data and will not be appropriate for single-cell datasets. These settings and additional improvements have also been tested subsequently and published in [Zhu2018] and [AhlmannEltze2020].

• Use test="LRT" for significance testing when working with single-cell data, over the Wald test. This has been observed across multiple single-cell benchmarks.
• Set the following :func:`DESeq` arguments to these values: useT=TRUE, minmu=1e-6, and minReplicatesForReplace=Inf. The default setting of minmu was benchmarked on bulk RNA-seq and is not appropriate for single cell data when the expected count is often much less than 1.
• The default size factors are not optimal for single cell count matrices, instead consider setting :attr:`.DESeqDataSet.sizeFactors` from scran::computeSumFactors.
• One important concern for single-cell data analysis is the size of the datasets and associated processing time. To address the speed concerns, DESeq2 provides an interface to glmGamPoi, which implements faster dispersion and parameter estimation routines for single-cell data [AhlmannEltze2020]. To use this feature, set fitType = "glmGamPoi". Alternatively, one can use glmGamPoi as a standalone package. This provides the additional option to process data on-disk if the full dataset does not fit in memory, a quasi-likelihood framework for differential testing, and the ability to form pseudobulk samples (more details how to use glmGamPoi are in its README).

Optionally, one can consider using the zinbwave package to directly model the zero inflation of the counts, and take account of these in the DESeq2 model. This allows for the DESeq2 inference to apply to the part of the data which is not due to zero inflation. Not all single cell datasets exhibit zero inflation, and instead may just reflect low conditional estimated counts (conditional on cell type or cell state). There is example code for combining zinbwave and DESeq2 package functions in the zinbwave vignette. We also have an example of ZINB-WaVE + DESeq2 integration using the splatter package for simulation at the zinbwave-deseq2 GitHub repository.

Approach to count outliers

RNA-seq data sometimes contain isolated instances of very large counts that are apparently unrelated to the experimental or study design, and which may be considered outliers. There are many reasons why outliers can arise, including rare technical or experimental artifacts, read mapping problems in the case of genetically differing samples, and genuine, but rare biological events. In many cases, users appear primarily interested in genes that show a consistent behavior, and this is the reason why by default, genes that are affected by such outliers are set aside by DESeq2, or if there are sufficient samples, outlier counts are replaced for model fitting. These two behaviors are described below.

The :func:`DESeq` function calculates, for every gene and for every sample, a diagnostic test for outliers called Cook's distance. Cook's distance is a measure of how much a single sample is influencing the fitted coefficients for a gene, and a large value of Cook's distance is intended to indicate an outlier count. The Cook's distances are stored as a matrix available in dds.layers["cooks"].

The :meth:`.DESeqDataSet.results` function automatically flags genes which contain a Cook's distance above a cutoff for samples which have 3 or more replicates. The p-values and adjusted p-values for these genes are set to NA. At least 3 replicates are required for flagging, as it is difficult to judge which sample might be an outlier with only 2 replicates. This filtering can be turned off with dds.results(cooksCutoff=FALSE).

With many degrees of freedom -- i.e. many more samples than number of parameters to be estimated -- it is undesirable to remove entire genes from the analysis just because their data include a single count outlier. When there are 7 or more replicates for a given sample, the :func:`DESeq` function will automatically replace counts with large Cook's distance with the trimmed mean over all samples, scaled up by the size factor or normalization factor for that sample. This approach is conservative, it will not lead to false positives, as it replaces the outlier value with the value predicted by the null hypothesis. This outlier replacement only occurs when there are 7 or more replicates, and can be turned off with DESeq(dds, minReplicatesForReplace=Inf).

The default Cook's distance cutoff for the two behaviors described above depends on the sample size and number of parameters to be estimated. The default is to use the 99% quantile of the F(p, m − p) distribution (with p the number of parameters including the intercept and m the number of samples). The default for gene flagging can be modified using the cooksCutoff argument to the :meth:`.DESeqDataSet.results` function. For outlier replacement, :func:`DESeq` preserves the original counts in dds.counts() saving the replacement counts as a matrix named "replaceCounts" in dds.layers. Note that with continuous variables in the design, outlier detection and replacement is not automatically performed, as our current methods involve a robust estimation of within-group variance which does not extend easily to continuous covariates. However, users can examine the Cook's distances in dds.layers["cooks"], in order to perform manual visualization and filtering if necessary.

Dispersion plot and fitting alternatives

Plotting the dispersion estimates is a useful diagnostic. The dispersion plot below is typical, with the final estimates shrunk from the gene-wise estimates towards the fitted estimates. Some gene-wise estimates are flagged as outliers and not shrunk towards the fitted value, (this outlier detection is described in the manual page for :func:`estimateDispersionsMAP`). The amount of shrinkage can be more or less than seen here, depending on the sample size, the number of coefficients, the row mean and the variability of the gene-wise estimates.

Independent filtering of results

The :meth:`.DESeqDataSet.results` function of the DESeq2 package performs independent filtering by default using the mean of normalized counts as a filter statistic. A threshold on the filter statistic is found which optimizes the number of adjusted p-values lower than a significance level alpha (we use the standard variable name for significance level, though it is unrelated to the dispersion parameter α). The theory behind independent filtering is discussed in greater detail :ref:`below<indfilttheory>`. The adjusted p-values for the genes which do not pass the filter threshold are set to NA.

The default independent filtering is performed using the :func:`~.results.filtered_p` function of the genefilter package, and all of the arguments of :func:`~.results.filtered_p` can be passed to the :meth:`.DESeqDataSet.results` function. The filter threshold value and the number of rejections at each quantile of the filter statistic are available as metadata of the object returned by :meth:`.DESeqDataSet.results`.

For example, we can visualize the optimization by plotting the filterNumRej attribute of the results object. The :meth:`.DESeqDataSet.results` function maximizes the number of rejections (adjusted p-value less than a significance level), over the quantiles of a filter statistic (the mean of normalized counts). The threshold chosen (vertical line) is the lowest quantile of the filter for which the number of rejections is within 1 residual standard deviation to the peak of a curve fit to the number of rejections over the filter quantiles:

Independent filtering can be turned off by setting independentFiltering to FALSE.

Tests of log2 fold change above or below a threshold

It is also possible to provide thresholds for constructing Wald tests of significance. Two arguments to the :meth:`.DESeqDataSet.results` function allow for threshold-based Wald tests: lfcThreshold, which takes a numeric of a non-negative threshold value, and altHypothesis, which specifies the kind of test. Note that the alternative hypothesis is specified by the user, i.e. those genes which the user is interested in finding, and the test provides p-values for the null hypothesis, the complement of the set defined by the alternative. The altHypothesis argument can take one of the following four values, where β is the log2 fold change specified by the name argument, and x is the lfcThreshold.

• greaterAbs - |β| > x - tests are two-tailed
• lessAbs - |β| < x - p-values are the maximum of the upper and lower tests
• greater - β > x
• less - β <  − x

The four possible values of altHypothesis are demonstrated in the following code and visually by MA-plots in the following figures.

Access to all calculated values

All row-wise calculated values (intermediate dispersion calculations, coefficients, standard errors, etc.) are stored in the :class:`~DESeqDataSet.DESeqDataSet` object, e.g. dds in this vignette. These values are accessible by inspecting the var attribute of dds. Descriptions of the columns are accessible through the description attribute:

The mean values μ_ij = s_jq_ij and the Cook's distances for each gene and sample are stored as matrices in the :attr:`layers` attribute:

The dispersions α_i can be accessed with the :attr:`.DESeqDataSet.dispersions` attribute:

The size factors s_j are accessible via the :attr:`.DESeqDataSet.sizeFactors` attribute:

For advanced users, we also include a convenience function :func:`coef` for extracting the matrix [β_ir] for all genes i and model coefficients r. This function can also return a matrix of standard errors, see the documentation of :func:`coef`. The columns of this matrix correspond to the effects returned by :meth:`.DESeqDataSet.resultsNames`. Note that the :meth:`.DESeqDataSet.results` function is best for building results tables with p-values and adjusted p-values:

The beta prior variance σ²_r is stored as an attribute of the :class:`~DESeqDataSet.DESeqDataSet`:

General information about the prior used for log fold change shrinkage is also stored in a slot of the :class:`.DESeqResults` object. This would also contain information about what other packages were used for log2 fold change shrinkage:

The dispersion prior variance σ²_d is stored as an attribute of the dispersion function:

Sample-/gene-dependent normalization factors

In some experiments, there might be gene-dependent dependencies which vary across samples. For instance, GC-content bias or length bias might vary across samples coming from different labs or processed at different times. We use the terms normalization factors for a gene x sample matrix, and size factors for a single number per sample. Incorporating normalization factors, the mean parameter μ_ij becomes:

with normalization factor matrix NF having the same dimensions as the counts matrix K. This matrix can be incorporated as shown below. We recommend providing a matrix with gene-wise geometric means of 1, so that the mean of normalized counts for a gene is close to the mean of the unnormalized counts. This can be accomplished by dividing out the current gene geometric means:

>>> normFactors = normFactors / np.exp(np.mean(np.log(normFactors), axis=0))
>>> dds.normalizationFactors = normFactors

These steps then replace :meth:`.DESeqDataSet.estimateSizeFactors` which occurs within the :func:`DESeq` function. The :func:`DESeq` function will look for pre-existing normalization factors and use these in the place of size factors (and a message will be printed confirming this).

"Model matrix not full rank"

While most experimental designs run easily using design formula, some design formulas can cause problems and result in the :func:`DESeq` function returning an error with the text: "the model matrix is not full rank, so the model cannot be fit as specified." There are two main reasons for this problem: either one or more columns in the model matrix are linear combinations of other columns, or there are levels of factors or combinations of levels of multiple factors which are missing samples. We address these two problems below and discuss possible solutions.

>>> dds.results(contrast=("grpY.cndB", "grpX.cndB"))

The DESeq2 model

The DESeq2 model and all the steps taken in the software are described in detail in our publication [Love2014], and we include the formula and descriptions in this section as well. The differential expression analysis in DESeq2 uses a generalized linear model of the form:

where counts K_ij for gene i, sample j are modeled using a negative binomial distribution with fitted mean μ_ij and a gene-specific dispersion parameter α_i. 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 β_i give the log2 fold changes for gene i for each column of the model matrix X. Note that the model can be generalized to use sample- and gene-dependent normalization factors s_ij.

The dispersion parameter α_i defines the relationship between the variance of the observed count and its mean value. In other words, how far do we expected the observed count will be from the mean value, which depends both on the size factor s_j and the covariate-dependent part q_ij as defined above.

An option in DESeq2 is to provide maximum a posteriori estimates of the log2 fold changes in β_i after incorporating a zero-centered Normal prior (betaPrior). While previously, these moderated, or shrunken, estimates were generated by :func:`DESeq` or :func:`nbinomWaldTest` functions, they are now produced by the :func:`lfcShrink` function. Dispersions are estimated using expected mean values from the maximum likelihood estimate of log2 fold changes, and optimizing the Cox-Reid adjusted profile likelihood, as first implemented for RNA-seq data in edgeR [@CR,edgeR_GLM]. The steps performed by the :func:`DESeq` function are documented in its manual page; briefly, they are:

1. estimation of size factors s_j by :meth:`.DESeqDataSet.estimateSizeFactors`
2. estimation of dispersion α_i by :meth:`.DESeqDataSet.estimateDispersions`
3. negative binomial GLM fitting for β_i and Wald statistics by :func:`nbinomWaldTest`

For access to all the values calculated during these steps, see the section :ref:`above<access>`.

Changes compared to R DESeq2

This module is a Python port of the R package DESeq2. The port is based on version 1.39.3, and changes to the R package since this version have not been reflected in the present module. Also note that the port is partial: not all features may have been ported yet. To help us prioritize our future focus, please open a "feature request" issue on inmoose GitHub repository.

The present page mirrors DESeq2 vignette, and despite our efforts, some code examples may not accurately reflect the state of the Python module.

The main changes in this module compared to the original DESeq2 package are as follows:

• :class:`anndata.AnnData` is used as the superclass for storage of input data, intermediate calculations and results.

Methods changes since the 2014 DESeq2 paper

• In version 1.18 (November 2017), we add two :ref:`alternative shrinkage estimators<shrink>`, which can be used via :func:`lfcShrink`: an estimator using a t prior from the apeglm packages, and an estimator with a fitted mixture of normals prior from the ashr package.
• In version 1.16 (November 2016), the log2 fold change shrinkage is no longer default for the :func:`DESeq` and :func:`nbinomWaldTest` functions, by setting the defaults of these to betaPrior=FALSE, and by introducing a separate function :func:`lfcShrink`, which performs log2 fold change shrinkage for visualization and ranking of genes. While for the majority of bulk RNA-seq experiments, the LFC shrinkage did not affect statistical testing, DESeq2 has become used as an inference engine by a wider community, and certain sequencing datasets show better performance with the testing separated from the use of the LFC prior. Also, the separation of LFC shrinkage to a separate function :func:`lfcShrink` allows for easier methods development of alternative effect size estimators.
• A small change to the independent filtering routine: instead of taking the quantile of the filter (the mean of normalized counts) which directly maximizes the number of rejections, the threshold 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. This change was introduced in version 1.10 (October 2015).
• For the calculation of the beta prior variance, instead of matching the empirical quantile to the quantile of a Normal distribution, DESeq2 now uses the weighted quantile function of the Hmisc package. The weighting is described in the manual page for :func:`nbinomWaldTest`. The weights are the inverse of the expected variance of log counts (as used in the diagonals of the matrix W in the GLM). The effect of the change is that the estimated prior variance is robust against noisy estimates of log fold change from genes with very small counts. This change was introduced in version 1.6 (October 2014).

Count outlier detection

DESeq2 relies on the negative binomial distribution to make estimates and perform statistical inference on differences. While the negative binomial is versatile in having a mean and dispersion parameter, extreme counts in individual samples might not fit well to the negative binomial. For this reason, we perform automatic detection of count outliers. We use Cook's distance, which is a measure of how much the fitted coefficients would change if an individual sample were removed [Cook1977]. For more on the implementation of Cook's distance see the manual page for the :meth:`.DESeqDataSet.results` function. Below we plot the maximum value of Cook's distance for each row over the rank of the test statistic to justify its use as a filtering criterion.

Contrasts

Contrasts can be calculated for a :class:`~DESeqDataSet.DESeqDataSet` object for which the GLM coefficients have already been fit using the Wald test steps (:func:`DESeq` with test="Wald" or using :func:`nbinomWaldTest`). The vector of coefficients β is left multiplied by the contrast vector c to form the numerator of the test statistic. The denominator is formed by multiplying the covariance matrix Σ for the coefficients on either side by the contrast vector c. The square root of this product is an estimate of the standard error for the contrast. The contrast statistic is then compared to a Normal distribution as are the Wald statistics for the DESeq2 package.

Expanded model matrices

For the specific combination of :func:`lfcShrink` with the type normal and using contrast, DESeq2 uses expanded model matrices to produce shrunken log2 fold change estimates where the shrinkage is independent of the choice of reference level. In all other cases, DESeq2 uses standard model matrices, as produced by :func:`patsy.dmatrix`. The expanded model matrices differ from the standard model matrices, in that they have an indicator column (and therefore a coefficient) for each level of factors in the design formula in addition to an intercept. This is described in the DESeq2 paper. Using type normal with :func:`coef` uses standard model matrices, as does the apeglm shrinkage estimator.

Independent filtering and multiple testing