Classification of Human Astrocytic Gliomas on the Basis of Gene Expression: A Correlated Group of Genes with Angiogenic Activity Emerges As a Strong Predictor of Subtypes¹²

Because of their diffusely infiltrating behavior, LGA⁴ (WHO grade II) cannot be resected completely and will usually recur. At relapse, progression to anaplastic astrocytoma (WHO grade III) or glioblastoma multiforme, the most malignant form of gliomas (WHO grade IV), is common. Most glioblastoma arise de novo without evidence of a less malignant precursor lesion and are termed PrGBM. Glioblastoma evolving from a previous lower grade astrocytoma are defined as ScGBM. Although PrGBM are indistinguishable from ScGBM by histology, the two types of tumors exhibit distinct genetic alterations and occur in different age groups. The mean age for PrGBM is ∼55 years, whereas ScGBM typically occur in younger patients (<45 years). Thus, PrGBM and ScGBM can be considered as two different diseases (1), despite a similarly grim outcome with a median survival of <1 year after diagnosis and no effective therapy. PrGBM are characterized by amplification/rearrangement and overexpression of the EGFR gene (in 40 and 60% of the patients, respectively) often in association with deletion of the INK4a/p14ARF gene locus (2). The hallmarks for ScGBM are TP53 mutations (60%) and overexpression of PDGF and PDGF receptor (3, 4). Development of mouse glioma models and developmental neurobiology have allowed for recent advances in the understanding of the molecular bases of these two distinct genetic pathways and their implication for tumor initiation and progression as well as the cell of origin (5). However, a substantial number of glial tumors cannot be characterized by either of the two pathways depicted above, suggesting additional not yet recognized pathogenetic pathways. The current knowledge of tumor genetics does not allow identifying clinically relevant factors predictive for outcome or response to therapy. More detailed knowledge of underlying mechanisms and their relevance for the cancer process will allow treating cancer specifically by targeting deregulated pathways, leading to rational design of future treatment modalities tailored according to the biology of the individual tumors (6).

An essential initial step toward this goal is the establishment of a taxonomy of tumors on the basis of their gene expression profiles. A search for alternative pathways must be based on identification of genes whose expression differs significantly between the various tumor classes. In particular, it is important to look for a group of (possibly) coregulated genes, some of which share some known biological function, and whose expression differentiates tumor classes. Identification of such groups leads to better understanding of the biological processes that underlie the distinction between the tumors and may provide clues for the roles of such genes in initiation and progression of cancer.

We aimed at identifying expression profiles that differentiate three groups of astrocytic glioma: (a) LGA; (b) their respective ScGBM; and (c) PrGBM. In a novel gene selection approach, we combined supervised statistical analysis with CTWC (7, 8, 9), an unsupervised method, to identify correlated groups of genes that distinguish between the various tumor subtypes. Here, we demonstrate that gliomas can be separated according to their gene expression profiles, with PrGBM exhibiting a much higher variation of expression profiles than LGA. A cluster of correlated genes was identified that separates PrGBM from the other tumors and contains genes that are known to be involved in angiogenesis, such as VEGF, but also IGFBP2, whose implication in angiogenesis is novel. By using expression data of the most informative separating gene clusters, we were able to construct an almost perfect tumor classifier. Our main findings, based on cDNA arrays, were validated on an independent set of glioma and by other methods, respectively.

Tumor biopsies, obtained from patients who underwent surgery at University Hospitals in Lausanne or Geneva or the Cantonal Hospital in Fribourg (Switzerland), were shock frozen and stored at −70°C. The use of biopsies and respective clinical data have been approved by the local ethics committee and the respective federal agency. The tumors were diagnosed according to the WHO classification 2000 (1). Twenty-one biopsies originated from 20 patients enrolled in a prospective pilot trial for newly diagnosed glioblastoma (10). From the tumor bank, an additional 32 gliomas from 22 patients were included, comprising 24 LGA and 8 respective high-grade recurrent tumors (comprising two astrocytoma WHO grade III and 6 grade IV). These tumors have been analyzed previously for TP53 mutations (11). The tumor samples were organized into two data sets according to their date of analysis: (a) subsequently used as a training set, comprising 14 PrGBM, 5 ScGBM, and 12 LGA; and (b) used as validation set, 4 PrGBM, 4 ScGBM, and 12 LGA. A summary of the information on the individual tumors, including TP53 status, age, and gender of the patients, and their organization into data sets, is available in Table 1. Normal brain tissue for RNA isolation was obtained from a lobectomy after brain edema, and additional samples of human normal brain total RNA were obtained from Clontech (human total RNA panel IV).

