Molecular Subtypes in Glioblastoma
Multiforme: Integrated Analysis Using
Agilent GeneSpring and Mass Profiler
Professional Multi-Omics Software
Application Note
Authors

Abstract

Durairaj Renu, Pritha Aggarwal,

Agilent GeneSpring and Mass Profiler Professional (MPP) software was used

Sunil C. Cherukuri, and Pramila Tata

to perform integrated multi-omics analysis of mRNA expression, miRNA

Strand Life Sciences

expression, protein expression, and copy number aberration (CNA) associatedBangalore, India

with glioblastoma multiforme (GBM), the most frequent and lethal of the central

Vadiraja Bhat, Carolina B. Livi, andMichael Rosenberg
Agilent Technologies Inc.
Santa Clara, California

nervous system tumors. The genomic data was obtained from The Cancer
Genome Atlas (TCGA) project. The molecular subtypes defined by TCGA for GBM
were based on mRNA expression. We applied GeneSpring metadata analysis
framework to explore the relationship between GBM subtypes and genomic CNA.
A multi-omics correlation analysis was performed between mRNAs involved

Mona Al-Gizawiy, Kathleen Schmainda,

in neurogenesis and their validated miRNAs highlighting the molecular eventsand Shama P. Mirza

in Proneural tumors that result in suppression of p53 signaling. Furthermore, aMedical College of Wisconsin

combined analysis of mRNA expression data from TCGA with an independent,

Milwaukee, Wisconsin

label-free GBM proteomics study revealed core subset of 54 signature genes andproteins highlighting the integrated biology analysis power of GeneSpring/MPP
software. The core genomics signature clearly differentiated the known tumor
subtypes in TCGA cohort. Intriguingly, while TCGA categorization of theproteomics samples were not experimentally determined, the core proteomics
signature separated this cohort in 2 or 3 groups, suggesting that different TCGA
subtypes may be present.

Introduction
Gliomas are the most frequent central
nervous system (CNS) tumors in adulthumans. Gliomas are classified into
astrocytomas, oligodendrogliomas,
and ependymomas based on histology.
Astrocytomas represent two thirds ofall diagnosed gliomas and are gradedfrom I to IV according to the World
Health Organization (WHO)1. GBM(WHO grade IV Astrocytoma) is themost frequent and lethal among thesetumors, with an estimated global
incidence rate of approximately 3–4 per
100,000 people per year (http://www.
pubcan.org/cancer/4938/glioblastoma/epidemiology). Subjects with GBM
have an overall median survival of only
15 months2,3.
The TCGA project has previously
identified genetic changes in the DNA
sequence and copy number, DNA
methylation, gene expression, andphysiological information for a set of
GBM tumors4. The results generated greatinterest in compiling and characterizingchanges in the genome, transcriptome,
and proteome of human GBM tumors,so that aberrantly functioning molecularpathways and tumor subtypes could
be identified. Phillips et al.5 andVerhaak et al.6, identified subtypes inprimary GBM tumors by clusteringmRNA expression datasets from TCGA.
Phillips et al.5, categorized GBM into
three subtypes: Proneural, Proliferative,
and Mesenchymal. Verhaak et al.6
integrated mRNA profiles frommultiple array platforms to define four
expression‑based subtypes called
Proneural, Classical, Mesenchymal,
and Neural. Noushmehr et al.7 profiled
promoter DNA methylation alterations
in GBM tumors to define a glioma-CpG
island methylator phenotype (G-CIMP),
which was preferentially enriched in theProneural subtype.
In this study, the metadata andcorrelation analysis framework in
GeneSpring/MPP multi-omics software
were used to understand the molecular
events associated with GBM subtypes.

Highlights
•

•

•

•

The molecular subtypes of GBM
defined by TCGA study6 based onmRNA expression were extended
to understand the relationshipbetween the subtypes and the
genomic copy number aberrations
using a metadata analysis
framework.
Correlation analysis performed
between mRNA and miRNAdatasets from TCGA identified
potential candidate miRNAs
regulating the known Proneuralsignature genes OLIG2, NKX2-2,
and SOX11. The molecular eventsin Proneural tumors leading tosuppression of p53 signaling werehighlighted.
Analysis was performed on an
independent discovery proteomics
study of GBM and control (brain
fragments from epilepsy surgery)subjects measured using alabel-free proteomics approach.
Sample-sample correlation and
PCA analysis resulted in clear
separation between control andtumor samples. In addition, thetumor samples were segregatedinto two to three subgroups.
Combined analysis of mRNA
expression data from TCGA study
and proteomics data from thelabel-free study identified a subsetof the GBM subtype signatures,
which clearly differentiate theknown subtypes in the larger TCGA
cohort.

Dataset and Methods
TCGA sample details

The mRNA expression dataset was
obtained from https://tcga-data.nci.nih.
gov/docs/publications/gbm_exp/ inthe form of unified expression values of
840 mRNA profiles for 173 GBM tumor
samples from multiple array platforms.
Expression data of 534 miRNAs studied
using the Agilent 8 × 15K array was
2

