Identification of prognostic biomarkers for glioblastomas using protein expression profiling

A set of proteins reflecting the prognosis of patients have clinical significance since they could be utilized as predictive biomarkers and/or potential therapeutic targets. With the aim of finding novel diagnostic and prognostic markers for glioblastoma (GBM), a tissue microarray (TMA) library consisting of 62 GBMs and 28 GBM-associated normal spots was constructed. Immunohistochemistry against 78 GBM-associated proteins was performed. Expression levels of each protein for each patient were analyzed using an image analysis program and converted to H-score [summation of the intensity grade of staining (0–3) multiplied by the percentage of positive cells corresponding to each grade]. Based on H-score and hierarchical clustering methods, we divided the GBMs into two groups (n=19 and 37) that had significantly different survival lengths (p<0.05). In the two groups, expression of nine proteins (survivin, cyclin E, DCC, TGF-β, CDC25B, histone H1, p-EGFR, p-VEGFR2/3, p16) was significantly changed (q<0.05). Prognosis-predicting potential of these proteins were validated with another independent library of 82 GBM TMAs and a public GBM DNA microarray dataset. In addition, we determined 32 aberrant or mislocalized subcellular protein expression patterns in GBMs compared with relatively normal brain tissues, which could be useful for diagnostic biomarkers of GBM. We therefore suggest that these proteins can be used as predictive biomarkers and/or potential therapeutic targets for GBM.


