In this manuscript, we address the above-mentioned issues by developing a novel approach to identify the signature of low oxygen tension (hypoxia) in a set of neuroblastoma cell lines. Oxygen is essential for aerobic metabolism in all mammalian cells. To maintain function and homeostasis, cells have to be able to sense and respond to inadequate oxygen levels. The O₂ levels within the neoplastic lesion are an important factor in determining the tumor phenotype [1] and hypoxia is associated with metastatic spread, resistance to radio- and chemotherapy and poor prognosis [1-3]. The cellular response to hypoxia is caused by changes in gene expression [4-6] through the activation of several transcription factors among which the hypoxia-inducible transcription factor-1α (HIF-1α) [1,7], and -2α (HIF-2α) [8] are those taken as indicators of a hypoxic status of the cell. HIFs transactivate the hypoxia-responsive element (HRE) present in the promoter or enhancer elements of many genes encoding angiogenic, metabolic and metastatic factors [3,9,10]. Although hypoxia responses are thought to be evolutionarily conserved in all mammalian cells [11,12] not every cell responds to hypoxia in an identical fashion. Although certain biochemical pathways are common hypoxia targets, the specific genes modulated by hypoxia within each pathways will depend heavily on the nature, type and genetic makeup of the responding cell [5,6,13]. In other words, hypoxia-induced common biochemical pathways may utilize different genes depending on the cell type.

Neuroblastoma is the most common pediatric solid tumor, deriving from immature or precursor cells of the ganglionic lineage of the sympathetic nervous system (SNS) [14,15]. Neuroblastoma shows notable heterogeneity, with regard to both histology and clinical behavior [16]. The outcome of the disease ranges from rapid progression and poor clinical outcome, to spontaneous regression into benign ganglioneuroma [17]. The heterogeneity of neuroblastoma (NB) cells is found also in the cell lines derived from the fresh tumors which manifest various degree of differentiation and chromosomal alteration. For example, the amplification and/or expression of MYCN oncogene is a relatively frequent event, that is indicative of poor prognosis in fresh tumors and is present in several cell lines which share an aggressive behavior [18]. Recent data of microarray analysis confirm the existence of different patterns of gene expression profile among different NB cell lines [13]. The heterogeneity of the NB cell transcriptome complicates the identification of specific gene expression signatures associated to defined biological responses such as environmental stimulation that, albeit biologically very important, may be overshadowed by major genetic alterations as those caused by oncogenes which impact on several aspects of cell physiology. This problem is of major concern when several different NB cell lines have to be compared in in vitro studies.

The problem of identifying a gene signature, namely a significant group of variables, is aggravated by the typical high dimensional nature of the data. Complexity grows even more when the heterogeneity of the cells must be factored in. Several feature selection techniques have been proposed to deal with these problems (for review see [19]). The number of data available for a single study is usually small with respect to the number of variables, and it is crucial to adopt sound methodologies and strict experimental protocols to ensure statistical robustness [20]. Cross validation loops are valid approaches to avoid selection bias [21] and to separate training and test phases. The standard categorization proposed by Blum et al. [22] groups variable selection techniques in three main classes: filters, wrappers and embedded. Filters [23-25] are mostly based on ranking criteria where the features are ordered and then selected or discarded according to a fixed threshold. These methods are broadly employed due to their simplicity and fast computation, despite the lack of guarantee that the selection is optimal with respect to the class discrimination. In wrapper methods [26-28] the relevance of a feature subset is determined according to prediction performance of the learning algorithm itself, though variable selection and training are two separate processes. In contrast, embedded methods [29-34] have the advantage of incorporating feature selection within the construction of the classifier or regression model, i.e. as part of the training phase. We applied the embedded feature selection technique l₁-l₂ regularization with double optimization to the analysis of gene expression profile. This technique is based on the optimization presented by Zou et al. [35]. Theoretical studies [36] and empirical experiments [37,38] showed that such technique fulfills all the desirable properties of a variable selection algorithm. Indeed, the use of regularization allows performing embedded feature selection in the supervised learning framework, since the particular type of penalty used in l₁-l₂ regularization forces the classifier or the regression model to depend on a small number of selected features. Another asset of l₁-l₂ regularization is that it is multivariate by design since its solution is a classification or regression model that takes into account the combined effect of multiple features, and the set of relevant features is selected while looking at all the features at the same time. A strong advantage of l₁-l₂ regularization over other embedded methods is also its ability to take into account correlation among variables. In other words, when one variable is considered relevant to the problem, its correlated variables are considered relevant as well. While most feature selection techniques are based on heuristics, l₁-l₂ regularization is asymptotically consistent from the statistical viewpoint, i.e. theoretical results [36] guarantee that the best possible estimator is found as the number of training samples increases. Finally, the use of the double optimization allows to identify the relevant genes and to provide accurate discrimination. This approach was successfully applied in different contexts ranging from computer vision [37] to computational biology [38].

