pseudoBulkDGE: Quickly perform pseudo-bulk DE analyses

In MarioniLab/scran: Methods for Single-Cell RNA-Seq Data Analysis

Quickly perform pseudo-bulk DE analyses

Description

A wrapper function around edgeR's quasi-likelihood methods to conveniently perform differential expression analyses on pseudo-bulk profiles, allowing detection of cell type-specific changes between conditions in replicated studies.

Usage

pseudoBulkDGE(x, ...)

## S4 method for signature 'ANY'
pseudoBulkDGE(
  x,
  col.data,
  label,
  design,
  coef,
  contrast = NULL,
  condition = NULL,
  lfc = 0,
  include.intermediates = TRUE,
  row.data = NULL,
  sorted = FALSE,
  method = c("edgeR", "voom"),
  qualities = TRUE,
  robust = TRUE,
  sample = NULL
)

## S4 method for signature 'SummarizedExperiment'
pseudoBulkDGE(x, col.data = colData(x), ..., assay.type = 1)

Arguments

x	A numeric matrix of counts where rows are genes and columns are pseudo-bulk profiles. Alternatively, a SummarizedExperiment object containing such a matrix in its assays.
...	For the generic, additional arguments to pass to individual methods. For the SummarizedExperiment method, additional arguments to pass to the ANY method.
col.data	A data.frame or DataFrame containing metadata for each column of x.
label	A vector of factor of length equal to ncol(x), specifying the cluster or cell type assignment for each column of x.
design	A formula to be used to construct a design matrix from variables in col.data. Alternatively, a function that accepts a data.frame with the same fields as col.data and returns a design matrix.
coef	String or character vector containing the coefficients to drop from the design matrix to form the null hypothesis. Can also be an integer scalar or vector specifying the indices of the relevant columns.
contrast	Numeric vector or matrix containing the contrast of interest. Alternatively, a character vector to be passed to makeContrasts to create this numeric vector/matrix. If specified, this takes precedence over coef.
condition	A vector or factor of length equal to ncol(x), specifying the experimental condition for each column of x. Only used for abundance-based filtering of genes.
lfc	Numeric scalar specifying the log-fold change threshold to use in glmTreat or treat.
include.intermediates	Logical scalar indicating whether the intermediate edgeR objects should be returned.
row.data	A DataFrame containing additional row metadata for each gene in x, to be included in each of the output DataFrames. This should have the same number and order of rows as x.
sorted	Logical scalar indicating whether the output tables should be sorted by p-value.
method	String specifying the DE analysis framework to use.
qualities	Logical scalar indicating whether quality weighting should be used when method="voom", see voomWithQualityWeights for more details.
robust	Logical scalar indicating whether robust empirical Bayes shrinkage should be performed.
sample	Deprecated.
assay.type	String or integer scalar specifying the assay to use from x.

Details

In replicated multi-condition scRNA-seq experiments, we often have clusters comprised of cells from different samples of different experimental conditions. It is often desirable to check for differential expression between conditions within each cluster, allowing us to identify cell-type-specific responses to the experimental perturbation.

Given a set of pseudo-bulk profiles (usually generated by sumCountsAcrossCells), this function loops over the labels and uses edgeR or voom to detect DE genes between conditions. The DE analysis for each label is largely the same as a standard analysis for bulk RNA-seq data, using design and coef or contrast as described in the edgeR or limma user guides. Generally speaking, edgeR handles low counts better via its count-based model but method="voom" supports variable sample precision when quality=TRUE.

Performing pseudo-bulk DGE enables us to reuse well-tested methods developed for bulk RNA-seq data analysis. Each pseudo-bulk profile can be treated as an in silico mimicry of a real bulk RNA-seq sample (though in practice, it tends to be much more variable due to the lower numbers of cells). This also models the relevant variability between experimental replicates (i.e., across samples) rather than that between cells in the same sample, without resorting to expensive mixed-effects models.

The DE analysis for each label is independent of that for any other label. This aims to minimize problems due to differences in abundance and variance between labels, at the cost of losing the ability to share information across labels.

In some cases, it will be impossible to perform a DE analysis for a label. The most obvious reason is if there are no residual degrees of freedom; other explanations include impossible contrasts or a failure to construct an appropriate design matrix (e.g., if a cell type only exists in one condition).

Note that we assume that x has already been filtered to remove unstable pseudo-bulk profiles generated from few cells.

Value

A List with one DataFrame of DE results per unique (non-failed) level of cluster. This contains columns from topTags if method="edgeR" or topTable if method="voom". All DataFrames have row names equal to rownames(x).

The metadata of the List contains failed, a character vector with the names of the labels for which the comparison could not be performed - see Details.

The metadata of the individual DataFrames contains design, the final design matrix for that label. If include.intermediates, the metadata will also contain y, the DGEList used for the analysis; and fit, the DGEGLM object after GLM fitting.

Comments on abundance filtering

For each label, abundance filtering is performed using filterByExpr prior to further analysis. Genes that are filtered out will still show up in the DataFrame for that label, but with all statistics set to NA. As this is done separately for each label, a different set of genes may be filtered out for each label, which is largely to be expected if there is any label-specific expression.

By default, the minimum group size for filterByExpr is determined using the design matrix. However, this may not be optimal if the design matrix contains additional terms (e.g., blocking factors) in which case it is not easy to determine the minimum size of the groups relevant to the comparison of interest. To overcome this, users can specify condition.field to specify the group to which each sample belongs, which is used by filterByExpr to obtain a more appropriate minimum group size.

Author(s)

Aaron Lun

References

Tung P-Y et al. (2017). Batch effects and the effective design of single-cell gene expression studies. Sci. Rep. 7, 39921

Lun ATL and Marioni JC (2017). Overcoming confounding plate effects in differential expression analyses of single-cell RNA-seq data. Biostatistics 18, 451-464

Crowell HL et al. (2019). On the discovery of population-specific state transitions from multi-sample multi-condition single-cell RNA sequencing data. biorXiv

See Also

sumCountsAcrossCells, to easily generate the pseudo-bulk count matrix.

decideTestsPerLabel, to generate a summary of the DE results across all labels.

pseudoBulkSpecific, to look for label-specific DE genes.

pbDS from the muscat package, which uses a similar approach.

Examples

set.seed(10000)
library(scuttle)
sce <- mockSCE(ncells=1000)
sce$samples <- gl(8, 125) # Pretending we have 8 samples.

# Making up some clusters.
sce <- logNormCounts(sce)
clusters <- kmeans(t(logcounts(sce)), centers=3)$cluster

# Creating a set of pseudo-bulk profiles:
info <- DataFrame(sample=sce$samples, cluster=clusters)
pseudo <- sumCountsAcrossCells(sce, info)

# Making up an experimental design for our 8 samples.
pseudo$DRUG <- gl(2,4)[pseudo$sample]

# DGE analysis:
out <- pseudoBulkDGE(pseudo, 
   label=pseudo$cluster,
   condition=pseudo$DRUG,
   design=~DRUG,
   coef="DRUG2"
)
out[[1]]
metadata(out[[1]])$design

MarioniLab/scran documentation built on Sept. 7, 2024, 6:25 a.m.