Fit a mean-dependent trend to the gene-specific variances in single-cell RNA-seq data.

1 2 3 4 5 6 |

`x` |
A numeric matrix of normalized expression values, where each column corresponds to a cell and each row corresponds to a spike-in transcript. Alternatively, a SCESet object that contains such values. |

`trend` |
A string indicating whether the trend should be fully loess-based, or a mixture of loess and parametric fitting. |

`span, family, degree` |
Arguments to pass to |

`start` |
A named list of numeric scalars, containing starting values for parametric fitting with |

`design` |
A numeric matrix describing the uninteresting factors contributing to expression in each cell. |

`subset.row` |
A logical, integer or character scalar indicating the rows of |

`...` |
Additional arguments to pass to |

`assay` |
A string specifying which assay values to use, e.g., |

`use.spikes` |
A logical scalar specifying whether the trend should be fitted to variances for spike-in transcripts or endogenous genes. |

The strategy is to fit an abundance-dependent trend to the variance of the log-normalized expression for the spike-in transcripts, using `trendVar`

.
For SCESet objects, these expression values can be computed by `normalize`

after setting the size factors, e.g., with `computeSpikeFactors`

.
Log-transformed values are used as these are more robust to genes/transcripts with strong expression in only one or two outlier cells.

The mean and variance of the normalized log-counts is calculated for each spike-in transcript, and a trend is fitted to the variance against the mean for all transcripts.
The fitted value of this trend represents technical variability due to sequencing, drop-outs during capture, etc.
Variance decomposition to biological and technical components for endogenous genes can then be performed later with `decomposeVar`

.

The design matrix can be set if there are factors that should be blocked, e.g., batch effects, known (and uninteresting) clusters. Otherwise, it will default to an all-ones matrix, effectively treating all cells as part of the same group.

A named list is returned, containing:

`mean`

:A numeric vector of mean log-CPMs for all spike-in transcripts.

`var`

:A numeric vector of the variances of log-CPMs for all spike-in transcripts.

`trend`

:A function that returns the fitted value of the trend at any mean log-CPM.

`design`

:A numeric matrix, containing the design matrix that was used.

By default, a robust loess curve is used for trend fitting via `loess`

.
This protects against genes with very large or very small variances.
Some experimentation with `span`

, `degree`

or `family`

may be required to obtain satisfactory results.
The fit is also dependent on the quality of the spike-ins – the fit will obviously be poor if the coverage of all spike-ins is low.

Alternatively, when `trend="semiloess"`

, a non-linear curve of the form

*y = ax/(x^n + b)*

is fitted to the variances against the means using `nls`

, and a loess curve is then fitted to the log-ratios of the observed to fitted values.
The parametric curve reduces the sharpness of the trend for easier loess fitting.
Conversely, the parametric form is not exact, so the loess curve models any remaining trends in the residuals.

In general, the semi-loess setting tends to give smoother curves than loess alone.
It is more robust to the uneven distribution of spike-in transcripts across the covariate range.
However, it tends to be susceptible to convergence issues, and may require some fiddling with the start values to converge properly.
Reasonably good starting values for `a`

, `n`

and `b`

are chosen automatically but these can be set manually with `start`

.

Spike-in transcripts can be selected in `trendVar,SCESet-method`

using the `use.spikes`

method.
By default, `use.spikes=TRUE`

which means that only rows labelled as spike-ins with `isSpike(x)`

will be used.

When spike-ins are not available, `trendVar`

can also be applied directly to the counts for endogenous genes by setting `use.spikes=FALSE`

(or by manually supplying a matrix of normalized expression for endogenous genes, for `trendVar,matrix-method`

).
This assumes that most genes exhibit technical variation and little biological variation, e.g., in a homogeneous population.

If `use.spikes=NA`

, every row will be used for trend fitting, regardless of whether it corresponds to a spike-in transcript or endogenous gene.
Users can also directly specify which rows to use with `subset.row`

.
This will override any setting of `use.spikes`

.

If `assay="exprs"`

, `trendVar,SCESet-method`

will attempt to determine if the expression values were computed from counts via `normalize`

.
If so, a warning will be issued if the size factors are not centred at unity.
This is because different size factors are typically used for endogenous genes and spike-in transcripts.
If these size factor sets are not centred at the same value, there will be systematic differences in abundance between these features.
This precludes the use of a spike-in fitted trend with abundances for endogenous genes in `decomposeVar`

.

For other expression values and in `trendVar,matrix-method`

, the onus is on the user to ensure that normalization preserves differences in abundance.
In other words, the scaling factors used to normalize each feature should have the same mean.
This ensures that spurious differences in abundance are not introduced by the normalization process.

Aaron Lun

Lun ATL, McCarthy DJ and Marioni JC (2016).
A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor.
*F1000Res.* 5:2122

`nls`

,
`loess`

,
`decomposeVar`

,
`computeSpikeFactors`

,
`computeSumFactors`

,
`normalize`

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | ```
example(computeSpikeFactors) # Using the mocked-up data 'y' from this example.
# Normalizing (gene-based factors for genes, spike-in factors for spike-ins)
y <- computeSumFactors(y)
y <- computeSpikeFactors(y, general.use=FALSE)
y <- normalize(y)
# Fitting a trend to the spike-ins.
fit <- trendVar(y)
plot(fit$mean, fit$var)
curve(fit$trend(x), col="red", lwd=2, add=TRUE)
# Fitting a trend to the endogenous genes.
fit.g <- trendVar(y, use.spikes=FALSE)
plot(fit.g$mean, fit.g$var)
curve(fit.g$trend(x), col="red", lwd=2, add=TRUE)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.