Gene-set tests#

Gene-set tests of one contrast of a linear model: limma’s camera and fry.

A gene set, such as a pathway, is tested as a unit for one contrast of the linear model fitted to every gene. camera() is competitive: it asks whether the genes of a set are more differentially expressed than the other genes, allowing for the correlation between the genes of the set (Wu and Smyth 2012). fry() is self-contained: it asks whether the genes of a set are differentially expressed at all, as the rotation test of Wu et al. (2010) with the mean statistic does, without drawing rotations (Giner and Smyth 2016). Results match limma 3.62.1 except where noted under “Deliberate differences”. The rotation tests, which draw rotations (limma’s roast, mroast and romer), are in cytepool.stats.rotation_gene_set_tests.

Effects:

Both tests work with the effects of each gene’s least squares fit (Wu et al. 2010, section 4.1). The design X, n samples by p coefficients, is reparametrised so that the contrast is its last coefficient: a contrast given as a column of the design moves that column last; a contrast given as weights c becomes the last column, X c / (c’c), after the columns X h_2, …, X h_p, where h_1, …, h_p are the columns of the Householder reflection H of c, so that the last coefficient is c’alpha (limma’s contrastAsCoef, found by probing). With the Householder QR decomposition (Golub and Van Loan 2013) of this design in R’s sign convention, a gene’s data y give the effects Q’y; the p-th is the contrast effect u_0, which divided by the residual standard deviation is the gene’s t-statistic, and the last d = n - p are its residual effects u_1, …, u_d, independent with variance sigma^2 under the model. With precision weights W each gene’s data and design are whitened by W^(1/2) first. Every sample must be observed, with a positive weight, in every gene, so that all genes have the same residual degrees of freedom d.

The residual effects are coordinates in a basis of the residual space, which the QR decomposition chooses. With observation weights each gene has a basis of its own, and the inter-gene correlation of camera and the statistics of fry depend on them; fry’s standardizations other than “none” depend on the basis even without weights, through the largest effect. The bases are limma’s: they depend on how the contrast is given, and a pivot of the decomposition that is zero up to rounding (below 1e-12 of its column), as in indicator designs with the samples sorted by group, counts as zero, which R’s qr() gives the positive sign. limma computes most such zeros exactly, but leaves rounding noise in a few percent of genes, a different few in camera and in fry, so that their bases, and limma’s results, change with the last bits of the data there; CytePool’s do not (found by probing).

camera:

Genes get moderated t-statistics, the contrast effect over the posterior standard deviation of squeeze_variances() on min(d + d0, G d) degrees of freedom (as eBayes caps them, for G genes), turned into z-scores by Hill’s (1970) approximation (Brophy 1987, equation 5), as limma’s zscoreT(approx = TRUE, method = "hill") does (found by probing). With the mean z-score of a set of m genes and the mean and variance of all G z-scores, the test statistic is the two-sample t-statistic of the set against the other G - m genes with the variance of the set’s mean inflated by VIF = 1 + (m - 1) rho (Wu and Smyth 2012). The inter-gene correlation rho is a value given for every set, or estimated for each set from the residual effects: each gene’s are divided by its residual standard deviation (its variance bounded below by 1e-8, found by probing), the m x d matrix is averaged over genes, and m/d times the sum of squares of that mean estimates the VIF. Unless negative correlations are allowed, rho enters the VIF as max(rho, 0); a given negative correlation that makes the VIF of a set negative, which no correlation matrix allows, gives the set a NaN p-value. The p-value is two-sided, from the t distribution on G - 2 degrees of freedom for a given correlation and on min(d, G - 2) for an estimated one (found by probing). The rank version replaces the t-statistic by the Wilcoxon-Mann-Whitney rank sum of the set’s z-scores among all genes, with the variance of the rank sum when the set’s values share a correlation rho (Wu and Smyth 2012, from Barry et al. 2008, equation 4.5), times Zar’s (1999, section 8.10) correction for ties, 1 - sum(t^3 - t) / (G^3 - G), and a continuity correction of 1/2; it is referred to the normal distribution for a given correlation and to the t distribution on min(d, G - 2) degrees of freedom for an estimated one (found by probing).

fry:

The effects of each gene are divided by a scale (see standardize), and the set statistic is the mean b over the set’s genes of these (d + 1)-vectors. Under the null hypothesis the d + 1 effects of a gene are independent with equal variances, so one rotation applied to the effects of every gene leaves the joint distribution of the data unchanged (Wu et al. 2010, section 3.4). Under a uniformly random rotation b turns uniformly over its sphere, so T = b_0 / (sum_k b_k^2 / d)^(1/2) has the t distribution on d degrees of freedom: the p-value the rotation test of the mean statistic would give with infinitely many rotations, exactly so when every gene has the same scale. The mixed (non-directional) statistic is the sum over genes of the squared standardized contrast effects. Under rotation it is r’(E’E)r for the set’s m x (d + 1) matrix E of standardized effects and r uniform on the unit sphere, so it equals sum_k lambda_k v_k with lambda_k the d + 1 eigenvalues of E’E and (v_k) Dirichlet(1/2, …, 1/2); scaled to [0, 1] by the smallest and largest eigenvalue, it is approximated by the beta distribution with its first two moments (found by probing to be limma’s approximation, for sets of more than d genes). With one gene the approximation is exact: the mixed p-value is the two-sided p-value of the gene’s t-statistic. For more genes it can be much too small in the tail, as limma’s is: 6.3e-4 where a simulation of the rotation test gives 3.2e-3.

The scales of standardize="residual_sd" and "posterior_sd" start from the mean square of a gene’s d + 1 effects without the largest, which leaves out the contrast effect of a strongly differentially expressed gene (found by probing). "posterior_sd" divides it by its expected value under the null hypothesis, which makes it unbiased for sigma^2, and moderates it with the prior of the residual variances as a variance on 0.92 d degrees of freedom, fewer than d since the largest effect is left out (limma’s rule, found by probing).

Deliberate differences from limma:

fry’s mixed p-value of a set of 2 to d genes. limma uses only the set’s m nonzero eigenvalues and treats them as an m-dimensional problem (found by probing), which can be far from the rotation p-value; here all d + 1 dimensions are used, as for larger sets. validation/gene_set_tests.py compares both with the rotation p-value by simulation.
fry’s standardize="posterior_sd". limma’s expected value of the trimmed mean square is about 7.3e-5 (d + 1) / d above the exact one used here (found by probing), which moves p-values by about 1e-4 of their value.
camera’s rank test with continuity correction can give both tails above 1/2; limma then reports twice the smaller, a p-value above 1, and here it is 1.
A gene listed twice in a set counts once; limma counts it twice.
A set without genes, or containing every gene in camera, has NaN p-values and no direction, and is left out of the adjusted p-values of the other sets. limma reports the direction “Up”, in camera’s rank test the p-value 2, which it adjusts with the others, and in fry a mixed p-value of 0 for an empty set.
camera takes robust: the moderated t-statistics of the robust empirical Bayes estimator (Phipson et al. 2016), each turned into a z-score on its own degrees of freedom. limma’s camera has no such option and ignores it.

References

Wu D, Smyth GK (2012). Camera: a competitive gene set test accounting for inter-gene correlation. Nucleic Acids Research 40, e133. doi:10.1093/nar/gks461

Wu D, Lim E, Vaillant F, Asselin-Labat ML, Visvader JE, Smyth GK (2010). ROAST: rotation gene set tests for complex microarray experiments. Bioinformatics 26, 2176-2182. doi:10.1093/bioinformatics/btq401

Giner G, Smyth GK (2016). FRY: a fast approximation to ROAST gene set test with mean aggregated set statistics [version 1; not peer reviewed]. F1000Research 5, 2605 (slides). doi:10.7490/f1000research.1113351.1

Goeman JJ, Buehlmann P (2007). Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics 23, 980-987. doi:10.1093/bioinformatics/btm051

Barry WT, Nobel AB, Wright FA (2008). A statistical framework for testing functional categories in microarray data. Annals of Applied Statistics 2, 286-315. doi:10.1214/07-AOAS146

Hill GW (1970). ACM Algorithm 395: Student’s t-distribution. Communications of the ACM 13, 617-619. doi:10.1145/355598.355599

Brophy AL (1987). Efficient estimation of probabilities in the t distribution. Behavior Research Methods, Instruments, & Computers 19, 462-466. doi:10.3758/BF03205616

Zar JH (1999). Biostatistical Analysis, 4th edition. Prentice Hall, Upper Saddle River.

Phipson B, Lee S, Majewski IJ, Alexander WS, Smyth GK (2016). Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Annals of Applied Statistics 10, 946-963. doi:10.1214/16-AOAS920