TAs have been constructed from archived paraffin blocs at the University Hospital in Lausanne (1984–2000) as described (12). All cases have been reviewed according to the WHO classification 2000 (1) by the neuropathologist (R. C. J.). The GBM-TA comprises 190 GBM, and the non-GBM glioma-TA includes 158 gliomas WHO I–III of different glioma subtypes.

Culture conditions for the glioblastoma cell line U87, LN229, LNZ308, and the TP53-inducible glioblastoma cell line 2024, which is derived from LNZ308 after introducing the Tet-On System, have been described before (13). Expression of wild-type TP53 was induced with 2 μg/ml doxocycline in cell line 2024 during 24 h before RNA isolation (2024+). For anoxia treatment, cells were cultured for 20 h in an anaerobic culture incubator (O₂ <1%; Scholzen, Microbiology Syst AG) filled with a mixture of N₂, H₂, and CO₂.

Before RNA isolation, a section of the frozen tumor biopsy was reevaluated after H&E staining by the neuropathologist (R. C. J.) to estimate the proportion of solid tumor, contaminating normal tissue, or infiltration zone. Pieces comprising >30% normal tissue were excluded. The setting of thresholds on the tumor fraction is treated differently by various groups, an issue that has been reviewed recently by Ramaswamy and Golub (14), who suggest to use “tumor cell-enriched material.” Fifty to 100 mg of frozen tumor tissue were homogenized in TRIzol solution (Life Technologies). RNA phase was purified in saturated phenol solution [60% phenol, 15% glycerol, and 0.1 m sodium acetate (pH 4)]. The quality of the RNA was evaluated on agarose gel. RNA from cell lines was isolated similarly.

In situ hybridization was performed according to the protocol supplied by Roche (Roche Applied Science) using cRNA probes labeled with digoxigenin during the in vitro transcription reaction. The plasmid pcDNA3 containing the 1400-bp human IGFBP2 cDNA (a kind gift from S. Babajko; Ref. 15) was linearized with EcoRI or HindIII and transcribed with T7 or Sp6 RNA polymerase, respectively, to obtain sense or antisense probes. The plasmid pBluescript-KS-M13+ containing the 650-bp human VEGF165 cDNA (pBspt-KS-VEGF165, a kind gift from K. Plate; Ref. 16) was linearized with EcoRI or BamHI and transcribed with T7 or T3 RNA polymerase, respectively, to obtain sense and antisense probes. Probes were further reduced to an average size of 100 bp by limited alkaline hydrolysis. Northern blot analysis was performed as described before using 10 μg of total RNA (13). The membrane was sequentially hybridized to the plasmid-derived probe for VEGF165 (EcoRI/BamHI-fragment of pBspt-KS-VEGF165) and the PCR-derived probes for IGFBP2, 3, and 5 (respective primer sequences provided by Clontech). Probes were radioactively labeled using the random primed DNA labeling kit (Boehringer) using [α-³²P]dCTP (3000 Ci/mmol, Amersham). Expression was quantified by phosphorimager (Fuji, BAS 1000).

Atlas Human Cancer 1.2 Array membranes (Clontech) were used for all experiments described. These nylon filters are spotted with 1185 genes, including reference genes, and 1176 genes related to cancer. Three to 5 μg of DNase-treated total RNA were used to prepare a labeled first strand cDNA using the Clontech kit, basically as recommended. Briefly, RNA in a volume of 2 μl was mixed with 1 μl of specific CDS primers and 1 μl of RNasin (Promega; 40 units/ml) and denatured at 70°C. Subsequently, 1 μl of Superscript II reverse transcriptase (Life Technologies, Inc.; 200 units/ml) and 3.5–5 μl of α-³²P-dATP (3000 Ci/mmol; Amersham) were added, and the reaction was performed at 48°C for 30 min. The probe was purified with Clontech Atlas Nucleospin extraction kit following the manufacturer’s recommendations. Before use, the membranes were boiled in 0.5% SDS solution for stripping and also at first use. Membranes were exposed to an imaging plate (BAS MS 2040; Fuji) for 1–8 days. The Atlas Cancer Arrays were used three times.

The tumors were screened for TP53 mutations using the yeast functional assay as described (17, 18), followed by direct sequencing, if the test was positive (Microsynth, Balgach, Switzerland).

Immunohistochemical determination for EGFR (Novo Castra; NCL-EGFR; dilution 1:40) and tenascin R [Santa Cruz Biotechnology; Tenascin-R (N20) sc-9874; dilution 1:1000] on paraffin sections was performed according to standard procedures using a high temperature epitope retrieval technique in citrate buffer (pH 6.0; pressure cooker, 3 min). Semiquantitative evaluation was performed independently by two researchers.

