pathVar: a new method for pathway-based interpretation of gene expression variability

Step 1: selecting a statistical measure to estimate gene expression variability

Variability is defined as the amount of dispersion in a given distribution (Larsen & Marx, 2017). Different statistical measures are available to estimate gene expression variability, and in genomics, the estimators that are most often employed are the standard deviation (SD) (Eq. (1)) (Hasegawa et al., 2015), the coefficient of variation (CV) (Eq. (2)) (Mar et al., 2011; Mason et al., 2014), and the median absolute deviation (MAD ( Eq. (3)) (Wijetunga et al., 2014). Conceptually, these statistics share similarities in their mathematical definition, and each one comes with their own advantages, as is often the case with any estimator that is applied to data.

Let x₁, …, x_n denote a univariate dataset, then (1)${\displaystyle \mathrm{SD}=\sqrt{\frac{1}{n-1}\sum _{i=1}^n{\left(x_i-\overline{x}\right)}^2}}$ (2)${\displaystyle \mathrm{CV}=\frac{\mathrm{SD}}{\overline{x}}}$ (3)${\displaystyle \mathrm{MAD}=median\left(|x_i-median\left(x\right)|\right).}$ A consensus on which estimator should be adopted for variability analysis remains unclear, and this is partly because performance of the estimators appears to be data-specific. The SD is often the preferred estimator for measuring gene expression variability since it is on the same scale as the average and therefore easy to interpret. The variance (i.e., SD²) is also characterized by the second central moment of a distribution, and hence the SD is directly linked to one of the fundamental metrics characterizing the probability distribution. A criticism of the SD, however, is that it may be dependent on average expression and therefore it is necessary to investigate the association between these two measures. To address this concern, the CV, which represents the ratio of the SD and the average, is often used; however it also has its own drawbacks. Most significantly, it can be affected by zero-inflation which occurs when near zero levels of average expression result in extremely large values of CV that do not necessarily reflect a large degree of overall variability. The MAD is a robust measurement of dispersion that behaves well in the presence of outlier data points. It has previously been used to study DNA methylation variability (Wijetunga et al., 2014).

Data simulations suggest that in general, the SD shows stronger performance as the variability estimator compared to MAD and CV in the pathVar method (see Text S1). The simulations were conducted to test the performance of each estimator under a wide variety of conditions. The choice of an estimator can also be motivated by the expectation that an ideal estimator of gene expression variability will be uncorrelated with average gene expression. If the variability estimator is highly correlated with the average expression, then trends observed for expression variability may simply be recapitulated by those observed for average expression. We therefore desire an estimator of variability that is the least correlated with average gene expression for an analysis of gene expression variability to be maximally informative. Our suggestion is to use the SD since overall this estimator displayed stronger performance in simulations or the estimator with the lowest correlation with average expression; however, ultimately, the final choice of the estimator is left to the user and can be easily specified in the software.

Step 2: identifying genes that belong to categories of high, medium or low levels of expression variability

Using a specified estimator, all genes are assigned to a discrete level of expression variability. In the one-group case, assignments are based on clustering the data using Normal mixture models via the mclust algorithm (Fraley & Raftery, 2002). The number of clusters or mixtures corresponds to the discrete levels of variability and is a parameter inferred by mclust. The mclust algorithm considers a finite range of values (starting with a minimum of one level to a maximum of four by default) and chooses the number that is most appropriate for the data using the Bayesian Information Criterion. The upper limit of four levels is recommended out of simplicity, where it is more useful to model a handful of variability levels, e.g., low, medium, high and very high, whereas for much larger numbers, the interpretation ceases to be as intuitive. In the pathVar package, the user is, however, free to use whatever upper limit is appropriate for the analysis.

For the two-group case, assignments are based instead on the 33rd and 66th percentiles that are computed from the combined gene expression variability distribution of all genes from both groups in the dataset. Low variability genes are defined as those with values falling between 0 and the 33rd percentile, medium variability genes are those between the 33rd percentile and the 66th percentile, and high variability genes being greater than the 66th percentile. The variability levels are defined by a fixed number of standardized percentiles, instead of inferred as in the one-group case, because it is possible that a different number of variability levels might be inferred for the two groups. To ensure a straightforward and balanced comparison between the two groups, the number of discrete categories are fixed and the boundaries for variability levels are based on percentiles calculated from all the data. Under both the one-group and two-group cases, the outcome of Step 2 is to identify the fixed boundaries that define each of the discrete levels of expression variability.

Step 3: testing pathways for aberrant gene expression variability signatures

The pathVar method decomposes each pathway into a set of counts corresponding to the number of genes in each discrete level of expression variability. By default, pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa & Goto, 2000) and REACTOME (Croft et al., 2014) are used; however, users may also import their own definitions. For a specific pathway, let O_i denote the observed count of genes annotated to this pathway with expression variability in the ith level, where i = 1, …, m discrete levels. The total number of genes n is defined as $n={\sum }_{i=1}^m{O}_i$. A statistical test is used to evaluate whether the set of counts (O₁, …, O_m) associated with a pathway deviates significantly from either a reference count distribution in the one-group case, or between the two phenotypes in the two-group case. In the one-group case, the reference distribution is obtained by counting the total number of genes in each level of gene expression variability (Fig. S1).