The first attempts to identify the hypoxia signature relied on three different unsupervised clustering analyses. These approaches detected major differences among the transcriptome of the cell lines driven by the characteristics of the cell lines themselves of which the cascade of events triggered by MYCN expression was a major component. However, the disturbance generated by these transcriptional patterns was such that the detection of more subtle changes induced by hypoxia was completely masked. This conclusion is supported by the results obtained with the supervised univariate analysis t-test which was able to identify a strong response associated to MYCN expression and, to a lesser extent to MYCN amplification, but not the hypoxia dependent response. The heterogeneity of neuroblastoma and neuroblastoma cell lines has been previously observed [43].

The impossibility to obtain a hypoxia signature by unsupervised approaches prompted us to consider different algorithms based on a robust supervised variable selection technique, capable of detecting the hypoxia-induced transcriptome even in the presence of the disturbance of a strong competing signal. The l₁-l₂ regularization allowed us to build a powerful discriminative rule and to define a signature of probesets also taking into account the presence of variables (probesets) correlated (collinear) with each other. The use of cross validation allows the selection protocol to generate an unbiased and objective output [21] beyond the theoretical results that guarantee the robustness of the core algorithm [36]. The strong discriminative power is proven by the 17% classification error that is a very low value when dealing with 18 samples and nearly 50 thousands variables.

We adopted a validation framework based on a double loop of leave-one-out cross-validation in order to extract unbiased estimates of the classification error. The outer loop produces 18 lists of relevant variables, from which we extract a common list by setting a frequency threshold. It can be appreciated from visual inspection of the frequency distribution that, for values lower than 30%, a large number of probesets is included, which are extremely unstable. For frequency above 70%, the number of selected probesets slowly decreases, and a plateau is present between 30 and 70%. Therefore, we set our frequency threshold to 50% that is the intermediate value of such a frequency plateau. The correlation parameter ε can be potentially tuned between 0 and +∞ in order to extract lists of probesets with different correlation degree. However, values of ε equal to or smaller than 1 provide the same minimal list which comprises 9 probesets. This list is minimal in that it does not include correlated probesets, and it can be viewed as the smallest set of variables needed to predict the hypoxic status without any prior information. Conversely, by increasing the correlation parameter ε we are able to expand the list to 11 probesets, which is obtained for ε ≥ 100. Since we are interested in all genes involved in the hypoxic condition, we define such a correlation-aware list as our hypoxia signature in neuroblastoma cell lines. This signature identifies also hypoxic status in astrocyte cell lines. We estimated a similar low cross validation error demonstrating that its application goes beyond the neuroblastoma lineage and providing an additional proof of its discriminatory power on out-of-set data. However, the neuroblastoma hypoxia signature is less efficient in discriminating hypoxic dendritic cells indicating the existence of a limited spectrum of hypoxic cell types that can be identified by this signature. Dendritic cells are terminally differentiated mononuclear phagocytes, biologically very far from the other cell types, deriving from hematopoietic precursors, short lived and programmed to serve as immunomodulatory cells [44]. This conclusion supports the concept of the heterogeneity of the response to hypoxia in different cell types [13,42,43]. The enrichment of the 11 probeset signature in the hypoxic phenotype of neuroblastoma cells and astrocytes provides biological validation of our approach by establishing a clear link between our signature and the hypoxic microenvironment.

