IDH mutation-specific radiomic signature in lower-grade gliomas

Demographic and clinical characteristics

A total of 158 patients diagnosed with LGG were enrolled as the training data set. Of these, 118 (74.7%) had IDH mutant (IDH^MUT) tumors and 40 (25.3%) had IDH wild type (IDH^WT) tumors. No significant differences were observed with respect to age, sex, WHO grade, and tumor location between the two groups. The IDH mutation rate in the validation data set was 75.5% (77 out of 102). The clinical and pathological characteristics of the training and the validation data sets are listed in Table 1.

Identification of LGGs with similar radiomic patterns

To assess the radiomic expression patterns, the quantitative radiomic features were extracted from the LGG patients in the training set. An unsupervised hierarchical clustering method with average linkage revealed two clusters of patients with similar radiomic expression patterns (Figure 1). By comparing the clinical parameters of the two clusters, we found that the second cluster was significantly associated with a high frequency of the IDH mutation (P = 0.0020, Fisher’s exact test, Figure 1), which indicates a tight association between IDH mutation status and quantitative radiomic features.

Identification of the IDH-mutation specific radiomic signature

Based on previous observations, our goal was to identify a set of radiomic features that would enable the pre diction of the IDH mutation status in LGGs. We first screened the differences in the radiomic features between the IDH^WT and IDH^MUT LGGs using the SAM algorithm. A total of 146 features, including 5 first order statistics features (energy, entropy, mean, median and root mean square), 1 shape and size based feature (surface to volume ratio), 8 textural features (GLCM: contrast, dissimilarity, energy, entropy, difference entropy, informational measure of correlation 1; GLRLM: short run emphasis, run percentage) and 132 wavelet features (Figure 2A), were found to be expressed differentially. All the 146 screened features were listed in Supplementary Table 1.

Then, we utilized the logistic regression algorithm to select the IDH mutation-specific signature. A series of ROC curves with relevant AUC was delineated using the given parameters. The results showed that only ten radiomic features were needed to divide 158 LGGs into an IDH^WT group and an IDH^MUT group with an AUC of 0.86, while the AUC was 0.92 for 20 features and 0.98 for 50 features (Figure 2B and 2C). Notably, using a classification model of 86 radiomic features, the enrolled patients were correctly classified into the IDH^WT and the IDH^MUT groups (AUC = 1.00, Figure 2B and 2C).

To further validate the classification of this radiomic signature, T2-weighted images from 102 patients and in vivo xenograft glioma models were subjected to feature extraction and logistic regression. Consistently, the results suggested that this signature, which was comprised of 86 radiomic features, could separate the 102 LGG patients into two groups with high sensitivity and specificity (AUC=0.99, Figure 2D). Fifty-four differentially expressed radiomic features were found in xenograft model, and 40 of them (such as Energy, Mean, Median and Root mean square) were shared between patients and xenograft model (Supplementary Table 2). Moreover, by using 40 of the differential features, we could divide the experimental mice into two groups in accordance with their IDH phenotype (Figure 2E and 2F). These in vivo experiments further indicated the robustness of radiomic features in differentiating IDH mutation status. Specifically, the immune associated features were also listed in Supplementary Table 2.

The potential molecular mechanism of IDH-mutation-specific radiomic features

CGGA and TCGA mRNA sequencing data were utilized to identify the biological processes and signaling pathways that differed significantly between the IDH^MUT and IDH^WT LGGs. Separate gene annotations of the CGGA data and the TCGA data consistently indicated that they primarily involved these biological processes: the immune response, cell adhesion, and vascular development (Supplementary Figure 1 and Supplementary Figure 2). In addition, the TCGA data suggested that the cell cycle phase and cell proliferation were also involved in the radiological manifestation of the IDH^WT LGGs (Supplementary Figure 2).

Meanwhile, to attempt to determine the exact biological significance of each radiomic feature, 48 LGGs samples with both transcriptomic and radiomics data were subjected to correlation analysis and gene annotation. The results indicated that the IDH-specific radiomic features, such as the surface to volume ratio (SVR) and entropy, were primarily associated with cell polarity, cell adhesion, cell growth, and immune processes, which were the same biological processes that were found to differ between the IDH^MUT and IDH^WT LGGs (Supplementary Figures 3–5). Specifically, LGGs with large SVRs were found to be positively associated with a high expression of oncogenes, such as OTX1, HES1, BAG5, and the top 10 genes that were positively associated with SVR were listed in Supplementary Table 3.

Identification of a prognostic-based radiomic signature

Accumulating evidence has revealed that the IDH mutation is a crucial predictor for LGG patient outcomes. The present study also showed that the IDH^MUT LGG patients had a longer OS than the IDH^WT patients (Figure 3B). To further explore the prognostic ability of the IDH-specific radiomic features, we extracted a compact signature consisting of features with P < 0.05 after a univariate Cox regression analysis (Table 2). The β value of each significant radiomic feature was used for risk evaluation. When they had an elevated risk score, the patients were prone to have a higher mortality and a greater frequency of IDH mutation (Figure 3A).

