Comparison of scores for bimodality of gene expression distributions and genome-wide evaluation of the prognostic relevance of high-scoring genes

Cluster Algorithms

Model-based clustering

The assumption for this approach is that the distribution Y of a gene with bimodal expression can be expressed as a mixture of two gaussian distributions Y₁ and Y₂ with parameters μ₁, and μ₂, , respectively.

(1)

where Δ ∈ {0, 1} with P(Δ = 1) = π. π is the probability of belonging to Y₁. Let ϕ _θ denote the normal density with parameters θ = (μ, σ²). Then the density of Y is given by

(2)

The parameters of the mixture model are determined using the EM-algorithm.

To fit the models to our data we make use of the R package mclust[14, 15]. Typically equality of variances is not assumed and we fit models with unequal variances. However, in cases with only one outlier in the expression distribution a model with unequal variances is not suitable since one component has variance 0. In this case we fit a two component model with equal variances.

Tests for bimodality

The first two scores for testing bimodality quantify whether the partitioning of the expression values into two groups is adequate. We want to quantify whether assuming two clusters is more appropriate than assuming one cluster. The tests are based on the decomposition of the Total Sum of Squares as

(3)

where BSS is the Between Cluster Sum of Squares and WSS the Within Cluster Sum of Squares.

Variance reduction score

The variance reduction score (VRS) is defined as the ratio of WSS and TSS:

(4)

VRS measures the proportion of variance reduction when splitting the data into two clusters. The value of this score lies in the interval [0; 1], and a low score indicates an informative split.

Weighted variance reduction score

The weighted variance reduction score (WVRS) measures the variance reduction independent of the cluster sizes. In the numerator we calculate the mean of the two within cluster variances. The denominator is the sample variance as above for the VRS. We define

(5)

In this case the variance of a cluster with few observations has the same influence as the variance of a large cluster. The value of this score can be larger than 1. Again, a low score reflects bimodality. WVRS has the ability to also identify splits into two clusters with extremely unequal sample sizes.

Kurtosis

Teschendorff et al. (2006) [8] proposed an approach to identify genes with bimodal density based on model-based clustering and kurtosis. The general approach is based on two steps. To select the optimal number of clusters and to find bimodal expression profiles a clustering algorithm is used together with a model selection criterion. Then the kurtosis can be used to rank the bimodal genes according to if they define major subgroups or outlier subgroups. Kurtosis is related to the fourth central moment and can be defined by

(6)

Given a gene's expression values x = (x₁,..., x _n ) an unbiased estimate for the kurtosis is given by

(7)

A gaussian distribution has kurtosis K = 0, whereas most non-gaussian distributions have either K > 0 or K < 0. Specifically, a mixture of two approximately equal mass normal distributions must have negative kurtosis since the two modes on either side of the center of mass effectively flatten out the distribution. A mixture of two normal distributions with highly unequal mass must have positive kurtosis since the smaller distribution lengthens the tail of the more dominant normal distribution [8]. If there is an 80%-20% split of the samples into two groups, then the kurtosis is close to 0 [9]. Therefore biologically interesting genes might be missed.

We calculate the kurtosis for all genes without a prior feature selection step.

Likelihood Ratio

Using the likelihood ratio of a normal model and a mixture normal model to identify bimodal distributions was suggested by Ertel and Tozeren (2008) [6]. We fit a normal model and a mixture normal model to the expression values for each gene using the mclust package and calculate the likelihood ratio

(8)

where L₁ is the likelihood of the normal model with one component and L₂ is the likelihood of a normal mixture model with two components with unequal variance.

Small ratios indicate that the distribution is unimodal, whereas large ratios suggest that the expression values are bimodally distributed.

Bimodality Index

The Bimodality Index introduced by Wang et al. (2009) [9] is a criterion to identify and rank bimodal signatures from gene expression data. It is assumed that the distribution of a gene with bimodal expression can be expressed as a mixture of two normal distributions with means μ₁ and μ₂ and equal standard deviation σ. The standardized distance δ between the two populations is given by

(9)

For identifying genes with bimodal distribution the null hypothesis is δ = 0. Then the bimodality index (BI) is defined by

(10)

where π is the proportion of observations in the first component.

For a given data set δ and π can be estimated using mclust. Then BI can be calculated using the estimated values. Larger values of BI correspond to bimodal distributions where the two components are easier to distinguish. Wang et al. recommend BI = 1.1 as a cutoff to select bimodally distributed genes.

