Abstract
Hypoxia is a condition of low oxygen tension occurring in the tumor and negatively correlated with the progression of the disease. We studied the gene expression profiles of nine neuroblastoma cell lines grown under hypoxic conditions to define gene signatures that characterize hypoxic neuroblastoma. The - regularization applied to the entire transcriptome identified a single signature of 11 probesets discriminating the hypoxic state. We demonstrate that new hypoxia signatures, with similar discriminatory power, can be generated by a prior knowledge-based filtering in which a much smaller number of probesets, characterizing hypoxia-related biochemical pathways, are analyzed. - regularization identified novel and robust hypoxia signatures within apoptosis, glycolysis, and oxidative phosphorylation Gene Ontology classes. We conclude that the filtering approach overcomes the noisy nature of the microarray data and allows generating robust signatures suitable for biomarker discovery and patients risk assessment in a fraction of computer time.
1. Background
Neuroblastoma is the most common pediatric solid tumor, deriving from immature or precursor cells of the ganglionic lineage of the sympathetic nervous system [1, 2] endowed with remarkable heterogeneity with regard to histology and clinical behavior [3, 4]. The neuroblastoma cell lines derived from the fresh tumors show various degrees of differentiation, chromosomal alterations, and morphology and consequently, a great variability in the gene expression profile. We studied the transcriptional response of neuroblastoma cell lines to hypoxia by microarray analysis [5].
Hypoxia is a condition of low oxygen tension that characterizes many pathological tissues and that is a critical determinant of tumor cell growth, susceptibility to apoptosis, and resistance to radio and chemotherapy [6–8]. The general response to hypoxia involves activation of biochemical pathways leading to alternative ways to generate energy that becomes scant in low oxygen [9]. Hypoxia modulates gene expression through the activation of several transcription factors, among which the hypoxia-inducible transcription factor-1 (HIF-1) [7, 10], and -2 (HIF-2) [11] are the most studied. Rapidly expanding neuroblastoma tumors present areas of hypoxia [12] and it has been reported that HIF-2α expression correlates with poor prognosis [13, 14] suggesting a central role of hypoxia in tumor progression. HIFs transactivate the hypoxia-responsive element (HRE) present in the promoter or enhancer elements of many genes encoding angiogenic, metabolic, and metastatic factors [8, 15, 16]. However, neuroblastoma cell lines respond differently to hypoxia and the nature of the modulated genes depends strongly on the type and genetic makeup of the cell [17]. Furthermore, amplification and/or overexpression of MYCN oncogene, occurring in poor prognosis tumors, influence the transcriptional response to hypoxia of neuroblastoma cell lines [5].
The identification of molecular markers capable of discriminating the hypoxic status of the tumor may result in the discovery of new risk factors for neuroblastoma patients’ stratification and potential targets for tumor therapy. To this end, we were interested in identifying hypoxia signatures that discriminate the hypoxia status of neuroblastoma cell lines. Unsupervised analysis of gene expression profile could not be applied to this system because the overwhelming effect of MYCN amplification on the transcriptome masked the response to hypoxia [5]. The application of a supervised approach represented by regularization with double optimization on microarray data, an embedded feature selection technique proposed by Zou and Hastie [18] and studied by De Mol et al. [19], identified 11 probesets capable of reliably subdividing hypoxic and normoxic cell lines [5]. These results raise the question as to whether this signature is the only possible outcome of the regularization algorithm, and hence the only source of neuroblastoma hypoxia markers, or whether additional signatures, with similar characteristics of performances and robustness can be derived from the experimental data set. Hypoxia induces massive transcriptional changes in the cell [20–22] and it is possible that additional signatures may be found by the algorithm under appropriate conditions.
The regularization algorithm has to deal with heterogeneity of the response of each cell line and with the background noise that is enhanced by the high dimensionality of the system composed by a low number of samples ( in this work) relative to the large number of the expression values for each sample (, 613). The scenario is a common issue in signal processing and machine learning [23, 24]. Furthermore, the strong response of each cell line to alteration of the genetic makeup (e.g., MYCN rearrangement) tends to overcame and mask the response to hypoxia. Here, we explore the possibility that feature selection algorithm may generate new hypoxia signatures following prior knowledge-based data filtering techniques as a preprocessing step to feature selection.
Most dimensionality reduction methods, such as PCA and other unsupervised learning methods [25], rely only on the input data and may be driven by strong concurrent signals which are unrelated with, and somehow hide, the problem under study. Alternative strategies of data filtering are based on some form of prior knowledge of the biology of the system. The information collected by Gene Ontology (GO), a project having the aim of classifying gene products in terms of their associated biological processes, cellular and molecular components [26] can help identifying the pathways related to hypoxia and restricts the analysis to smaller sets of data.
In this paper, we demonstrate that regularization applied separately to probesets representing genes belonging to selected GO ontologies has the capability to generate robust hypoxia signatures, equivalent to that generated by the whole data set yielding biologically relevant information in a fraction of computer time.
2. Materials and Methods
2.1. Microarray Experiments
Microarray data were downloaded from the Gene Expression Omnibus public repository at National Center for Biotechnology Information database (accession number GSE15583). These data represent the gene expression profile of nine neuroblastoma cell lines cultured under normoxic (20% O2) or hypoxic (1% O2) conditions for 18 hours as detailed in [5], to obtain a total of 18 samples. Affymetrix HG-U133 Plus 2.0 GeneChip (Affymetrix, SantaClara, CA) were used for this study. Gene expressions were then extracted from CEL files and normalized using the Robust Multichip Average (RMA) method [27] by running an R script using the Bioconductor [28] package affy.
Comparative analysis of hypoxic relative to normoxic expression profiles for each cell line was conducted on GeneSpring 7.3 software (Agilent Technologies). Gene expression data were normalized using “per chip normalization” and “per gene normalization” algorithms implemented in GeneSpring. First, each signal was normalized based upon the median signal in that chip (“per chip normalization”). We then performed a median centering using “per gene normalization” function by which each normalized value is corrected based upon the median of the measurements for that gene in all samples. Finally, only genes that were modulated by at least 2-fold between hypoxic and normoxic cells were considered differentially expressed.
2.2. Gene Ontology
The biological groups were obtained from the literature and they were divided into three main categories depending on the biological characteristics of our experimental system (see Table 2): (i) hypoxia related groups [9, 17, 29]; (ii) MYCN related groups [30–32]; (iii) neuroblastoma related groups [31, 33, 34]. The selected functional groups were then filtered to avoid overlapped or duplicated categories and were defined according to predetermined pathways and functional categories annotated by the Gene Ontology project [26].
2.3. Supervised Methods for Gene Selection: Regularization
The core of our approach is the regularization originally presented in [18] and further developed and studied in [35, 36]. To describe such method we first fix some notation in the learning framework. Assume we are given a collection of n examples/subjects, each represented by a p-dimensional vector x of gene expressions. Each sample is associated with a binary label Y, assigning it to a class (e.g., patient or control). The dataset is therefore represented by a matrix X, where and Y is the n-dimensional labels vector. We consider a linear model . Note that 1, …, is a vector of weight coefficients and each probeset is associated to one coefficient. A classification rule can be then defined taking sign . If β is sparse, that is some of its entries are zero, then some genes will not contribute in building the estimator. The estimator defined by regularization solves the following optimization problem: where the least square error is penalized with the and norm of the coefficient vector. The least square term ensures fitting of the data whereas adding the two penalties allows avoiding overfitting. The relative weight of the two terms is controlled by the parameter ε. The role of the two penalties is different, the term (sum of absolute values) enforces the solution to be sparse while the term (sum of the squares) preserves correlation among the genes. This approach guarantees consistency of the estimator [19] and enforces the sparsity of the solution by the term, while preserving correlation among input variables with the term. Differently to [18] we follow the approach proposed in [36], where the solution , computed through the simple iterative soft thresholding, is followed by regularized least squares (RLSs) to estimate the classifier on the selected features. The parameter in the regularization is fixed a priori and governs the amount of correlation. By tuning in (0, ) we obtain a one-parameter family of solutions which are all equivalent in terms of prediction accuracy, but differ on the degree of correlation among the selected features. In practice, has an upper bound, , such that for selection does not change, because all correlated features were already selected with . By setting , the maximal value, the maximal gene list, which is correlation aware, is obtained. Conversely, the minimal list is obtained for values of equal to or lower than 1.
The training for selection and classification requires the choice of the regularization parameters for both regularization and RLS denoted with and , respectively. Hence, statistical significance and model selection is performed within double-selection bias-free cross-validation loops (see [37] for details). The classification performance of the system is measured by the leave-one-out error that is the percentage of misclassified samples. In other words, leave-one-out error is equal to one minus accuracy. In order to assess a common list of probesets, it is necessary to choose an appropriate criterion [38]. We based ours on the frequency, that is, we decided to promote as relevant variables the most stable probesets across the lists. The complete validation framework comprising the regularization is implemented in MATLAB code.
2.4. Correlation Analysis
The correlation among the probesets selected by the algorithm was performed as previously described in [39]. Briefly, we build blocks of correlated probesets using a variation of well-known agglomerative clustering techniques based on Pearson distance. We first examine the minimal list, which genes are clustered via hierarchical clustering with correlation distance and average linkage. Since no objective algorithm, other than heuristics, is available for establishing the number of clusters, for each GO class the cut of the hierarchical graph determining the number of clusters is chosen following visual examination of the graph. In particular we set the cut at 0.75 of the maximum linkage value in the dendrogram. For each GO class the cut of the hierarchical graph determining the number of clusters is chosen following visual examination of the graph. Each probeset in the maximal list is then assigned to the cluster which average correlation with the given probeset is the highest. In this way we populate the clusters built from the minimal list with all the probesets coming from the maximal list. The correlation analysis was performed using MATLAB Statistic Toolbox.
2.5. HRE Analysis
We mapped the HRE elements in the promoter regions of the genes represented in the Affymetrix HG-U133 Plus 2.0 GeneChip. We downloaded the annotation file for the HG-U133 Plus 2.0 from NetAffx Analysis Center (http://www.affymetrix.com/) and the dataset was restricted to the known mRNA sequences listed in the Ensembl database V56 [40]. The regulatory regions were retrieved from Ensembl database using Ensembl Perl APIs. We operationally defined as ‘‘promoter’’ the first 2,000 base pairs upstream the transcription initiation site and generated a dataset containing the promoters of the genes coding for the mRNAs spotted on the chip. The HRE matrix has been obtained from 69 experimental validated human HRE sequences [41] with MatDefine tool (Genomatix Software GmbH). HRE consensus elements [(GCAT) (CGTA) (GCAT) (TGCA) (AG) (CGTG) (CGTA) (GCAT) (GCTA) (CGTA)] were searched in the promoter sequences with MatInspector software (Genomatix) with core and optimized matrix similarity. About 33% of the promoters contain at least one HRE consensus element. χ2 was used to evaluate the significance of the HRE frequency in the promoter regions of genes belonging to the different signatures. was considered significant
3. Results and Discussion
We studied nine neuroblastoma cell lines [2] heterogeneous with respect to MYCN amplification and morphology (Table 1). The cell lines were cultured under normoxic and hypoxic conditions for 18 hours and the total RNA was tested for gene expression profiling using the Affymetrix HG-U133 Plus 2.0 platform. The response to hypoxia of each individual cell line was first analyzed by measuring the fold change as the ratio of the expression level between hypoxic and normoxic samples. We found that the response of each neuroblastoma cell line to hypoxia is characterized by a high number of modulated genes ranging from 855 to 1609 for the upregulated and from 758 to 1317 for the downregulated probesets (Table 1). However, the modulated genes changed from cell line to cell line (data not shown) and only the application of a strong feature selection technique, represented by the regularization, allowed to identify a single signature of 11 probesets (All-chip signature) discriminating the normoxic and the hypoxic status [5].
The large amount of hypoxia-modulated genes suggests that additional hypoxia signatures may be identified if we reduce the background noise of the system. To this end, we applied a data filtering strategy based on prior knowledge. We restricted our analysis to the genes known to be involved in the hypoxic response on the bases of our reading of the literature and comprised in the biological processes according to the Gene Ontology (GO) classification [13]. The selection of the GO classes was based on the reports of hypoxia modulated genes without attempting to distinguish the various cell types under investigation. 13 biological processes that are involved in hypoxia response were selected (Table 2). We reasoned that this approach would restrict the analysis to the probesets that have a high impact on the hypoxic response while potentially eliminating the noisy features. To explore the potential interference from MYCN status in the classification process, we selected and tested 7 biological processes involved in MYCN activity (Table 2). Finally, we selected a third group of GO processes related to the neuroblastoma biology as a control. For each of the 38 classes shown in Table 2, the algorithm selected a list of hypoxia discriminating probesets and calculated the corresponding classification leave-one-out error. The output of the regularization algorithm depends on the parameter that governs the amount of correlation allowed among the probesets. We set , the maximal value, to obtain the most comprehensive signature maximizing the number of correlated probesets to be included in the output [5].
The validation has been performed by leave-one-out cross-validation on the 18 samples. The 18 cross-validation loops produced 18 lists of probesets. Then, a unique list is obtained as the union of the probesets included in the 18 lists, with a frequency score calculated as the frequency of each probeset in the 18 lists generated by the cross validation loops. Stable probesets were defined as those characterized by a frequency score equal to, or higher than, 50% as previously reported in [5]. The use of cross validation allows the selection protocol to generate an unbiased and objective output [42] beyond the theoretical results that guarantee the robustness of the core algorithm [19]. The discriminatory power of the probeset lists is represented by the classification performance. A leave-one-out error of 20% was chosen as the cutoff level for the classification performance. The leave-one-out error of the All-chip signature is 17% [5].
The only classes characterized by a list of selected probesets capable of generating a leave one-out error lower than the 20% cutoff were apoptosis (17%), glycolysis (11%), and oxidative phosphorylation (11%) (Table 2), all of them belonging to the hypoxia biological group. These results demonstrate that, within each of the above classes, there is a list of probesets capable of discriminating the condition of the cell lines thereby defining three new neuroblastoma hypoxia signatures, named apoptosis signature, glycolisis signature, and oxidative phosphorylation signature. As expected, there were no GO classes belonging to the MYCN or neuroblastoma biological groups that generated hypoxia signatures, supporting the validity of our choice of hypoxia-related GO functional classes. Although MYCN represents a strong signal that drives major transcriptome difference in neuroblastoma cell lines [5], our results show that there are no enough discriminatory genes in the MYCN-related processes. These results demonstrate that the feature selection method applied is capable of revealing the differences occurring among hypoxic and normoxic neuroblastoma cell lines by filtering out strong competing signals, such as MYCN amplification status.
The list of the probesets comprising the 11 probesets (All-chip signature) [5] and the newly identified signatures is shown in Table 3 and consists of 10 probesets for apoptosis signature, 3 for glycolysis-signature, and 32 for the oxidative phosphorylation signature. The new signatures highlight 41 probesets that were not previously included in the All-chip signature and contribute to the discrimination of the hypoxic status. Furthermore, the 32 probesets of the oxidative phosphorylation signature does not overlap with the All-chip signature, demonstrating that the increased resolution generated by data filtering allows the identification of previously discarded relevant GO processes. The hypoxia signatures present in the literature show different sizes and gene composition [9, 43–46]. Since different cell types respond heterogeneously to hypoxia by modulating different set of genes, we decided to compare our results with the published hypoxic gene signatures obtained from neuroblastoma cell lines [47] (Table 4). In order to make the comparison feasible, the probesets constituting our signatures have been collapsed to gene symbol. The overlapping genes are underlined in bold in Table 4. While important differences among the signatures exist, the comparison highlights a general consistency. In fact, an overlap is present in All-chip (3/8 genes), apoptosis (2/4 genes), and glycolysis (2/2 genes) signatures. Interestingly, there is no overlap (0/24 genes) among the results published by Jögi et al. [47] and the oxidative phosphorylation signature.
About 33% of the genes spotted on the chip present a HRE sequence in the promoter region. We investigated whether there was enrichment in HRE containing promoter in the genes composing our signatures. We found that all the signatures are significantly enriched () in genes containing HRE (Table 4). In particular, all the genes included in All-chip, apoptosis, and glycolysis signatures contain at least one HRE, while HRE containing genes constitute 91% of the oxidative phosphorylation signature. These results support the idea that our signatures are associated with the hypoxia status.
The whole signature, rather than individual genes, is important for discriminating the hypoxic status. For example, VEGF is a gene whose expression is strongly related to hypoxia [45] and is part of the apoptosis and angiogenesis classes, both of which are part of the hypoxia biological group. However, the contribution of VEGF probesets is not sufficient to reach the discriminatory power required to generate a significant signature out of the angiogenesis class as opposed to the apoptosis class.
The strong discriminatory power of the signatures can be visualized by a 3-dimensional representation of the probesets projected on their 3 principal components. algorithm produces a multigene model but the multidimensional representation can be well approximated by the tridimensional picture when the number of probesets is not too large. Figure 1 depicts the separation of probesets belonging to the glycolysis and shows that the two classes of normoxic and hypoxic cell lines are clearly separated in the multidimensional space.
In conclusion, we demonstrate that, upon data reduction the algorithm can identify new hypoxia signatures that have equivalent discriminatory power relative to that obtained by the analysis of the whole transcriptome. From the computational stand point, this process allows to reduce the computer time by approximately 10 times, from days to hours with an average machine, facilitating the analysis of the data. Finally, it is important to highlight the possibility of applying our method to different experimental settings by choosing appropriate selection of GO processes.
Our prior knowledge-based method produces nested lists of relevant probesets but does not highlight the correlation among them [39], and it should be completed by a postprocessing step depicting the correlation structure. The correlation within the oxidative phosphorylation signature is shown in Figure 2. We computed a distance matrix based on the expression values of the probesets and subdivided it into 8 modules by hierarchical clustering. These modules represent subgroups of correlated probesets that are positively or negatively associated to the hypoxic status. This information is important to pick the correct probesets in order to assess the expression of these markers in the in vivo setting. Furthermore, these data lend themselves to the tuning of the parameter that is part of the algorithm.
The output of the regularization algorithm depends on the free parameter that governs the amount of correlation allowed among the probesets and selects the amount of probesets to be included in the signature. By setting , the maximal value, we can obtain a comprehensive signature more descriptive of the biology of the system. By setting , we can obtain an equally discriminating signature with fewer genes thereby more effective in identifying critical biomarkers for diagnostic applications [5]. We analyzed the effects of tuning on the oxidative phosphorylation signature. The results are shown in Table 5, where the probesets selected by regularization with both and for oxidative phosphorylation are listed. The results demonstrated that the reduction in is associated with a smaller signature (from 32 to 16 probesets) as expected by the fact that correlated probesets tend to be discarded.
4. Conclusions
The identification of signatures discriminating the hypoxic status of the tumor cell may be important for our understanding of the biology of neuroblastoma tumors and for the stratification of the patients in risk groups. One way to generate a robust and reliable hypoxic signature is the application of a supervised approach represented by regularization that generates an 11 probesets signature discriminating the hypoxic status of our panel of nine neuroblastoma cell lines.
Here, we demonstrate that feature selection algorithm generates new and robust hypoxia signatures following prior knowledge-based data filtering techniques as a preprocessing to feature selection. These new signatures have the same discriminatory power as that generated by the whole data set and yield biologically relevant information in a fraction of computer time.
The data filtering is based upon the use of the prior information contained in GO and the literature, and it allows restricting the analysis to smaller data sets. This process filters out not only many noisy probesets but also the probesets selected from the all-chip analysis whose strong relation with hypoxia hid some weaker but important genes. regularization algorithm following data filtering selects probesets that were not the first chosen when all the probesets were considered. The prior knowledge utilized in setting up the filter, comes from the current literature from which we derived the molecular pathways that are important for the response of the cell to the hypoxic environment. These pathways were gathered in the hypoxia biological group. Interestingly, the new signatures were found only in this group and not in other collections of GO pathways like those related to the effects of the MYCN oncogene or to the neuroblastoma biology. In general, the identification of the GO classes related to the phenomenon under investigation may be an empirical, but effective way to target the potential source of signatures to be fed to the regularization. We speculate that this approach could be used to address questions that go beyond the hypoxic status and may find signatures characterizing other pathophysiological situations provided that there is a relevant cellular model and there are sufficient insights in the underlying molecular mechanisms.
The nested structure of the selected gene lists allows the choice of the desired level of complexity, which is the magnitude of signature, maintaining all the information extracted from the data. For example, the minimal list may be preferable when interested in finding biomarkers to be used on large-scale diagnostic tests due to potential constrains on time, cost, and resources.
Finally, working on a limited number of probesets has a major impact on the computational time required for the analysis that changes from days to hours, thereby allowing more leeway to the study of the dataset.
Acknowledgment
This paper was supported by Fondazione Italiana per la Lotta al Neuroblastoma, the Italian Association for Cancer Research (AIRC), the Italian Health Ministry, the EU Integrated Project Health-e-Child IST-2004-027749 and Compagnia di San Paolo Project 4998- ID/CV 2007.0887.