Introduction
Glioblastomas (GBMs), the most common primary brain tumor in adults, show a median survival of <12 months due to its resistance to current medical treatments. Several genetic aberrations have been shown in GBMs (1). However, the molecular markers that correlate with the clinical outcome of GBMs are still required to establish the comprehensive molecular fingerprint. Such molecular profiling may eventually lead to diagnostic biomarkers and/or targeted therapeutic approaches that can improve the clinical outcome (2).
The Cancer Genome Atlas Research Network provided an integrative analysis of mRNA and DNA data (3). The landmark study was of particular importance because they integrated a large number of datasets and presented an unbiased and systematic cancer genome analysis (4). Yet, what is still lacking is a robust protein and clinical dataset, which could provide functional level understanding and thereby complement genomic and transcriptomic data. Therefore, we applied the tissue microarray (TMA) analysis to acquire the protein expression data (5,6).
Previous reports demonstrated brain tumor profiling using TMA is feasible (2,(7)(8)(9)(10). However, the conventional profiling has inherent limitations because they only provide categorical values. Moreover, acquisition of TMA data is often plagued with high subjectivity originating from intra-and inter-observer variations (11). Now that image analysis technique has improved, it would benefit us to use it to reduce systematic errors and human bias (12,13). Therefore, here, we used image analysis tool to more objectively quantify the values.
We applied the TMA technology for an exhaustive immunoprofiling of GBMs with integrating clinical variables in an attempt to identify new biomarkers for diagnostic and/or prognostic utility.

Materials and methods
Patients and tissue collection. Sixty-two GBM samples were obtained at Samsung Medical Center (Seoul, Korea) with written informed consent in accordance with the institutional review board between January 2004 and December 2006 (original set). Patients were managed according to established diagnostic and therapeutic protocols, including surgical resection and subsequent chemoradiotherapy. A macroscopic total resection was performed in 47 of 62 patients (75.8%), a partial resection in 14 of 62 patients (22.6%), and a biopsy only in one of 62 patients (1.6%). Tumor samples were re-evaluated by neuropathologists to confirm the diagnosis according to World Health Organization criteria. All patients underwent subsequent radiotherapy (60 Gy in 2 Gy fractions) after surgical resection. Fifty of 62 patients received temozolomide concurrent chemoradiotherapy with a median of 4 cycles (range, 1-9 cycles). However, the 12 others (20%) did not receive chemotherapy because of clinical deterioration during radiotherapy. There was no follow-up loss.
Another independent set consisting of 82 GBM patients diagnosed between January 2004 and December 2007 from the same institute was used for validation (validation set). The 82 GBMs of the validation set were also managed according to established diagnostic and therapeutic protocols. A macroscopic total resection was performed in 59 of 82 patients (72.0%), a partial resection in 19 of 82 patients (23.2%), and a biopsy in 4 of 82 patients (4.9%). All patients underwent subsequent radiotherapy after surgical resection. Fifty of 82 patients had temozolomide concurrent chemoradiotherapy with a median of 5 cycles (range, 1-17 cycles). There were 6 follow-up losses.
Of the 82 patients, 39 were included in both the original and validation datasets owing to the difficulty of acquiring adequate number of surgical samples or public TMA data of GBMs. Instead, the surgical samples of the 39 GBMs were newly processed for TMA and expression of proteins was re-examined. Detailed patients' clinical data of the original dataset was also reported in Kong et al (14).
Construction of tissue microarray. Surgical samples were fixed in 10% formalin solution (Sigma-Aldrich) for 24 h at 4˚C within 24 h after surgery and then paraffin-embedded. A representative area of each GBM was marked on an H&E section of each patient's paraffin block avoiding necrosis and extensively vascularized area. Corresponding tissue core of 2 mm diameter was extracted from the original donor block using an arraying machine (MTA-1, Beecher Instruments). The cores were fit into a vertical hole that was bored in a recipient paraffin block. Recipient blocks were incubated at 58˚C for 5 min, pressed on a hot plate for 3 min, and cooled in ice water to enable tissue cores to integrate into the recipient block. Sections of 4 µm thickness were cut from each array block.
Immunohistochemistry and generation of protein expression values. Immunohistochemical staining against 78 tumor-associated gene products was done using a standard procedure (14). Briefly, sections of TMAs were prepared on the slide, baked at 55˚C for 30 min, deparaffinized in xylene and rehydrated in graded concentrations of ethanol. Antigens were retrieved in 10 mM citrate buffer (pH 6.0, Dako) for 5 min in MicroMED T/T Mega microwave (Milestone). Endogenous peroxidase activity was blocked by incubation in 0.3% hydrogen peroxide in methanol. Primary antibodies (overnight at 4˚C, Table I), biotinylated secondary antibodies (1:200, 1 h at room temperature, Vector) and an ABC kit (1 h at room temperature, Vector) were applied sequentially. Diaminobenzidine tetrahydrochloride (DAB) was used as the enzyme substrate. Specificity of primary antibodies was validated by i) no immunoreactivity in GBM sections which were reacted without primary antiserum as negative controls, and ii) comparing staining pattern of each antibody with previous studies that used the same antibody.
One hundred and eight subcellular location-specific proteins were extracted based on previously reported cellular localizations of 78 proteins (e.g., Myc_Cytoplasm). Immunostained TMA slides were scanned by Aperio Scan Scope CS System and converted into image files. Staining intensity of each protein (grade 0, negative; 1, weak but detectable above control; 2, moderate; and 3, strong) and percentage of positively stained cells for each intensity were automatically analyzed in the whole area of each spot by TissueMine (Bioimagene Co.) (15).
TissueMine determines positive versus negative staining based on the colorimetric differences between stained and unstained subcellular location (i.e. nuclei; positive nuclei demonstrating brown staining consisted of pixels with red greater than blue). These algorithms then used the gray scale (0-255) to quantitate intensity of staining. The grade of intensity is determined by percentage of positive cells. For example, the grade 1 means that the percentage of positive cells is between 10-25%.
For each observed tissue component, a summary value we refer to as H-Score was calculated. This consists of a sum of the percentages of positively stained cells multiplied by a weighted intensity of staining: where P i is the percentage of stained cells in each intensity category, and i is the intensity for i = 0, 1, 2, 3 (16)(17)(18)(19)(20). Kruskal-Wallis test and plotting were performed to verify the utility of H-score and to confirm the correlation of it with the manual grade. The validity of this automated process was confirmed by clinical experts (Fig. 1).