Outlier-Sum Statistic

Another approach to group genes by expression values is based on the calculation of outlier sums as proposed by Tibshirani and Hastie (2007) [12]. This method is able to detect genes with unusually high or low expression in some but not all samples. Tibshirani and Hastie assume that the samples fall into a reference and a disease group. However, as we only consider tumor samples we modified the method for our purpose.

Let med _i and mad _i be the median and the median absolute deviation of the expression values of gene i. First, the expression values x _ij of each gene are standardized as follows:

(11)

Let q _r (i) be the r th percentile of the values and IQR(i) = q₇₅(i) - q₂₅(i) the interquartile range. All values smaller than IQR(i) - q₂₅(i) or greater than IQR(i) + q₇₅(i) are defined to be outliers. The outlier-sum statistic for positive outliers is defined as

(12)

W _i is large if there are many outliers in the data or few extreme outliers with high values, and W _i is zero if there are no outliers. Analogously the outlier-statistic for negative outliers is defined as

(13)

The outlier-sum statistic is the maximum of W _i and W' _i in absolute value.

For the survival analysis we need to split patients into two groups. We consider the outliers that are used to calculate the outlier sum as one group and all other observations as the other group.

Testing for prognostic relevance

For each gene we compare the hazard rates of two patient subgroups obtained by the different methods. The goal is to determine whether the survival in the two groups differs. We test the hypothesis H₀: h₁(t) = h₂(t), for all t ≤ τ, against the alternative that the hazard rates differ for some t ≤ τ. τ is the largest time at which both groups have at least one subject at risk. We make use of the logrank test which is a nonparametric test for right-censored data first proposed by Mantel (1966) [17].

Let t ₁ < t₂ < ... < t _D be the distinct death times in the pooled sample. Then d _ij is the number of events in the j th sample out of Y _ij individuals at risk at time t _i , j = 1, 2, i = 1,..., D. Let d _i = d_{i 1}+ d_{i 2}and Y _i = Y_{i 1}+ Y_{i 2}be the number of events and the number at risk in the combined sample at time t _i . The test statistic Z of the logrank test is given by

(14)

When the null hypothesis is true and the sample sizes in the two groups are similar, asymptotically Z has a standard normal distribution. In the case of highly unequal sample sizes, especially if one group contains less than five patients, the asymptotic distribution under the null hypothesis is not appropriate anymore [18, 19]. Thus we do not make use of the normal distribution for obtaining significance values. Instead, we use a permutation test. We calculate the test statistic for 100.000 random permutations of the patients and determine the percentage of values more extreme than the observed value of the test statistic. Using k-means we always get two groups of patients. However, using the model-based clustering approach and the outlier sum approach it is possible that for some genes all patients are in the same group so that the logrank test can not be applied. As these genes rather have unimodal expression profiles we decided to set the corresponding p-value 1.

Adjustment for multiple testing

We use the logrank test described in the previous Section to test for prognostic relevance for each of the genes. Here, we have to take the multiple testing problem into account. The significance level α of a hypothesis test controls the Type I error, i.e. the probability of rejecting a true null hypothesis. If m independent true hypotheses are tested simultaneously one would expect m · α p-values p < α. Therefore one has to adjust for multiple testing.

There are different correction methods. The Bonferroni-Holm method[20] controls the Familywise Error Rate (FWER), i.e. the probability of rejecting at least one true null hypothesis.

The method suggested by Benjamini and Hochberg (1995) [21] controls the False Discovery Rate (FDR), which is the expected proportion of falsely rejected hypotheses among all rejected hypotheses. This procedure is less conservative than the Bonferroni-Holm method and often more appropriate for high-dimensional microarray data.

Let p₍₁₎, . . ., p_(m)denote the ordered p-values of m hypothesis tests, i.e. p_(i)≤ p_(i+1). Then the adjusted p-values, also referred to as q-values, can be calculated iteratively as follows

(15)

where q_(m)= p_(m).

Analyzing gene lists for enrichment with prognostic genes

To determine whether there is an enrichment of the top-scoring bimodality genes with prognostic genes, we rank the genes based on the different bimodality scores. We then select the 250 genes with smallest p-values of the logrank test. These genes differ for the different approaches (i.e. k-means, model-based clustering, outlier sum). In this context, we call these genes prognostic genes.

