1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
|
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/pseudoBulkDGE.R
\name{pseudoBulkDGE}
\alias{pseudoBulkDGE}
\alias{pseudoBulkDGE,ANY-method}
\alias{pseudoBulkDGE,SummarizedExperiment-method}
\title{Quickly perform pseudo-bulk DE analyses}
\usage{
pseudoBulkDGE(x, ...)
\S4method{pseudoBulkDGE}{ANY}(
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
)
\S4method{pseudoBulkDGE}{SummarizedExperiment}(x, col.data = colData(x), ..., assay.type = 1)
}
\arguments{
\item{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.}
\item{...}{For the generic, additional arguments to pass to individual methods.
For the SummarizedExperiment method, additional arguments to pass to the ANY method.}
\item{col.data}{A data.frame or \linkS4class{DataFrame} containing metadata for each column of \code{x}.}
\item{label}{A vector of factor of length equal to \code{ncol(x)},
specifying the cluster or cell type assignment for each column of \code{x}.}
\item{design}{A formula to be used to construct a design matrix from variables in \code{col.data}.
Alternatively, a function that accepts a data.frame with the same fields as \code{col.data} and returns a design matrix.}
\item{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.}
\item{contrast}{Numeric vector or matrix containing the contrast of interest.
Alternatively, a character vector to be passed to \code{\link{makeContrasts}} to create this numeric vector/matrix.
If specified, this takes precedence over \code{coef}.}
\item{condition}{A vector or factor of length equal to \code{ncol(x)},
specifying the experimental condition for each column of \code{x}.
Only used for abundance-based filtering of genes.}
\item{lfc}{Numeric scalar specifying the log-fold change threshold to use in \code{\link{glmTreat}} or \code{\link{treat}}.}
\item{include.intermediates}{Logical scalar indicating whether the intermediate \pkg{edgeR} objects should be returned.}
\item{row.data}{A \linkS4class{DataFrame} containing additional row metadata for each gene in \code{x},
to be included in each of the output DataFrames.
This should have the same number and order of rows as \code{x}.}
\item{sorted}{Logical scalar indicating whether the output tables should be sorted by p-value.}
\item{method}{String specifying the DE analysis framework to use.}
\item{qualities}{Logical scalar indicating whether quality weighting should be used when \code{method="voom"},
see \code{\link{voomWithQualityWeights}} for more details.}
\item{robust}{Logical scalar indicating whether robust empirical Bayes shrinkage should be performed.}
\item{sample}{Deprecated.}
\item{assay.type}{String or integer scalar specifying the assay to use from \code{x}.}
}
\value{
A \linkS4class{List} with one \linkS4class{DataFrame} of DE results per unique (non-failed) level of \code{cluster}.
This contains columns from \code{\link{topTags}} if \code{method="edgeR"} or \code{\link{topTable}} if \code{method="voom"}.
All DataFrames have row names equal to \code{rownames(x)}.
The \code{\link{metadata}} of the List contains \code{failed},
a character vector with the names of the labels for which the comparison could not be performed - see Details.
The \code{\link{metadata}} of the individual DataFrames contains \code{design}, the final design matrix for that label.
If \code{include.intermediates}, the \code{\link{metadata}} will also contain
\code{y}, the DGEList used for the analysis; and \code{fit}, the DGEGLM object after GLM fitting.
}
\description{
A wrapper function around \pkg{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.
}
\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 \code{\link{sumCountsAcrossCells}}),
this function loops over the labels and uses \pkg{edgeR} or \code{\link{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 \code{design} and \code{coef} or \code{contrast} as described in the \pkg{edgeR} or \pkg{limma} user guides.
Generally speaking, \pkg{edgeR} handles low counts better via its count-based model
but \code{method="voom"} supports variable sample precision when \code{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 \emph{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 \code{x} has already been filtered to remove unstable pseudo-bulk profiles generated from few cells.
}
\section{Comments on abundance filtering}{
For each label, abundance filtering is performed using \code{\link{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 \code{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 \code{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 \code{condition.field} to specify the group to which each sample belongs,
which is used by \code{filterByExpr} to obtain a more appropriate minimum group size.
}
\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
}
\references{
Tung P-Y et al. (2017).
Batch effects and the effective design of single-cell gene expression studies.
\emph{Sci. Rep.} 7, 39921
Lun ATL and Marioni JC (2017).
Overcoming confounding plate effects in differential expression analyses of single-cell RNA-seq data.
\emph{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.
\emph{biorXiv}
}
\seealso{
\code{\link{sumCountsAcrossCells}}, to easily generate the pseudo-bulk count matrix.
\code{\link{decideTestsPerLabel}}, to generate a summary of the DE results across all labels.
\code{\link{pseudoBulkSpecific}}, to look for label-specific DE genes.
\code{pbDS} from the \pkg{muscat} package, which uses a similar approach.
}
\author{
Aaron Lun
}
|
