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
|
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/fitTrendCV2.R
\name{fitTrendCV2}
\alias{fitTrendCV2}
\title{Fit a trend to the CV2}
\usage{
fitTrendCV2(
means,
cv2,
ncells,
min.mean = 0.1,
nls.args = list(),
simplified = TRUE,
nmads = 6,
max.iter = 50
)
}
\arguments{
\item{means}{A numeric vector containing mean normalized expression values for all genes.}
\item{cv2}{A numeric vector containing the squared coefficient of variation computed from normalized expression values for all genes.}
\item{ncells}{Integer scalar specifying the number of cells used to compute \code{cv2} and \code{means}.}
\item{min.mean}{Numeric scalar specifying the minimum mean to use for trend fitting.}
\item{nls.args}{A list of parameters to pass to \code{\link{nls}}.}
\item{simplified}{Logical scalar indicating whether the function can automatically use a simpler trend if errors are encountered for the usual paramterization.}
\item{nmads}{Numeric scalar specifying the number of MADs to use to compute the tricube bandwidth during robustification.}
\item{max.iter}{Integer scalar specifying the maximum number of robustness iterations to perform.}
}
\value{
A named list is returned containing:
\describe{
\item{\code{trend}:}{A function that returns the fitted value of the trend at any value of the mean.}
\item{\code{std.dev}:}{A numeric scalar containing the robust standard deviation of the ratio of \code{var} to the fitted value of the trend across all features used for trend fitting.}
}
}
\description{
Fit a mean-dependent trend to the squared coefficient of variation,
computed from count data after size factor normalization.
}
\details{
This function fits a mean-dependent trend to the CV2 of normalized expression values for the selected features.
Specifically, it fits a trend of the form
\deqn{y = A + \frac{B}{x}}{y = A + B/x}
using an iteratively reweighted least-squares approach implemented via \code{\link{nls}}.
This trend is based on a similar formulation from \pkg{DESeq2} and generally captures the mean-CV2 trend well.
Trend fitting is performed after weighting each observation according to the inverse of the density of observations at the same mean.
This avoids problems with differences in the distribution of means that would otherwise favor good fits in highly dense intervals at the expense of sparser intervals.
Low-abundance genes with means below \code{min.mean} are also removed prior to fitting, to avoid problems with discreteness and the upper bound on the CV2 at low counts.
Robustness iterations are also performed to protect against outliers.
An initial fit is performed and each observation is weighted using tricube-transformed standardized residuals (in addition to the existing inverse-density weights).
The bandwidth of the tricube scheme is defined as \code{nmads} multiplied by the median standardized residual.
Iterations are performed until convergence or \code{max.iters} is reached.
Occasionally, there are not enough high-abundance points to uniquely determine the \eqn{A} parameter.
In such cases, the function collapses back to fitting a simpler trend
\deqn{y = \frac{B}{x}}{y = B/x}
to avoid errors about singular gradients in \code{\link{nls}}.
If \code{simplified=FALSE}, this simplification is not allowed and the error is directly reported.
}
\examples{
library(scuttle)
sce <- mockSCE()
normcounts <- normalizeCounts(sce, log=FALSE)
# Fitting a trend:
library(DelayedMatrixStats)
means <- rowMeans(normcounts)
cv2 <- rowVars(normcounts)/means^2
fit <- fitTrendCV2(means, cv2, ncol(sce))
# Examining the trend fit:
plot(means, cv2, pch=16, cex=0.5,
xlab="Mean", ylab="CV2", log="xy")
curve(fit$trend(x), add=TRUE, col="dodgerblue", lwd=3)
}
\references{
Brennecke P, Anders S, Kim JK et al. (2013).
Accounting for technical noise in single-cell RNA-seq experiments.
\emph{Nat. Methods} 10:1093-95
}
\seealso{
\code{\link{modelGeneCV2}} and \code{\link{modelGeneCV2WithSpikes}}, where this function is used.
}
\author{
Aaron Lun
}
|