An enrichment of the top-scoring bimodality genes with prognostic genes corresponds to small ranks of prognostic genes in the ranked list. We test the null hypothesis that the ranks of the prognostic genes are uniformly distributed among the bimodality-ranked genes with a one-sample Kolmogorov-Smirnov test, which is a nonparametric test of equality with a given one-dimensional reference function. The Kolmogorov-Smirnov statistic quantifies a distance between the empirical distribution function of the sample and the cumulative distribution function of the reference distribution.

We perform the one-sided Kolmogorov-Smirnov test for testing the hypothesis that the cumulative distribution function of the ranks of the prognostic genes is not greater than the distribution function of a uniform distribution. Let F _n denote the empirical distribution function of the sample and F the distribution function of the reference distribution, then the test statistic is given by

(16)

The null distribution of this statistic is calculated under the null hypothesis that the sample is drawn from the reference distribution. If F is continuous the distribution of does not depend on F but just on the sample size n. For n > 40 approximative p-values can be derived from the asymptotical distribution of [22].

Patients and gene array analysis

The Mainz study cohort consists of 194 node-negative breast cancer patients who were treated at the Department of Obstetrics and Gynecology of the Johannes Gutenberg University Mainz between the years 1988 and 1998 [1]. All patients underwent surgery and did not receive any systemic therapy in the adjuvant setting. Gene expression profiling of the patients' RNA was performed using the Affymetrix HG-U133A array, containing 22283 probe sets, and the GeneChip System [23]. The normalization of the raw data was done using RMA [24] from the Bioconductor [25] package affy[26]. The raw .cel files are deposited at the NCBI GEO data repository [27] with accession number GSE11121.

Correlation between bimodality scores

We calculated the seven scores for quantifying bimodality for all 22283 genes measured in the Mainz study. First, we analyze the global correlations between the measures. The results show that indeed a clear group structure exists which enables a classification of the measures for our application.

Figure 1 shows dendrograms generated using 1-C as distance matrix, where C is the correlation matrix of the seven scores. Both figures clearly show two groups. On the one side of the dendrogram there are the dip statistic, the bimodality index and -(kurtosis), i.e. the kurtosis score starting with largest values. Overall there is no large Pearson correlation between these scores, the distance is in most cases close to 1. The other group contains the highly correlated measures log-likelihood ratio, WVRS, kurtosis, and the outlier-sum statistic, as well as VRS which has a large distance to the other scores.

In summary, the outlier-sum statistic and the log-likelihood ratio have the highest correlation (Pearson 0.863 and Spearman 0.786). In both dendrograms these scores are thus very close. WVRS and kurtosis also belong to this group, whereas VRS, dip statistic, and bimodality index seem to measure different distributional features, at least according to the overall correlation.

Figure 2 shows smoothed scatterplots for the comparison of selected scores. The strong positive correlation between the log-likelihood ratio and the outlier-sum statistic is clearly visible (Figure 2(a)). The areas of highest local density are at low values for the both scores. The correlation between the log-likelihood ratio and the kurtosis is smaller, but still obvious. However, at higher values of these scores there is a large dispersion; there are some genes with small log-likelihood ratio but very large positive kurtosis (Figure 2(b)). WVRS and the kurtosis are also positively correlated. Some of the genes with largest positive kurtosis have medium WVRS. The highest local density is at WVRS values between 0 and 1.

Distributional shapes of genes with extreme bimodality scores

In order to obtain an impression of types of distributions detected by the different bimodality scores we present histograms of the expression values for the top 6 genes identified by each score, respectively (Figure 3, Figure 4, Figure 5). The most striking differences can be observed between the shapes detected by the correlated measures log-likelihood ratio, WVRS, kurtosis and outlier-sum statistic on the one side and all other scores on the other side. The former scores all identify genes with a main unimodal expression distribution and few additional outliers. The other scores detect mainly genes whose expression values can be grouped into two separated groups.

This difference can be exemplified with the clustering approaches resulting in the VRS and WVRS scores. The histograms for the top genes of VRS all indicate a large group and a smaller group (cf. Figure 3(a)). Except for FOXA1, the smaller group has larger expression values. On the other hand, the expression values of the top genes according to WVRS show a unimodal density with one or two outliers. Typically, most expression values of these genes are low, often in a range that could be regarded as noise, whereas the outliers have high expression values. This different behavior between VRS and WVRS is due to the larger impact of a subset with few samples when using WVRS, since then the variance reduction is calculated by averaging the variances of the two resulting subsets.