downloaded from the TCGA portal. The
mRNA, miRNA, and protein expression
datasets were imported as a generic
single-color expression experiment in
GeneSpring 13.0 (http://genespringsupport.com/).
Copy number aberrations were obtained
from the sample-centric table (https://tcga-data.nci.nih.gov/docs/publications/gbm_2013/) generated by Brennan et al.8
The table provided gene-wise CNA forthe samples in the TCGA cohort. CNA
is characterized by the measured copy
number (expressed as log2 ratio versus
normal euploid genome), and by theextent of the change in the genome.
Aberration thresholds determined by
Brennan et al.8 were:•Deletion (x ≤ –1)•

Partial deletion (–1 < x ≤ –0.7)•

Loss (x < –0.2)•

Gain (0.2 ≤ x < 1)•

Intermediate amplification(1 ≤ x < 2)•

Amplification (x ≥ 2)

Proteomics sample details

Both GBM and epilepsy specimens
were obtained from temporal, frontal,
and occipital regions of the brain. GBM
tumor tissue and epileptic brain tissue
specimens were collected at the time oftissue resection and snap-frozen in liquid
nitrogen within 30 minutes of removal.
The samples were stored at –80 °C in
the Brain and Spine Tissue Bank at the
Medical College of Wisconsin (MCW)
until used for sample analysis by mass
spectrometry (MS).
A comprehensive analysis of 10 tumor
tissue specimens of GBM and 10 epilepsy
controls was performed using ahigh‑throughput label-free quantification
approach by MS. All the tissue samples
were processed, and trypsin digestion
was performed as described earlier9. The
homogenized tumor tissue proteome
was fractionated through SDS-PAGE.
Internal DNA markers were added during
PAGE to facilitate consistent gel-based
fractionation. The in-gel tryptic digested

The unified expression values of
840 mRNA profiles for 173 GBM samples
are shown in Figure 1A. As expected,
the hierarchical clustering analysis
revealed the four predominant clusters
that coincide with the GBM subtypes.
Additionally, the clustering on entities
resulted in nine clusters, defined by theirexpression characteristics in the GBM
subtype (Figure 1A); see miRNA section
for discussion.

DOWN
Proneural

UP ClassicalUP
Mesenchymal

DOWN Neural

UP Classical

DOWN
Mesenchymal

To exclude the possibility that theobserved tumor subtyping was due tobatch effects, sample quality control
parameters such as batch, sample
purity (% tumor nuclei), and sample
quality (% necrosis) were visuallycorrelated with the samples using themetadata framework in GeneSpring.
Figure 1B summarizes the sample quality
parameters, aligned against the cluster
tree for the 173 GBM samples, as a heat
map.

A

B

SubtypeClassical

Batch info
Batch 1

Batch 4

Mesenchymal

Batch 2Batch 5

Gender-7
0
Female% Necrosis
Male

NeuralBatch 3Batch 60

Batch 710

7Subtype
Batch info
Gender% Necrosis% Nuclei

Sample quality parameters

Principal Component Analysis (PCA) wasused to confirm the GBM subtypes fromTCGA sample data, and to identify sample
groups in the proteomics study. The data
was mean-centered and scaled prior toperforming the PCA.
Gene ontology (GO) analysis was
performed based on the GO terms
provided by the GO Consortium(http://geneontology.org). The clusters
obtained from hierarchical clustering
were translated into an Agilent wholegenome array experiment in GeneSpring
to identify significant GO terms.
Hypergeometric distribution was used
to identify the GO terms enriched in
the clusters with respect to all known
annotated genes.

UP Proneural

All analyses reported in this study were
performed using GeneSpring/MPP 13.0
multi-omic bioinformatics software. GBM
subtypes were recreated by hierarchicalclustering of the 840 mRNA profiles
using Euclidean Distance metric and
Ward's Linkage rule. Genomic aberration
information, that is, copy number values
and other sample parameters obtained
from TCGA were imported as metadatain GeneSpring. The GeneSpring metadata
framework was used to assess batch
effects and sample quality parameters
commonly known to interfere with gene
expression analysis. The metadata
framework was also used to visuallycorrelate the expression subtypes withcopy number events.

mRNA expression-based molecularsubtypes in GBM

Correlation was performed between
mRNA profiles and miRNA profiles for
GBM using 159 GBM samples for which
both mRNA and miRNA measurements
were available. The targeting miRNA forgenes involved in neurogenesis wereidentified by querying the miRWalk
database (http://www.umm.uniheidelberg.de/apps/zmf/mirwalk/).

UP Neural

Data analysis

Results and Discussion

The GeneSpring correlation analysis
framework supports correlation between
molecules studied by a single -omics
platform or between two different -omics
platforms. It also supports correlation on
samples measured using a single -omics
technology. In this study, all sample andentity correlation analysis was performed
using the Pearson similarity metric.

DOWN Classicalpeptides were analyzed in triplicate indata-dependent mode using an Agilent
HPLC-Chip/MS system interfaced toan Agilent 6550 iFunnel Q-TOF MS.
The Q-TOF data was first processed by
Agilent Spectrum Mill Software against
the UniProt human database, and search
results were auto-validated at 1 % FDR.
Spectrum Mill has the unique ability toexport protein abundances from databasesearch results to the MPP module of
GeneSpring software for statisticalanalysis.

40% Tumor nuclei

Proneural

10

100

Figure 1. A) Molecular subtypes in GBM defined by Verhaak et al. validated in Agilent GeneSpring. The
entities could be grouped into nine different clusters. B) Sample quality metadata visualized as a heat map.
63

Batch effects are a notorious problem
in high-throughput experiments; insome cases they can mask biologicalphenomena. We did not observe batch
effects or gender bias in the GBM
samples. There is also no concordance
between sample purity and themolecular subtypes defined for GBM.
The metadata analysis confirms that themolecular subtypes defined for GBM
based on mRNA expression are truly
biological in nature and are not causedby experimental artifacts, which arecommonly observed in genomic analysis.
Pearson correlation analysis performedon the 173 GBM samples, using the same
mRNA expression profiles, confirmed
the subtypes observed in the clustering
analysis. The sample correlation analysis
indicates the similarity in the expressionof samples within a group and condition
if it exists between the sample groups.
As shown in Figure 2A, most samples
within the same subtype were positively
correlated to each other, indicating thesimilarity in their expression profiles.
The inter-sample variability was foundto be minimal for Proneural (Figure 1A
and Figure 2). The Proneural subtype
is known to mimic the expression ofa differentiated neuron and has a low
rate of proliferation5. This could explain
the observation that Proneural samples
cluster more tightly compared to othersubtypes. Proneural tumors are also
known to be less proliferative thanMesenchymal and Classical subtypes.
Consistent with this fact, we observedthat almost all of the Proneural samples
are negatively correlated to Mesenchymal
and Classical samples, indicating the
striking contrast in their expression
profiles. The largest sample‑sample
variability was observed in the Neuralsubtype.
PCA performed based on the 840 mRNA
profiles further confirmed the presenceof GBM subtypes. The scores plot for
the PCA is shown in Figure 2B. In PCA,
principle component 1 (PC1) is the
Eigen vector that captures the primaryvariation in the dataset. In our analysis,though the sample groups are separated
by PC1, maximum variation across PC1

is observed between the Proneural andMesenchymal tumors. PC2 captured
the variation between the Classical and
Neural sample groups.

Copy number events correlatedwith molecular subtypes in GBM

To gain insights into the genomic events
differentiating the expression subtypes,
we examined their relationship withgenomic CNA. Figure 3A provides the
copy number amplification and deletionof EGFR, CDKN2A, PDGFRA, TP53,
NF1, PTEN, and FGFR2 genes as a heat
map along with the subtype heat map.
Ordering the samples by their molecularsubtype reveals certain interesting
relationships. For example, the GBM
samples of Classical subtype do notharbor any copy number changes in TP53
or NF1. Instead, deletion of CDKN2A andamplification of EGFR is more prominent
in the Classical subtype. These findings
confirm observations reported by
Verhaak et al.6

B

A

PC2

Proneural

NeuralMesenchymal
PC1Classical-10Classical

NeuralMesenchymal

Proneural

1

Figure 2. Presence of four subtypes in the TCGA cohort confirmed by using sample-sample correlation and PCAs of 173 GBM samples in Agilent GeneSpring(A) Pearson correlation between GBM samples and (B) PCA.

4

PC3

Sorting the heat map by copy number
values of the key marker genes canfurther clarify their distribution between
GBM subtypes. It is evident from
Figure 3B, (heat map sorted on NF1)

that, in addition to Classical and Neural,
Proneural tumors also typically do nothave NF1 mutations. Instead, as reported
by Verhaak et al.6, NF1 hemizygous
deletion is a distinct characteristic of the
Mesenchymal subtype.

GeneSpring metadata framework reveals
other notable events occurring in GBM
subtypes:•

EGFR and CDKN2A mutations are
the most common, occurring ina large number of GBM subjects.
While more prominent in Classical,
these mutations are not confined
to any one subtype, but are present
across all subtypes (Figure 3A).

•

TP53 deletion is conspicuously
absent in the Classical subtype andis observed predominantly in theProneural subtype, as reported.

•

High level amplification of PDGFRA
has been reported to be mostly
restricted to Proneural6, but wealso observe that intermediatelevels of amplification of PDGFRA
are not restricted to any single
subtype (data not shown, can beobserved by sorting on PDGFRA
heat map).

SubtypeClassical

NeuralMesenchymal

Proneural-7
Copy number aberration
-2

0

2

7

Sorted by

Mesenchymal
NeuralProneural

Copy number aberration-20Sorted by

2

SubtypeClassical

MGMT status
Methylated

Mesenchymal

Unmethylated

NeuralNA-70Sorted by

7

Subtype
MGMT status

SubtypeClassical

Subtype
EGFR
CDKN2APDGFRA
TP53
NF1
PTEN
FGFR2

C

Subtype
EGFR
CDKN2APDGFRA
TP53
NF1
PTEN
FGFR2

B

0

Subtype
EGFR
CDKN2APDGFRA
TP53
NF1
PTEN
FGFR2

Copy number aberrationA

Proneural

Figure 3. A) Copy number aberrations of some of the key genes- EGFR, CDKN2A, PDGFRA, TP53, NF1, PTEN, and FGFR2 in relation to the GBM samples sortedby their molecular subtypes. B) Genomic copy number aberrations sorted based on NF1 and FGFR2, respectively, as indicated by the arrow (C) GBM samples
sorted based on their MGMT gene methylation status. Note the total loss of clustering pattern observed in Figure 3A.

