Rotation gene-set tests

Rotation gene-set tests#

Rotation gene-set tests of one contrast of a linear model: limma’s roast, mroast and romer.

A rotation test (Langsrud 2005) draws random rotations of the residual space of the linear model fitted to every gene, which leave the joint distribution of the data unchanged under the null hypothesis, and compares a set statistic with its values under the rotations. roast() and mroast() are self-contained tests of one gene set and of many (Wu et al. 2010): whether the genes of a set are differentially expressed for the contrast at all. romer() is a competitive test of many sets (Majewski et al. 2010): whether the genes of a set rank higher than the other genes. fry() gives the limit of mroast with the mean statistic over infinitely many rotations without drawing any. Results match limma 3.62.1 within Monte Carlo error, and exactly when CytePool is given limma’s rotations, except where noted under “Deliberate differences”.

Effects and rotations:

Each gene’s contrast effect u_0 and residual effects u_1, …, u_d are those of cytepool.stats.gene_set_tests: coordinates of its data, whitened by its precision weights, in the basis of the Householder QR decomposition of the design reparametrised so that the contrast is its last coefficient. Under the null hypothesis they are independent with variance sigma^2, so one rotation applied to the (d + 1)-vectors of every gene leaves the joint distribution of the data, and the correlation between genes, unchanged (Wu et al. 2010, section 4.2). A rotation is given by its first row r, a point drawn uniformly from the unit sphere in d + 1 dimensions (a vector of standard normal draws divided by its length): a gene’s rotated contrast effect is u* = r’u and its rotated residual sum of squares |u|^2 - u*^2. The rotated moderated t-statistic of a gene uses the prior variance s_0^2 and prior df d_0 estimated once from the observed residual variances of all genes (Wu et al. 2010, section 2.4): t* = u* / s*, with s*^2 = (d_0 s_0^2 + |u|^2 - u*^2) / (d_0 + d). The rotations come from rng.standard_normal, 4,096 at a time, which bounds the memory they take; numpy gives the same stream however its draws are split, so a seeded generator gives the same p-values every time. limma draws its own rotations from R’s generator, so the two programs agree within Monte Carlo error.

roast and mroast:

Each gene’s moderated t becomes a z-score by Bailey’s (1980) normalizing transformation, limma’s default, z = (8 nu + 1) / (8 nu + 9) (nu + 19/12)^(1/2) log(1 + t^2 / (nu + 1/12))^(1/2) with the sign of t (limma’s zscoreT(approx = TRUE) to 1e-14, found by probing), on nu = min(d_0 + d, 10^4) degrees of freedom, limma’s limit (found by probing). With an infinite prior df the exact z-score is t itself, which the limit lowers by about (t^2 + 1) / (4 10^4) of its value (1.2e-4 at t = 2, 2.5e-3 at t = 10); the same transformation applies to the observed and the rotated statistics, so the test remains exact. With contribution weights a_g (1 unless given; a gene of weight zero leaves the set), A = sum |a_g| and m genes in the set, the set statistics of Wu et al. (2010, section 3.3) are, for the up, down and mixed alternatives:

  • "mean": sum a_g z_g / A, its negative, and sum |a_g z_g| / A.

  • "floormean": sum max(a_g z_g, 0) / A, sum max(-a_g z_g, 0) / A, and sum |a_g| max(|z_g|, c) / A, with c = 0.674, the square root of the median of the chi-squared distribution on 1 df.

  • "mean50": the mean of the h largest a_g z_g, of the h largest -a_g z_g, and of the h largest |a_g z_g|, with h = ceiling(m / 2), the smallest integer at least half the number of genes, as the paper’s text and limma have it (its formula, ceiling((m + 1) / 2), is larger for even m).

  • "msq": sum |a_g| z_g^2 / A over the genes with a_g z_g > 0, over those with a_g z_g < 0, and over all genes.

The two-sided directional statistic is the larger of the up and down statistics. A rotated statistic at least as large as the observed one counts as at least as extreme, and the p-value is (b + 1) / (B + 1) for b such rotations out of B (Phipson and Smyth 2010). For the up and down alternatives each rotation r also gives -r, which negates every rotated z-score, so b counts 2B rotated statistics and the p-value is (b + 1) / (2B + 1), as in limma. A gene is active for the up alternative when sign(a_g) z_g > sqrt(2), for the down alternative when sign(a_g) z_g < -sqrt(2), and for the mixed alternative when |z_g| > sqrt(2) (Wu et al. 2010, section 3.5); the two-sided alternative reports the larger of the up and down proportions, as limma does. mroast’s direction is “up” when the up p-value is smaller than the down p-value and “down” otherwise (found by probing), and its adjusted p-values are Benjamini and Hochberg’s, applied to the mid-p-values p - 1 / (2 (B + 1)) (Routledge 1994) and then raised to at least the p-value, as limma does by default. All sets share the same rotations.