Similarly, the density of the genes with the largest positive kurtosis is unimodal with a small outlier group with high expression values (cf. Figure 4(a)), whereas the top 6 genes according to negative kurtosis show bimodal densities with almost equal size groups or a flat density (SCGB1D2). For the top genes according to the likelihood ratio one can see a large group and an outlier group with large variance. Also the outlier-sum statistic detects genes with densities with heavy tails in the higher intensities. For the dip statistic the top genes have clearly visible bimodal densities (Figure 3(c)). The genes with the largest bimodality index show bimodal densities, where the two groups can have unequal sizes (Figure 5(a)). There is not much overlap between the top six genes of all scores, except for the gene HLA-DQA1 with a clear bimodal histogram. The gene is among the top list for the scores detecting bimodal shapes, i.e. VRS, dip, negative kurtosis, and bimodality index, and additionally for WVRS.

In conclusion, all considered scores are applicable for the identification of genes with characteristic distributions. However, the distribution curves differ particularly with respect to the number of patients (samples) constituting the high expression group.

Prognostic relevance of genes with extreme bimodality scores

For the top genes presented in the previous section, we plot in Figure 6, Figure 7 and Figure 8 Kaplan-Meier curves for the top-scoring genes and corresponding p-values of the logrank test when discriminating patients according to low or high expression values. As described in the methodology section we use different methods to assign patients to a low and a high expression group. For VRS and WVRS we use the k-means algorithm, for the kurtosis based approach, the likelihood ratio, and the bimodality index we use the results of the model-based clustering, and for the outlier-sum statistic we use the main group and the outlier group. We argue that the k-means approach is also appropriate for the results of the dip test as we see that the k-means algorithm is able to separate the groups adequately at least for the top ranked genes.

For three of the top genes detected by VRS the survival in the two groups differs considerably. For MAGEA2 and ZIC1, the survival prediction is worse for the group with high expression values, for FOXA1 for the group with low expression values. The only patient in the high value group for UGT1A@ and CTNNA2 drops out at an early point, resulting in a small p-value (p = 0.057, see Figure 6(b)). This patient is the same for both genes. Among the top 6 genes for the dip statistic there are two with p-values smaller than 0.05, namely TFAP2B and ERAP2. For both genes the survival in the group with lower expression values is worse. For two of the 6 genes with largest positive kurtosis, UGT1A@ and MYH7, the survival in the patient subgroups differs significantly. Patients with high expression values have worse prognosis (Figure 7(a)). Among the top 6 genes for negative kurtosis, there is just one (TFAB2B) that splits the samples in two groups with significantly different prognosis. MAGEA2 and MAGEA6 are the two genes identified as bimodal by the likelihood ratio score whose p-values of the logrank test are small. Here, the p-value for MAGEA2 is highly significant (p < 0.001). Among the 6 top-scoring genes based on the bimodality index the two genes FOXA1 and ZIC1 lead to a prognostic split. Among the top 6 genes according to the outlier-sum statistic only VGLL1 generates a p-value smaller than 0.05.

Prognostic relevance of bimodality scores on a genome-wide scale

For the genome-wide evaluation, we first sorted the genes based on the bimodality scores. For VRS and WVRS we generated gene lists with smallest scores at the top. For the dip statistic, the bimodality index, the likelihood ratio and the outlier-sum statistic we generated decreasing lists. For the kurtosis we look at both the increasing as well as the decreasing list. We then selected the 250 genes with smallest p-values for survival differences. These gene lists are different for the different approaches.

For WVRS, among the 1000 top-scoring genes there are just 13 of the 250 prognostic genes (p = 0.333). Apparently there is no enrichment with prognostic genes on top of the list and the global Kolmogorov-Smirnov test can not reject (p = 0.547). However, many of the prognostic genes have large WVRS values (see Additional file 1). These genes have expression distributions with a main unimodal distribution and an additional outlier group with large variance. For the list based on the dip statistic just two of the prognostic genes are among the top 200 and 9 among the top 1000 genes. Here, the null hypothesis of the Kolmogorov-Smirnov test can not be rejected (p = 0.749).

Among the first 1000 genes on the list with decreasing kurtosis values there are 17 highly prognostic genes (p = 0.059), whereas on the increasing list there are 8 genes (p = 0.878). The result of the Kolmogorov-Smirnov test for the decreasing kurtosis list is borderline significant (p = 0.003). The prognostic genes are overrepresented among the first half of the ranked gene list (see Additional file 1). For the likelihood ratio the ranks of the prognostic genes are significantly not uniformly distributed (p < 10^-12), and among the top 1000 genes we find 21 prognostic genes (p = 0.005). For the bimodality index the global enrichment is also highly significant (p < 10^-12). There is a significant enrichment for the top-ranked genes. 22 of the 250 prognostic genes are among the first 1000 genes according to the bimodality index (p = 0.002).