Using the calculated risk score, we divided the LGG patients into high-risk and low-risk groups. Patients with a high-risk score had a worse prognosis (Figure 3C, HR = 0.3207, P = 0.0017). Intriguingly, IDH^MUT patients could be further categorized into IDH^MUT-high risk and IDH^MUT-low risk groups with significantly different OSs (Figure 3D, HR = 0.4232, P = 0.0483). A multivariate Cox analysis also demonstrated that the radiomic risk score could serve as a prognostic indicator for LGGs (Table 3). Moreover, the OS was not significantly different between the IDH^WT–high risk groups and the IDH^MUT-high risk groups (Figure 3E, P = 0.199), further emphasizing the prognostic value of the IDH-specific radiomic signature.

A radiogenomic analysis found that a high-risk score was positively associated with genes that included ERCC1, G6PD, SOX9, and EGLN2, which are primarily enriched during the regulation of programmed cell death, cell growth, and metabolic processes. This could partially account for the radiological malignancy of the high-risk group (Figure 4). Representative samples of T2-weighted images with relevant radiomics and clinical features are presented in Figure 5. The first case was a 39-year old male patient with IDH mutant LGG. This patient was classified into the IDH^MUT group with a relatively low radiomics risk score. Case 2 was a 46-year-old male with the IDH wildtype LGG who was correctly classified into the IDH^WT group with a high-risk score.

Additionally, the current study revealed that the radiomic risk score was significantly positively correlated with glioma stem cell markers such as TWIST1 (R = 0.503, P < 0.001) and CD133 (R = 0.346, P = 0.016) (Supplementary Figure 6).

Patients and samples

A total of 158 LGG patients who underwent surgical treatment between March 2006 and December 2012 were retrospectively selected as a training set. An additional 102 LGG patients who underwent surgical treatment from January 2016 to June 2016 were prospectively included as a validation set. The inclusion criteria were as follows: (1) pathologically confirmed lower grade glioma; (2) available IDH status; (3) available high-resolution preoperative T2-weighted MR images; (4) available clinical characteristics. The basic clinical characteristics of the patients, including age, gender, WHO grade, and tumor location are summarized in Table 1. The study was approved by the institutional review board of Beijing Tiantan Hospital. The design of the present study is illustrated in Figure 6.

Tumor masking and normalization

