Modelling p-value distributions to improve theme-driven survival analysis of cancer transcriptome datasets

In clinical oncology, the discovery of knowledge from gene expression microarray experiments is often based on what has been described a "top-down" or hypothesis-free approach, where a prognostic model is derived from global tumour gene expression data in a way that does not require a priori knowledge of biological function [1]. A subsequent goal is to extract one or more higher-level biological "themes" from this primary model. In practice the task of inferring themes regarding biological function from a prognostic geneset, which may consist of a relatively small number of genes representing very distinct biological processes, has proved to be a major analytical bottleneck.

A contrasting approach is a "bottom-up" or hypothesis-driven method where sets of biologically-related genes (e.g. all transcripts coding for proteins of a particular signalling pathway) are defined first, and then these sets of genes are analysed to determine if any set acts as a prognostic indicator. Since the biological functions of the genesets are already known, sets with prognostic power can then play a more informative role in clinical decision-making [1]. We will refer to this particular approach as "theme-driven", where the biological function of a geneset becomes a theme, so as to distinguish it from other geneset-based methods in which genes do not share common functionality. Although useful in principle, a theme-driven approach in clinical oncology has not been widely used until recently, mainly due to the difficulty of defining sets of genes of common cellular function. However, the emergence of function-oriented classifications of genes, such as those developed by the Biocarta [2], Gene Ontology (GO) [3], Gene Map Annotator and Pathway Profiler (GenMAPP) [4], and Kyoto Encyclopedia of Genes and Genomes (KEGG) [5] consortia, in addition to those commercially available (e.g. Ingenuity [6], GeneGO [7]), has facilitated the comprehensive analysis of collections of biologically-related genes.

Common methods of theme-driven cancer survival analysis can involve clustering to segregate cancer samples into two groups according to the gene expression levels of the geneset of interest, and a subsequent survival test to determine whether a difference in prognosis exists between both groups [8, 9]. It is important to clarify what it is being tested in such studies. This type of experiments has two inherent null hypotheses. Goeman and Bühlmann [10] have discussed these two very different hypotheses, which result from analysing genesets in terms of subject-sampling and gene-sampling models. Although discussed in a different context, the implications are relevant to any geneset analysis of gene expression microarray data. According to the authors, in a subject-sampling model the subjects (i.e. tissue samples or patients) are treated as the sampling units, and it is natural to formulate a self-contained null hypothesis. In the authors' example, this null hypothesis states that a geneset is not differentially expressed. A significant p-value against this self-contained null hypothesis would mean that the same association found for the subjects analysed would also be found for a new set of subjects with high certainty. Yet, this p-value does not make any statement about how the geneset relates to others. On the other hand, in a gene-sampling model the genes become the sampling units, and it is straightforward to formulate a competitive null hypothesis. This hypothesis states that a particular geneset is differentially expressed as often as other genesets. In this case, a significant p-value would mean that the geneset would still be differentially expressed for the same subjects when new sets of genes are included on a new microarray. However, this p-value does not imply that the same association would be found for a new set of subjects. Throughout this paper, the expressions p-value₁ and p-value₂ will refer to p-values obtained from the first (self-contained) and the second (competitive) hypothesis tests, respectively. The expression p-value will be used in general terms.

In theme-driven survival studies, the two null hypotheses (with an exception mentioned below) must be tested. The first null hypothesis, which is self-contained and based on subject-sampling, claims that a particular geneset, independent of its composition (i.e. made up of either biologically-related or random genes), has no predictive power of survival time. A survival test, such as a log-rank test, is the usual choice for testing this hypothesis. A significant p-value₁ for a geneset would give one confidence that the association between the geneset and survival would also be found for new patients, but it does not demonstrate that the theme of the geneset is actually relevant. This is addressed with a competitive null hypothesis based on gene-sampling. This second null hypothesis states that a set of biologically-related genes is, at most, as correlated with survival as are random sets of unrelated genes. This hypothesis can be tested by an empirical test that compares the correlation to survival of a given geneset with that of random sets of unrelated genes. A significant p-value₂ against this hypothesis would mean that the theme of the geneset is associated with survival time, but only for the patients analysed in the dataset. Therefore, both null hypotheses must be rejected in order to claim predictive power of a theme.