Data analysis
Exclusion criteria and normalization. We excluded samples according to the exclusion criteria: i) GBMs missing >15% of the protein expression values, ii) recurrent GBMs, iii) GBMs without clinical data. To remove non-biological origin of variation between arrays, we performed the quantile-normalization on the assumption that the distribution of protein expressions for each patient would be the same, similar to what is already observed for patient's gene expression distribution (21). In the quantile-normalization, a qqplot was utilized to compare distribution of datasets from patients. Projecting each value onto the unit diagonal in the qqplot makes distributions identical. To adapt the method for TMA data, we transposed TMA data matrix to set columns and rows to samples and proteins and used normalize.quantiles function in affy package of R program.
Hierarchical clustering and survival analysis. Hierarchical clustering (distance, standard Euclidean distance; criterion, complete-link) was performed using the R program package (22). With clustering results, we performed Kaplan-Meier survival  Selection of biomarker candidates. Two clusters differentiated by the hierarchical clustering from the original dataset were compared by Student's t-test which is then used for finding differentially expressed proteins. We performed multiple testing corrections using Benjamini & Hochberg false discovery rate to correct for the occurrence of false positives (23).
Supervised analysis using the classification methods. We used supervised methods [support vector machine (SVM), random forest (RF), and general linear model (GLM)] to confirm the utility of whether the set of biomarker protein candidates distinguished the groups appropriately into survival status. Among the 56 cases, two-thirds of the cases from each group were randomly assigned to the training set, and the remaining one-third was assigned to the test set. A receiver operating characteristic (ROC) curve was constructed, and the area under the curve (AUC) was calculated. The class, defined as the prognosis outcome for any patient's protein expression data using the output computed by the model, was predicted by R program package (22). Validation of the biomarkers. We validated our ten prognostic markers against another independent GBM TMA dataset. Same preprocessing works (immunohistochemistry, generation of protein expression values, normalization, hierarchical cluster, survival analysis) were performed. We re-evaluated the 10 biomarkers using a publically available DNA microarray data [GSE4271 from the gene expression omnibus (24,25)]. We normalized the gene expression CEL files using Robust Multichip Averaging procedure, and PM-MM difference model was used to obtain the expression values. With nine candidate genes which corresponded to our ten prognostic markers, we calculated the Euclidean distance between the patients and constructed a corresponding distance matrix. The resulting nine-dimensional data was rescaled to a two-dimensional map using the multidimensional scaling (MDS) method. We performed the k-means clustering on the MDS result with the number of clusters set to two (k=2) to correspond to the clustering performed in TMA analysis. Survivals of the two groups were compared by the Kaplan-Meier survival analysis and log-rank test.

Results
Expression profiling of 108 subcellular location-specific proteins. The expressions of 78 proteins (Table I) were studied by immunohistochemistry on TMAs containing 62 GBM and 28 tumor-associated normal (normal brain tissue extracted around tumor tissue) spots proteins. For each protein, grade value (0-4) of each spot was made by pathological reading. H-scores (a sum of the percentages of positively stained cells multiplied by a weighted intensity of staining: grade 0, negative; 1, weak but detectable above control; 2, moderate; and 3, strong) were also generated by an image analysis program (12). In the generation of H-score, each protein was read in each subcellular localization, i.e., nucleus, cytoplasm and membrane, based on previously reported cellular localizations of 78 proteins (e.g., Myc_Cytoplasm). Consequently, we obtained 108 subcellular location-specific proteins (Table I). According to criteria for exclusion, 56 tumor spots and 26 normal spots were included in further analysis (male, 34; female, 22).
To verify the utility of the H-score, we compared manually constructed numerical values against automatically generated H-scores in 20 randomly selected proteins. We concluded that using H-score was reasonable for profiling protein expressions as those values were well-correlated (p<0.05, Kruskal-Wallis test, Fig. 1).

