|Home | About | Journals | Submit | Contact Us | Français|
Individual MRI characteristics (e.g., volume) are routinely used to identify survival-associated phenotypes for glioblastoma (GBM). This study investigated whether combinations of MRI features can also stratify survival. Furthermore, the molecular differences between phenotype-induced groups were investigated.
Ninety-two patients with imaging, molecular, and survival data from the TCGA (The Cancer Genome Atlas)-GBM collection were included in this study. For combinatorial phenotype analysis, hierarchical clustering was used. Groups were defined based on a cutpoint obtained via tree-based partitioning. Furthermore, differential expression analysis of microRNA (miRNA) and mRNA expression data was performed using GenePattern Suite. Functional analysis of the resulting genes and miRNAs was performed using Ingenuity Pathway Analysis. Pathway analysis was performed using Gene Set Enrichment Analysis.
Clustering analysis reveals that image-based grouping of the patients is driven by 3 features: volume-class, hemorrhage, and T1/FLAIR-envelope ratio. A combination of these features stratifies survival in a statistically significant manner. A cutpoint analysis yields a significant survival difference in the training set (median survival difference: 12 months, p = 0.004) as well as a validation set (p = 0.0001). Specifically, a low value for any of these 3 features indicates favorable survival characteristics. Differential expression analysis between cutpoint-induced groups suggests that several immune-associated (natural killer cell activity, T-cell lymphocyte differentiation) and metabolism-associated (mitochondrial activity, oxidative phosphorylation) pathways underlie the transition of this phenotype. Integrating data for mRNA and miRNA suggests the roles of several genes regulating proliferation and invasion.
A 3-way combination of MRI phenotypes may be capable of stratifying survival in GBM. Examination of molecular processes associated with groups created by this combinatorial phenotype suggests the role of biological processes associated with growth and invasion characteristics.
Despite aggressive treatment, the prognosis for patients diagnosed with glioblastoma (GBM) remains poor, with a median survival time of 12–15 months.32 Surgery never completely removes the tumor due to widespread infiltration of tumor cells, leading to recurrence. Tumor samples are used to establish diagnosis, but the rest of the tumor is unavailable for pathologic or genomic interrogation. The nonuniformity of GBM may explain why it is so difficult to control with conventional therapies.
Imaging provides the only nondestructive means to evaluate the morphology of the entire tumor. Imaging features potentially represent a phenotypic characterization of the genetics and biology of the tumor.4 This study seeks combinations of imaging features that are predictive of tumor behavior and therapeutic response.
Identifying noninvasive surrogates of clinical outcome and associated molecular differences might suggest the etiology of the underlying disease state.4,29 Several studies13,24,38 indicate that MRI phenotypes might potentially serve as noninvasive proxies of cancer-associated molecular processes. Although the standard of care for treating GBM is relatively nonspecific (i.e., resection, radiotherapy, and temozolomide33), personalized treatment may have different efficacies among GBM subtypes.25 Better patient stratification using clinical imaging methodology may offer a more patient-specific treatment. Molecular information concerning the GBM is only available from surgical specimens. Thus, a noninvasive radiographic method to determine molecular and clinical characteristics may provide prognostic data to help guide the frequency of surveillance imaging for tumor progression or recurrence, possibly leading to earlier and more effective therapy.
The Cancer Genome Atlas (TCGA) contains comprehensive genomic data (e.g., expression, copy number, methylation, and microRNA [miRNA]). The Cancer Imaging Archive (TCIA) is a publicly available repository of patient-derived images across multiple modalities (SPECT, CT, and MRI) and tumor types. The TCIA images are from patients who have corresponding molecular data in TCGA. Using genomic and imaging data from the TCGA-GBM archive, we performed a computational analysis to identify radiological feature combinations that might stratify survival. Two phenotype groups were defined based on these feature combinations. Differential analysis of mRNA and miRNA expression was performed to identify molecular differences between these groups.
We showed that a combinatorial phenotype comprising 3 radiological characteristics (volume-class, T1/FLAIR ratio, and hemorrhage) can stratify survival. The median survival difference between the less aggressive phenotype and the more aggressive phenotype was 8 months. Molecular examination of these 2 groups suggests that genes and miRNAs related to growth, invasion, and proliferation underlie this phenotypic difference.
Patients were included in the study based on the availability of MRI data from TCIA and were selected to ensure consistency in glioma history and treatment profile. Specifically, patients with no prior glioma and those who received standard of care (concurrent chemotherapy and radiotherapy after resection) were selected, for a total of 92 cases. The clinical data for these untreated, primary GBM cases were obtained from the cBioPortal for Cancer Genomics (Memorial Sloan Kettering Cancer Center; http://www.cbioportal.org/public-portal/). To assess stratification of survival by the combinatorial phenotype, we divided data into a training set (44 patients; 14 females, 30 males) and a validation set (48 patients; 18 females, 30 males), based on median survival and vital (alive vs dead) status. We also verified that patients in these 2 sets had similar demographic and functional (i.e., Karnofsky Performance Scale score) characteristics.
Radiology annotations for the 92 patients with GBM and their preoperative images (from TCIA) were obtained from the TCGA Glioma Phenotype Research Group (https://wiki.cancerimagingarchive.net/display/Public/TCGA+Glioma+Phenotype+Research+Group) based on the VASARI (Visually AcceSAble Rembrandt Images) standardized feature set.9 This contains annotations produced via consensus reads by a team of radiologists affiliated with the TCGA Glioma Phenotype Research Group. See details at https://wiki.nci.nih.gov/display/CIP/VASARI and https://wiki.cancerimagingarchive.net/display/Public/VASARI+Research+Project. Gutman et al.9 discuss details of specific imaging protocols and the semantic scoring system used by the radiologists.
VASARI characteristics studied were volume-class, enhancement quality, proportion contrast-enhancing tumor, proportion noncontrast-enhancing tumor, proportion necrosis, proportion edema, T1/FLAIR-envelope ratio, enhancing margin thickness, distribution, and hemorrhage. This subset corresponded to the imaging-specific rather than anatomical (location) subset of the VASARI features. Several of these imaging characteristics were described as independently prognostic in previous studies10,18,26 and most are qualitative (categorical data), whereas some (i.e., major and minor axis lengths) are quantitative (continuous-valued). For consistency, all variables were made categorical.
Instead of using the length variables separately, we computed an “approximate” 3D volume using the formula (4π*a2*b/3)/8, with “a” being the major axis length and “b” being the minor axis length. These measurements were derived from the region of T2 signal abnormality. For consistency (i.e., to make all variables categorical), we converted the volume into an ordinal scale based on its magnitude with respect to the median volume across all cases (1 = ≤ median, 2 = > median). This dichotomized version of the volume is referred to as “volume-class.”
The T1/FLAIR ratio attribute compares the relative size of the abnormality in the T1-weighted images to the size in the T2-FLAIR images. It can have 3 categories: expansive (size of precontrast T1 abnormality region approximately equal to FLAIR envelope), mixed (T1 abnormality less than FLAIR abnormality), and infiltrative (T1 region significantly smaller than FLAIR region). Hemorrhage represents intrinsic hemorrhage in the tumor matrix (observed in T1-weighted images or T2-weighted images).
Level 3 expression data for both mRNA (Affymetrix U133Av2 BI Platform: Affymetrix HT HG U133A) and miRNA (UNC Agilent Human miRNA 8 × 15 K) were downloaded from the TCGA data portal. These data are normalized and were used for differential expression analysis.
Patient-specific survival data came from the cBioPortal for Cancer Genomics. Clinical data are available at http://www.cbioportal.org/public-portal/study.do?cancer_study_id=gbm_tcga. We extracted “Overall Survival (months)” and “Overall Survival Status” for each patient. Further demographic information (sex, age, prior glioma status, and treatment history) was also obtained from the cBioPortal (or from the TCGA portal [http://cancergenome.nih.gov/]).
To examine whether these image features stratify the cases in any clinically relevant manner, we hierarchically clustered them (using R package “cluster”) using Gower’s similarity metric with complete linkage. Gower’s distance is applicable for mixed (combination of discrete, interval-scaled, and categorical) data7 and is defined as
where dijk is the distance between individuals i and j based on variable k. wijk is a weight given to the ijk comparison, taking values of “0” for an invalid comparison and “1” for a valid comparison.
Viewing the cluster plot allowed identification of imaging variables with a range consistent within a cluster that also changed across the cluster boundary. This led to “cluster blocks,” within which the variable is approximately constant and between which there is a change. This allowed creation of a combinatorial phenotype via element-wise multiplication of the individual imaging variables. For example, a case with “high” volume (Category 2), “low” T1/FLAIR ratio (Category 1), and “no” hemorrhage (Category 1) will take the value 2 (2*1*1) for its combinatorial phenotype.
Using survival data from the 92 patients, we used a k-adaptive partitioning scheme (R package “kaps”) to estimate a cutpoint on the combinatorial phenotype that induces a statistically significant difference in survival between the 2 groups. This cutpoint is initially estimated on the training set (44 cases) and validated on a validation set (remaining 48 cases absent from the training set). The survival difference between the 2 cutpoint-induced groups was estimated using a log-rank test. The split value obtained in the training set was also used to compare differences in progression-free survival (PFS) via the log-rank test. The performance of the combinatorial phenotype relative to individual variables in the combination is assessed via area under the survival receiver operating characteristic (ROC) curve. Area under the curve (AUC) can take a value between 0 (poor) and 1 (perfect discrimination). An AUC > 0.5 suggests predictive ability better than random chance.
Cutpoint analysis (described above) identified a split value on the combinatorial radiophenotype that partitioned the data into 2 distinct survival groups. Next, we looked for molecular differences (differential expression of mRNA and miRNA) between these 2 groups, using the Comparative Marker Selection module within the GenePattern Suite (Broad Institute).
Normalized Level 3 expression data for both mRNA and miRNA were downloaded from the TCGA data portal. These data were used for differential expression analysis via the Comparative Marker Selection module within the publicly available GenePattern platform (http://genepattern.broadinstitute.org/). The algorithm uses a 2-sided t-test to identify genes/miRNAs differentially expressed between the 2 phenotype classes. These classes are induced by the split value (“2”) obtained from the partitioning algorithm used above. Cases with “volume-class:T1/FLAIR:hemorrhage” combination values greater than “2” are designated Group 1. The rest are designated Group 0.
Functional analysis (i.e., pathway analysis) of the differentially expressed genes and miRNAs was done using Gene Set Enrichment Analysis (GSEA; Broad Institute). We used the GSEA desktop version (http://www.broadinstitute.org/gsea/) with the settings “t-test” metric and “equalized and balanced” randomization during 1000 permutations. The t-test is used for consistency with the metric used for differential expression analysis (from GenePattern). Multiple testing correction for significance was done using false discovery rate (FDR) computation.1
For relating differentially expressed miRNAs with differentially expressed mRNAs, the miRNA target filter feature was used within Ingenuity Pathway Analysis (IPA) software (Ingenuity Systems) to find miRNAs targeting the differentially expressed genes. This integrative miRNA:mRNA analysis looks for target relationships between miRNAs and the genes derived from differential expression analysis. Also, the miRNA target filter in IPA looks for concordant changes in expression (i.e., anticorrelated expression changes in miRNA:mRNA abundance).
Core analysis and functional analysis were performed on the gene and miRNA lists. We explored “Direct Interactions” using “Experimentally Observed or High Confidence Predictions” in the IPA Knowledge Base (“Genes only”).
Study of the hierarchical clustering plot (based on Gower’s distance and complete linkage; Fig. 1) suggested 3 categorical variables were influencing clustering structure: volume-class, T1/FLAIR ratio, and hemorrhage. Next, the combinatorial phenotype was derived as an element-wise product of the categorical values taken by each of the 3 variables (volume-class, T1/FLAIR ratio, and hemorrhage). For example, if volume-class is Category 1, T1/FLAIR ratio is Category 2, and hemorrhage is Category 2, then the combined phenotype takes a value equal to 1*2*2 (i.e., 4).
The association of this combinatorial phenotype with survival was examined via cutpoint analysis on the 92 cases that included complete survival data. Using a k-adaptive partitioning scheme (R package “kaps”), we found that a split value of 2 induced a statistically significant survival difference on the training cohort (44 cases, chosen randomly) (see Kaplan-Meier plot, Fig. 2 upper). The low survival group had a median survival of approximately 7 months and the high survival group had a median survival of approximately 19 months. This survival difference was confirmed in a validation set (the remaining 48 of 92 cases) (Fig. 2 lower). This split value partitioned the data into 2 groups (cases with volume-class:T1/FLAIR:hemorrhage combination ≤ 2, and those with volume-class:T1/FLAIR:hemorrhage combination > 2).
The group defined by a split value of “2” was also examined for differences in PFS. The p value of the log-rank test was 0.0577, suggesting a trend toward PFS stratification by this combinatorial phenotype.
The predictive performance (measured using area under the ROC curve, AUC) of the combinatorial phenotype is superior relative to the individual variables (volume-class, T1/FLAIR ratio, and hemorrhage). The AUCs for the combinatorial phenotypes and the individual variables, for 4 different survival times (12, 15, 24, and 36 months), are represented graphically in Fig. 3.
Application of the Comparative Marker Selection procedure to the mRNA and miRNA data yielded 384 genes and 23 miRNAs showing differential expression between the 2 phenotypic classes. Tables 1–3 show the ontological analysis of these 384 genes (using the IPA tool). The key pathways activated in this set are related to cell-to-cell signaling and interaction, cellular assembly and organization, mitochondrial dysfunction, oxidative phosphorylation (OXPHOS), and so on.
To find pathways that were overrepresented between these 2 radiophenotype classes, we used GSEA. Several pathways (gene sets) associated with metabolism and immune function were found to be significant at an FDR cutoff of 0.1 (Table 4).
Using the IPA tool, we combined the miRNA expression results with the mRNA expression results and found that at least 8 differentially expressed miRNAs (miR-214–3p, miR-499–5p, miR-129–5p, miR-433–3p, miR-485–5p, miR-614, miR-617, and miR-637) target several differentially expressed genes while exhibiting concordant expression correlation (lowering of mRNA expression and elevation of miRNA expression; Table 5).
Several potentially interesting networks were indicated. The top-scoring canonical pathways pertained to immunological signaling (iCOS-iCOSL signaling in T helper cells, RXR activation, and PI3K signaling in B lymphocytes), mitochondrial dysfunction, and OXPHOS. The primary biological processes (based on Fisher exact tests for overrepresentation analysis) were “cell-to-cell signaling and interaction” and “cellular assembly and organization” (Table 1).
In the present study, we show that the combinatorial radiographic phenotype (element-wise product), comprising volume-class, T1/FLAIR-envelope ratio, and hemorrhage, stratifies survival in GBM in a statistically significant manner. It also shows a trend toward stratifying PFS in a statistically significant manner.
To summarize, the volume is derived based on T2 abnormality. The T1/FLAIR ratio attribute compares the relative size of the abnormality in the T1-weighted images to the size in the T2/FLAIR images. Hemorrhage represents intrinsic hemorrhage in the tumor matrix (observed high signal in T1-weighted images or low signal in T2-weighted images). Prior work in this domain has focused on the identification of single image variables (i.e., contrast-enhancing tumor,9 nonenhancing tumor17) as prognostic factors and identified molecular processes that might underlie such disease-associated phenotypes (spanning processes such as invasion38 and hypoxia4). Independently, each of these imaging characteristics may be relevant to patient outcome, but together they become a much more powerful prognostic indicator. Patients with a particular imaging profile may be at risk for earlier recurrence of disease and poorer survival. Because of their clinical relevance, these imaging characteristics should be more objectively defined. For example, volumetric quantification is rarely, if ever, reported beyond a general description of size. Although the role of volume in stratifying survival is expected, we observed that combining this with T1/FLAIR ratio and hemorrhage has predictive value. The T1/FLAIR ratio is related to extent of tumor infiltration. Hemorrhage is potentially a histological manifestation (presumably related to extent of necrosis or microvascular proliferation). Examination of this combinatorial phenotype’s association with survival leads to observations that might have clinical consequences, e.g., low-volume, infiltrating disease portends better outcome than high-volume, moderately infiltrating disease; or low-volume, moderately infiltrating disease is worse than low-volume disease with hemorrhage. Based on the molecular and phenotypic correlates of this combined radiophenotype, we hypothesize that it might correlate with the aggressiveness of the tumor.
We sought to obtain a molecular characterization of the groups induced by the combinatorial phenotype, rather than based on survival characteristics alone. Examination of the miRNAs and mRNAs in Table 5 suggests the functional role of several of these molecules in glioma biology (spanning proliferation, invasion, and growth characteristics). Among the miRNAs, miR-499 has been shown to promote tumor metastasis and cellular invasion.20 miR-199a is known to regulate cellular proliferation and survival.31 miR-429 inhibits cellular growth, invasion, and epithelial-to-mesenchymal transition (EMT)34 and is seen to be upregulated in the less aggressive phenotype.
Furthermore, hsa-miR-146b and hsa-miR-193a are associated with survival in GBM.32 hsa-miR-146b was shown to be involved in glioma cell migration and invasion.36 Also, miR-499 regulates apoptosis, mitochondrial dynamics, and represses CNTNAP2, PPP3CB, and SYT1 (invasion-associated genes upregulated in the more aggressive phenotype). Upregulation of miR-146b is observed in the less aggressive phenotype (Group 0), concordant with the observation that it reduces migration and invasiveness in glioma.15,36 miR-125 and miR-129 were associated with tumor progression and glioma proliferation.19 Although other genes/miRNAs in Table 5 indicate phenotypic modulators of interest, their specific roles in GBM remain to be elucidated. Low expression of miR-200 family members has been linked to the epithelial-to-mesenchymal phenotype and subsequent metastatic disease.11,16,35
Several of the pathways obtained via GSEA point to interesting aspects of GBM biology (Table 4). An examination of pathways enriched at the top of the differentially expressed gene lists indicates the potential involvement of several immune-modulatory pathways underlying this phenotypic transition: natural killer cell activity, interferon gamma response, and T-cell lymphocyte differentiation. Additionally, EMT is associated with the more aggressive phenotype, an association further strengthened by the enrichment of the mesenchymal subtype signature in the same group.37 Molecules such as PAR1 and HOXC6 are implicated in tumor invasiveness.12 Furthermore, 20q13 region amplification is known to be relevant to glioma risk and GBM.2,17,23,27 However, an in-depth validation in a reliable in vivo model is crucial before drawing any conclusions pertinent to this specific combinatorial phenotype.
Similarly, several pathways, mostly related to cellular metabolism, are overrepresented in the differential enrichment list, e.g., PGC1alpha activity is associated with OXPHOS,5 which is involved in brain tumor progression.8 These sources all point to the potential association of metabolism-related activity with the transition of this combinatorial phenotype. Moreover, OXPHOS has been suggested to confer tumorigenic potential to GBM cancer stem cells.14 This finding is also concordant with that of altered metabolic activity underlying tumor-invasive characteristics in GBM.8
IPA indicated that the coordinated activity of several miRNAs and mRNAs might underlie this phenotype. Specifically, we observed the presence of several immune-modulatory genes and key signaling molecules. However, determining the roles of these molecules in modulating this phenotypic transition will require very careful study.
Molecular analysis of this combinatorial radiophenotype thus suggests the role of proliferation/invasion characteristics underlying its transition. Recent literature suggests that targeting invasion/proliferation mechanisms may hold therapeutic promise.3,6,21,22,28 Thus, this combinatorial radiophenotype might also be used as a screening tool for assessing response to drugs that target glioma cell proliferation and invasion characteristics. A potential avenue might also be the surveillance of this combinatorial phenotype to assess treatment response, thereby serving as a complement to existing response evaluation criteria (such as MacDonald, Response Evaluation Criteria In Solid Tumors [RECIST] criteria). In addition, because this phenotype has a survival correlate that is presumably related to the aggressiveness of the disease, modulating this phenotype might have consequences for treatment decision making. Indeed, there are several approaches targeting invasion,37 proliferation, EMT, and metabolic31 characteristics for GBM. The role of immune-modulatory processes underlying this combinatorial phenotype also suggests the possibility of targeting disease through immune modulation25 and for assessing response to immune-targeted therapeutics using this phenotype.
Our work has several limitations, providing abundant scope for further investigation. An unbiased search for other feature combinations (within the VASARI set) that can stratify survival might reveal other interesting characteristics. Beyond indicating interesting combinations of features, an integrative molecular analysis comprising mRNA, miRNA, methylation, and copy number from these phenotype-induced groups might identify other molecular processes relevant to glioma biology. This study focused exclusively on the association between combinatorial radiographic phenotypes and survival. In the future, assessment of the phenotype’s clinical value will require a multivariate model including clinical factors (e.g., age, Karnofsky Performance Scale score) and known molecular predictors such as MGMT methylation and molecular subtype.
Our results suggest that combining radiographic phenotypes provides additional variables that stratify survival. Specifically, we showed that combination of multiple MR image–derived features has more predictive value for survival than some individual ones. Because volume-class is important in this combinatorial predictor, exploring precise 3D volume measurement based on a 3D segmentation markup of the lesion data may be useful. Currently, only major and minor axis lengths are available. Further, it would be interesting to explore quantitative assessments of the hemorrhage component’s volume and T1/FLAIR ratio as additional predictors for stratifying survival, rather than the current categorical measurements.
This project has been funded in whole or in part with federal funds from the National Cancer Institute, National Institutes of Health, under Contract No. HHSN261200800001E. The content of this publication does not necessarily reflect the views or policies of the Department of Health and Human Services, nor does mention of trade names, commercial products, or organizations imply endorsement by the U.S. Government. We would like to thank Dr. Chad Holder from the TCGA Phenotype Research Group for guidance with data interpretation. We are grateful to Dr. Keith Baggerly and Dr. Rehan Akbani for discussions about the TCGA methylation and expression data. We also thank Dr. Rehan Akbani for helpful discussions about clustering categorical data. We are also grateful to Ms. Priyanka Raja and Mr. Suvindra Vijayan for assistance with coding the radiological attributes. In addition, Ms. Anisha Luthra assisted with the miRNA:mRNA analysis. Support of the TCGA Glioma Phenotype Research Group is gratefully acknowledged. We would also like to thank Ms. Audria Patrick and Mr. David Wildrick for editorial assistance. We gratefully acknowledge the support of institutional start-up funds, an MD Anderson Cancer Center institutional research grant (IRG), support from MD Anderson Cancer Center Support Grant P30 CA016672 (the Bioinformatics Shared Resource), and a Career Development Award from the Brain Tumor SPORE (granted to A.R.) for conducting this research.
Mr. Kirby has direct stock ownership in Myriad Genetics.
Author ContributionsConception and design: A Rao, G Rao. Acquisition of data: A Rao, G Rao, Gutman, Flanders, Hwang, Rubin, Colen, Zinn, Jain, Wintermark. Analysis and interpretation of data: A Rao, G Rao. Drafting the article: A Rao, G Rao. Critically revising the article: A Rao, G Rao, Flanders, Hwang, Rubin, Colen, Zinn, Jain, Wintermark, Kirby, Freymann, Jaffe. Reviewed submitted version of manuscript: all authors. Approved the final version of the manuscript on behalf of all authors: A Rao. Statistical analysis: A Rao. Administrative/technical/material support: A Rao, G Rao, Kirby, Freymann, Jaffe. Study supervision: A Rao, G Rao.
Dr. Wintermark: Department of Radiology, Neuroradiology Section, Stanford Medical Center, Stanford, CA.