Expression was quantified by phosphorimager (Fuji BAS1000) and analyzed with Atlasimage 1.5 software (Clontech). After background subtraction signals were normalized using the “sum method” comparing the sum of the intensity of all genes in the experiment to the sum calculated from the reference sample, yielding the coefficient of normalization, “c.” The reference used in all experiments represents the calculated average of five independent expression profiles derived from “normal” brain. The ratio (R) was calculated for every gene as follows: r = [(s − s0) × c + K]/[(r − r0) + K], where K = b × [s0 + (c × r0)]; s, gene intensity in sample; s0, background in sample; c, coefficient of normalization obtained using the sum method; r, gene intensity in reference experiment; r0, background in reference experiment; and b = 2. For s ≫ s0 and r ≫ r0, this formula generates normalized ratios R shrunk toward 1 for low expressed genes in both experiments, sample and reference, respectively. Furthermore, it circumvents the problem of losing valuable information if the experiment or the reference display no expression for a given gene (division by 0).

This method projects the data points from the high dimension in which it is embedded to a low (two or three) dimensional space, in a way that best preserves their relative distances. Euclidian distances in the space of logarithms of the ratios R were used. MDS was performed with the implementation for the classical metric scaling (also known as principal coordinate analysis; Ref. 19) available in the mva package for R available online at the Comprehensive R Archive Network.⁵

A variation filter was applied; only those genes were kept, for which the ratio of the maximal and minimal R values (obtained for the 36 experiments) exceeded 2. For each gene (row), the log of the ratio R was mean-centered (subtracting the average) and normalized. Euclidean distances, measured between all pairs of genes and between all pairs of tumors, served as the input to our clustering procedure. CTWC has been described elsewhere; see Getz et al. (7) for full details and comparisons with other methods; (8) for its applications to leukemia, colon, and breast cancer data analysis; and (9) describing the use of the publicly available CTWC–Server.⁵ CTWC starts with clustering all genes on the basis of the data from all tumors and clustering all tumors, using data from all genes. In both resulting dendrograms, we identify stable (i.e., statistically significant) gene and sample clusters; these are denoted, respectively, as GX or SY (where X,Y are running indices). Note that G1 represents the set of all genes and S1 that of all samples. In the second iterative step, each stable gene cluster GX is used to characterize and cluster the members of every sample cluster SY and vice versa.

To use CTWC, we must be able to identify stable clusters. One of the few algorithms that provide a stability index to each cluster is Super Paramagnetic Clustering, a physics-based method that has been described in full detail in Blatt et al. (20). We register a cluster C as stable only if it exceeds a certain size and when Stab(C), its stability index defined in terms of the range of resolution parameters, T, through which cluster C “lives” (21), exceeds a certain threshold.

For binary class comparisons, two standard statistical tests were used: (a) the two-sample t test (with unknown but equal variances); and (b) the Wilcoxon rank-sum test. To address contamination with false positive genes associated with multiple comparisons, we use the method of Benjamini and Hochberg (22) that bounds the average FDR. The outcome of this method is a list of differentiating genes; the expected fraction of false positive genes is at most q.

The class discrimination power of the selected sets of genes was validated by training a k-nearest neighbor (k-NN) classifier on one set of tumor tissue samples and testing the class prediction on an independent validation set. Computations were performed with the class package for R and with k = 3. Euclidian distance in the space of the ratios R was used, classification was decided by majority vote, with ties broken at random. When ties occurred, classification was repeated 100 times, and Ps were averaged. Ps for the significance of the deviation from independence between true and predicted labels were obtained with Fisher’s exact test (ctest package for R) and with the alternative hypothesis set to “greater” in the 2 by 2 case.

To characterize and classify gliomas by their gene expression profiles, RNA isolated from frozen gliomas and three glioblastoma cell lines was analyzed using cDNA arrays comprising 1185 genes.

The configuration of the 51 astrocytic gliomas in Euclidian space of overall gene expression, as visualized in Fig. 1 using MDS, clearly suggests that gene expression profiles contain information discriminating the three classes of astrocytic gliomas. The LGA and ScGBM show a higher degree of spatial intermingling, whereas the best separation is between LGA and PrGBM. This correctly reflects other biological characteristics of the ScGBM, in that they progress from LGA but share the malignancy grade (WHO IV) with PrGBM and are indistinguishable from them histopathologically. Interestingly, the pairwise distances in space of all genes (1185) are highest among PrGBM (mean, 10.06; SD, 2.1) and shortest between LGA (mean, 6.09; SD, 1.18), whereas ScGBM are intermediate (mean, 7.27; SD, 1.18). Using 358 genes that pass a variation filter (see CTWC, below) yields a nearly similar MDS picture (data not shown). Thus, MDS of gene expression suggests more biological heterogeneity of PrGBM, without firm evidence for obvious subclasses, reflecting well the heterogeneous morphological and clinical characteristics of glioblastoma multiforme.

