Precision weights for counts

Precision weights for counts#

Precision weights for log counts from their mean-variance trend, as in limma’s voom.

For counts r_gi of feature g in sample i and effective library sizes R_i (library sizes times normalization factors), voom (Law, Chen, Shi and Smyth 2014) works with the log2 counts per million

y_gi = log2((r_gi + 0.5) / (R_i + 1) * 1e6).

The design is fitted to every feature by ordinary least squares, which gives fitted values mu_gi and residual standard deviations s_g. The square roots s_g^(1/2) are smoothed by LOWESS against the average log2 count of each feature,

x_g = mean_i(y_gi) + mean_i(log2(R_i + 1)) - log2(1e6),

and the smoothed curve lo() is interpolated linearly between its points and held constant beyond them. Each observation is then weighted by the inverse of its predicted variance, w_gi = lo(lambda_gi)^(-4), where

lambda_gi = log2(2^mu_gi (R_i + 1) / 1e6)

is its fitted log2 count. A weighted least squares fit of the log2 counts per million with these weights, followed by empirical Bayes moderation, is limma-voom.

With prior weights p_gi, for samples, observations or features, the design is fitted by weighted least squares with them, which gives the residual standard deviations and fitted values of the trend, and the weights are w_gi p_gi. limma’s voom returns w_gi alone, which its help page says modify the prior weights in a multiplicative fashion: limma’s weights times the prior weights are CytePool’s, the weights a fit needs.

voom_with_quality_weights() combines the weights with sample quality weights (Liu et al. 2015), as limma’s voomWithQualityWeights does (found by probing): voom; sample quality weights a from its log2 counts per million and weights, by sample_quality_weights(); voom again, with a as prior weights; sample quality weights b from the new trend’s weights w’; and the weights w’_gi b_i.

Results match limma 3.62.1, called with adaptive.span = TRUE where the span is chosen. Conventions that the paper leaves open were settled from limma’s documentation (its help page and NEWS) or by running limma as a black box; comments mark the latter “found by probing”:

  • Features with no counts in any sample are left out of the trend (found by probing); their weights are read from it at their fitted log2 counts, as for any other feature.

  • The LOWESS is R’s lowess with its defaults other than the span: three robustness iterations, and interpolation within 1% of the range of x (found by probing).

  • By default the span is chosen from the number of features n by limma’s rule, choose_lowess_span(): 1 up to 50 features, and 0.3 + 0.7 (50 / n)^(1/3) beyond. limma’s voom chooses it so with adaptive.span = TRUE, an option since limma 3.60 and the default since 3.66 (limma’s NEWS); its default span was 0.5 before. n counts every feature, also those without counts that the trend leaves out, as limma’s nrow(counts) does (its help page; confirmed by probing). A span given as a number is used as it is, as with adaptive.span = FALSE.

  • The interpolation averages the curve over tied values of x (limma’s NEWS); LOWESS gives tied points the same value, so each distinct x is one point of the curve. The x are computed as R computes them, to the last bit, so that their ties are limma’s: LOWESS treats tied points differently from points a rounding error apart.

  • x averages every log2 count per million of a feature, those of zero prior weight included, as limma’s Amean does (found by probing).

  • A design that leaves no residual degrees of freedom gives every observation weight one, times its prior weight, with a warning (limma’s NEWS).

  • There must be at least two features, as limma requires, and at least two with counts (and residual degrees of freedom under the prior weights) to fit a trend.

Deliberate deviations, where limma returns values without meaning or fails:

  • A sample whose library size is zero is an error. limma gives it the log2 counts per million of a count of 0.5 in a library of one read.

  • A trend that is not positive where weights are read from it is an error. limma takes its inverse fourth power anyway, which is positive for a negative trend.

  • When every feature in the trend has the same x, the curve is constant; limma fails.

  • A design that is not of full column rank is rejected by fit_linear_model(); limma drops the aliased columns.

  • span must lie in (0, 1], as limma documents; limma accepts larger values.

  • A feature whose prior weights leave no residual degrees of freedom is left out of the trend. limma passes its missing standard deviation to R’s lowess, which then returns a curve without meaning for every feature.

  • Where zero prior weights leave coefficients of a feature that its other samples cannot estimate, its fitted values are those of the least squares fit over the samples of positive weight, and its observations of zero prior weight get weight zero. limma returns missing weights for every observation of such a feature.

