modelGeneVarByPoisson: Model the per-gene variance with Poisson noise
In scran: Methods for Single-Cell RNA-Seq Data Analysis

Description Usage Arguments Details Value Computing p-values Handling uninteresting factors Author(s) See Also Examples

Model the variance of the log-expression profiles for each gene, decomposing it into technical and biological components based on a mean-variance trend corresponding to Poisson noise.

modelGeneVarByPoisson(x, ...)

## S4 method for signature 'ANY'
modelGeneVarByPoisson(
  x,
  size.factors = NULL,
  block = NULL,
  design = NULL,
  subset.row = NULL,
  npts = 1000,
  dispersion = 0,
  pseudo.count = 1,
  ...,
  equiweight = TRUE,
  method = "fisher",
  BPPARAM = SerialParam()
)

## S4 method for signature 'SummarizedExperiment'
modelGeneVarByPoisson(x, ..., assay.type = "counts")

## S4 method for signature 'SingleCellExperiment'
modelGeneVarByPoisson(x, size.factors = sizeFactors(x), ...)

`x`	A numeric matrix of counts where rows are (usually endogenous) genes and columns are cells. Alternatively, a SummarizedExperiment or SingleCellExperiment containing such a matrix.
`...`	For the generic, further arguments to pass to each method. For the ANY method, further arguments to pass to `fitTrendVar`. For the SummarizedExperiment method, further arguments to pass to the ANY method. For the SingleCellExperiment method, further arguments to pass to the SummarizedExperiment method.
`size.factors`	A numeric vector of size factors for each cell in `x`, to be used for scaling gene expression.
`block`	A factor specifying the blocking levels for each cell in `x`. If specified, variance modelling is performed separately within each block and statistics are combined across blocks.
`design`	A numeric matrix containing blocking terms for uninteresting factors of variation.
`subset.row`	See `?"scran-gene-selection"`, specifying the rows for which to model the variance. Defaults to all genes in `x`.
`npts`	An integer scalar specifying the number of interpolation points to use.
`dispersion`	A numeric scalar specifying the dispersion for the NB distribution. If zero, a Poisson distribution is used.
`pseudo.count`	Numeric scalar specifying the pseudo-count to add prior to log-transformation.
`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 `block` is specified.
`method`	String specifying how p-values should be combined when `block` is specified, see `combinePValues`.
`BPPARAM`	A BiocParallelParam object indicating whether parallelization should be performed across genes.
`assay.type`	String or integer scalar specifying the assay containing the log-expression values.

For each gene, we compute the variance and mean of the log-expression values. A trend is fitted to the variance against the mean for simulated Poisson counts as described in fitTrendPoisson. The technical component for each gene is defined as the value of the trend at that gene's mean abundance. The biological component is then defined as the residual from the trend.

This function is similar to modelGeneVarWithSpikes, with the only difference being that the trend is fitted on simulated Poisson count-derived variances rather than spike-ins. The assumption is that the technical component is Poisson-distributed, or at least negative binomial-distributed with a known constant dispersion. This is useful for UMI count data sets that do not have spike-ins and are too heterogeneous to assume that most genes exhibit negligible biological variability.

If no size factors are supplied, they are automatically computed depending on the input type:

If size.factors=NULL for the ANY method, the sum of counts for each cell in x is used to compute a size factor via the librarySizeFactors function.
If size.factors=NULL for the SingleCellExperiment method, sizeFactors(x) is used if available. Otherwise, it defaults to library size-derived size factors.

If size.factors are supplied, they will override any size factors present in x.

A DataFrame is returned where each row corresponds to a gene in x (or in subset.row, if specified). This contains the numeric fields:

mean:: Mean normalized log-expression per gene.
total:: Variance of the normalized log-expression per gene.
bio:: Biological component of the variance.
tech:: Technical component of the variance.
p.value, FDR:: Raw and adjusted p-values for the test against the null hypothesis that bio<=0.

If block is not specified, the metadata of the DataFrame contains the output of running fitTrendVar on the simulated counts, along with the mean and var used to fit the trend.

If block is specified, the output contains another 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 metadata contains the output of fitTrendVar for the cells of that block.

The p-value for each gene is computed by assuming that the variance estimates are normally distributed around the trend, and that the standard deviation of the variance distribution is proportional to the value of the trend. This is used to construct a one-sided test for each gene based on its bio, under the null hypothesis that the biological component is equal to zero. The proportionality constant for the standard deviation is set to the std.dev returned by fitTrendVar. This is estimated from the spread of variance estimates for the simulated Poisson-distributed counts, so the null hypothesis effectively becomes “is this gene more variable than a hypothetical gene with only Poisson noise?”

Setting block will estimate the mean and variance of each gene for cells in each level of 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:

For means and variance components, this is done by averaging values across levels. If 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.
Per-level p-values are combined using combinePValues according to 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 equiweight depends on the chosen method.
Blocks with fewer than 2 cells are completely ignored and do not contribute to the combined mean, variance component or p-value.

Use of block is the recommended approach for accounting for any uninteresting categorical factor of variation. In addition to accounting for systematic differences in expression between levels of the blocking factor, it also accommodates differences in the mean-variance relationships.

Alternatively, uninteresting factors can be used to construct a design matrix to pass to the function via design. In this case, a linear model is fitted to the expression profile for each gene and the residual variance is calculated. This approach is useful for covariates or additive models that cannot be expressed as a one-way layout for use in block. However, it assumes that the error is normally distributed with equal variance for all observations of a given gene.

Use of block and design together is currently not supported and will lead to an error.

Aaron Lun

fitTrendVar, for the trend fitting options.

modelGeneVar, for modelling variance without spike-in controls.

library(scuttle)
sce <- mockSCE()

# Using spike-ins.
pois <- modelGeneVarByPoisson(sce)
pois

plot(pois$mean, pois$total, ylim=c(0, 10))
points(metadata(pois)$mean, metadata(pois)$var, col="red", pch=16)
curve(metadata(pois)$trend(x), add=TRUE, col="dodgerblue")

# With blocking.
block <- sample(LETTERS[1:2], ncol(sce), replace=TRUE)
blk <- modelGeneVarByPoisson(sce, 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, ylim=c(0, 10))
    points(metadata(current)$mean, metadata(current)$var, col="red", pch=16)
    curve(metadata(current)$trend(x), add=TRUE, col="dodgerblue")
}