We applied t test and Wilcoxon’s rank-sum test to look for genes that differentiate two known classes, using all 1185 genes and the training set of 31 astrocytic gliomas. The threshold for the FDR was set at q = 0.05, i.e., the expected number of false positives was kept at 5%. The results for three two-way comparisons are summarized in Table 2. The lists of differentiating genes determined by the two tests are very similar. Because of the low number of the ScGBM in the training set, comparisons involving ScGBM, ScGBM versus PrGBM, and ScGBM versus LGA yielded only very few separating genes (1 and 3, respectively). Full results of these supervised analyses, including a list of the respective genes selected by these binary comparisons, can be found in Table S1 of the supplementary information.

Our efforts to find genes correlated with TP53 status failed, whether we used all tumors or only subsets thereof. Neither did we find any gene cluster separating groups of samples on the basis of their TP53 status (see TP53 status of biopsies in Table 1). Our results are in line with a previous report on a series of LGA, which also did not find evidence for genes correlated with TP53 status (23). This result may not be surprising in consideration of the fact that the TP53 pathway has been found to be inactivated, in >70% of all astrocytic gliomas, regardless of tumor grade, by mutational alteration of TP53 or alternative mechanisms, such as overexpression of mdm2, a negative regulator of TP53, or inactivation of p14ARF, a negative regulator of mdm2, by deletion or promoter methylation of the gene (24).

The data derived from 36 experiments (Table 1; of 33 gliomas and three experiments with glioblastoma cell lines) were subjected to the variation filter described in “Materials and Methods,” reducing the number of genes (rows) to 358. This filtering reduces the “noise” generated by genes that do not vary significantly over the samples. Changing the filtering parameters from 2 to 1.5 (or 3) changed the number of genes that passed to 601 (or 190), but the changes induced in the observed gene clusters were comparable with the fluctuations observed in repeated runs (with 358 genes). The first operation, G1(S1), clustered the set of all 358 genes (G1), using the data from the set of all 36 experiments (S1); the second, S1(G1), clustered the tumors S1 using their expression profiles measured for all of the genes G1. The two resulting dendrograms are shown in Fig. 2,A together with the correspondingly reordered expression matrix. G1(S1) identified 15 stable gene clusters (Stability > 8, corresponding to P < 0.01) that are indicated in the gene dendrogram of Fig. 2 A and denoted G2 to G16. The operation S1(G1) yields a dendrogram with two stable clusters, S2 and S3. Details of the clusters obtained in the CTWC can be viewed and searched online at the respective Web page.⁵

We used all of the 15 stable gene clusters found above, one by one, to recluster S1. Whenever a stable partition of the samples was found, we checked whether it divided the tumors according to some biologically relevant attribute. In particular, we looked for gene clusters that partition the tumors into GBM versus LGA.

In the first level of CTWC, we identified G5 as a stable cluster of nine genes (Fig. 2, A and B). This cluster comprises hallmarks of angiogenesis, such as VEGF, FLT1 (fms-related tyrosine kinase 1; also called VEGFR1), and pleiotrophin (25). Some of the other genes have also been related to some aspects of angiogenesis before (Table 3), namely, IGFBP2 (insulin-like growth factor binding protein 2), which has been suggested to be activated by hypoxia-inducible factor-1, although in an indirect manner because of the absence of a defined hypoxia-response element (26). When the expression levels of only these genes were used to characterize the tumors, a large and stable cluster, S11, of 21 tumors emerged (Fig. 2,B). This cluster contained all of the 12 LGA and 4 of 5 ScGBM. Of the remaining 5 samples of S11, three were cell lines, whose expression profiles were consistently different. Hence, the expression levels of G5 gave rise to a nearly perfect separation of clinically defined tumor entities, namely PrGBM from LGA and ScGBM. These genes were significantly up-regulated in PrGBM (Fig. 2 B) as compared with LGA and ScGBM.

The gene cluster G12 contains 13 genes that are almost exclusively related to the immune system and inflammation (Table 3). Four of which are known to be regulated by γ-IFN (FCGR1A, CYBB, IFITM1, and MX1), and some of the Fc γ receptors and DAP12 represent markers for natural killer cells, although their expression is not exclusive. The dendrogram obtained by clustering the tumors using these genes is presented in Fig. 2 C. Of its two large stable clusters, S39 contains nine tumors, eight of which are PrGBM; the expression level of the G12 genes in these tumors is high. The other large cluster, S38, contains all but one of the LGA and ScGBM samples, together with four PrGBM. Three of these four tumors are from patients whose survival was relatively long (>18 months; Ref. 10). Hence, the immune system-related genes of G12 appear to have lower expression levels in LGA and ScGBM and in those PrGBM that have longer survival.