References

Law CW, Chen Y, Shi W, Smyth GK (2014). voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology 15, R29. doi:10.1186/gb-2014-15-2-r29

Liu R, Holik AZ, Su S, Jansz N, Chen K, Leong HS, Blewitt ME, Asselin-Labat ML, Smyth GK, Ritchie ME (2015). Why weight? Modelling sample and observational level variability improves power in RNA-seq analyses. Nucleic Acids Research 43, e97. doi:10.1093/nar/gkv412

class cytepool.stats.voom.VoomResult(log_cpm, weights, effective_library_sizes, mean_log_count, sqrt_sd, sqrt_sd_trend, span, sample_weights=None)#

Log2 counts per million and their precision weights, from voom().

Fit a linear model with fit_linear_model(result.log_cpm, design, weights=result.weights), or, for the weights of voom_with_dream_weights(), a mixed model with fit_mixed_model. Arrays with a sample axis have one row per sample; the others have one value per feature, in the order of the columns of the counts. The arrays are read-only views; the constructor does not copy them.

Variables:

  • log_cpm – (n_samples, n_features) log2 counts per million, log2((r + 0.5) / (R + 1) * 1e6) with R the effective library size.

  • weights – (n_samples, n_features) precision weights: the inverse of each observation’s variance predicted by the trend, times the prior weights given to voom() (each feature’s divided by their mean, for voom_with_dream_weights) or the sample_weights.

  • effective_library_sizes – (n_samples,) library sizes times normalization factors. Pass them back as library_sizes without normalization factors to repeat the computation.

  • mean_log_count – (n_features,) average log2 count of each feature, the x of the mean-variance trend. NaN for features with no counts, which are left out of the trend.

  • sqrt_sd – (n_features,) square root of each feature’s residual standard deviation, the y of the trend. NaN for the features left out of the trend: those with no counts, those whose prior weights leave no residual degrees of freedom, and every feature when the design leaves none.

  • sqrt_sd_trend – (n_features,) the LOWESS trend of sqrt_sd at mean_log_count; NaN where sqrt_sd is. Sort by mean_log_count to draw it as a curve.

  • span – The span of the LOWESS trend, the fraction of the features in the trend that each local fit uses, as chosen or given. NaN when the design leaves no residual degrees of freedom, so that no trend is fitted.

  • sample_weights – (n_samples,) sample quality weights that weights include, from voom_with_quality_weights(); None otherwise.

Raises:

  • TypeError – If an array is a numpy masked array, or span is not a real number.

  • ValueError – If the shapes of the arrays disagree, span is neither in (0, 1] nor NaN, or a sample weight is not positive and finite.

Explained in: Precision weights

cytepool.stats.voom.voom(counts, design, *, library_sizes=None, norm_factors=None, span=None, weights=None)#

Log2 counts per million with precision weights from the mean-variance trend.

limma’s voom (Law et al. 2014); the module docstring gives the method. The design is fitted by fit_linear_model(), which checks it. Filter out features with very low counts first, for example with filter_by_expression(), as limma advises.

Continue with fit_linear_model(result.log_cpm, design, weights=result.weights) and moderated_statistics(). With observation weights, limma’s contrasts.fit approximates the standard errors of contrasts of correlated coefficients; contrasts_fit(method="limma") reproduces it, and the default computes them exactly.