The tumor masking was conducted as previously described [44]. Most of the MR images were obtained using a Trio 3.0T scanner (Siemens, Erlangen, Germany), and the remaining clinical structural images were acquired on a Magnetom Verio 3T (Siemens AG, Erlangen, Germany). The T2-weighted image parameters were as follows: repetition time = 5800 ms; echo time = 110 ms; flip angle = 150 degrees; 24 slices; field of view = 240×188 mm²; voxel size = 0.6×0.6×5.0 mm³; matrix = 384×300. Tumors were traced directly using MRIcron (http://www.mccauslandcenter.sc.edu/mricro/mricron). Masks of the brain tumors were drawn on each patient’s T2-weighted images in native space by two board-certified neuroradiologists (K.W. and X.C.), who were blind to the patients’ clinical information. Areas that produced abnormal hyperintense signals on the T2-weighted images were identified as LGG tumor areas. When a greater than 5% discrepancy existed between these two masks, the masks utilized were determined by a senior neuroradiologist (S.L.). The intensities of the voxels in each tumor mask were normalized to the z distribution ([gray value – mean gray value] / SD), on the individual level, to ensure consistency in the distribution of the gray values among the cases in the cohort.

Quantitative radiological features extraction and selection

The extraction of the radiomic features was performed as previously reported [26]. A total of 431 image features were extracted from the tumor masks. The features were categorized into four groups. Group 1 quantitatively described the distribution of voxel intensities in the MR image using 14 descriptors. Group 2 described 8 three-dimensional features based on the shape of the tumor regions. Group 3 described textual features for quantifying intra-tumor heterogeneity, which were calculated from gray level co-occurrence (GLCM, 22 descriptors) and gray level run-length (GLRLM, 11 descriptors) texture matrices. Group 4 calculated the intensity and textural features from wavelet decompositions of the original image. All the algorithms were implemented in MATLAB (2014a).

Transcriptomic comparison between IDH^WT and IDH^MUT LGG

All available whole genome mRNA sequencing data of LGG patients and clinical information were acquired from the Chinese Glioma Genome Atlas database (http://www.cgga.org.cn) and the Cancer Genome Atlas database (http://cancergenome.nih.gov/). The differentially expressed genes were selected by a significance analysis of the microarray (SAM) algorithm using R programming language (http://cran.r-project.org), and with the criteria of fold change > 20% and false discovery rate (FDR) < 0.05 (Benjamini-Hochberg). Finally, the genes significantly overexpressed in the IDH^WT group (CGGA, 1509; TCGA, 1262) or in the IDH^MUT group (CGGA, 914; TCGA, 582) were processed using Gene Cluster and Gene Treeview software to construct a heatmap, and the online Database for Annotation, Visualization, and Integrated Discovery (DAVID, http://david.ncifcrf.gov/) program [45] for the gene ontology (GO).

Xenograft model of glioma and radiomic analysis

Lentiviral vectors carrying IDH1 wildtype or IDH1 R132H mutant cDNA sequences were transduced into the U87MG cells with polybrene (Sigma), as previously described [42]. Stably transduced cells were selected during three days of puromycin (Sigma) treatment. The Western Blot was conducted to verify the expression of the mutant IDH1 enzyme using a cell lysate and the IDH1^R132H antibody (1:200, DIA-H05, Dianova). Then, IDH1 wild type or R132H mutant U87MG cells (5 × 10⁵ cells per mouse in 5 µL) were intracranially injected into 5 to 6 weeks old female nude mice (Beijing Vital River Laboratory Animal Technology), as described earlier [46]. After 40 days, the tumors were measured using a 7T MR Image System (Bioclinscan, Bruker). The relevant T2-weighted MR images were subjected to radiomic feature extraction and analysis. After tumor segmentation, a total of 431 radiomic feature were extracted from the T2-weighted MR images from each mouse. The differentially expressed radiomic features were selected using SAM algorithm, with the criteria of false discovery rate (FDR) < 0.05 (Benjamini-Hochberg).

IDH phenotype classification

First, similar to the radiomic analysis in xenograft model, SAM algorithm was conducted on the 431 radiomic features to select the differentially expressed radiomic features between IDH^WT and IDH^MUT tumors, with the criteria of false discovery rate (FDR) < 0.05 (Benjamini-Hochberg). Next, the selected differentially expressed radiomic features were utilized for IDH phenotype classification. The classification process was conducted using a logistic regression model (Y = expr _{feature 1} × β _{feature 1} + expr _{feature 2} × β _{feature 2} + … + expr _{feature n} × β _{feature n} + ε) based on the MATLAB (2014a) software. In this model, ‘Y’ represents the IDH status (1 indicates mutation, 0 indicates wildtype), while ‘expr’ represents the expression value of each radiomic feature. β represents the model parameter to be estimated, and ε is the estimated residual. The prediction results were further interpreted using the receiver operating characteristic (ROC) curve. To select the most predictive parameters, the logistic regression algorithm was applied repeatedly. Using the dimensionality reduction principle, the radiomic feature with the highest P value in predicting IDH status was excluded from the model each time until features with the best predictive effect were identified. The signature derived from the training set was subsequently applied to the validation set.

Identification of the prognosis-based signature

The prognostic values for each IDH-specific feature in patients with LGG were calculated using a univariate Cox regression after elimination of the patients with an OS of less than 30 days. The significant radiomic features with P values < 0.05 were selected to develop a prognosis-based signature. The risk score was calculated using the linear combination of the selected features weighted by the regression coefficient derived from the univariate Cox regression analysis (β), referring to previous studies [47, 48]. The risk score for the OS of each individual was calculated as follows: Risk score = expr _{feature 1} × β _{feature 1} + expr _{feature 2} × β _{feature 2} + … + expr _{feature n} × β _{feature n}. ‘Expr’ represents the expression value of each radiomic feature. The ‘expr’ and β used in methods 2.7 have no relationship with those used in methods 2.6.

We next divided the patients into high-risk and low-risk groups based on the risk score. The cutoff value was determined when the P value was the smallest in log-rank test in the training set, and then the cut-off value was fixed and applied to all Kaplan-Meier curve analyses.

Radiogenomic analysis

Forty-eight LGG patients involved in the radiomic analysis were subjected to an Agilent Whole Human Genome Array analysis according to the manufacturer’s instructions [49]. The data were acquired using the Agilent G2565BA Microarray Scanner System and Agilent Feature Extraction Software (version 9.1). The probe intensities were normalized using GeneSpring GX 11.0. The IDH mutations of the training data set were assessed by pyrosequencing, while the IDH status for the validation data set was identified by immunohistochemistry (IDH1^R132H, DIA-H05, Dianova).

Radiogenomic analysis was further performed. The Pearson correlation coefficients were calculated between the genes and IDH-associated radiomic features. The association was identified to be statistically significant when the absolute value of Pearson correlation coefficient was > 0.4 and the P value was < 0.05. Gene ontology (GO) analysis was conducted to investigate the underlying biological processes of the radiomic features based on the DAVID Bioinformatics Resources (http://david.ncifcrf.gov/). The top 200 positive/negative genes that were significantly associated with each feature were subjected to GO analysis to reveal the underlying biological processes of each feature. Using this method, the underlying biological processes of the risk score and the differentially expressed features in xenograft model were also investigated. In addition, the relationship between the radiomic risk score and glioma stem cell signatures was also assessed with using the Pearson correlation analysis.

Statistics

The significant differences between the two groups were estimated using a Student’s t-test. Chi-square and Fisher’s exact tests were used to compare the frequencies between the groups. The OS curves were plotted according to the Kaplan–Meier method, with the log-rank test applied for comparison. A Cox regression was used to determine the prognostic value of each radiomic feature for the OS in LGG patients. All the differences were considered statistically significant at the two-sided P < 0.05 level.