In a second step, we used the results of the supervised analysis described above to identify those clusters in the dendrogram yielded by G1(S1) that are rich in discriminating genes. Of the 205 genes that were found to differentiate PrGBM from ScGBM and LGA (Table 2; Table S1 of supplementary information), 91 passed the variation filter and were included in G1. As depicted in Fig. 3,A, most separating genes belonged to one of four marked clusters. The same four clusters were identified also as rich in genes that partition samples according to the other supervised binary comparisons mentioned above. Two of the identified clusters were indeed G5 (six of nine genes are among the 91 separating ones) and G12 (5 of 13). However, two additional clusters were identified as rich in separating genes; these are marked in Fig. 3,A as nodes N505 (28 of 29) and N566 (6 of 7). These sets of genes (Table 4) were not selected as stable clusters by CTWC because they were either too small or they “leaked” too many genes as the control parameter T was increased. However, as clearly shown in the distance matrix (Fig. 3 B), the genes of both of these nodes had relatively highly correlated expression profiles.

These additional clusters contain some interesting genes; members of the Wnt- and Notch-pathway, genes involved in adhesion/migration, and apoptosis-related genes in node N505 (Fig. 3,C; Table 4) and in node N566 genes like TGFβ2 and IGFBP3 that have been reported before to be overexpressed in high-grade gliomas (Fig. 3,D; Table 4).

Next, we used these two additional clusters to partition all of the samples S1. As shown in Fig. 3,C, node N505 gives rise to a nearly perfect separation of the samples into three clusters: (a) “Non-PrGBM”; of 17 samples, comprising 12 of 12 LGA, 3 of 5 ScGBM, and 1 PrGBM and the OAIII, both of which (1308 and 1357) joined the non-PrGBM cluster in S1(G5); (b) mainly PrGBM, 13 of 14 PrGBM, the RecGBM (1497), 2 of 5 ScGBM, one of which (1342) had joined the PrGBM in S1(G5), 1 cell line (U87); and (c) a cluster with the remaining cell line under two conditions (wild-type TP53 on and off; 2024+, 2024−). Hence, node N505 gives a clear separation into LGA versus PrGBM, with the ScGBM split evenly. As we also see in Fig. 3 D, node N566 identifies a cluster of 20 samples comprising 12 of 12 LGA; thus, this node also separates the samples largely into LGA versus PrGBM, with the ScGBM distributed evenly.

A gene expression based tumor classifier was built to validate our findings. A successful test of such a classifier based on the selected genes and applied to an independent validation set of 20 samples indicates that the separation found previously is robust and was not caused by random fluctuations. For each separating cluster of genes (G5, G12, node N505, and node N566) and their combination, we generated a k-nearest neighbor classifier (k = 3) and tested its ability to perform class partitions. In all cases, the separation power was found to be significant (Table 5). In addition, we tested whether a three-way separation of the tumor types can be performed based on the combination of the four gene clusters. This gave rise to only one to two errors out of the 20 gliomas, thus validating our claim that the tumor types have different gene expression profiles that can be used to distinguish between them (P = 2.49E-06). When using the gene sets separately, the distinction of PrGBM versus ScGBM was the least successful. Most classification errors disappeared when all four selected clusters were used. This suggests that distinct sets discriminate between different aspects of the classes. In fact, genes of the same cluster were generally highly correlated across the samples of the same class as a result of the selection procedure, whereas genes of different clusters were not. Although the best performance was obtained by combining all genes, the G5 set alone provided already a good separation power in all cases, including three-way classification where only one to two (8%) errors were made (Table 5).

Genes found most significantly up-regulated in these glioma as compared with normal brain comprise genes involved in the cell cycle, angiogenesis, migration, and immune response, such as vimentin, fibronectin, EGFR, VEGF, CDK4, PDGFRA, SPARC, HLA-G, and HLA-DRA. Similar lists have been published before (23, 27, 28, 29). An overview of the most differentially expressed genes is available in Table S2 of supplementary information.

To further validate our findings obtained by gene expression profiling, we confirmed a gene known to be overexpressed, EGFR, on respective paraffin sections by immunohistochemistry. Furthermore, we validated differential expression of TNR on a larger panel of tumors using immunohistochemistry on tissue arrays.

A good correlation between EGFR overexpression as measured on cDNA arrays and by immunohistochemical detection of the EGFR on respective paraffin sections (no/low expression versus high expression) was observed using a subset of 26 gliomas. The mean expression levels were 265 (SD 101) for samples with no/low expression as determined by immunohistochemistry and 1052 (SD 1369), respectively, for biopsies with high immunopositivity (unpaired t test, P = 0.035). Please note that expression of EGFR as determined by immunohistochemistry can be focal.