romer:

Every gene’s observed moderated t-statistic is ranked among all genes (ties get their average rank), and a set’s statistic is the mean rank of its genes: the ranks of t for the up alternative, of -t for the down alternative and of |t| for the mixed alternative ("mean"); the same after setting negative values to zero, and for the mixed alternative |t| below 1 to 1 ("floormean", the floor of 1 found by probing); or the mean of the ceiling(m / 2) largest ranks of the set ("mean50"). Each rotation ranks the rotated t-statistics of all genes again. The p-value of each alternative is (b + 1) / (B + 1), where b counts the rotations whose statistic is at least the observed one.

Before rotating, each gene’s contrast effect is shrunk, so that the differential expression of the gene does not inflate its rotated statistics (limma’s shrink.resid; formula found by probing): u_0 becomes u_0 / (1 + P_g v_0 / v)^(1/2). P_g is the posterior probability that the gene is differentially expressed and v_0 the prior variance of the contrast for such genes, from the log-odds (B-statistics) of moderated_statistics() with the t-statistics on d_0 + d degrees of freedom (their limit for an infinite prior df) and 1 - pi_0 differentially expressed genes, and v is the unscaled variance of the contrast. The proportion of true null hypotheses pi_0 is (2 / (n + 1)) sum min(p_(i), i / n) over the n ordered p-values p_(i) of the moderated t-statistics (limma’s propTrueNull(method = "lfdr"), found by probing). The rotated residual sum of squares remains |u|^2 - u*^2 with |u| the length of the gene’s unshrunk effects, as in limma (found by probing). romer multiplies the residual effects by the sign of the last diagonal element of R, as limma’s romer does and its roast does not (found by probing); a uniformly random rotation does not depend on the basis of the residual space, but limma’s rotations then mean the same here. The bounds of v_0 are those of moderated_statistics (stdev_coef_lim), in log-expression units, so, as in limma, romer’s p-values can change when the data are rescaled.

Deliberate differences from limma:

  • roast and mroast count a rotated statistic equal to the observed one as at least as extreme, as Wu et al. (2010) and Phipson and Smyth (2010) define the p-value (and as limma’s romer does); limma’s roast counts only larger ones. The two differ only where the observed statistic is the smallest the statistic can take, which happens with positive probability for "floormean" and "msq": a set with no gene on one side of zero has a one-sided p-value of 1 here and of about 1/2 in limma, and with "floormean" a set whose |z| are all below 0.674 has a mixed p-value of 1 here.

  • "msq" with contribution weights follows Wu et al. (2010). With weights of both signs, limma’s up, down and two-sided msq p-values differ from it, and with more than 1,000 rotations (limma draws them 1,000 at a time) its msq p-values change when every weight is multiplied by the same number (found by probing). Its mixed msq p-values, and those of sets whose weights have one sign, agree with CytePool’s for up to 1,000 rotations.

  • "floormean"’s mixed statistic weighs each gene by |a_g|, as the paper’s mixed statistics do; limma weighs it by a_g, so a negative weight turns a gene with a large |z| into evidence against the set (a one-gene set’s mixed p-value goes from 0.003 to 0.997 when its weight changes sign).

  • A gene of contribution weight zero leaves the set: it counts neither in n_genes nor among the values that "mean50" averages. limma keeps it in both (and leaves it out of the active proportions, as CytePool does).

  • romer’s posterior probability P_g is 1 / (1 + exp(-lods)). limma computes exp(lods) / (1 + exp(lods)), which is NaN beyond lods = 709.8 (|t| above about 38 with an infinite prior df); the gene’s shrunk contrast effect is then NaN and the p-values of every set change (found by probing).

  • A set without genes, or whose genes all have weight zero, has NaN p-values and no direction and is left out of the adjustment of the other sets; limma reports “down” and NA, or the smallest possible p-value when all weights are zero. A gene listed twice counts once; limma counts it twice.

  • mroast tests every set with the same rotations; limma draws new rotations for each set. Each set’s p-value is a valid rotation p-value either way.

  • A z-score is computed without overflow for |t| above 1e154, where limma’s is infinite.