Parameters:

  • counts – (n_samples, n_features) counts: finite and non-negative, integer or not (estimated counts are fine). Do not pass values already scaled by library size, such as CPM or TPM: voom needs the library sizes.

  • design – (n_samples, n_coefficients) design matrix of full column rank, or anything else that fit_linear_model accepts, such as a Design.

  • library_sizes – Library size of each sample; defaults to the total counts. Must be positive.

  • norm_factors – Normalization factor of each sample, such as TMM factors from tmm_factors(), multiplied into the library sizes; defaults to ones. Must be positive.

  • span – Fraction of the features in the trend (those with counts) that each local fit of the LOWESS uses, in (0, 1]. By default it is chosen from the number of features, those without counts included, by choose_lowess_span(): 1 up to 50 features, 0.65 at 400 and 0.395 at 20,000, as limma’s voom chooses it by default since limma 3.66 (adaptive.span = TRUE). A number fixes it, as adaptive.span = FALSE, span = <number> does in limma, whose default before 3.66 was 0.5.

  • weights – Prior weights, as limma’s voom(weights = ...) takes them: (n_samples,) sample weights, such as sample quality weights, (n_samples, n_features) observation weights, or (1, n_features) feature weights. Finite and non-negative; a zero weight leaves an observation out of the fit of the trend and gives it weight zero. The returned weights include them. A vector is always sample weights, also when there are as many features as samples, where limma takes it for feature weights.

Returns:

The log2 counts per million, their weights (finite and non-negative; positive where the prior weights are), the effective library sizes, and the mean-variance trend with its span.

Warns:

RuntimeWarning – If the design, without the observations of zero prior weight, leaves no residual degrees of freedom, so that no trend can be fitted; every weight is then its prior weight, or one.

Raises:

  • TypeError – If an argument has the wrong type, or is or holds a numpy masked array.

  • ValueError – If counts is not a 2-D array of finite, non-negative values with at least two features, the design does not have one row per sample or is not of full column rank, a library size or normalization factor is not positive, span is not in (0, 1], the weights have another shape or a negative or non-finite value, fewer than two features have counts (with residual degrees of freedom), or the trend is not positive where the weights are read from it.

Used in: Quick start

Explained in: Choosing a method, Precision weights

cytepool.stats.voom.voom_with_quality_weights(counts, design, *, library_sizes=None, norm_factors=None, span=None, variance_design=None, variance_groups=None, method='gene_by_gene', prior_n=10.0, max_iterations=50, tolerance=1e-5)#

voom’s precision weights combined with sample quality weights.

limma’s voomWithQualityWeights (Liu et al. 2015): the module docstring gives the steps, and sample_quality_weights() the sample weights. Samples that follow the model less closely, such as pseudobulk samples pooled from few cells or from an outlying donor, get lower weights instead of being removed. Continue as after voom(), with fit_linear_model(result.log_cpm, design, weights=result.weights).

Parameters:

  • counts – (n_samples, n_features) counts, as for voom().

  • design – The design, as for voom().

  • library_sizes – Library size of each sample, as for voom().

  • norm_factors – Normalization factor of each sample, as for voom().

  • span – The span of the LOWESS trend, as for voom(): by default chosen from the number of features, as limma’s voomWithQualityWeights does with adaptive.span = TRUE.

  • variance_design – Design of the variance model, as for sample_quality_weights; defaults to a weight for each sample.

  • variance_groups – One label per sample, for a weight per group of samples; give at most one of the two.

  • method – "gene_by_gene", limma’s default here, or "reml".

  • prior_n – Prior number of genes at equal weights, as for sample_quality_weights; limma’s voomWithQualityWeights uses 10.

  • max_iterations – The most REML iterations for each estimate of the weights.

  • tolerance – Convergence tolerance of REML.

Returns:

The log2 counts per million, their weights (the second trend’s times the sample weights), the effective library sizes, the second trend, its span, and the sample_weights.

Warns:

RuntimeWarning – If the design leaves no residual degrees of freedom (every weight is one), fewer than two (every sample weight is one), or the sample weights warn as sample_quality_weights does.

Raises:

  • TypeError – If an argument has the wrong type, or is or holds a numpy masked array.

  • ValueError – If an argument is invalid as for voom() or sample_quality_weights, or the trend is not positive where weights are read from it.

Explained in: Precision weights