The null hypothesis that pathVar evaluates is:

• H₀₍₁₎: the expression variability counts observed for a specific pathway were generated under the same distribution as the reference counts for the one-group case.

• H₀₍₂₎: the variability-based counts for both groups were drawn from the same underlying distribution, for the two-group case.

To assess the deviation observed between the variability count distribution and an expected distribution for a specific pathway, we provide two statistical tests as options available for this analysis. Option 1 is the multinomial exact test, and Option 2 is the Chi-squared test.

For Option 1, the exact test models the counts as a multinomial random variable (O₁, …, O_m) ∼ Multi(n, p) where p = (p₁, …, p_m) and p_i is the probability that a gene belongs to the ith variability level for i = 1, …, m in the reference distribution. Under the exact test, the P-value is obtained by summing over all possible events that are less likely than the set of counts observed. If the P-value falls below a specified significance threshold (e.g., P-value < 0.05), sufficient evidence exists to reject the null hypothesis and we conclude that the pathway has an aberrant distribution of expression variability counts. In the one-group case, this means that the variability count distribution of the specific pathway deviates from the reference distribution i.e., all genes surveyed. In the two-group case, this result means that the variability is not identically distributed between the two contrasting phenotypes.

While the exact test is attractive because it calculates the P-value exactly, from a practical perspective, these kinds of tests can often be time and memory-intensive for genomic data, especially as the number of levels m and the number of genes in the pathway n grows. For example, consider a pathway with 30 genes, where the number of possible sets of counts to consider with three variability levels for the calculation of the P-value is 496. If the size of the pathway increases to 100, then the number of possibilities to consider grows to 5,151. From this basic example, it can be seen how increasing n or m can lead to some very extensive calculations.

Option 2 overcomes this limitation, as it is less computationally intensive, and tests the same null hypotheses using the Chi-squared test as an alternative to the exact test. The test statistic $X^2={\sum }_{i=1}^m\frac{{\left(O_i-E_i\right)}^2}{E_i}$ follows a Chi-squared distribution $X^2\sim {\chi }_{m-1}^2$ with m − 1 degrees of freedom: $X^2\sim {\chi }_{m-1}^2$. The expected counts E_i = n∗p_i are the expected number of genes in each level of expression variability within a specific pathway. The Chi-squared test achieves its computational efficiency because it is based on an asymptotic approximation, where the test becomes more accurate as n, the total number of genes in a pathway increases.

Both the exact test and the Chi-squared test assess whether a pathway has a significant change in gene expression variability. The resulting P-values from all pathways tested are adjusted using the Benjamini–Hochberg method which is based on controlling the false discovery rate (Hochberg & Benjamini, 1990). Finally, a pathway of interest can be investigated further using a Binomial test to assess within each level of variability, whether a significant deviation exists between either the pathway and the reference in the one-group case, or between the two phenotypic groups. This pairwise difference for each independent level provides a means to pinpoint the subset of genes within a pathway showing deviation in expression variability.

Power calculations indicate that pathVar has greater power to detect changes for pathways with skewed expression variability distributions than symmetric ones

Simulations were designed to investigate how experimental parameters influence the power of pathVar (Text S1). The results demonstrate that for a fixed effect size, the power of the Chi-squared test increases when the number of genes in a pathway also increases. Similarly, for a pathway of fixed size, the power of the test is higher for a larger effect size (Text S1, Fig. S1.5). The simulations also showed that the Chi-squared test had more power to detect differences between the reference distribution and a pathway that has a skewed variability distribution compared to a symmetric one. As an example, consider pathway p₁ where an effect size of 0.24 results in 80% power for approximately 150 genes in the pathway. For the symmetric distribution in pathway p₂, in order to obtain the same level of power, a pathway size of at least 200 genes is required (Text S1, Fig. S1.6).

The influence of sample size was also investigated on the power of the Chi-squared test for fixed pathway sizes (Text S1, Fig. S1.7). Higher levels of power were obtained for larger effect sizes, and this trend was more pronounced with smaller numbers of genes. As expected, the asymptote to 100% power occurred with smaller numbers of samples when the effect size was larger (0.5 or 0.9 versus 0.1). Overall, these calculations support the recommendation that experimental designs with larger samples will yield higher powered analyses, especially when the effect size is large.

Application of pathVar to human embryonic stem cell datasets identify significant pathways with distinct profiles of gene expression variability

