With treatments for solid cancers transitioning towards personalized therapy, mutational profiling for identifying target gene mutation status or gene amplification status using tissue biopsies is becoming status-quo for most cancers. However, a significant challenge in reliable assessment of gene mutation or amplification status, is the underlying genomic heterogeneity which makes it challenging to identify the “true” mutational status based on random tissue sampling . Multiple studies have shown that certain gene mutations (e.g. MGMT promoter methylation, EGFR) have varying expressions across different parts of the tumor or between primary or secondary metastatic sites [14, 13]. There is hence a need for developing “virtual biopsy” techniques on imaging that can comprehensively capture the gene mutation heterogeneity of solid tumors, and potentially assist in surgical navigation to identify sampling sites for biopsy targeting.
The field of radiogenomics has provided a surrogate mechanism to predict gene mutational status on routine imaging (i.e. MRI) by training machine-learning models. Most of the existing radiogenomic models fall in two categories: (1) deep learning[2, 7, 8]/radiomics [9, 11, 4], and (2) atlas-based probabilistic approaches [3, 1]. In the absence of co-localized biopsy sites on MRI, deep learning/radiomic approaches employ features from the entire tumor to predict the gene mutational status. In contrast, atlases-based approaches obtain the likelihood of the mutational status of driver genes such as MGMT and EGFR at different spatial locations by creating probabilistic radiographic atlases obtained from a large population. These population-based approaches however do not leverage any tumor-specific information in the model.
In this work, we present the first-attempt at creating “virtual biopsy” radiogenomic maps for predicting gene mutational status on MRI, by combining two complementary attributes that capture mutational heterogeneity at: (1) population-level via spatial-priors for presence or absence of mutation status () using probabilistic atlases from a retrospective cohort, and (2) local tumor-level by incorporating context-priors that capture mutational heterogeneity via radiomic attributes obtained from a stereotactically co-localized biopsy site within the tumor. The spatial and context priors are combined within a Least Absolute Shrinkage and Selection Operator (LASSO) regression model to obtain a per-voxel probability of the likelihood of increased expression of the gene mutation () at that location. The prediction probabilities obtained for every voxel are further improved using probabilistic pairwise Markov models. In this work, we evaluate these Spatial and Context Aware (SpACe) maps in the context of two problems in Glioblastoma (GBM): (1) predicting EGFR status (amplified (), non-amplified ()), and (2) predicting MGMT status (methylated (), non-methylated ()), from routine MRI scans. The pipeline of the entire workflow is illustrated in Figure 1.
We deﬁne an image scene as , where is a spatial grid of voxels , in a 3D space, . Each voxel, is associated with an intensity value . represents the co-localized biopsy location on MRI scans, such that . denotes the feature set obtained for every . For gene , defines mutated/methylated, while defines non-mutated/unmethylated.
2.2 Computing context-aware mutational heterogeneity from stereotactic biopsy locations ()
We define “context” as local heterogeneity attributes computed from the co-localized biopsy site on imaging, using radiomic features including Haralick features (capture image heterogeneity ), Gabor features (capture structural details at different orientations and scales), Laws (capture spots and ripples-like patterns), and CoLlAGe features (capture localized gradient orientation changes ). Specifically, for every , we extract a set of 3D radiomic features (i.e. Haralick, Gabor, Laws, CoLlAGe). We define , where is the type of feature family (e.g. Haralick, Gabor), and , where is the number of feature attributes for every feature family. Feature pruning is then conducted on the extracted features using Spearman’s correlation metric, to eliminate redundant features. The pruned “context-aware” features (152 for EGFR cohort, 149 for MGMT cohort) are finally aggregated into one feature descriptor .
2.3 Computing spatially-aware priors () for likelihood of gene mutation status () using probabilistic atlases
Using the lesion segmentation obtained for every patient in the training set, two different population atlases for gene are constructed using subjects that belong to either or . This is done to quantify the frequency of occurrence of every voxel across and , and compute voxel-wise probability values, , (). All scans need to first be registered to an isotropic reference atlas (i.e. MNI152; Montreal Neurological Institute). The intensity values are then averaged across across all the annotated binary images of all patients involved in the study. This means that for , two probability values from these two atlases could be obtained, that characterize the probability of a voxel being or . The 2 probability values () for every voxel are finally aggregated in the spatial feature descriptor .
2.4 Creating SpACe maps for predicting voxel-wise mutational heterogeneity in the tumor
In order to obtain a voxel-wise prediction of the gene mutation status, the context descriptor (), spatial descriptor (), age (), and gender () of every patient in the training set, are incorporated within a LASSO model 
. LASSO model is selected to obtain the probability score using a parsimonious feature set by utilizing its capability in reducing variance when shrinking features, while simulteneously not increasing the bias. We designed the LASSO model to perform regularization of feature parameters as follows:  = , where is the shrinked set of coefficients obtained after regularization, , and is the penalty term. The voxel-wise probability is then computed as the weighted sum of the selected features for the set of coefficients , as follows: , where , and
is the number of features selected by LASSO. After obtaining the probabilistic map for every, we incorporate probabilistic pairwise Markov models (PPMMs) to improve voxel-wise gene mutation prediction 
. PPMMs are adopted from Markov Random Fields, through formulating priors in terms of probability density functions, hence allowing the creation of more robust prediction models. The input to this model is the voxel-wise probability values obtained from LASSO model. Interaction between neighboring sites is then modelled, to improve voxel-wise probability scores, and to finally obtainmaps.
2.5 Applying SpACe maps on testing sets for predicting voxel-wise mutational heterogeneity within the tumor
The top features selected on the training set are applied to the entire tumor on the test set, for obtaining voxel-wise probabilities for predicting the mutation status. For the purpose of computing accuracy of our model, we predict the mutation status [, ] based on pooled probability values for an already known biopsy site, and compare the prediction with the known mutation status. As an additional qualitative analysis, we threshold the probability values obtained from the entire tumor (threshold obtained empirically), followed by connected component analysis and PPMM, to obtain 2-3 hot-spots of high probability mutation sites. These hot-spots prospectively could be used to drive surgical navigation as potential sites for biopsy localization.
3 Experimental Design
3.1 Data description and preprocessing
We employed a unique retrospective dataset of a total of 100 GBM patients who underwent CT-guided biopsy for disease confirmation, since surgical resection was not feasible (due to location or other clinical reasons) for these patients. Segmentation of the enhancing lesion was conducted by an experienced radiologist on the MR scans. The biopsy site was co-localized by co-registering CT images with the MRI scans, followed by expert evaluation for confirmation. All scans were then registered to an MNI152 atlas and then bias-corrected using N4 bias correction . These studies were then divided into two cohorts: (a) : EGFR amplified () versus non-amplified () studies, and (2) : MGMT methlated () versus unmethylated (). For , we had a total of 91 subjects of which 70 were used for training (35 amp, 35 non-amp), and the remaining 21 (6 , 15 ) were used for validation. For , of a total of 81 subjects, 60 (28 , 32 ) were used for training , while 21 subjects were used for validation (5 , 16 ).
3.2 Implementation details
Two experiments were set-up using cohorts (Experiment 1: versus ) and (Experiment 2: versus ), respectively. For both experiments, we extracted a total of 3D context-features for every (where
was a 1-cm diameter sphere in our case), including 1 raw feature, 8 gray features, 13 gradient features, 26 Haralick features, 64 Gabor features, 152 Laws features, and 52 CoLlAGe features, extracted using 2 window sizesand . These features are pruned as detailed in Section 2.2 to obtain . In addition, population atlases were constructed to quantify the frequency of occurrence of versus and versus as detailed in Section 2.3. Similarly, for both experiments, was created following 10 runs of 10-fold cross validation, as detailed in Section 2.4. The median value of the probabilities across voxels of all subjects was used as a threshold to determine for every voxel. Finally, majority voting was used to obtain the mutation status for every biopsy site.
3.3 Comparative strategies
In order to evaluate the efficacy of SpACe model, we compared our results with two state-of-the-art methods, radiomic-based, and deep-learning-based that employ features from the entire tumor to predict amplification/methylation status. For the radiomic- experiment (
, a total of 316 radiomic features were extracted from the entire enhancing lesion of every subject (same attributes that were extracted from biopsy sites), and the feature vector was constructed from the 4 statistics: median, variance, skewness, and kurtosis values that were computed for every feature across all voxels for this patient, for total of 1264 features. After feature pruning, 283 features were fed to the LASSO model to compute patient-wise scores that determined their gene status.
For the DL approach to predict the mutation status, we used a deep residual neural network as described in. ResNet has previously been used to predict EGFR and MGMT mutation in GBM and other cancers [7, 17]
. Specifically, patches of size 128x128 were sampled from the center of the selected MRI slices and augmented using horizontal flips and random rotations to enlarge the limited training data. Following patch sampling, we trained separate deep Res-Net networks with 18 layers (ResNet-18) for the two experiments. In order to train the networks on MRI scans, we used pre-trained model on ImageNet and performed transfer learning using the sampled patches from MRI scans. We selected ResNet-18 because it removes the vanishing gradient problem and the network has several layers containing composite function of operations such as batch normalization (BN), convolution (Conv), rectified linear units (ReLU) and pooling for non-linear transformation of the input. We trained each model for 25 epochs with dropout of 0.2 to avoid overfitting. Models with minimum loss were locked down to test the patches obtained from the test set.
4 Results and Discussion
4.1 Experiment 1: Determining EGFR amplification status
Using alone, training and testing accuracies were reported as 61.43% and 66.67% respectively. Combination of features with features into LASSO model yielded training and testing accuracies of 80% and 90.48% respectively, using 8 (1 raw, 1 gray, 2 gradient, 1 Haralick, 3 Gabor) and 2 features. This implies that incorporating improved the model’s performance, rather than using alone. Next, we evaluated the efficacy of including , , and clinical features (, ) into our model to predict the mutation status. Clinical features did not improve accuracy of the model. PPMMs were then employed, and successfully corrected the amplification status for 7 subjects from the training set, yielding final training and testing accuracies of 90% (1 amp, 6 non-amp subjects were misclassified out of 70) and 90.48% (2 non-amp subjects are misclassified out of 21) respectively. Results on 2 different patients are illustrated in Figure 2.
Using radiomic features from the entire tumor to predict mutation status yielded training and testing accuracies of 80% and 71.43% respectively. Further, the Res-Net model to predict EGFR status yielded training and testing accuracies of 74.28% and 65.52%, significantly underperforming in comparison to the SpACe model.
4.2 Experiment 2: Determining MGMT methylation status
Using alone, the model achieved training and testing accuracies of 76.67% and 57.14% respectively. When combining features with features into LASSO model, we obtained training and testing accuracies of 81.67% and 61.9% respectively, which implies that incorporating improved the model’s performance, rather than using alone. Next, using , , and clinical features, the model picked a set of 12 features that included 8 features; 1 gray, 3 Haralick, and 4 Gabor features, in addition to , , , and . This model yielded training and testing accuracies of 83.3% and 66.67% respectively. Applying PPMMs for predicting methylation status on these results corrected the mutation status for 3 training subjects as well as 1 testing subject, with final accuracies of 88.3% and 71.5% respectively. Results on 2 different patients are illustrated in Figure 3.
When using radiomic features from the entire tumor to predict methylation status, training and testing accuracies were 76.67% and 52.38%. In addition, the DL model that was trained to predict methylation status gave training and testing accuracies of 79.37% and 68.40%, suggesting that results obtained using SpACe maps outperformed both comparative approaches.
5 Concluding Remarks
In this work, we presented the first-attempt at creating “virtual biopsy” radiogenomic maps for predicting gene mutational status on MRI, by combining two complementary attributes: (1) spatial-priors for presence or absence of mutation status via probabilistic atlases from a retrospective cohort, and (2) context-priors to capture mutational heterogeneity using radiomic attributes obtained from a stereotactically co-localized biopsy site within the tumor. These spatial-and-context aware (SpACe) maps were evaluated in the context of two experiments: predicting (1) EGFR amplification status, and (2) MGMT amplification status, on Glioblastoma. Our results demonstrated that SpACe outperformed state-of-the-art radiomic and deep learning approaches that were performed on the entire tumor, instead of learning features from the co-localized biopsy site. The virtual biopsy maps created using SpACe could not only improve prediction of gene mutation status of the tumor, but could also serve as surgical navigation to guide potential biopsy sites for specific gene mutations.
-  (2016) Population-based mri atlases of spatial distribution are specific to patient and tumor characteristics in glioblastoma. NeuroImage: Clinical 12, pp. 34–40. Cited by: §1.
-  (2018) . American Journal of Neuroradiology 39 (7), pp. 1201–1207. Cited by: §1.
-  (2013) Probabilistic radiographic atlas of glioblastoma phenotypes. American Journal of neuroradiology 34 (3), pp. 533–540. Cited by: §1.
-  (2019) Defining egfr amplification status for clinical trial inclusion. Neuro-oncology 21 (10), pp. 1263–1272. Cited by: §1.
-  (1973) Textural features for image classification. IEEE Transactions on systems, man, and cybernetics (6), pp. 610–621. Cited by: §2.2.
-  (2002) Discriminating vital tumor from necrotic tissue in human glioblastoma tissue samples by raman spectroscopy. Laboratory Investigation 82 (10), pp. 1265–1277. Cited by: §1.
-  (2017) Residual deep convolutional neural network predicts mgmt methylation status. Journal of digital imaging 30 (5), pp. 622–628. Cited by: §1, §3.3.
-  (2017) Deep learning based radiomics (dlr) and its usage in noninvasive idh1 prediction for low grade glioma. Scientific reports 7 (1), pp. 1–11. Cited by: §1.
-  (2018) Multiregional radiomics features from multiparametric mri for prediction of mgmt methylation status in glioblastoma multiforme: a multicentre study. European radiology 28 (9), pp. 3640–3650. Cited by: §1.
-  (2010) High-throughput detection of prostate cancer in histological sections using probabilistic pairwise markov models. Medical image analysis 14 (4), pp. 617–629. Cited by: §2.4.
-  (2016) Intratumoral heterogeneity identified at the epigenetic, genetic and transcriptional level in glioblastoma. Scientific reports 6, pp. 22477. Cited by: §1.
-  (2016) Co-occurrence of local anisotropic gradient orientations (collage): a new radiomics descriptor. Scientific reports 6, pp. 37241. Cited by: §2.2.
-  (2012) MGMT expression and promoter methylation status may depend on the site of surgical sample collection within glioblastoma: a possible pitfall in stratification of patients?. Journal of neuro-oncology 106 (1), pp. 33–41. Cited by: §1.
-  (2017) Intratumoral heterogeneity: pathways to treatment resistance and relapse in human glioblastoma. Annals of Oncology 28 (7), pp. 1448–1456. Cited by: §1.
-  (1996) Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological) 58 (1), pp. 267–288. Cited by: §2.4.
-  (2010) N4ITK: improved n3 bias correction. IEEE transactions on medical imaging 29 (6), pp. 1310–1320. Cited by: §3.1.
-  (2019) Implementation strategy of a cnn model affects the performance of ct assessment of egfr mutation status in lung cancer patients. IEEE Access 7, pp. 64583–64591. Cited by: §3.3.