Hierarchical clustering and survival analysis.
The H-scores were normalized with quantile normalization to make the data statistically comparable. We clustered the normalized data with hierarchical clustering method. Tumor and relativelynormal spots were well-divided (p=9.9904x10 -16 , Fisher's exact test), confirming that this method was suitable for separating spots into groups having characteristic protein expressions. Using only tumor spots, we performed quantile normalization ( Fig. 2A) and clustering analysis. The combined protein expression patterns defined 56 GBMs into cluster 1 (n=19) and cluster 2 (n=37) (Fig. 2B).
We performed survival analysis to determine whether the two clusters represent distinct prognostic subgroups. Indeed, cluster 2 had a significantly better survival than cluster 1 [median survival (months) of cluster 1, 11.7 (range: 0.7-20); cluster 2, 13.2 (range: 1.4-30.4), p<0.05, log-rank test, Fig. 2C]. The clinical characteristics of the two groups showed no significant difference (Table II). On multivariate survival analysis, molecular classification (clusters 1 and 2) and age (<70 and ≥70 years) were the most significant factors to distinguish patients by prognosis (p<0.01 and p<0.01, respectively, Table III). Based on these significant results, we determined that two groups from clustering analysis were clinically distinct.
Biomarker candidates and supervised analysis. Student's t-tests were utilized to identify proteins whose expression significantly differed between the two groups with Benjamini & Hochberg False Discovery Rate. Ten subcellular location-specific proteins were identified (Table IV, q<0.05, CDC25B_nuclear, cyclinE_ nuclear, p16_nuclear, p16_cytoplasm, TGF-β_cytoplasm, histone_nuclear, p.VEGFR2.3_cytoplasm, DCC_nuclear, survivin_nuclear, p.EGFR_cytoplasm). Of these, four proteins were overexpressed while the remaining six were underexpressed  The log-rank test shows that the difference between two curves is significant (p<0.05).
A B C in cluster 1 (Fig. 3). With a more stringent statistical criterion (q<0.01), CDC25B_nuclear, cyclinE_nuclear, p16_nuclear and p16_cytoplasm were significantly different. Functional categorization of these proteins using the PANTHER (26) ontology database showed that they are significantly related with p53 pathway, angiogenesis, cell cycle, Ras pathway and VEGF signaling pathway (p<0.05). Three supervised (classification) models (SVM, RF and GLM) were used to confirm that these proteins were available for classification for prognosis of GBMs. SVM shows the best performance in our data. The prediction accuracy of SVM for patient survival was 87.5% with 0.05 adjusted p-value (AP) and 80.4% with 0.01 AP. The sensitivity for better prognosis group of classifier for patient survival was 100% (specificity: 63.2%) with 0.05 AP and 94.6% (specificity: 52.6%) with 0.01 AP. ROC curves for 0.05 AP and 0.01 AP showed that classification model using 0.05 AP performs better than 0.01 AP (data not shown).
Interestingly, survivin_nuclear was the most significant factor for 0.05 AP and cyclinE_nuclear was the most significant factor for 0.01 AP for each classification model in GLM. Contrary to our expectations, the same protein was not the most important factor in both models. This discrepancy was, however, rectified when considering that these two proteins both belong to the chemo-radiation-resistant group (poorer prognosis group). In addition, further analysis by protein clustering showed that they belonged to the same cluster, implying that these two proteins had a similar pattern. Consequently, with this analysis, we could confirm that the result of hierarchical clustering, unsupervised method was well explained by supervised methods.
Validation of selected biomarker candidates. In order to validate the clinical significance of the 10 markers, an independent TMA consisting of 82 GBMs was analyzed. We obtained expression values of the 10 subcellular location-specific proteins and calculated H-scores. Two groups were separated (cluster 1A and 2B, n=59 and 23, respectively). These two groups also showed statistically significant difference in the survival [median survival (days) of cluster 1A, 391 (range, 42-730); cluster 2B, 415 (range, 22-929), p<0.05, log-rank test, Fig. 4A]. The clinical characteristics of the two groups showed no significant difference (Table V). These results showed that the biomarkers obtained from the original set are valid to another independent data set.
Protein expression data are arguably more useful than transcriptional data because it provides the data at the functional level of cells, and thus it is not uncommon to observe differences in the expression of proteins and mRNAs. Nevertheless,   Fig. 4B). Survivin and cyclin E were also the most significant contributors for discriminating between two prognostic groups (GLM). According to the TMA analysis, these proteins were overexpressed in the chemo-radiationresistant group. Interestingly, the same observation was made for their corresponding genes from the DNA microarray analysis, confirming other previous reports (27,28) and validating our data. The results here show that our TMA analysis is indeed a valid approach in differentiating prognostic groups.

