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
|
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/fitTrendVar.R
\name{fitTrendVar}
\alias{fitTrendVar}
\title{Fit a trend to the variances of log-counts}
\usage{
fitTrendVar(
means,
vars,
min.mean = 0.1,
parametric = TRUE,
lowess = TRUE,
density.weights = TRUE,
nls.args = list(),
...
)
}
\arguments{
\item{means}{A numeric vector containing the mean log-expression value for each gene.}
\item{vars}{A numeric vector containing the variance of log-expression values for each gene.}
\item{min.mean}{A numeric scalar specifying the minimum mean to use for trend fitting.}
\item{parametric}{A logical scalar indicating whether a parametric fit should be attempted.}
\item{lowess}{A logical scalar indicating whether a LOWESS fit should be attempted.}
\item{density.weights}{A logical scalar indicating whether to use inverse density weights.}
\item{nls.args}{A list of parameters to pass to \code{\link{nls}} if \code{parametric=TRUE}.}
\item{...}{Further arguments to pass to \code{\link{weightedLowess}} for LOWESS fitting.}
}
\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 variances of the log-normalized expression values derived from count data.
}
\details{
This function fits a mean-dependent trend to the variance of the log-normalized expression for the selected features.
The fitted trend can then be used to decompose the variance of each gene into biological and technical components,
as done in \code{\link{modelGeneVar}} and \code{\link{modelGeneVarWithSpikes}}.
The default fitting process follows a two-step procedure when \code{parametric=TRUE} and \code{lowess=TRUE}:
\enumerate{
\item A non-linear curve of the form
\deqn{y = \frac{ax}{x^n + b}}{y = ax/(x^n + b)}
is fitted to the variances against the means using \code{\link{nls}}.
Starting values and the number of iterations are automatically set if not explicitly specified in \code{nls.args}.
\item \code{\link{weightedLowess}} is applied to the log-ratios of the variance of each gene to the corresponding fitted value from the non-linear curve.
The final trend is defined as the product of the fitted values from the non-linear curve and the exponential of the LOWESS fitted value.
If any tuning is necessary, the most important parameter here is \code{span}, which can be passed in the \code{...} arguments.
}
The aim is to use the parametric curve to reduce the sharpness of the expected mean-variance relationship for easier smoothing.
Conversely, the parametric form is not exact, so the smoothers will model any remaining trends in the residuals.
If \code{parametric=FALSE}, LOWESS smoothing is performed directly on the log-variances using \code{\link{weightedLowess}}.
This may be helpful in situations where the data does not follow the expected parametric curve,
e.g., for transformed data that spans negative values where the expression is not defined.
(Indeed, the expression above is purely empirical, chosen simply as it matched the shape of the observed trend in many datasets.)
If \code{lowess=FALSE}, the LOWESS smoothing step is omitted and the parametric fit is used directly.
This may be necessary in situations where the LOWESS overfits, e.g., because of very few points at high abundances.
}
\section{Filtering by mean}{
Genes with mean log-expression below \code{min.mean} are not used in trend fitting.
This aims to remove the majority of low-abundance genes and preserve the sensitivity of span-based smoothers at moderate-to-high abundances.
It also protects against discreteness, which can interfere with estimation of the variability of the variance estimates and accurate scaling of the trend.
Filtering is applied on the mean log-expression to avoid introducing spurious trends at the filter boundary.
The default threshold is chosen based on the point at which discreteness is observed in variance estimates from Poisson-distributed counts.
For heterogeneous droplet data, a lower threshold of 0.001-0.01 may be more appropriate, though this usually does not matter all too much.
When extrapolating to values below the smallest observed mean (or \code{min.mean}), the output function will approach zero as the mean approaches zero.
This reflects the fact that the variance should be zero at a log-expression of zero (assuming a pseudo-count of 1 was used).
When extrapolating to values above the largest observed mean, the output function will be set to the fitted value of the trend at the largest mean.
}
\section{Weighting by density}{
All fitting (with \code{\link{nls}} and \code{\link{weightedLowess}}) is performed by weighting each observation according to the inverse of the density of observations at the same mean.
This ensures that parts of the curve with few points are fitted as well as parts of the trend with many points.
Otherwise, differences in the distribution of means would favor good fits in highly dense intervals at the expense of sparser intervals.
(Note that these densities are computed after filtering on \code{min.mean}, so the high density of points at zero has no effect.)
That being said, the density weighting can give inappropriate weight to very sparse intervals, especially those at high abundance.
This results in overfitting where the trend is compelled to pass through each point at these intervals.
For most part, this is harmless as most high-abundance genes are not highly variable so an overfitted trend is actually appropriate.
However, if high-abundance genes are variable, it may be better to set \code{density.weights=FALSE} to avoid this overfitting effect.
}
\examples{
library(scuttle)
sce <- mockSCE()
sce <- logNormCounts(sce)
# Fitting a trend:
library(DelayedMatrixStats)
means <- rowMeans(logcounts(sce))
vars <- rowVars(logcounts(sce))
fit <- fitTrendVar(means, vars)
# Comparing the two trend fits:
plot(means, vars, pch=16, cex=0.5, xlab="Mean", ylab="Variance")
curve(fit$trend(x), add=TRUE, col="dodgerblue", lwd=3)
}
\seealso{
\code{\link{modelGeneVar}} and \code{\link{modelGeneVarWithSpikes}}, where this function is used.
}
\author{
Aaron Lun
}
|
