pseudoBulkDGE {scran} | R Documentation

## 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,
      sample,
      label,
      design,
      coef = ncol(design),
      contrast = NULL,
      condition = NULL,
      lfc = 0
    )

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

### Arguments

`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 `ncol(x)`, specifying the sample of origin 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 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 `samples`.
`coef` | Integer scalar or vector indicating the coefficients to drop from `design` to form the null hypothesis.
`contrast` | Numeric vector or matrix containing the contrast of interest. Takes precedence over `coef`.
`condition` | A vector or factor of length equal to `nrow(design)`, specifying the experimental condition for each sample (i.e., row of `design`). Only used for filtering.
`lfc` | Numeric scalar specifying the log-fold change threshold to use in `glmTreat`.
`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 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.

### Value

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.

### Comment 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` 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.

### Note about the design matrix

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.

### 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.

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

### Examples

    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
    )

[Package scran version 1.16.0 Index]