Relative overexpression of TNR in LGA as compared with PrGBM and ScGBM was indicated by supervised analysis and its presence in the separating gene cluster N505. We chose to validate TNR, a multifunctional extracellular matrix glycoprotein that is of interest for further evaluations because of its role in adhesion, migration, and differentiation during central nervous system morphogenesis, in particular, differentiation of oligodendrocytes (30). This member of the tenascin family is exclusively expressed in the central nervous system in particular by oligodendrocytes. TNR expression was validated on a large panel of gliomas by immunohistochemistry on our GBM-TA and non-GBM glioma-TA, respectively. The criteria for semiquantitative analysis were no/low expression as opposed to high expression. Twenty-nine “pure” LGA were available for analysis on the array applying the same pathology criteria as for the LGA included for gene expression profiling. LGA exhibited significantly higher expression of TNR (83%, 24 of 29) than GBM (63%, 84 of 133 available for analysis; Fisher’s exact test, P = 0.03).

Most of the genes in G5 are involved in angiogenesis in a broader sense as detailed above. Within this cluster, IGFBP2 expression displayed the highest correlation (r = 0.85) with VEGF expression as shown in Fig. 2,B. Both genes exhibited significant differential expression between any two tumor subclasses as visualized in a box-plot in Fig. 4, using all 51 astrocytic gliomas of the test and the validation set. To further understand the biology of IGFBP2, we addressed the question whether the two genes are involved in the same process. Therefore, in situ hybridization was performed on serial frozen sections of glioblastomas for IGFBP2 and VEGF, respectively. Expression of IFGBP2 and VEGF was most prominent in pseudopalisading cells surrounding tumor necrosis of GBM (Fig. 5; VEGF in situ hybridization; data not shown). Thus, the IGFBP2 expression pattern in the glioblastoma tissue is superposable to the reported pattern for VEGF (16).

Subsequently, we investigated the regulation of IGFBP2 and two other family members, IGFBP3 and 5, under anoxic conditions. In PrGBM, increased expression of IGFBP2 was often paralleled by overexpression of either IGFBP3 (comprised in N566; Fig. 3,D) or IGFBP5 (comprised in the list of separating genes, see Table S1 of supplementary information). Glioblastoma cell lines LNZ308 and LN229 displayed IGFBP2 up-regulation concomitant with VEGF after 20 h of anoxia (Fig. 6). Similarly, IGFBP5 was induced in these two cell lines, whereas U87 that is wild-type for TP53 exerted most prominent induction of IGFBP3, a TP53 target gene. Thus, IGFBP3 and 5 might also be induced in the same process in gliomas. Induction of different members of the family in these cell lines is in line with a certain redundancy of these genes observed in IGFBP2 knockout mice. These show increased levels of other IGFBPs and display a much less dramatic phenotype than the one initially predicted based on the fetal IGFBP2 expression pattern (31).

Our characterization of astrocytic gliomas by their gene expression profiles revealed consistent inherent differences between LGA, ScGBM, and PrGBM. Thus, the previous histological and clinical recognition of these three glioma entities can be strongly supported by cDNA-array data as reflected by MDS of overall gene expression (Fig. 1). The LGA were found to be the most closely related group and well separated from PrGBM that displayed the most heterogeneous expression profiles. The ScGBM share some of the features of both subgroups, reflecting the common pathogenetic pathway with LGA and the malignant behavior with PrGBM.

Analysis of the expression profiles for correlated genes separating tumor subtypes allowed identification of sets of genes implicating particular biological features that are associated with tumor-specific subtypes. Most strikingly, a cluster (G5) of correlated genes suggests that inherent angiogenic activity, manifested in overexpression of known angiogenic factors, including the most potent, VEGF, and FLT1 and PTN, distinguishes PrGBM from the other two groups. In fact, this feature alone was sufficient to classify most astrocytic gliomas (92%) correctly according to their subtype. Glioblastoma are among the most vascularized tumors with characteristic, structurally and functionally abrogated blood vessels of multilayered glomeruloid pattern. Angiogenesis, an essential requisite for tumor growth, is a tightly controlled, complex process involving various pro and antiangiogenic factors that need to act in a concerted fashion to yield functional blood vessels (32). The angiogenic factors VEGF, angiopoietins, and Ephrins act specifically on endothelial cells, because their respective receptors FLT1 and KDR, Tie receptors, and Ephrin receptors are exclusively expressed on these cells. In contrast, other angiogenic factors highly abundant in glioma, such as FGF2, PDGF, and TGF-β, act also on other cell types (33).