5

•MGMT methylation is a predictive
marker for treatment response in
the Classical subtype8. However,
we did not find any patternthat would associate MGMT
methylation with GBM subtypes,
suggesting that yet unidentified
factors may exist (Figure 3C).

•

shows that about one third of thesamples have deletion of FGFR2.
The deletion pattern coincides withthe deletion of PTEN.
••

A comparison of the two sorted
heat maps in Figure 3B showsthat the FGFR2 deletion patternclosely mirrors the deletion patternof PTEN in a substantial numberof samples, indicating that thesesamples might have had a singledeletion event leading to theloss of both PTEN and FGFR2
located on 10q23.3 and 10q26.1,
respectively.

FGFR2 has been suggested tofunction as a tumor suppressor
gene in Glioma10,11. However,
PathCards does not include FGFR2
in Glioma or the GBM pathway(http://pathcards.genecards.
org/pathway/2538). Figure 3B(heat map sorted on FGFR2),

Klinz et al.12 had summarized thefrequency of copy number eventsin GBM subtypes (Table 1). The
analysis performed using themetadata framework in GeneSpring
was found to be broadly in sync
with the copy number eventssummarized by Klinz et al.12

miRNA regulation insubtype‑specific expression

Figure 4 illustrates the mRNA expression
pattern of the nine clusters noted inhierarchical clustering (Figure 1A). GO
functional enrichment analysis for theclusters identified the 840 signature
genes to be primarily involved in cellular%