References

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

Majewski IJ, Ritchie ME, Phipson B, Corbin J, Pakusch M, Ebert A, Busslinger M, Koseki H, Hu Y, Smyth GK, Alexander WS, Hilton DJ, Blewitt ME (2010). Opposing roles of polycomb repressive complexes in hematopoietic stem and progenitor cells. Blood 116, 731-739. doi:10.1182/blood-2009-12-260760

Langsrud O (2005). Rotation tests. Statistics and Computing 15, 53-60. doi:10.1007/s11222-005-4789-5

Phipson B, Smyth GK (2010). Permutation p-values should never be zero: calculating exact p-values when permutations are randomly drawn. Statistical Applications in Genetics and Molecular Biology 9, article 39. doi:10.2202/1544-6115.1585

Bailey BJR (1980). Accurate normalizing transformations of a Student’s t variate. Journal of the Royal Statistical Society C 29, 304-306. doi:10.2307/2346906

Routledge RD (1994). Practicing safe statistics with the mid-p. Canadian Journal of Statistics 22, 103-110. doi:10.2307/3315826

Smyth GK (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology 3, article 3. doi:10.2202/1544-6115.1027

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

Hoare CAR (1961). Algorithm 65: find. Communications of the ACM 4, 321-322. doi:10.1145/366622.366647

Wirth N (1976). Algorithms + Data Structures = Programs. Prentice-Hall, Englewood Cliffs.

cytepool.stats.rotation_gene_set_tests.mroast(y, design, contrast, gene_sets, *, gene_ids=None, weights=None, gene_weights=None, set_statistic='mean', n_rotations=1999, trend=False, robust=False, mid_p=True, rng=None, n_threads=None)#

Self-contained rotation tests of many gene sets for one contrast.

limma’s mroast: roast() for each set, all sets with the same rotations, and with their p-values adjusted for testing many sets. The module docstring gives the method and the deliberate differences from limma.

Parameters:

  • y – Log-expression values or a VoomResult, as for roast().

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

  • contrast – The contrast to test, as for roast().

  • gene_sets – Mapping of set name to its genes, each as gene_set of roast(): gene ids, positions or a boolean mask, or a mapping of gene to its contribution weight. 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, unique, to name genes by id.

  • weights – Precision weights, as for roast().

  • gene_weights – (n_genes,) contribution weights, as for roast(); not with sets that give their own weights.

  • set_statistic – "mean", "floormean", "mean50" or "msq".

  • n_rotations – Number of rotations B, at least 1.

  • trend – Whether the prior variance follows the mean log-expression.

  • robust – Whether to estimate the prior robustly.

  • mid_p – Whether the adjustment uses mid-p-values, p - 1 / (2 (B + 1)), as limma does by default (limma’s midp), which suit false discovery rates; the adjusted p-values are then at least the p-values.

  • rng – A numpy.random.Generator, an integer seed, or None.

  • n_threads – The most threads to use for the rotations.

Returns:

gene_set; n_genes; proportion_down and proportion_up, the proportions of the set’s genes active down and up; direction, "up" or "down" (None where the p-value is NaN); p_value, the two-sided directional p-value; adjusted_p_value, Benjamini-Hochberg over the sets (Benjamini and Hochberg 1995), limma’s default (its other adjust.method choices are not offered); mixed_p_value, the p-value of the mixed alternative; and adjusted_mixed_p_value. A set without genes has NaN proportions and 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 – As roast(), or if gene_sets is empty or no gene id of any set is among gene_ids.

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

cytepool.stats.rotation_gene_set_tests.roast(y, design, contrast, gene_set, *, gene_ids=None, weights=None, gene_weights=None, set_statistic='mean', n_rotations=1999, trend=False, robust=False, rng=None, n_threads=None)#

Self-contained rotation test of one gene set for one contrast.

limma’s roast (Wu et al. 2010): tests whether the genes of the set are differentially expressed for contrast, for each alternative: down, up, either (two-sided) and mixed (in either direction, gene by gene), by comparing a set statistic of the genes’ moderated z-scores with its values under random rotations of the residual space. The module docstring gives the method and the deliberate differences from limma. With many sets, mroast() reports one row per set.

Parameters:

  • y – (n_samples, n_genes) log-expression values, finite, or a VoomResult, whose log-CPM values are tested with its weights. All genes are used to estimate the prior of the residual variances.

  • 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, as for camera(): the name of a named contrast of a Design, the name or 0-based position of a design column, or weights on the design columns.

  • gene_set – The set’s genes: gene ids (matched to gene_ids; ids not among them are left out), 0-based positions of columns of y when gene_ids is not given, or a boolean mask over the genes; or a mapping of gene (id or position) to its contribution weight. A gene listed twice counts once, and a gene of weight zero leaves the set. One set of a GeneSets, gene_sets[name], is a tuple of gene ids.

  • gene_ids – One id per gene (column of y), unique, to name genes 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.

  • gene_weights – (n_genes,) contribution weights, finite, one per gene: a negative weight marks a gene expected to change in the opposite direction to the set, and larger weights count more in the set statistic. Not with a set that gives its own weights.

  • set_statistic – "mean" (limma’s default), "floormean", "mean50" or "msq"; see the module docstring. The mean detects sets in which most genes change in the same direction; the others detect sets in which fewer genes change, msq the fewest.

  • n_rotations – Number of rotations B, at least 1. The smallest p-value is 1 / (B + 1), and 1 / (2B + 1) for the up and down alternatives.

  • trend – Whether the prior variance follows the mean log-expression of the genes (limma’s trend).

  • robust – Whether to estimate the prior robustly (limma’s robust).

  • rng – A numpy.random.Generator to draw the rotations from, or an integer seed for one, or None for a fresh, unseeded generator.

  • n_threads – The most threads to use for the rotations; None uses numba’s current setting. The result does not depend on it.

Returns:

A table with one row per alternative, "down", "up", "up_or_down" (two-sided) and "mixed": alternative; n_genes, the number of the set’s genes among those tested; active_proportion, the proportion of the set’s genes active for that alternative; and p_value. A set without genes has NaN proportions and p-values.

Raises:

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

  • ValueError – If y is not finite, 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 the set is invalid, a weight is not positive and finite, a contribution weight is not finite, an option is invalid, or no gene id of the set is among gene_ids.

Explained in: Choosing a method, Gene-set tests

cytepool.stats.rotation_gene_set_tests.romer(y, design, contrast, gene_sets, *, gene_ids=None, weights=None, set_statistic='mean', n_rotations=9999, shrink_contrast_effects=True, rng=None, n_threads=None)#

Competitive rotation test of gene sets by the ranks of their genes.

limma’s romer (Majewski et al. 2010): tests whether the genes of each set rank higher than the other genes for contrast, by their moderated t-statistics, up, down or in either direction gene by gene (mixed), comparing the set’s mean rank with its values when every gene’s rotated statistics are ranked again. Like gene set enrichment analysis it pits the sets against one another, but it uses rotations instead of permuting samples, so it works for any contrast of a linear model with residual degrees of freedom. The module docstring gives the method.

Parameters:

  • y – (n_samples, n_genes) log-expression values, finite. All genes are ranked.

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

  • contrast – The contrast to test, as for roast().

  • gene_sets – Mapping of set name to its genes: gene ids (matched to gene_ids; ids not among them are left out), 0-based positions, or a boolean mask over the genes. A gene listed twice counts once. A GeneSets is such a mapping, of gene ids.

  • gene_ids – One id per gene, unique, to name genes by id.

  • weights – (n_samples,) sample weights, positive and finite (limma’s array.weights). romer has no observation weights.

  • set_statistic – "mean" (limma’s default), "floormean" or "mean50"; see the module docstring.

  • n_rotations – Number of rotations B, at least 1; limma’s default is 9999. The smallest p-value is 1 / (B + 1).

  • shrink_contrast_effects – Whether each gene’s contrast effect is shrunk before rotating, as limma does by default (limma’s shrink.resid).

  • rng – A numpy.random.Generator to draw the rotations from, or an integer seed for one, or None for a fresh, unseeded generator.

  • n_threads – The most threads to use for the rotations; None uses numba’s current setting. The result does not depend on it.

Returns:

gene_set; n_genes; and up_p_value, down_p_value and mixed_p_value, the p-values of the three alternatives, which, as in limma, are not adjusted for testing many sets (adjust_p_values() adjusts them). A set without genes 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 – As mroast(), or if y is a VoomResult or weights are observation weights, which romer does not take.

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