If the distribution of p-values₁ for random sets of unrelated genes is uniform, which in turn implies that most sets of unrelated genes are non-informative with respect to survival, testing for the second null hypothesis becomes redundant. To illustrate this point with an example, assume that a biologically-related geneset g has prognostic power as assessed by rejecting the first null hypothesis at a p-value₁ p = 0.004 below a significance threshold α. If this geneset is systematically compared to n = 1000 non-informative (random) genesets (i.e. testing the second null hypothesis), it will outperform them approximately (1 - p) × n = 996 times (assuming a uniform distribution of p-values₁ for non-informative genesets). It is expected that 4 random genesets will be at least as significant as geneset g, which means that the empirical p-value₂ for g will also be equal to 4/1000, that is, 0.004 (see Methods). Since p-value₁ and p-value₂ are equal, it would be enough to only test the first null hypothesis. However, it has been reported that the distribution of p-values is often not uniform in gene expression experiments [11], which questions the adequacy of current theme-driven survival methods that assume uniformity. In our study, we are therefore interested in the distribution of p-values₁ for random sets of unrelated genes, since by incorrectly assuming a uniform distribution of these p-values₁, a researcher may falsely conclude that the theme of a geneset has prognostic power.

One situation that may lead to false positive conclusions in theme-driven survival analysis is when there are widespread changes in gene expression that lead to the patient samples clustering in a prognostically informative manner. Take as an extreme example a dataset consisting of two tumour types: benign early lesions and aggressive metastatic disease. Assume that the degree of transcriptome changes is so great that the data for any randomly selected group of genes clusters in such a way as to separate the patients into benign and aggressive categories. In such a scenario any biologically-based geneset (say all the genes within a particular GO category) might likewise cause an informative clustering of the data, and thus have prognostic power. The erroneous conclusion would be that this particular GO geneset has prognostic power greater than a random geneset, that is, the biological process (i.e. the theme of this geneset) is associated with prognosis. Formally speaking, the error has been caused by assumptions regarding the distribution of the p-values₁ for random sets of unrelated genes.

Although this hypothetical experiment is somewhat extreme, it illustrates the caution that must be taken when interpreting the p-values from simple applications of theme-driven survival techniques. The key to a more rigorous approach is to model the distribution of p-values₁ for random genesets, so as to then also provide accurate empirical p-values₂ for the biologically-related genesets. Although some hypothesis-free studies [12, 13] and some geneset-based experiments focusing on differential expression patterns [14] have implicitly dealt with the possible non-uniform distribution of p-values₁ for random genesets, theme-driven cancer survival studies have not, to date, acknowledged this problem.

Therefore, in this study we set out to investigate the distribution of p-values₁ for random sets of unrelated genes and its effects on theme-driven cancer survival methods. Using two publicly available datasets we demonstrate that such distribution is indeed non-uniform. We then describe a method that empirically approximates this distribution using a permutation approach based on 100,000 random sets of unrelated genes. Finally, we apply this methodology to these two datasets and compare our results with those of a previous theme-driven study [8]. This study has reported a significant association of a gene expression signature of fibroblasts in response to serum (core serum response, CSR, signature as called by the authors) with patient survival for many cancer datasets, including the breast and lung cancer sets used in our study. Our data question the validity of some of those associations, since the method employed in such study assumes that p-values₁ for random genesets are uniformly distributed. However, our data reveal that such assumption is frequently incorrect, which may result in erroneous conclusions, and we provide an empirical methodology to address this problem.

Unlike hypothesis-free studies in survival analysis, theme-driven experiments aim to test whether a set of biologically-related genes is associated with survival. In other words, researchers seek an association between the theme of the geneset and cancer prognosis. This implies that such predefined genesets must be significant predictors as determined by a survival test, which constitutes the first null hypothesis, and also must perform "better" than most sets of unrelated genes, which constitutes the second null hypothesis. If the distribution of p-values₁ for random genesets is uniform, which implies that most genes in gene expression datasets are uninformative with respect to survival time, then the second null hypothesis need not to be tested. However this is rarely the case; p-values₁ are usually non-uniformly distributed.. Therefore, any attempt to determine significance of the theme of a geneset based on the first hypothesis test is rendered meaningless [11]. Our results illustrate this (Figure 2), and raise the issue that simple methods in theme-driven survival studies that assume uniformity of p-values₁ for random genesets may yield false conclusions. Even when an empirical distribution is used to calculate empirical p-values₂ for the original themes of interest, it must be noted that such p-values₂ do not claim anything about the predictive power of the themes for new sets of patients. Theoretically, this claim has to be assessed by a prior (or posterior) independent hypothesis test that verifies an association between a geneset and survival time. In our study, such hypothesis testing has been done beforehand when a p-value₁ was calculated for each geneset with a log-rank test. It is arguable that such a simple test can have predictive power for new patients, but it is a common procedure in survival analysis and a thorough examination of this topic is outside the scope of this paper.