Table 1. Frequency of copy number alterations in Glioblastoma subtypes (Klinz et al.12).
100
90
80
70
60
70
40
60
20
10
0

Classical
Mesenchymal
Proneural
NeuralAmplification events

EGFR
CDK6
MET
(7p11.2)(%) (7q21.2)(%) (7q31.2)(%)

100
95
54
96

92
89
46
96

86
91
54
92

Cluster 00

Deletion events(7q34)(%)

PDGFRA(4q12)(%)

86
91
52
92

5
9
35
13

NF1
PTEN
CDKN2A/C
RB1(17q11.2)(%) (10q23)(%) (9p21.3)(%) (13q14)(%)

Cluster 01

NeuralProneural ClassicalMesenchymal

NeuralCluster 03

Proneural ClassicalMesenchymal

NeuralProneural ClassicalMesenchymal

NeuralProneural ClassicalMesenchymal

Mesenchymal

Mesenchymal

NeuralProneural ClassicalMesenchymal

Proneural ClassicalMesenchymal

Cluster 08
4
2
0
-2

4
2
0
-2
Proneural Classical

Proneural Classical

4
2
0
-2

Cluster 07

4
2
0
-2
NeuralNeuralCluster 05

4
2
0
-2

Cluster 06

16
53
52
46

Cluster 02

Cluster 04

4
2
0
-2

95
67
56
71

4
2
0
-2

4
2
0
-2

4
2
0
-2

100
87
69
96

5
38
6
17

NeuralProneural ClassicalMesenchymal

Figure 4. Box plot visualization of the nine clusters showing subtype specific expression patterns.

6

Neuralphysiology events such as extracellularmatrix organization, response towound healing, response to stress,
nervous system development, glial
cell differentiation, actin cytoskeleton
organization, and chromatin organization.
Additional analysis performed in
GeneSpring confirmed that there is no
significant overlap between the genes
known to participate in pathways
implicated for GBM (such as RTK,
AKT, RB, P53, and PI3K) and the genesignature that defines the molecularsubtypes (data not shown).

miR-29a, miR-27a, miR-155, and miR‑21
displayed negative correlation to the
25 mRNAs (Figure 6). The miRNAs
up‑regulated in Mesenchymal subtype
by 1.5 fold in comparison to Proneuralare significantly enriched in the clusterthat displayed negative correlation to thegenes in nervous system development.
In Figure 6, the enriched miRNAs arehighlighted by red bars. Most of themiRNA-mRNA pairs in this cluster had
negative correlation values in the rangeof –0.3 to –0.5.

