Unsupervised clustering of gene expression data points at hypoxia as possible trigger for metabolic syndrome

Differences in gene expression between the FOREL clusters

Next, we identified genes that were most differentially expressed between cluster 3 (metabolically healthy) and cluster 6 (high prevalence of DM2, low insulin sensitivity). To select such genes we applied caGEDA tools [21] on unabbreviated list of genes (see Methods). We found that clusters occupying the extremes ends of the line connecting the centroids varied with respect to the expression of genes encoding ribosomal proteins. Other significant function categories included actin, myosin, hemo- and myoglobins, enzymes involved in glucose and fructose metabolism, lipid metabolism and genes encoding proteins involved in oxidative phosphorylation (see Additional file 2 in supporting materials [22]).

We also performed functional annotation of these same genes at the extremes of clusters 3 and 6 using EASE [23]. We found several categories of genes that were significantly overrepresented as compared to their expected occurrence (see Additional file 3 in supporting materials [22]). The most prominent functional categories included ribosomal genes, genes involved in skeletal muscle contraction, actin, cytoskeleton and protein biosynthesis.

The dataset produced at the East Carolina University [24] has a different purpose, it is focused on skeletal muscle fatty acid metabolism in connection with obesity, but utilized the Affymetrix GeneChip (U133). This chip is different from those used in experiments of Mootha et al. (U95) and Patti et al. (HS6800). Given that obesity and insulin resistance are clearly linked [25, 26] it is interesting to compare this study with our re-analyses of the Mootha et al. dataset. The 24 samples were collected from 3 groups of patients with normal weight (BMI 23.8 ± 0.58), obese (BMI 30.9 ± 0.81) and morbidly obese (BMI 53.8 ± 3.5). Surprisingly, the FOREL clustering results from the Hulver dataset presented on Figure 3 are highly similar to the clustering results from the Mootha et al. dataset. Four large clusters, produced by FOREL analysis form a line in the gene expression space, one end of the line is occupied by a cluster composed almost entirely by the most metabolically sound specimens, with the opposite cluster again containing the least sound samples. Analysis of the genes most differentially expressed between the most distant clusters along this line reveals a remarkably similar list of genes to the first analysis (see Additional file 5 in supporting materials [22]).