According to our analysis, PrGBM appear to have higher angiogenic activity mediated by the genes comprised in G5 than LGA but also higher activity than ScGBM. The difference between PrGBM and ScGBM is unexpected, because they have the same malignancy grade (WHO IV) and cannot be discriminated by classical histology. This may reflect the fact that PrGBM are rapidly growing tumors that cannot keep up with their increasing need for blood supply and thus may suffer from more severe hypoxic conditions triggering angiogenic activity. In contrast, ScGBM may progress over years from LGA and thus may use different pathways for sustained angiogenesis. As compared with PrGBM, LGA display relative overexpression of the angiogenic factor FGF2 that emerged in the separating gene node N505 (Fig. 3 C; all LGA versus all PrGBM, P < 0.01, unpaired t test).

An important role for FGF2 in LGA is also supported by findings from expression profiling efforts reported from small sets of pediatric gliomas (6 LGA versus 7 high-grade astrocytoma; Ref. 34) revealing significantly higher FGF2 expression in LGA. Differential expression of VEGF, however, did not reach statistical significance, possibly because of the small numbers of samples. Interestingly, in this study (34) that appeared during revision of this manuscript, the authors used the opposite strategy to analyze their gene expression data obtained from these pediatric gliomas. A knowledge-based list of 133 genes related in some way to angiogenesis was constructed and successfully used to separate LGA from GBM according to hierarchical clustering and PCA. This supervised gene reduction yielded similar separation power as obtained using their full set of 9198 genes present on the chip, thus again emphasizing the important role of angiogenesis in tumor progression.

The biological differences between mechanisms of angiogenesis may have important implications for response to antiangiogenic therapy, which is about to enter the clinics (35).

An important implication of the CTWC method concerns investigation of the genes that belong to identified clusters, here in particular G5, for their biological function, assuming that coexpressed genes might be coregulated and therefore might be involved in the same biological process. The fact that IGFBP2 is strongly correlated with the angiogenesis genes came as a surprise that will be worth detailed additional studies. The presented association supported by our finding of coexpression of IGFBP2 and VEGF in pseudopalisading cells surrounding necrosis in GBM (Fig. 5) provides strong evidence that IGFBP2 overexpression in astrocytic gliomas may reflect its implication in angiogenesis. These correlations add to the fact the IGFBP2 can be induced under anoxic conditions in glioblastoma cell lines (Fig. 6) and mouse embryonic stem cells (26).

Overexpression of IGFBP2 in GBM has first been discovered by Fuller et al. (36) using gene expression profiling. This finding was confirmed (23, 27, 34, 37) and extended to other tumor types. Subsequent immunohistochemical investigations for IGFBP2, a secreted protein, suggested a significant association with the histopathological malignancy grade in glioma (37) with more pronounced immunopositivity in macrophage/microglial and glioma cells around focal necrosis (38). The IGFBPs were first identified and characterized based on their role to bind and modulate the in vivo bioactivity of the mitogenic and anabolic peptides IGF-I and IGF-II. More recently, it has become clear that IGFBPs are multifunctional proteins with IGF-dependent and -independent functions in controlling growth and metabolism. The intrinsic bioactivity of IGFBPs, in a positive or negative fashion, depends not only on the cell type but also on their differentiation state and physiological/pathophysiologic condition (39). Furthermore, additional genes of cluster G5 have been linked to the IGF/IGFBP system (SOCS2 and 3, ZFP36L1; Refs. 40 and 41). For IGF-I itself, a role in angiogenesis has been established (26, 42). Thus, our observations support the idea to develop a therapeutic approach targeting the IGF-I/IGFBP system in cancer.

Another separating gene cluster emerging from CTWC analysis associated with known biological function is G12 (Fig. 2; Table 3). This gene cluster partitioned a subgroup of PrGBM with increased expression from LGA and comprises only genes related to the immune system. Half of the genes are known to be IFN inducible and are likely expressed by lymphoid cells and/or macrophages and microglia; others represent markers of natural killer cells. These results confirm the observation by Fathallah-Shaykh et al. (29) who reported elevated expression of MHC class I and II genes. However, some of the genes are potentially expressed by tumor cells according to our preliminary results. It is worth to note that studies trying to relate immune response in gliomas with survival have not yielded conclusive results. The notion that gliomas produce an immuno-compromised environment by secreting cytokines, such as TGFβ-2 (43, 44), is also supported by this study; TGFβ-2 emerged in node N566 that separates PrGBM from LGA (Fig. 3 D). Along the same lines, we found HLA-G to be overexpressed (over normal brain, mean 2.9) in PrGBM. This nonclassical class I antigen is mainly expressed at the materno-fetal interface during pregnancy where it plays an important role in maternal-fetal tolerance. HLA-G expression has been associated with tumor progression in several tumor types and has been suggested to contribute to escape from immune surveillance (45, 46). Furthermore, functional assays with human glioma cell lines indicated that few HLA-G-positive cells within a population of HLA-G-negative tumor cells in vitro are sufficient to exert a significant immune inhibitory effect (46).