GO analysis identified nervous systemdevelopment (that is, neurogenesis)as the primary biological processenriched in the cluster of genes thatwere up-regulated in Proneural anddown-regulated in Mesenchymal tumors(Figure 1A). Twenty-five genes in thecluster were functionally annotated to be
involved in nervous system development(Figure 5). We selected these 25 genes asrepresentative Proneural genes. Selected
genes included previously reported
Proneural signature genes OLIG2,
NKX2‑2, SOX11, and ASCL15,6.
MicroRNAs contribute to tumorigenesis
by modulating both on oncogenic andtumor suppressor pathways in Glioma13,14.
We identified the miRNAs which havebeen reported to be targeting genes
involved in nervous system development,by querying the miRWalk database15;
miRNAs tagged as validated in themiRWalk database were considered fordownstream analysis.
Pearson pair-wise correlation, for all
subtypes, was performed between the
25 mRNAs involved in nervous systemdevelopment and the 38 validated
miRNAs, using the multi-omic correlation
analysis framework in GeneSpring.
A group of 159 GBM samples fromTCGA for which both mRNA and miRNA
measurements were available were used
for the analysis. Of the 38 miRNAs that
had significant differential expression(moderated t-test p ≤ 0.05) between theProneural and Mesenchymal subtypes,
miR-222, miR-221, miR-223, miR-34a,

To further understand the mRNA-miRNA
regulatory component within a given
GBM subtype, we selected the miRNAs
that displayed negative correlation in
Figure 6, and four reported Proneuralmarkers OLIG2, NKX2-2, SOX11, andASCL1. The correlation of the selected
mRNAs and miRNAs was performed
based on their expression similarity
within the Proneural and Mesenchymal
samples, respectively. The results ofthe correlation performed are shown in
Figure 7 and Table 2.

Proneural

Mesenchymal

Mesenchymal

Proneural
Condition
Subtype
PAFAH1B3
ZEB2
NKX2-2
MTSS1
CRMP1
MYT1
SOX11
DBN1
DPYSL4TNRC4CDK5R1
STMN1
EPHB1
VAX2
NRXN1
RUFY3
NRXN2
MAPT
OLIG2
CSPG5
ASCL1
CKB
BEX1
MPPED2
SH3GL2

-4.5

0

4.5

Figure 5. GO analysis of the genes up-regulated in Proneural and down-regulated in Mesenchymal tumors,against the entire set of genes, identified 25 genes that are involved in nervous system development.

7

0

1

hsa-miR-155
hsa-miR-27a
hsa-miR-29a
hsa-miR-34a
hsa-miR-223
hsa-miR-221
hsa-miR-222
hsa-miR-29c
hsa-miR-185
hsa-miR-184
hsa-miR-452
hsa-miR-10a
hsa-miR-143
hsa-miR-145
hsa-miR-31
hsa-miR-202
hsa-miR-205
hsa-miR-21
hsa-miR-187
hsa-miR-214
hsa-miR-23b
hsa-miR-27b
hsa-let-7a
hsa-miR-96
hsa-miR-182
hsa-miR-183
hsa-miR-9
hsa-miR-9*
hsa-miR-138
hsa-miR-153
hsa-miR-504
hsa-miR-181a
hsa-miR-181a*
hsa-miR-19a
hsa-miR-93
hsa-miR-19b
hsa-miR-106a
hsa-miR-20a

-0.7

ZEB2
SH3GL2
MPPED2
CSPG5
CKB
MAPT
ASCL1
TNRC4
CDK5R1
NRXN2
BEX1
NRXN1
RUFY3
EPHB1
MTSS1
NKX2-2
VAX2
STMN1
MYT1
DPYSL4
CRMP1
OLIG2
DBN1
PAFAH1B
SOX11

Figure 6. Pearson pair-wise correlation between the 25 genes involved in nervous system development and 38 miRNAs that had significant differential
expression (moderated t-test p ≤ 0.05) between the Proneural and Mesenchymal subtypes. Expanded view shows the second quadrant with miRNA on x-axis
and mRNAs on y-axis. Similarity is calculated based on the profile of the mRNAs and miRNAs across the four subtypes.

hsa-miR-31

hsa-miR-34a

hsa-miR-29a

ASCL1
hsa-miR-223

ASCL1
hsa-miR-27a

OLIG2

hsa-miR-222

OLIG2

hsa-miR-221

SOX11

hsa-miR-214

SOX11

hsa-miR-21

NKX2-2hsa-miR-155

NKX2-2hsa-miR-34a

hsa-miR-31

hsa-miR-29a

hsa-miR-223

hsa-miR-27a

hsa-miR-222

hsa-miR-221

hsa-miR-21

B

hsa-miR-214

hsa-miR-155

A

Figure 7. Correlation of OLIG2, NKX2-2, SOX11, and ASCL1 with miRNAs that showed up-regulation in
Mesenchymal tumors. Correlation in Proneural samples (A) and Mesenchymal samples (B).
Table 2. Observed correlation coefficients between the Proneural markers NKX2-2, SOX11, and OLIG2 and someof the miRNAs down-regulated in proneural samples.
OLIG2