To demonstrate the utility of pathVar in practice, the method was applied to three gene expression datasets of human ESCs (Text S2). The Bock dataset (Bock et al., 2011) had twenty ESC samples that were generated using microarray profiling, and the Yan dataset (Yan et al., 2013) had hESC samples profiled using single cell RNA-sequencing (RNA-seq) for eight cells at passage 0 (p0), and 26 cells at passage 10 (p10). pathVar was run independently on the three datasets, and pathways with a statistically significant deviation in their gene expression variability profile relative to the reference distribution were detected (Table S1, P-value < 0.01). Significant KEGG pathways reflected aberrant gene expression counts in ribosomes, metabolism (oxidative phosphorylation), the spliceosome and neurodegenerative pathways (Alzheimer, Parkinson and Huntington) (Tables S2A– S4A). Significant REACTOME pathways fell into three main classes representing cell cycle, metabolism and infectious disease (Tables S2B–S4B). Considerable overlap was observed in the significant results obtained between the three datasets (Fig. S2), where cell cycle was the most highly represented REACTOME category. Similar percentage distributions were observed (Fig. S2), demonstrating consistency in the results obtained from pathVar despite differences in technology platforms, ESC lines and passage number.

We next inspect the variability count distributions for individual pathways of interest. As an example, consider the Bock dataset, where the significant KEGG pathways for spliceosome, oxidative phosphorylation and ECM-receptor interaction each cover different aspects of hESC regulation. For both the spliceosome and oxidative phosphorylation pathways, greater transcriptional stability was manifested through a significantly higher number of low variability genes, in addition to significantly fewer medium variability genes compared to the reference distribution (Figs. 3A–3C). A significant reduction in genes with high variability in the ESCs was also observed for the oxidative phosphorylation pathway compared to the reference (Fig. 3C). The opposite trend was observed for genes in the ECM-receptor interaction where there was a significant increase in genes with medium, high and very high levels of expression variability (Fig. 3D). This pathway also had a concurrent reduction in low variability genes compared to the reference distribution.

Highlighting differences between ESC and iPSC usage of global gene expression programs using pathVar

Human iPSCs were also profiled by Bock using microarrays, and we use this data to investigate how ESCs and iPSCs differ with respect to expression variability. pathVar identified five KEGG and thirty REACTOME statistically significant pathways between the ESCs and iPSCs (Table S5, P-value < 0.01). The significant KEGG pathways reflected aberrant gene expression activity in ribosome, oxidative phosphorylation, DNA replication and disease processes (Huntington’s disease and Parkinson’s disease, Table S6). The significant REACTOME terms were associated with cell cycle, splicing, and metabolic processes as well as DNA replication and repair, namely homologous recombination pathways (Table S6).

We then compared the variability count distributions in ESCs versus iPSCs for two significant KEGG pathways, oxidative phosphorylation and DNA replication (Fig. 4). Both pathways showed increased gene expression stability in ESCs compared to iPSCs. For the oxidative phosphorylation pathway, a significantly higher number of low variability genes were in ESCs versus iPSCs. The same trend was observed for the DNA replication pathway where there was a significant reduction in the number of highly variable genes in ESCs compared to iPSCs.

Benchmarking pathVar results using gene expression variability versus average expression and two different GSEA approaches

pathVar results were benchmarked against those obtained using an average expression statistic to investigate the utility of expression variability analyses. Under this average-based setting, genes were first classified into discrete levels using average expression, corresponding to low, medium, and high levels of absolute expression. This average-based implementation allows for direct comparison with the results that are obtained when studying gene expression variability by substituting the sample mean directly in place of the variability statistic. In this way, we aim to distinguish the pathways that are identified as having changes in expression variability that are not detected by changes in average expression. Four statistically significant pathways (P-value < 0.01, Table S7) were identified as having differences in average expression between iPSCs and ESCs (three REACTOME terms: Heme biosynthesis, Sphingolipid metabolism, Metabolism of porphyrins; and one KEGG pathway African trypanosomiasis). These results do not seem very informative or relevant to stem cell regulation. This could be due to the over-simplified nature in which the average-based comparison was performed and suggested that this benchmarking approach could be improved.