Aberrant protein expressions in GBM.
Aberrant or mislocalized protein expression patterns could be useful for diagnosis of GBM. Accordingly, we further compared expressions of all 234 subcellular location-specific proteins (78 proteins) of GBMs with those of relatively-normal tissues using the original dataset. H-scores of 16 subcellular location-specific proteins were decreased and 16 ones increased significantly in GBMs (data not shown). Of those, cyclin E expression in the cytoplasm was increased whereas the expression in the nuclear was   for each group from clustering analysis on public DNA microarray data were compared with overall survival. Clustering analysis was performed with nine candidate genes (corresponding to proteins of ten pairs based on their cellular localization). Differences in the survival were tested for statistical significance by the log-rank test.
A B decreased in GBMs. Therefore, this particular protein could be one of mislocalized candidates specific for GBMs. Previously, cyclin E was validated as a prognostic biomarker for GBM (27), which could underline the significance of cyclin E in GBMs.

Discussion
This study shows that molecular classification of GBMs can be accomplished based on the immunohistochemical profiles of biomarkers using the TMA methods. Recent studies have also supported that clinical application of the TMA methods to the molecular classifications of various cancers (29,30). To provide a more precise data analysis, this study utilized scaled continuous protein expression values (H-score) and could find more significant differences using parametric statistic tests with continuous variables. The overall protein expression patterns using 108 subcellular location-specific proteins defined two survival-associated clusters. However, 108 expression values were too many for practical use. To reduce the number of markers, we performed a feature selection using Student's t-test. Ten and four subcellular location-specific proteins were selected with 0.05 and 0.01 adjusted p-value cut, respectively. Survivin_nuclear was the most significant factor for 0.05 AP and cyclinE_nuclear for 0.01 AP. Survivin is known to be an inhibitor of apoptosis that acts via a pathway independent of Bcl-2. Previous studies showed a correlation between increased protein/mRNA levels of surviving and adverse prognosis in various cancers (28,(31)(32)(33)(34). Cyclin E mediates the initiation of DNA synthesis in the late G1 phase by activating cyclin-dependent kinases 2. Abnormal expressions of cyclin E (27,35,36) and other biomarkers (37)(38)(39)(40)(41)(42) we selected have frequently been found in cancer cells.
In this study, phosphorylation-specific antibodies against p70 S6 kinase, Akt, PDGFR-α, PDGFR-β, VEGFR2, VEGFR2/3 and EGFR were included since they are the most important targets of newly-developing targeting agents. In contrast to our expectation, phosphorylated EGFR and VEGFR2/3 were inversely correlated with worse prognosis. Although EGFR and VEGFR are heavily activated in GBMs, agents targeting them showed disappointing results in clinical trials (43,44). Therefore, our results could reflect those results and indicate that other specific targets that have prognostic significance need to be elucidated.
In this study, we could not obtain non-neoplastic human brains as controls. Alternatively, tumor-associated normal tissues that were adjacent to the tumors were utilized. These tissues may harbor GBM cells, therefore, we selected the outermost region of the surgical samples. In addition, the tissues were clearly separated from the tumor cores by hierarchical clustering. This result indicates that there was minimal GBM cell-contamination, and, though not the most ideal, they are valid as controls.
It could be argued that due to the immense heterogeneity of cancer specimens, the normalization method we adopted may not seem pertinent to our study. Also, the quantile-normalization has an admitted limitation in that it removes outliers which could be of meaningful value. Our H-score method scaled the range of expression of all the cells found in the entire area of each spot of TMA, therefore, it would represent the heterogeneity well. Moreover, since the scaled range cannot have outliers, we circumvented the limitations of normalization by our scoring method.
Our 10 biomarkers are proven powerful to predict prognosis of GBM by analyzing three independent datasets. Moreover, the suggested biomarkers are optimal for practical use in pathology laboratories with respect to cost and time of prognostic evaluations. They also can provide the basis for developing new personalized approach as well as drug discovery. In addition we expect that our findings at the protein level can complement the transcriptomic and genomic data of GBMs (3,45) and would improve the molecular understanding of GBM.
For clinical application of the results, further analysis with a larger sample set and more detailed validation are still necessary. Nevertheless, the approach used in this study suggests that subjective interpretation of TMA profiling can be minimized and that larger scaled TMA analysis can be simultaneously performed. To our knowledge, this is the first TMA data analysis study based on quantitative values of protein expression using an image analysis tool. By applying bioinformatics techniques with large-scale data, it is now possible to perform comprehensive analysis and identify the cluster set of protein biomarkers for GBM.