NKX2-2

SOX11

Proneural

Mesenchymal

Proneural

Mesenchymal

Proneural

Mesenchymal

hsa-miR-21–0.278–0.0360.18–0.2720.168–0.322hsa-miR-214

0.193–0.2380.218–0.1060.2–0.014hsa-miR-221–0.124

0.12

0.139

–0.17

0.36

–0.34

hsa-miR-222–0.112

0.078

0.155

–0.113

0.276

–0.31

hsa-miR-223–0.149–0.236–0.162–0.178–0.127–0.23

hsa-miR-27a–0.1750.072–0.007–0.208–0.145–0.317hsa-miR-29a

–0.12

0.18

–0.093–0.216–0.267–0.41

hsa-miR-31–0.113–0.1950.101–0.1260.067–0.239hsa-miR-34a–0.3120.177–0.024–0.0690.091–0.288

8

2-T
8-T
9-T
6-T
10-T
7-T
1-T
3-T
4-T
9-E
8-E
2-E
6-E
10-E
5-E
1-E
7-E
3-E
4-E

B
12
Normalized intensity value

2
1
0
-1
-2

11
10
9
8
7

Mesenchymal

Proneural

Mesenchymal

Proneural

Figure 8. Box plot representation of (A) OLIG2 and (B) miR-34a expression in Proneural and Mesenchymal
samples.

Proteomics analysis

mRNA expression does not always
translate directly to protein abundancelevels. To validate genomic observations
in GBM at the protein level, we generated
an independent proteomics dataset.
The motivation was to identify GBM
signatures relevant at both RNA andprotein levels, and to confirm TCGA
observations in an independent cohort
using an independent technology.

Label-free proteomics analysis

Initial profiling analysis was performed
in MPP software on a label-freeproteomic dataset for 10 GBM tumorsand 10 controls (epileptic brain tissue)
samples. The protein database search
results identified a total of 14,187 uniqueprotein groups across the entire data set.

Expression profiling showed that oneof the tumor samples was close to thecontrol group samples, and therefore,that tumor sample was excluded from
subsequent analysis (data not shown).
Statistical analysis of the nine GBM
and 10 control samples identified 587
differentially expressed proteins (T-test,
corrected p-value ≤ 0.05 and FC ≥ 2.0).
These 587 proteins clearly discriminatedthe tumor from control samples. Thiswas confirmed by the sample-sample
correlation and PCA analysis (Figure 9A
and 9B). The PCA and correlation analysis
reflected clear differences between theprotein expression profile of the control
and tumor samples. Our analysis also
identified potential patterns of tumor
subgroups within the limited number oftumor samples which were studied.

B
20
Component 2 (9.06 %)

A

A

2-T
8-T
9-T
6-T
10-T
7-T
1-T
3-T
4-T
9-E
8-E
2-E
6-E
10-E
5-E
1-E
7-E
3-E
4-E

The oligodendrocytic development
regulator OLIG2 displayed negative
correlation to miR-34a, miR-21, andmiR-155 in the Proneural samples,but either positive or no correlationin the Mesenchymal samples. Earlier
reports identified TP53 mutationsand LOH events as characteristic to
Proneural subtype. Here we observed
increased expression of oligodendrocytic
development regulator OLIG2 anddecreased expression of miR‑34a in theProneural cohort (Figure 8). OLIG2 plays acentral role in directing cell fate decisions
and controlling cell proliferation and it
has been reported to repress the p53
tumor‑suppressor pathway effector p2118.
miR-34a is a tumor suppressor regulating
the p53 signaling pathway19,20. Thissuggests the possibility that functionalloss of TP53 in Proneural samples iscaused either by deletion of the gene orby its repression due to elevated levelsof OLIG2 and decreased levels of miR-34a.

3
Normalized intensity value

miR-21, known to induce proliferation andinvasiveness16, is positively correlated
to transcriptional activators NKX2‑2and SOX11 in Proneural samples
but displayed negative correlation in
Mesenchymal samples. Literature on
targets for miR-2117, and our query in themiRWalk database do not show NKX2-2
and SOX11 as targets for miR-21. Our
study indicates the possibility of potential
regulation of NKX2-2 and SOX11 by
miR‑21 in GBM tumors.

Control
Tumor-1
Tumor-2
Tumor-3

10
0
-10
-20
-40

-0.2

0.5

1

-20

0
20
Component 1 (54.4 %)40

Figure 9. Separation between control (-E) and tumor (-T) samples based on 587 proteins differentially expressed between control and tumor samples.
A) Sample-sample correlation; B) PCA.

9

The GeneSpring metadata framework wasused to verify if the observed subgrouppattern can be explained by available
metadata for the tumor samples. Figure
10 shows observed tumor subgroupsaligned to sample metadata such as age,sex, tumor location, and methylation
status of MGMT (known prognostic
marker for GBM). The metadata did
not clearly explain the sub-grouping
observed within the tumors based on the
differentially expressed proteins.

