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
192
193
194
195
196
197
198
|
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/modelGeneCV2WithSpikes.R
\name{modelGeneCV2WithSpikes}
\alias{modelGeneCV2WithSpikes}
\alias{modelGeneCV2WithSpikes,ANY-method}
\alias{modelGeneCV2WithSpikes,SummarizedExperiment-method}
\alias{modelGeneCV2WithSpikes,SingleCellExperiment-method}
\title{Model the per-gene CV2 with spike-ins}
\usage{
modelGeneCV2WithSpikes(x, ...)
\S4method{modelGeneCV2WithSpikes}{ANY}(
x,
spikes,
size.factors = NULL,
spike.size.factors = NULL,
block = NULL,
subset.row = NULL,
...,
equiweight = TRUE,
method = "fisher",
BPPARAM = SerialParam()
)
\S4method{modelGeneCV2WithSpikes}{SummarizedExperiment}(x, ..., assay.type = "counts")
\S4method{modelGeneCV2WithSpikes}{SingleCellExperiment}(
x,
spikes,
size.factors = NULL,
spike.size.factors = NULL,
...,
assay.type = "counts"
)
}
\arguments{
\item{x}{A numeric matrix of counts where rows are (usually endogenous) genes and columns are cells.
Alternatively, a \linkS4class{SummarizedExperiment} or \linkS4class{SingleCellExperiment} containing such a matrix.}
\item{...}{For the generic, further arguments to pass to each method.
For the ANY method, further arguments to pass to \code{\link{fitTrendCV2}}.
For the \linkS4class{SummarizedExperiment} method, further arguments to pass to the ANY method.
For the \linkS4class{SingleCellExperiment} method, further arguments to pass to the SummarizedExperiment method.}
\item{spikes}{A numeric matrix of counts where each row corresponds to a spike-in transcript.
This should have the same number of columns as \code{x}.
Alternatively, for the SingleCellExperiment method,
this can be a string or an integer scalar specifying the \code{\link{altExp}} containing the spike-in count matrix.}
\item{size.factors}{A numeric vector of size factors for each cell in \code{x}, to be used for scaling gene expression.}
\item{spike.size.factors}{A numeric vector of size factors for each cell in \code{spikes}, to be used for scaling spike-ins.}
\item{block}{A factor specifying the blocking levels for each cell in \code{x}.
If specified, variance modelling is performed separately within each block and statistics are combined across blocks.}
\item{subset.row}{See \code{?"\link{scran-gene-selection}"}, specifying the rows for which to model the variance.
Defaults to all genes in \code{x}.}
\item{equiweight}{A logical scalar indicating whether statistics from each block should be given equal weight.
Otherwise, each block is weighted according to its number of cells.
Only used if \code{block} is specified.}
\item{method}{String specifying how p-values should be combined when \code{block} is specified, see \code{\link{combineParallelPValues}}.}
\item{BPPARAM}{A \linkS4class{BiocParallelParam} object indicating whether parallelization should be performed across genes.}
\item{assay.type}{String or integer scalar specifying the assay containing the counts.
For the SingleCellExperiment method, this is used to retrieve both the endogenous and spike-in counts.}
}
\value{
A \linkS4class{DataFrame} is returned where each row corresponds to a gene in \code{x} (or in \code{subset.row}, if specified).
This contains the numeric fields:
\describe{
\item{\code{mean}:}{Mean normalized expression per gene.}
\item{\code{total}:}{Squared coefficient of variation of the normalized expression per gene.}
\item{\code{trend}:}{Fitted value of the trend.}
\item{\code{ratio}:}{Ratio of \code{total} to \code{trend}.}
\item{\code{p.value, FDR}:}{Raw and adjusted p-values for the test against the null hypothesis that \code{ratio<=1}.}
}
If \code{block} is not specified,
the \code{metadata} of the DataFrame contains the output of running \code{\link{fitTrendCV2}} on the spike-in transcripts,
along with the \code{mean} and \code{cv2} used to fit the trend.
If \code{block} is specified, the output contains another \code{per.block} field.
This field is itself a DataFrame of DataFrames, where each internal DataFrame contains statistics for the variance modelling within each block and has the same format as described above.
Each internal DataFrame's \code{metadata} contains the output of \code{\link{fitTrendCV2}} for the cells of that block.
}
\description{
Model the squared coefficient of variation (CV2) of the normalized expression profiles for each gene,
using spike-ins to estimate the baseline level of technical noise at each abundance.
}
\details{
For each gene and spike-in transcript, we compute the variance and CV2 of the normalized expression values.
A trend is fitted to the CV2 against the mean for spike-in transcripts using \code{\link{fitTrendCV2}}.
The value of the trend at the abundance of each gene is used to define the variation attributable to technical noise.
The ratio to the trend is then used to define overdispersion corresponding to interesting biological heterogeneity.
This function is almost the same as \code{\link{modelGeneCV2}}, with the only theoretical difference being that the trend is fitted on spike-in CV2 rather than using the means and CV2 of endogenous genes.
This is because there are certain requirements for how normalization is performed when comparing spike-in transcripts with endogenous genes - see comments in \dQuote{Explanation for size factor rescaling}.
We enforce this by centering the size factors for both sets of features and recomputing normalized expression values.
}
\section{Computing p-values}{
The p-value for each gene is computed by assuming that the CV2 estimates are normally distributed around the trend, and that the standard deviation of the CV2 distribution is proportional to the value of the trend.
This is used to construct a one-sided test for each gene based on its \code{ratio}, under the null hypothesis that the ratio is equal to 1.
The proportionality constant for the standard deviation is set to the \code{std.dev} returned by \code{\link{fitTrendCV2}}.
This is estimated from the spread of CV2 values for spike-in transcripts, so the null hypothesis effectively becomes \dQuote{is this gene \emph{more} variable than spike-in transcripts of the same abundance?}
}
\section{Default size factor choices}{
If no size factors are supplied, they are automatically computed depending on the input type:
\itemize{
\item If \code{size.factors=NULL} for the ANY method, the sum of counts for each cell in \code{x} is used to compute a size factor via the \code{\link{librarySizeFactors}} function.
\item If \code{spike.size.factors=NULL} for the ANY method, the sum of counts for each cell in \code{spikes} is used to compute a size factor via the \code{\link{librarySizeFactors}} function.
\item If \code{size.factors=NULL} for the \linkS4class{SingleCellExperiment} method, \code{\link{sizeFactors}(x)} is used if available.
Otherwise, it defaults to library size-derived size factors.
\item If \code{spike.size.factors=NULL} for the \linkS4class{SingleCellExperiment} method and \code{spikes} is not a matrix, \code{\link{sizeFactors}(\link{altExp}(x, spikes)} is used if available.
Otherwise, it defaults to library size-derived size factors.
}
If \code{size.factors} or \code{spike.size.factors} are supplied, they will override any size factors present in \code{x}.
}
\section{Explanation for size factor rescaling}{
The use of a spike-in-derived trend makes several assumptions.
The first is that a constant amount of spike-in RNA was added to each cell, such that any differences in observed expression of the spike-in transcripts can be wholly attributed to technical noise.
The second is that endogenous genes and spike-in transcripts follow the same mean-variance relationship, i.e., a spike-in transcript captures the technical noise of an endogenous gene at the same mean count.
Here, the spike-in size factors across all cells are scaled so that their mean is equal to that of the gene-based size factors for the same set of cells.
This ensures that the average normalized abundances of the spike-in transcripts are comparable to those of the endogenous genes,
allowing the trend fitted to the former to be used to determine the biological component of the latter.
Otherwise, differences in scaling of the size factors would shift the normalized expression values of the former away from the latter, violating the second assumption.
If \code{block} is specified, rescaling is performed separately for all cells in each block.
This aims to avoid problems from (frequent) violations of the first assumption where there are differences in the quantity of spike-in RNA added to each batch.
Without scaling, these differences would lead to systematic shifts in the spike-in abundances from the endogenous genes when fitting a batch-specific trend
(even if there is no global difference in scaling across all batches).
}
\section{Handling uninteresting factors}{
Setting \code{block} will estimate the mean and variance of each gene for cells in each level of \code{block} separately.
The trend is fitted separately for each level, and the variance decomposition is also performed separately.
Per-level statistics are then combined to obtain a single value per gene:
\itemize{
\item For means and CV2 values, this is done by taking the geometric mean across blocking levels.
If \code{equiweight=FALSE}, a weighted average is used where the value for each level is weighted by the number of cells.
By default, all levels are equally weighted when combining statistics.
\item Per-level p-values are combined using \code{\link{combineParallelPValues}} according to \code{method}.
By default, Fisher's method is used to identify genes that are highly variable in any batch.
Whether or not this is responsive to \code{equiweight} depends on the chosen method.
\item Blocks with fewer than 2 cells are completely ignored and do not contribute to the combined mean, variance component or p-value.
}
}
\examples{
library(scuttle)
sce <- mockSCE()
# Using spike-ins.
spk <- modelGeneCV2WithSpikes(sce, "Spikes")
spk
plot(spk$mean, spk$total)
points(metadata(spk)$mean, metadata(spk)$var, col="red", pch=16)
curve(metadata(spk)$trend(x), add=TRUE, col="dodgerblue")
# With blocking (and spike-ins).
block <- sample(LETTERS[1:2], ncol(sce), replace=TRUE)
blk <- modelGeneCV2WithSpikes(sce, "Spikes", block=block)
blk
par(mfrow=c(1,2))
for (i in colnames(blk$per.block)) {
current <- blk$per.block[[i]]
plot(current$mean, current$total)
points(metadata(current)$mean, metadata(current)$var, col="red", pch=16)
curve(metadata(current)$trend(x), add=TRUE, col="dodgerblue")
}
}
\seealso{
\code{\link{fitTrendCV2}}, for the trend fitting options.
\code{\link{modelGeneCV2}}, for modelling variance without spike-in controls.
}
\author{
Aaron Lun
}
|
