pseudoBulkDGE {scran} | R Documentation |
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.
pseudoBulkDGE(x, ...) ## S4 method for signature 'ANY' pseudoBulkDGE( x, sample, label, design, coef = ncol(design), contrast = NULL, condition = NULL, lfc = 0 ) ## S4 method for signature 'SummarizedExperiment' pseudoBulkDGE(x, ..., assay.type = 1)
x |
For the ANY method, a numeric matrix of counts where rows are genes and columns are pseudo-bulk profiles. For the SummarizedExperiment method, 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. |
sample |
A vector or factor of length equal to |
label |
A vector of factor of length equal to |
design |
A numeric matrix containing the experimental design for the multi-sample comparison.
The number of rows should be equal to the total number of samples and the row names should be unique levels of |
coef |
Integer scalar or vector indicating the coefficients to drop from |
contrast |
Numeric vector or matrix containing the contrast of interest.
Takes precedence over |
condition |
A vector or factor of length equal to |
lfc |
Numeric scalar specifying the log-fold change threshold to use in |
assay.type |
String or integer scalar specifying the assay to use from |
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 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 user's guide.
The analysis for each label is independent, i.e., the GLM fitting is performed separately,
to minimize problems due to differences in abundance and variance between labels.
Performing pseudo-bulk DGE enables us to re-use 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 relatively higher numbers of cells). Pseudo-bulk analysis also allows direct modelling of variability between experimental replicates (i.e., across samples) rather than that between cells in the same sample. The former is more relevant to any statistical analysis that aims to obtain reproducible results.
In some cases, it will be impossible to perform an edgeR analysis as there are no residual degrees of freedom.
This will be represented by a DataFrame with log-fold changes but NA
p-values and FDRs.
In other cases, all statistics in the DataFrame will be NA
if the contrast cannot be performed,
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.
A List with one DataFrame of DE results per unique level of cluster
.
This will contain at least the fields "LogCPM"
, "PValue"
and "FDR"
,
and usually "logFC"
depending on whether an ANOVA-like contrast is requested in coef
or contrast
.
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,
e.g., due to lack of residual d.f.
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
to specify the group to which each sample belongs,
which is passed to filterByExpr
via its group
argument to obtain a more appropriate minimum group size.
We require the user to supply a design matrix and contrast for safety's sake. Technically, we could accept a formula and dynamically construct the design matrix for each label, which would be more convenient by avoiding the need for up-front construction. We do not do this as it can silently give incorrect results if the dimensions of the design matrix change between labels, especially if different labels have different numbers of pseudo-bulk profiles after QC filtering. The most obvious example is when the number of blocking levels change between labels, altering the identity of the column corresponding to the contrast of interest.
Aaron Lun
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
sumCountsAcrossCells
, to easily generate the pseudo-bulk count matrix.
decideTestsPerLabel
, to generate a summary of the DE results across all labels.
pbDS
from the muscat package, which uses a similar approach.
set.seed(10000) library(scater) 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) # Determining the experimental design for our 8 samples. DRUG <- gl(2,4) design <- model.matrix(~DRUG) rownames(design) <- seq_len(8) # DGE analysis: out <- pseudoBulkDGE(pseudo, sample=pseudo$sample, label=pseudo$cluster, design=design )