Combined analysis of proteomics
and gene expression

In the previous section, we described
separation between tumor and normal
samples and identification of the tumor
subgroups based on the abundanceof the 587 differentially expressed
proteins. Figures 12A and 12B show thesample‑sample correlation heat map
and PCA plot based on abundance ofthe 54 proteins. The same separation
between tumor and normal samples aswell as the tumor subgroups are clearlyvisible in these figures.

To explain the subgroups observed
within GBM tumors, we conducted
integrated analysis of the differentiallyexpressed proteins identified in our
proteomics study and the differentiallyexpressed genes identified by TCGA. Theintersection between the 587 differentiallyexpressed proteins (see previous section)and the 840 subtype‑specific genes
reported by Verhaak et al.6 contains
54 genes and the 54 proteins they encode(Figure 11). Using these 54 molecules,
we were able to confirm some of the key
findings of both studies.
-8.2

0

8.2

48

62.4

7-T
1-T
3-T
4-T
2-T
8-T
9-T
6-T

Group
Control
Tumor-1
Tumor-2

Sex
Female
Male

MGMT
None

Location
Left frontal

Right frontal

Right parietal

Methylated

Left temporal

Right media parietal

Right temporal

Unmethylated

None

Right occipetal

Right temporal

Tumor-3

Age
25.6

Survival (days)

Group
Age
Sex
Survival
MGMT
Location

10-T

576

0

Figure 10. Sample groups observed by clustering of tumor samples aligned to sample metadata: age, sex, survival in days, MGMT gene methylation status, andbrain location.

655 of the 840 subtype-specificgenes from TCGA measured inproteomics experiment587 proteins differentiallyexpressed from profiling
analysis. T-test p-value ~ 0.05
and FC ¡ 2.0

533

54

601

Figure 11. Venn diagram displays the overlap between 587 proteins from discovery differential analysis
and 840 subtype specific genes from TCGA; 655 of the 840 genes were measured in the proteomics
experiment.

10

To validate if the subgroups observedin this proteomics study are related
to the molecular subtypes reported
in GBM, genes corresponding to the
54 proteins were clustered based on
TCGA mRNA expression data. Figure 13
shows a gene‑sample clustering analysis
performed on the mRNA expression
profiles of the 54 genes. Figure 13
confirms that expression of the 54 genes
classify GBM tumors into four reported
molecular subtypes with a high degree ofcertainty.

9-E
8-E
1-E
5-E
7-E
3-E
2-E
6-E
4-E
10-E
3-T
7-T
1-T
4-T
2-T
6-T
10-T
8-T
9-T

Conclusions

Comparing data from several -omics
technologies revealed previously
unknown information about the biology
of GBM tumors. For example, theTCGA project identified 840 genes
capable of classifying GBM tumors infour known subtypes. Our proteomics
data identified 587 proteins whichwere differentially expressed in tumor
versus normal tissues and indicated anability to classify tumors based only on
their protein expression. By comparing
genomics and proteomics datasets in the
GeneSpring/MPP correlation framework,
we identified a core signature comprised

-0.2

B

0.5

1

8

Control
Tumor-1
Tumor-2
Tumor-3

6
Component 2 (10.1 %)GeneSpring/MPP software offers a suiteof powerful data analysis tools enablinga thorough investigation of complex
biological datasets. In this case study,
we performed multidimensional analysis
of GBM tumors leveraging publishedand newly generated genomics andproteomics datasets. Using innovative
metadata and correlation frameworks
now available in GeneSpring/MPP
software, we reproduced and expanded
key findings of a large-scale TCGA study.
For example, in a single GeneSpringexperiment we were able to confirm thegenomics subtype classification of GBM
tumors reported by TCGA. Furthermore,while somatic deletions of key genes in
GBM were reported by TCGA, GeneSpring
metadata framework revealed a higher
than expected rate of simultaneousdeletions of FGFR2 and PTEN in thesame samples. This important finding
suggests that these two genes may
become aberrant due to a single
genomic rearrangement rather than two
independent deletion events.

9-E
8-E
1-E
5-E
7-E
3-E
2-E
6-E
4-E
10-E
3-T
7-T
1-T
4-T
2-T
6-T
10-T
8-T
9-T

A

4
2
0
-2
-4
-6
-8

-15

-10

-5
0
5
Component 1 (62.38 %)10

15

Figure 12. Separation between control (-E) and tumor (-T) samples based on the 54 proteins
identified using a translation of gene signatures from the TCGA mRNA expression experiment inAgilent GeneSpring. A) Sample-sample correlation; B) PCA.

Classical

NeuralMesenchymal

Proneural-70Figure 13. Reported molecular subtypes in GBM established by clustering 173 TCGA tumor samples
based on the expression of 54 mRNA signatures identified using the proteomics study.

117