The tumor classifier using expression data of all four identified discriminatory clusters yielded correct tumor prediction in 93% of an independent set of astrocytic gliomas. This remarkable success has to be appreciated in light of the inherent inaccuracy of the initial tumor diagnosis for gliomas, because of sampling errors and subjective histological criteria that lead to striking interobserver differences among neuropathologists (47). In this context, three tumors originally designated as PrGBM attracted our attention (1357, 1297, and 1308), because they appear evenly in clusters characteristic for PrGBM and LGA, respectively, depending on the analysis (Figs. 2 and 3). Central pathology review revised the first (1357) as OAIII, whereas the other two cases were confirmed. Tumor 1297, clinically defined as a PrGBM, presented with genetic characteristics of a ScGBM, namely TP53 mutation and most prominent overexpression of PDGF receptor. Finally, glioblastoma 1308 exhibited overexpression of EGFR, characteristic for PrGBM. However, the biopsy used for gene expression profiling comprised a part of the tumor displaying only features of WHO grade II (the most malignant part of tumor defines diagnosis). This demonstrates that regions of a glioblastoma with the appearance of a lower grade astrocytoma exert already features only found in PrGBM, exemplified here as overexpression of the EGFR gene and the genes comprised in G12 and N566, characteristic for PrGBM. In contrast, after reclustering S1 (all samples) according to G5, comprising genes indicating high angiogenic activity, this tumor (PrGBM 1308) partitions with LGA. This is biologically convincing, because lower grade astrocytoma have no necrosis and thus are expected to suffer from less hypoxic stress. This observation might be of clinical relevance, knowing the difficulties in tumor sampling at surgery in particular stereotactic biopsies. Thus, the true tumor grade may be underestimated. It follows, the identification of markers specific for glioblastoma that are detectable in apparent lower grade areas of such tumors would allow to improve the diagnostic reliability and ultimately the choice of therapy. Interestingly, in MDS, both 1297 and 1308 were located at the interface of ScGBM and PrGBM.

Combined unsupervised and supervised analysis is a novel approach of gene selection that allows identification of clusters rich in genes informative (by supervised analysis) for tumor classification. CTWC was particularly suitable for this task, because it is designed to go beyond clustering of all genes on the basis of the data from all tumors, and clustering of all tumors, using data from all genes. This is of importance because most of the genes for which expression levels have been measured are irrelevant for the partition sought. CTWC proposes identification of correlated groups (clusters) of genes and using only data from one such group at a time to recluster the tumors, or vice versa.

Our gene selection approach has some interesting features, and allowed us: (a) to reduce the many discriminatory genes obtained by binary class comparisons; (b) enabled an almost correct classification of tumor entities, namely discrimination of LGA from ScGBM, and ScGBM from PrGBM, although pairwise comparisons yielded only few separating genes (partially because of small numbers of ScGBM in the training set, n = 5; threshold of FDR of q < 0.05); (c) by using the signal of a group of correlated genes, the noise of the individual measurements averages out and is reduced; and most importantly (d) the identified gene clusters yielded information on the biological context of coexpressed and possibly coregulated genes. Such insight may give some indication on biological processes determining tumor entities. This fact is highlighted in this study in particular by identification of G5, with high discriminatory power on its own, featuring angiogenesis-related genes.

Analyzing gene expression profiles has yielded biologically relevant information using a limited set of genes (1185). This approach is a first step toward molecular diagnostics for astrocytic gliomas that may improve tumor diagnoses in the future by adding objective criteria. Gene expression profiles may ultimately provide a tool for the identification of patients who are most likely to benefit from targeted therapy. It follows that such approaches for outcome prediction need to be established in clinical trials that in turn will allow rational design of future therapies.

Here, we discussed mostly findings demonstrating similarities within and differences between tumor subtypes for classification purposes. Numerous interesting clusters emerging from this analysis are awaiting attention to gain additional insights into the biology of astrocytic gliomas, e.g., in regard to differentiation and migration (e.g., TNR, node 505), and the involvement of the wnt/notch pathways as indicated in N505, and other clusters that have not been discussed in this study [found in G1(S2), see full CTWC analysis on the web⁵ ].

We thank our colleagues who made this study possible by participating actively in the clinical trial and providing primary tumor samples, namely Drs. J-G. Villemure, F. Porchet, O. Vernet, P. Otten, A. Reverdin, and B. Rilliet. We also thank Dr. P. Reymond for helpful discussions; Drs. S. Ostermann, M. Albertoni, and G. Pizzolato for their collaboration; and Dr. P. Walker for critical reading of this manuscript.