The 11 probesets represents 8 genes all of which are known from the literature to be modulated by hypoxia in different cell types and to be part of key biological processes associated with the response to hypoxia, indicating, once more, the biological roots of our signature. ALDOC and PDK1 belong to the glycolytic pathway and are known to be up regulated by hypoxia in neuroblastoma [45]. Potentiating of the oxygen-independent glycolytic pathway to comply with the energy demand, is one of the major cellular response to hypoxia [46]. The energetic balance in cell metabolism has to be controlled by different mechanisms. For example, AK3L1 is known to be modulated by hypoxia and catalyzes the interconversion of adenine nucleotides playing an important role in cellular energy homeostasis in mitochondria [47]. DDIT4, induced by hypoxic stimulus, is an inhibitor of the mTOR signalling pathway, that results in inhibition of protein synthesis which, in turn, may affect the cellular tolerance to hypoxia by promoting energy homeostasis [48]. VEGF, a direct target of HIF-1, is secreted by a large variety of different hypoxic cells and promotes angiogenesis thereby favouring tumor growth and metastasis [4]. BNIP3 and E2IG5 are two genes promoting hypoxia-induced apoptosis observed mainly at very low oxygen concentrations [49]. E2IG5 is localized to mitochondria and facilitates apoptotic cell death via permeability transition, cytochrome c release, and caspase 9 activation [50]. Hypoxia is also known to increase histone H3 methylation through histone methyltransferase G9a [51]. WDR5B encodes for a protein that is the core component of histone methylation complexes, which are essential for histone H3 methylation. Thus, hypoxia might regulate chromatin organization and gene transcription by modulating WDR5B. Finally, the GenBank entry W57613 is part of the signature it is associated with a transcribed hypothetical protein FLJ11267 and was not previously known to be induced by hypoxia.

The novelty of our work is to introduce a rigorous and robust feature selection technique that can be exported to other experimental models and that is able to identify discriminative genes even in an adverse setting where the cell lines express great heterogeneity. The hypoxia signatures present in the literature show different sizes and composition [4,42,52-54]. The MSigDB [55] represents a valuable source of gene sets associated to the response to hypoxia. A first attempt to discuss our 11 selected probesets within the contest of 9 hypoxia signatures contained in the MSigDB is based on the analysis of the overlap among signatures (Table ​(Table4).4). One limitation of this comparison is that it must be based on gene names rather than probesets, because of the heterogeneity of the platforms. All the genes of our hypoxia signature, but one, are represented at least once in the 9 signatures. DDIT4 is the only hypoxia inducible gene that is included only in our signature. This comparison lends further support to the conclusion that the l₁-l₂ algorithm selected biologically relevant genes that have been included in other hypoxia signatures. There are at least three major reasons for the variability among the hypoxia signatures. The first is the diversity of the cell types as shown by Chi et al. [13] and Mense et al. [42] and supported by the observations on the heterogeneity among neuroblastoma cell lines by Fredlund et al. [43] and ourselves in this paper. Each cell responds to hypoxia on the bases of its own genetic make up, epigenetic constrains and differentiation stage. In fact, we show that our signature does not apply to dendritic cells, that are biologically very different from astrocytes and neuroblastoma. The second reason is the difference in the experimental setting and gene expression platforms. The need to collapse the microarray probes to gene names for comparisons is a direct consequence of this problem. The third, and more important issue, is the criterion used for assembling the signature. The majority of the signatures described so far, are based upon the representation of hypoxia associated biochemical pathways or the inclusion of differentially expressed genes rather than the essentiality and the discriminating power that we have chosen. Having defined the hypoxia signature as the minimal number of probesets capable of distinguishing normoxic and hypoxic gene expression profiles, our list is relatively short, not specific for a biochemical pathway, not relying on prior biological knowledge, but endowed with high discriminating power.

The classification performance of our signature is evident in representations that indicate as first approximation the up-regulation of the signature in hypoxic condition. However, the multidimensional visualization is needed to fully appreciate the strong discriminative power of our signature because it takes into account its multivariate nature. In fact, when projecting over the individual probesets of the signature, the two classes are only approximately separated, since they appear either partially overlapping or very close. Indeed, since the l₁-l₂ regularization is a multivariate method, there is no need to expect a single probeset to have perfect discriminatory power on the classes, but one has to take into account the 11-dimensional model. While the normoxic cell lines are highly grouped and close to low expression values, the hypoxic lines are well spread over the multidimensional space, though well separated by the normoxic ones. Again, this behavior can be detected only by means of a multivariate analysis, since the analysis of individually regulated genes allows detecting only those probesets which multidimensional representation would see the hypoxic cell lines very well lumped together.