Lastly, for the outlier-sum statistic a highly significant overrepresentation of the prognostic genes among the top-scored genes can be observed (p < 10^-12). Among the first 1000 genes are 26 highly prognostic genes (p < 10^-4).

This shows that on a global level the likelihood ratio, the bimodality index and the outlier-sum statistic work best for the identification of prognostic genes with bimodal density, many prognostic genes have densities with heavy tails. The clustering based approaches VRS is better suited to detect prognostic genes than the dip statistic.

It is important to note that we found that the threshold 250 for determining the prognostic genes is not critical. We tried various other cutoffs. Using 50 prognostic genes, no significant results are obtained. In the range from 100 to 1000 genes we found that the results for VRS, bimodality index, likelihood ratio and the outlier-sum statistic are similar.

Validation on other data sets

To validate our results we repeated our analysis on an other free available data set. The raw .cel files and clinical parameters of the Rotterdam cohort [28, 29] were downloaded from the NCBI GEO data repository with accession number GSE2034 (n = 286). Here, the outlier-sum statistic and the likelihood ratio also have the smallest p-values of the logrank test (see Additional file 2: Supplemental Figure 2). For the bimodality index there is also a significant enrichment (p = 0.01). However, we do not observe a significant overrepresentation of the prognostic genes among the top-scoring genes according to VRS. For the gene list based on decreasing kurtosis the we also find a significant enrichment (p < 10^-7).

To determine whether our results also hold for known prognostic subgroups we used a pooled cohort of 766 patients from different free available data sets: GSE11121 (n = 200), GSE2034 (n = 286), GSE7390 (n = 177) [30] and GSE6532 (n = 103) [31, 32]. We look at three subgroups defined by the expression of the two genes ESR1 and erbB2, namely ESR1+/erbB2- (n = 519), erbB2+ (n = 107) and ESR1-/erbB2- (n = 140). In this analysis the variable for assessing prognostic power of the genes is the distant metastasis free interval (MFI).

In the ESR1+/erbB2- subgroup of the pooled cohort the bimodality index has the smallest p-value of the logrank test (p < 10^-14, see Additional file 3: Supplemental Figure 3). A significant enrichment of the top scoring genes with prognostic genes can also be observed for the likelihood ratio (p < 10^-8) and for the outlier-sum statistic (p < 10^-5). For the kurtosis score the prognostic genes have mainly negative values. There are also many genes with large positive kurtosis but these are not the top genes for positive kurtosis. For VRS the result of the Kolmogorov-Smirnov test is borderline significant (p = 0.018).

In the erbB2+ subgroup of the pooled cohort the Kolmogorov-Smirnov test significantly rejects the null hypothesis for various measures (see Additional file 4: Supplemental Figure 4). Some of the prognostic genes are at the top of the ranked gene lists based on VRS, WVRS, increasing kurtosis, the bimodality index and the likelihood ratio. The smallest p-value can be observed for the negative kurtosis (p < 10^-9). Here, for the outlier-sum statistic no significant result is obtained (p = 0.625). In the ESR1-/erbB2- subgroup the only considerable enrichment can be observed for negative kurtosis (p < 10^-4, see Additional file 5: Supplemental Figure 5).

It is important to note that the pooling of data sets can be problematical since in some cases bimodal expression distributions are an artifact of the pooling of experiments (see Additional file 6: Supplemental Figure 6).

Analysis of established genes with bimodal distribution

Techniques for the identification of bimodally expressed genes should be able to identify genes, whose bimodal distribution is already known. The best established genes with bimodal frequency distributions in breast cancer are the estrogen receptor (ESR1) and erbB2 [1, 2]. Histograms of the expression profiles of these genes in our data set are shown in Figure 9 and clearly indicate subgroups.

ErbB2 also has the smallest rank when using the bimodality index (338) and has small ranks for VRS, the outlier-sum statistic, and the likelihood ratio. The only extreme difference can be observed for WVRS, where erbB2 has rank 21505. The large value of WVRS is due to the large variance in the smaller group. Dip finds groups with almost equal sizes and thus fails to identify the groups with unequal sizes generated by ESR1 and erbB2 in our data set.