We next evaluated the performance of pathVar against two versions of GSEA. The first method was based on generating P-values from limma (Ritchie et al., 2015) to determine differential expression between the two groups being compared. The P-values for genes belonging to the same pathway were then used as input to a Kolmogorov–Smirnov (KS) test which assessed whether the distribution of this set of P-values was identical between the two groups. Since the KS test is affected by tied values, we used the bootstrap version of the KS test (R function ks.boot). This version of GSEA is similar to the implementation of (Mootha et al. (2003). The second GSEA method is a stand-alone tool by Oron, Jiang & Gentleman (2008) which uses a linear model to assess the differential usage of a pathway or gene set and calculates permutation-based P-values to determine the significance of each pathway. The resulting set of P-values for all pathways was adjusted for multiple testing using the Benjamini–Hochberg method (Benjamini & Hochberg, 1995).

In general, the two GSEA approaches detected significant pathways that had a moderate to high degree of overlap with the significant pathways that were identified by pathVar for all three of the cancer-based comparisons, and the mouse hippocampus versus striatum comparison. This suggests that there are significant pathways common to both GSEA and pathVar, but the latter also identifies distinct pathways since the highest overlap for all six comparisons for GSEAlm was 61.1% (KEGG), 96.2% (REACTOME), and for GSEA using limma P-values was 22.2% (KEGG), 33.3% (REACTOME). There are a considerable number of pathways identified uniquely by pathVar, suggesting that there is value in using multiple approaches to determine differential pathway usage or expression.

Regulatory insights from other datasets confirm utility of looking at pathways with changes in both gene expression variability and average gene expression using pathVar and GSEA

Examples using stem cells illustrate how pathVar works in practice; however, the method can be applied to virtually any gene expression dataset. To highlight the generalizability of pathVar, we selected ten other datasets that cover a variety of biological and experimental variables. Collectively, these ten datasets were generated from multiple technology platforms that featured samples from human, mouse and parasite which represent a range of different disease phenotypes (see Text S3).

Three cancer RNA-seq datasets from the Cancer Genome Atlas (TCGA) were selected; these were the ovarian serous cystadenocarcinoma (OVC) (Cancer Genome Atlas Research, 2011), acute myeloid leukemia (AML) (Cancer Genome Atlas Research, 2013), and glioblastoma multiforme (GBM) (Brennan et al., 2013) cohorts. An infectious disease was included where transcriptomes from patients infected with cerebral malaria were profiled using microarrays (Daily et al., 2007), as well as the Plasmodium falciparum parasites that the patients were infected with (Feintuch et al., 2016). A genetic disorder was featured where patient-derived iPSCs were collected from Down syndrome (DS) donors and profiled using microarrays, with a set of matched controls from healthy subjects (Briggs et al., 2013). A microarray dataset from a normal human population via the Geuvadis study (1000 Genomes Project) was used (Genomes Project et al., 2010), as well as two mouse datasets that profiled tissues from two different regions of the brain, the hippocampus and the striatum using microarrays (Park et al., 2011). pathVar identified statistically significant pathways from KEGG and REACTOME pathway terms for the ten different one-group analyses (Table S8), and five independent two-group comparisons (Table S9) of all ten datasets.

The results from the different analyses were used to investigate the uniqueness of analyses based on GSEA versus gene expression variability. The number of significant pathways that had changes in both mean and variability via GSEA and pathVar respectively showed that for all cases, the amount of overlap in significant pathways differed depending on the datasets that were used. This suggests that average and variability-based statistics reflect different ways in which cells may use their transcription programs depending on the biological context (Tables S8 and S9). It is interesting to note that for the KEGG pathways, pathVar results for the DS versus WT iPSC comparison, and the mouse hippocampus versus striatum comparison both had zero overlap between average-based and variability-based significant pathways (Table S9A). In fact, both comparisons also yielded no significant pathways with a difference in average, whereas pathways were found to have a difference in gene expression variability. These two comparisons are extreme examples where the analyses of gene expression variability identify changes in the transcriptional program, whereas average-based analyses do not yield significant results.

Overall, it was apparent that the transcriptional features responsible for distinguishing one phenotype from another are exerted through changes in average expression or variability in expression for key pathways. To further investigate the relevance of these different modes, we focused only on the ten most significant pathways from the pathVar results obtained for the three cancer versus normal comparisons (Tables S10 and S11). For all three cancers, the KEGG DNA replication pathway and REACTOME “DNA strand elongation term” had significant changes in both average and variability of expression. Other terms with changes in both average and variability were related to DNA damage response pathways, such as “base excision repair” (Table S10) for AML versus normal, and “non-homologous end-joining” for OVC versus normal (Table S11).

Five KEGG pathways had changes in variability only, that were consistent in all three cancer comparisons; these were the pathways involved in Epstein-Bar virus infection, cell cycle, Fanconi anemia, lysosome and apoptosis (Table S10). The Epstein-Bar virus is associated with certain kinds of cancer like lymphoma or carcinoma. Apoptosis is also an important pathway for tumors because its inactivation is central in the development of cancer. Similarly, for the REACTOME terms, those unique to changes in variability were related to DNA repair and replication (SLBP dependent processing of replication-dependent histone pre-mRNAs) for the AML and GBM comparisons. For OVC, several terms were related to the cell cycle, e.g., G1 phase, cyclin D associated events in G1, cyclin A/B1 associated events during G2/M transition (Table S11C).