In the analysed breast cancer dataset, many genesets were correlated to cancer prognosis when considering p-values₁ from the first hypothesis test, both for OS and for RFS. After testing the second null hypothesis, most of these genesets did not reach statistical significance. This is consistent with the skew towards low p-values₁ of the empirical distributions for random genesets (Figure 2A). Conversely, the lung cancer dataset identified several genesets that were significant for both the first and second hypothesis tests. This is in turn consistent with the lesser skew towards low p-values₁ (Figure 2B).

Chang et al. have reported a correlation of their core serum response (CSR) gene expression signature to survival time in many cancer datasets, including the same breast and lung cancer sets used in our study. As with other theme-driven approaches, these authors did not consider the impact of non-uniformly distributed p-values₁ for random genesets, so they tested only the first null hypothesis. Our findings suggest that this may have biased their results, leading to false correlations of the CSR signature to cancer survival. While it remains possible that this geneset is a significant predictor of survival in new patients for breast cancer (rejection of the first null hypothesis), it cannot be claimed that the biology of that geneset is relevant, given that many random genesets achieve at least the same predictive power (acceptance of the second null hypothesis). Our empirical approach improves upon their methodology and provides a technique to avoid this pitfall.

Although we have illustrated this issue using some commonly available ontologies, it is important to note that current biological databases still offer a somewhat limited perspective with respect to the annotation of gene function. Although biological ontologies (e.g. KEGG, Biocarta, GO) have successfully grouped genes in terms of biological function, such lists do not account for relationships and dependencies among genes. Moreover, different databases have overlapping annotations, hence analysing genesets from various ontologies may result in redundant tests due to duplication or multiplication of genesets. In principle our methodology can be applied equally well to other geneset collections, including proprietary sets and those available through commercial providers.

In this study we used an FDR method to control for multiple testing, which assumes that the test statistics are only weakly correlated [15]. This is conservative, because it does not account for the strong dependencies and correlations in our genesets. Filtering by use of a non-stringent p-value (or FDR value) has been suggested, since screening studies with small sample sizes will have inherently unstable feature rankings in any statistical test; thus, reproducibility of gene expression microarray experiments will be increased by more permissive thresholds [16], which is particularly important for theme-driven cancer survival studies. Many genesets that were significant for both hypothesis tests and also after correction for multiple testing are biologically plausible in the light of long-established biological knowledge and recently published studies. Our results show a correlation of the GO term "fatty acid metabolism" to breast cancer overall survival. This mirrors a recent experimental finding that reports an association between the failure to respond to tamoxifen treatment and the expression patterns of genes involved in cholesterol and fatty acid metabolism for estrogen receptor positive (ER+) breast tumours [17]. Our data also revealed an association between the GO term "receptor mediated endocytosis" and overall survival time in lung cancer. Again, the importance of proteins involved in receptor mediated endocytosis in lung cancer development, progression, and metastasis has been discussed in the literature [18]. Our results show additionally a relationship between the KEGG pathway "MAPK signaling" and overall survival in lung cancer. This is a confirmation of the well-known oncogenic properties of the MAPK signaling pathway, whose role in cancer development and progression have been extensively addressed in the literature for lung and many other cancers [19–21]. Lastly, our data suggests a correlation of the GO term "apical plasma membrane" with lung cancer overall survival. Researchers have reported that lipid trafficking to the apical plasma membrane of pulmonary epithelial cells is thought to be a protective stress response [22]. It can be speculated that lung cancer cells may develop by deregulation of this mechanism.

Our results show that p-values₁ for random genesets in theme-driven survival studies are frequently not uniformly distributed. This impairs common-practice methods, and results in false positive associations between the theme of a geneset and survival. Our proposed empirical approach correctly avoids such pitfalls, and the plausible biological observations obtained suggest that a theme-driven method that correctly tests the two inherent null hypotheses can contribute to further understanding of the biology of cancer (and other clinical survival datasets), and help bridge the gap between gene expression predictors and biological knowledge.