Benjamini Y, Hochberg Y (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B 57, 289-300. doi:10.1111/j.2517-6161.1995.tb02031.x

Golub GH, Van Loan CF (2013). Matrix Computations, 4th edition. Johns Hopkins University Press, Baltimore.

cytepool.stats.gene_set_tests.camera(y, design, contrast, gene_sets, *, gene_ids=None, weights=None, inter_gene_correlation=0.01, allow_negative_correlation=False, use_ranks=False, trend=False, robust=False)#

Competitive test of gene sets for one contrast, allowing for their correlation.

limma’s camera (Wu and Smyth 2012): tests whether the genes of each set are more differentially expressed for contrast than the other genes, by a two-sample test of their moderated z-scores against those of the other genes, with the variance of the set’s mean inflated by the correlation between its genes. The module docstring gives the method. The default correlation of 0.01 is limma’s (since limma 3.29.6): it ranks co-regulated, interpretable sets higher than a correlation estimated for each set, whose error rate control is strict.

Parameters:

y – (n_samples, n_genes) log-expression values, finite, or a VoomResult, whose log-CPM values are tested with its weights.
design – (n_samples, n_coefficients) design matrix of full column rank, or anything else that fit_linear_model accepts, such as a Design. It needs residual degrees of freedom.
contrast – The contrast to test: the name of a named contrast of a Design, the name or 0-based position of a design column (a coefficient; negative positions count from the last column), or weights on the design columns: a vector in their order, a mapping of column name to weight (columns left out weigh 0), a pandas Series indexed by them, or a Contrasts holding one contrast.
gene_sets – Mapping of set name to its genes: gene ids (matched to gene_ids; ids not among them are left out, as limma’s ids2indices does), 0-based positions of columns of y when gene_ids is not given, or a boolean mask over the genes. A gene listed twice counts once. A GeneSets, such as one read from a GMT file by read_gmt(), is such a mapping, of gene ids.
gene_ids – One id per gene (column of y), unique, to name genes in gene_sets by id.
weights – Precision weights, positive and finite: (n_samples,) for sample weights or (n_samples, n_genes) for observation weights. Not with a VoomResult, which carries its own. limma’s per-gene weights (a vector of one weight per gene) have no counterpart.
inter_gene_correlation – The correlation between the genes of every set, in (-1, 1), or None to estimate it for each set from its residuals.
allow_negative_correlation – Whether a negative correlation, given or estimated, may shrink the variance of a set’s mean; otherwise it counts as zero.
use_ranks – Whether to test the ranks of the z-scores (a correlation-adjusted Wilcoxon-Mann-Whitney test) rather than their values.
trend – Whether the prior variance of the moderation follows the mean log-expression of the genes (limma’s trend.var).
robust – Whether to moderate the variances with the robust estimator, which limma’s camera does not offer.

Returns:

gene_set; n_genes, the number of the set’s genes among those tested; correlation, the estimated inter-gene correlation (NaN for a set of one gene) or the given one; direction, "up" or "down" (None where the p-value is NaN); p_value, two-sided; and adjusted_p_value, Benjamini-Hochberg over the sets (Benjamini and Hochberg 1995). A set with no genes, or with every gene, has NaN p-values.

Return type:

A table with one row per gene set, in the order of gene_sets

Raises:

TypeError – If an argument has the wrong type, or is or holds a numpy masked array.
ValueError – If y is not finite or has fewer than three genes, the design does not have one row per sample, is not of full column rank or leaves no residual degrees of freedom, the contrast or a gene set is invalid, a weight is not positive and finite, or no gene id of any set is among gene_ids.

Explained in: Choosing a method, Gene-set tests, Multiple testing

cytepool.stats.gene_set_tests.fry(y, design, contrast, gene_sets, *, gene_ids=None, weights=None, standardize='posterior_sd', trend=False, robust=False)#

Self-contained test of gene sets for one contrast, without rotations.

limma’s fry (Giner and Smyth 2016): tests whether the genes of each set are differentially expressed for contrast, in the same direction (the directional p-value) or in either (the mixed p-value), by the limit of the rotation test of Wu et al. (2010) with the mean statistic over infinitely many rotations. The module docstring gives the method and the deliberate differences from limma. limma’s gene contribution weights (gene.weights, or weights inside a set), which give genes signs and sizes in the set statistic, are not implemented.