The other data set used in this study (Patti et al.) is much smaller; it contains 15 samples of which 10 represent normal and 5 represent DM2. Again using the FOREL clustering process, we found that all samples fell into 4 clusters. The largest cluster (#4) contains 5 samples, of which 3 are DM2, Cluster #2 is entirely composed of samples from normoglycemic subjects, clusters #1 and #3 include mostly normoglycemic subjects but each includes one DM2 subject. It is reasonable to suppose that clusters #2 and #4 represent the most contrast groups of samples in the data and the line that stretches between centroids of these clusters would be analogous to the trend we observe in the Mootha et al. data set. The results of our re-analysis of the Patti et al. data can be found in supporting materials [22] (see Additional file 4, supplementary figure 7 in Additional file 1).

When comparing these three studies, the samples are absolutely independent and differ in many ways (age, gender, type of microarray used, feature extraction algorithm etc.), yet the genes that are most differentially expressed between the most and least metabolically sound are again remarkably similar to those found in our re-analyses of the Mootha and Hulver datasets. (see supporting materials [22], Additional files 2, 3, 4, 5) These discriminant genes include GAPDH, myoglobins, titin, FHL1 and a score of ribosomal proteins. Although genes involved in oxidative phosphorylation are among those differentially expressed in all three datasets, they are not at the top of the list when ranked by the fold change between clusters. Thus these analyses identified common genes and functional categories that were unanticipated based on the earlier analyses.

After these efforts to cluster the subjects, we next clustered the genes independent of the subjects. Cluster analysis of the genes throughout the data set revealed a number of highly correlated clusters. Again using the FOREL technique with the Mootha dataset, we found several interesting clusters of genes. Based on a functional annotation analysis [23] of these genes, we found that one of the clusters contains a group of genes related to lipid transport and metabolism. However, the largest cluster includes 221 genes and among these genes are the same ribosomal, glucose and fructose metabolism genes that discriminated between the most and least metabolically sound subjects.

Natural classification of subjects based on their transcriptome separates subjects into categories recognizable by their clinical characteristics

We next asked whether the natural classification of the molecular phenotypes using the FOREL analysis corresponded to the clinical phenotype as measured by BMI, insulin sensitivity, etc. The clinical data, provided by Mootha et al. was used to compare the 'metabolic soundness' of the patients across cluster membership from the FOREL analysis (Figure 4). Patients from clusters 2, 5, 6 and 7 have a higher body mass index (BMI) and lower insulin sensitivity (M-value) when compared to the 'normal' cluster. This visual analysis was confirmed by ANOVA (supplementary data, supplementary Figures 1–4). From this and the other results (Figures 1 and 2) the subjects represented in cluster 3 are in exceptional physical form for their age. They form a tight group with a normal BMI and normal insulin sensitivity. Interestingly, maximal aerobic capacity (VO₂max) does not clearly discriminate this group from the least metabolically sound cluster with a high proportion of DM2. The two singletons, both NGT have either normal (#8) or almost normal (#1) BMI and high insulin sensitivity (M-value). A few NGT samples are also found in clusters 2, 4, 5, even less are found in clusters 6, and 7 that are dominated by DM2 and IGT samples. Thus in spite of the high BMI and low insulin sensitivity these individuals manage to remain normo-glycemic. This is consistent with the natural history of DM2 where beta cell failure is a late event with a variable onset but invariably preceded by insulin resistance. We speculate that the metabolically unsound individuals are at risk for the development of DM2 as beta cell apoptosis ensues and islet neogenesis fails to compensate [27]. These 'natural' clusters coincided with the classic markers of diabetes risk, namely obesity and insulin resistance, but did not follow the accepted clinical diagnosis of DM2 as defined by the presence or absence of hyperglycemia. Unfortunately, there are other important parameters missing in all three data sets, such as time of onset of hyperglycemia and estimation of physical activity. Since all samples were collected from the patients of general hospitals we can speculate that the level of physical activity and its distribution among sample donors is typical for the urban population of US East Coast and on average not very high. It is possible that molecular mechanisms of insulin resistance discussed below are significantly influenced by the level of patients' physical activity.

Relationship of the FOREL classifications with clinical characteristics

In the next step of analysis we superimposed the clusters in the space of first principal component. For all three datasets, the centroids of the clusters were situated along a single line or plane. Clusters with the most metabolically sound subjects with NGT and the least sound subjects (with DM2) occupy the extreme ends of this line, while most IGT samples are found in the mixed clusters between the extremes. In the Mootha dataset, the healthiest group is represented by a single cluster, which has 3 non-diabetic and 1 IGT subjects. On the opposite end, two dominant clusters contain only one non-diabetic sample each, while all other members are either diabetic or insulin resistant. This picture suggests a simple interpretation: instead of categorical classes of normal and diabetic samples we observe a continuous trend. On one end of the trend are the samples from the most metabolically sound individuals with a low prevalence of DM2. Prevalence of DM2 increases along the line, but can be mitigated by a number of factors. Many of these factors are beyond the scope of this research and probably relate to β-cell apoptosis and/or failure of beta cell neogenesis/replication [27, 29, 30]. As expected, unsupervised clustering revealed a few distinct classes. However, we did not anticipate that these clusters would be aligned along a single principal component. Rather than multiple 'subtypes' these analyses suggest a single underlying disorder. Importantly, these natural classes differ in average clinical characteristics and in the number of diabetic and non-diabetic members.

Discriminant genes

Once the direction of the trend is established it is relatively easy to select and annotate the genes that are most differentially expressed between the extremities of the trend. The list of such genes does include those involved in the energy metabolism and oxidative phosphorylation, as reported by Mootha et al. and Patti et al. The expression level of such genes changes dramatically between clusters on the extreme ends of the trend. Oxidative phosphorylation genes possibly play an important role in the etiology of DM2 although this role is not decisive and is insufficient for discrimination between classes. This also explains the difficulties Mootha et al. had selecting the genes that discriminate between DM2 and NGT samples.

These discriminant genes fall into several interesting gene ontogeny categories. First, we observed an increase in expression of genes encoding proteins required for glycolytic rather than oxidative metabolism such as GAPDH. This is consistent with prior data demonstrating a reduction in oxidative capacity of skeletal muscle from insulin resistant diabetics [31].

Second, we discovered that myoglobin, lactate dehydrogenase A, GAPDH, aldolase A, and pyruvate kinase, genes that are classically upregulated when oxygen tension falls [32] were upregulated in the least metabolically sound samples. Human skeletal muscle adaptation to altitude hypoxia has been well characterized and includes a seemingly paradoxical reduction in oxidative capacity, an increase in myoglobin, and a shift to glycolytic metabolism [32]. All of these are dominant features of the gene expression data reported herein. Altitude hypoxia reduces the overall mitochondrial volume by 20%, but with a preferential reduction in sub-sarcolemmal mitochondria (-43%) [33] compared to a modest change in the intra-myofibrilar mitochondria (-13%), a finding strikingly similar to the mitochondrial defects observed in type 2 diabetes [6]. Reduced oxidative enzymes are also observed in S. Cerevisiae as hypoxia ensues [34]. Hypoxia/anoxia also shifts substrate metabolism from aerobic to glycolytic in diverse species from S. Cerevisiae to man [32, 35]. Multiple transcription factors are involved in the transcriptional response to hypoxia in S. Cerevisiae including Rox1 [34], Hap1–4, Sut1, UPC2, and Mox 1–3, (reviewed in [35–37]) and in man HIF-1a [38]. The role of these transcription factors in DM2 is unknown and warrants exploration. If hypoxia is involved in the observed molecular and cellular phenotype it does not cause a compensatory effect on the capillary density. The capillary density and maximum oxygen capacity (VO2max) as reported in the supplementary materials (Supplementary Table 6) by Mootha et al. shows no significant difference across FOREL clusters. This again, is consistent with the literature on structural adaptations to hypoxia where no change in capillary density are seen unless intense exercise was employed as part of the experimental paradigm [33, 39, 40]. This finding does not necessarily present a contradiction: a shift from oxidative to glycolytic pathway metabolism in skeletal muscle tissue may reduce demand for oxygen as an alternative to capillary growth. Alternately, insulin is capable of activating the bundle of genes that are similar to those activated by HIF-1a. In vivo, insulin upregulates aldolase A, and GAPDH, both of which are classic targets of HIF-1a [38, 41]. There is some controversy whether these effects are dependent on changes HIF-1 mRNA [42], occur through the stabilization of HIF-1a protein or through an unknown HIF response element binding protein [43]. It is established, however, that insulin upregulation of these genes requires phosphoinositol 3-kinase [42, 44]. In summary, the unsupervised FOREL analysis identified upregulation of a cassette of genes strikingly similar to hypoxia inducible genes in the obese, insulin resistant state.

Lastly, we discovered a set of genes encoding ribosomal proteins separated metabolically sound from unsound subjects with DM2 and IGT. On the surface, this is an anomalous finding. However, microarray studies of petite, mitochondria deficient (ρ-) yeast data demonstrate an upregulation of genes encoding ribosomal proteins as compared to mitochondria replete (ρ+) yeast [45]. Similarly, as yeast become anoxic and shift from an oxidative to ethanol producing metabolism (the diauxic shift) many of these same genes encoding ribosomal proteins are up/downregulated. This is congruous with prior data showing a downregulation of oxidative phosphorylation genes and disrupted mitochondria in DM2 [31]. An alternate explanation for the increased expression of ribosomal and other genes associated with cell growth reflects alterations upstream of the mammalian Target of Rapamycin (mTOR) pathway. mTOR is a major regulator of metabolic processes in cell. It has been reported to stimulate ribosomal biogenesis and repress nutrient turnover [46, 47]. Prior studies have implicated mTOR in nutrient sensing and diabetes [46, 48]. Acute insulin infusion upregulates the expression of many ribosomal genes [49], providing an alternate pathway for the observed changes in ribosomal gene expression.

In contrast to the subtle changes in gene expression of the oxidative phosphorylation genes identified by Mootha et al., we found several genes that discriminanted lean healthy people and individuals with DM2 or IGT. For example, Aldolase A was 3.3 times overexpressed in obesity/insulin resistance and GAPDH 4 times higher. Similarly, hemoglobin was upregulated in obesity/insulin resistance; by 13, 12.1 and 11.8 fold for subunits alpha 2, beta and alpha 1 respectively. The novel transcription factor FHL1 was ~3.5 fold upregulated. The second and third data sets (Hulver et al. and Patti et al.) although limited by sample size and without complete clinical phenotypes, produced almost identical cluster layout with a single group of the most metabolically sound subjects and several mixed groups with differing prevalence of DM2. Similarly, the genes most differentially expressed between the two most and least metabolically sound clusters were remarkably similar across the three datasets. This suggests all three datasets identified similar underlying pathophysiological mechanisms leading to DM2.

Further study will be required to determine which factor plays the primary role, whether impaired oxidative phosphorylation causes a shift to glycolytic metabolism or oxidative phosphorylation is repressed as a result of unknown upstream transcriptional events [50]. Regardless, the observed upregulation of ribosomal gene expression might be a marker of the risk to develop type 2 diabetes.

Data sets

We downloaded the data published by Mootha et al [5], Patti et al [4] and Hulver et al. [24]. The Mootha et al. data set consists of 43 age-matched samples representing 17 normal, 8 with impaired glucose tolerance and 18 with DM2. Affymetrix U133 biochips were used to measure expression of over 22,000 genes. This data, along with the clinical information, was downloaded from the author's data portal at the Broad Institute [51]. The other data set [22] consists of 10 normal non-diabetic skeletal muscle tissue samples and 5 samples from the patients diagnosed with DM2. Five replicate microarrays are present in this dataset. In our analysis these samples are always found in one cluster close to each other. We have excluded these samples from the final data set.

Clustering algorithm

As a starting point for the algorithm development we took the heuristic concept proposed by Zagoruiko et al. [52, 53], which included the original idea of a limited radius hypersphere, moving stepwise to the mass center of captured objects. This idea represents a departure from widely used k-means and other hypersphere-based algorithms. This algorithm has been also extensively discussed in a reference book for applied statistics in economics [54]. Since their first introduction algorithms of the FOREL family have been widely applied in taxonomic analysis of biomedical data, pattern recognition in geology and image analysis [55–58]. The original algorithm underwent a significant development the Pennington Center to accommodate for extreme dimensionality of gene expression data. The results of this development and case study in classification of molecular subtypes of lung cancer has been recently published [22]. Our algorithm is based on dynamic amalgamation of objects (for example, expression profiles) in vicinity of an artificially introduced object (FORmal ELement). The vicinity is defined by equal distance from a point in all directions by selected inter-object distance metric (such as Euclidean, correlation, binary, etc.). Although theoretically the vicinity could be defined as any geometrical shape around the given point, only hyperspheric vicinity has been implemented and used in this study. FOREL clustering is based on the perception of the data set O as:

$O={\displaystyle \underset{i=1}{\overset{k}{\cup }}O(i);}$

where O(i) is a cluster of n _i elements. Clusters are extracted from the general population in order of their statistical fitness (see Cluster validation). This perception is fundamentally different from the popular k-means algorithm, which shares certain similarity with FOREL, but in k-means concept the whole data set is a sum of distinct classes rather then a union. FOREL clusters can partially or completely overlap in space or even share the same centroid, but can be separated as long as they differ in other statistical characteristics, for example density. In a nutshell, a white and a yolk of an egg would be separate classes by FOREL, while inseparable by k-means. Other hypersphere-based algorithms such as k-means imply Gaussian distribution of objects (phenotypes) in clusters [59]. FOREL is more flexible and does not require such assumptions. FOREL effectively combines the best features of k-means and hierarchical clustering approaches for the price of increased computation complexity. However, the performance of our C++ implementation is acceptable; up to a few hundreds of microarrays can be clustered on a Pentium IV PC within one hour. The algorithm starts with positioning of a hypersphere with a radius R₀ and a center C₀ in a certain coordinate, which can be one of the objects or a centroid of pre-defined cluster or any other point of interest in the expression space. Position of the "formal element" element is calculated as a center of mass of all objects, for which the distance ρ _i (C _i ) ≤ R _i . After the mass center of all captured objects is calculated, its center is moved to the new position. If new objects are found within the radius from the new position, they are added to the provisional cluster and the mass center is recalculated. This process is repeated until the no more objects can be added on the current step of the algorithm and the hypersphere stops.

Cluster validation

The version of FOREL developed at Pennington Biomedical Research Center (PBRC) consists of alternating steps of cluster isolation and cluster validation. Each completed walks of a hypersphere with R _n and a center C _m produces a provisional cluster O(R _n C _m ), which is temporarily stored. We perform an exhaustive coverage of the data trying each element of the original data set as potential starting point. For each starting point we perform series of clustering steps with different hypersphere radius, ranging between minimum R = min(D(C _i , C _j )) + μ and maximum R = min(D(C _i , C _j )) - μ. Here D(C _i , C _j ) is a distance (for example, Euclidean) between any two objects in the data set and μ is a margin, introduced to reduce computation time. The step of R increment is also a parameter. The resulting provisional clusters are fuzzy subsets, each captured by a hypersphere with specific radius as it moves gradually from the starting to the resting point. The validity of the provisional cluster can be verified by a statistical utility measure based on density, variance, sum of inter-cluster distances, etc (see [60] for review). If the cluster meets the selection criteria, it is removed from the original data set and the process reiterated until no more statistically valid (according to the chosen metric) unclassified objects left or the best provisional cluster does not satisfy the minimal fitness requirement. The PBRC implementation accommodates a few different metrics for computational cluster validation, but only one density-based metric has been applied in this study:

This metric is a reasonable compromise between precision and performance, which has proven to be effective in analysis of microarray data [22]. The "brute force" approach to computational cluster validation provides more reliable results compared to re-sampling, but results in considerably longer execution times. Depending on the validation metric applied and the parameter settings complete clustering of a large data set, such as Mootha et al. data can take up to a few hours on a typical computer. The demand for computational power is significantly mitigated by effective C++ implementation and generally affordable, considering the time required to collect such data.

Software implementation

The implementation developed by A. Ptitsyn [22] and available for anonymous download [61] employs a complete test of every object as a possible cluster seed or hypersphere starting point. By default the current version of the program implements an iterative solution: all possible radii are tested with a certain step within a limited range. The step (or precision) is derived from the analysis of variation of distances within the whole data set. The range is defined by the minimal and maximal distance found within the whole dataset. These are extreme values, with a radius less then a minimal distance, the algorithms can produce only singletons, and on the other hand with radius equal to the maximal distance, all objects are guaranteed to fall into one large cluster. The best radius is one that produces a provisional cluster with the best quality. Cutting percentile margins from the possible radius range can reduce the computation demands of the program. By default 20% of the range is cut from both minimal and maximal radius values.

The implementation of FOREL clustering algorithm made at the PBRC runs on the standard PC under MS Windows (win32 console application) or Linux. FOREL execution time is practically unaffected by the dimensionality of the data, but can be sensitive to the number of objects (samples) in the set. Large data sets of over 1000 microarrays would require a high-performance computer.

Identification of informative genes

We have selected genes differentially expressed between most distant clusters using J5 algorithm [21] with Jackknife option set to 4. We used the software implemented at the University of Pittsburg Cancer Institute with access through the Web interface [62].

Functional Annotation and Interpretation of the results

The consequent hyper-cluster analysis was performed in order to reveal relations between clusters and their juxtaposition in the gene expression space. Centroids of the original FOREL clusters have been used to produce a hierarchical classification (Euclidean distance, UPGMA) with the help of Spotfire Decisionsite for Molecuar Genetics (Spotfire Inc.). Functional annotation of the discriminating genes has been performed with EASE and DAVID tools [23]. Only a small portion of the results can be included within this publication. The additional supplementary materials can be found at [22].