of only 54 genes and proteins. The coresignature demonstrated classification
power comparable to both the 840 genes
and the 587 proteins, highlighting thepower of multi-omics data analysis in
GeneSpring. Our observation suggeststhat heterogeneity of GBM tumors may
be best characterized by both mRNAand protein abundance of the signature
genes, opening a possibility for follow‑up
studies with potential significance beyond
the research scope.

Acknowledgement
The authors wish to thank members ofthe Human Toxome Mapping Consortiumfor sharing data, feedback, and helpful
discussions during development of
the Agilent GeneSpring metadata
and correlation frameworks. MR
was supported in part by NIH grant
1R01ES020750-01 “MAPPING THE
HUMAN TOXOME BY SYSTEMSTOXICOLOGY” (PI T. Hartung). SPM was
supported by the James C. Benjamin
Fund for Brain Tumor Research, theFroedtert Foundation Grant, and theAdvancing a Healthier Wisconsin Clinical and Translational Science
Institute Pilot Research Grant.

References
1. Kleihues, P; Sobin, L..H; World HealthOrganization classification of tumors.
2000, Volume 88, Issue 12, p 2887.
2. Brandes, A. A; et al. Temozolomide
as a second-line systemic regimen
in recurrent high-grade glioma: A
phase II study. Ann. Oncol. 2001,12,
pp 255‑257.
3. Martinez, R; Rohde, V; Schackert, G.
Different molecular patterns inglioblastoma multiforme subtypes
upon recurrence. J. Neurooncol. 201096, pp 321-329.

4. Cancer Genome Atlas Research
Network. Comprehensive genomic
characterization defines humanglioblastoma genes and core
pathways. Nature 2008, 455,
pp 1061‑1068.
5. Phillips, H. S; et al. Molecularsubclasses of high‑grade glioma
predict prognosis, delineate a patternof disease progression, and resemble
stages in neurogenesis. Cancer Cell
2006, 9, pp 157‑173.

13. Ventura, A; Jacks, T. MicroRNAs andcancer: short RNAs go a long way. Cell
2009, 136(4), pp 586–591.
14. Calin, G. A; Croce, C. M. MicroRNAsignatures in human cancers. Nat. Rev.
Cancer 2006, 6(11), pp 857‑866.
15. Dweep, H; et al. miRWalk-database:prediction of possible miRNA bindingsites by “walking” the genes of 3
genomes. Journal of Biomedical
Informatics 2011, 44, pp 839-47.

6. Verhaak, R. G; et al. Integrated
genomic analysis identifies clinically
relevant subtypes of glioblastoma
characterized by abnormalities in
PDGFRA, IDH1, EGFR, and NF1. Cancer
Cell 2010, 17, pp 98-110.

16. Møller, H. G; et al. A Systematic
Review of MicroRNA in Glioblastoma
Multiforme: Micro-modulators in the
Mesenchymal Mode of Migration andInvasion. Mol. Neurobiol. 2013, 47,
pp 131-144.

7. Noushmehr, H; et al. Identification of a
CpG island methylator phenotype thatdefines a distinct subgroup of glioma.
Cancer Cell 2010, 17, pp 510–522.

17. Buscaglia, L. E; Li, Y. Apoptosis andthe target genes of microRNA-21.
Chin. J. Cancer. Jun 2011, 30(6),pp 371‑380.

8. Brennan, C. W; et al. The Somatic
Genomic Landscape of Glioblastoma.
Cell Oct 10 2013, 155(2), pp 462-77.

18. Ligon, K. L; et al. Olig2-regulated
lineage‑restricted pathway controls
replication competence in neural stem
cells and malignant glioma. Neuron
2007, 53, pp 503-517.

9. Heroux, M. S; et al. Comprehensive
characterization of glioblastoma tumor
tissues for biomarker identification
using mass spectrometry-basedlabel-free quantitative proteomics.
Physiological Genomics 2014. 46,
pp 467-481.
10. Toedt, G; et al. Molecular signaturesclassify astrocytic gliomas by IDH1
mutation status. Int. J. Cancer. 2011,
128(5), pp 1095-1103.
11. PhosphoSitePlus: http://www.
phosphosite.org/proteinAction.
do?id=1262&showAllSites=false
12. Klinz, F. J; et al. In Genetic Diversity of
Glioblastoma Multiforme: Impact on
Future Therapies, Glioma - Exploring
Its Biology and Practical Relevance;
Dr. Anirban Ghosh, Ed.; 2011;
ISBN: 978-953-307-379-8.

19. Hermeking, H. p53 Enters theMicroRNA World. 2007, Volume 12,
Issue 5, pp 414-418.
20. Silber, J; et al. miR-34a repression
in proneural malignant gliomas
upregulates expression of its target
PDGFRA and promotes tumorigenesis.
PLoS One 2012, 7(3):e33844.

www.agilent.com/chem
This information is subject to change without notice.
For Research Use Only.
Not for use in diagnostic procedures.
PR7000-0097© Agilent Technologies, Inc., 2016
Published in the USA, January 4, 2016
5991-5505EN