Comprehensive Analysis of Radiomic Datasets by RadAR

Radiomics refers to the process of extracting and analyzing data from diagnostic images such as CT, MRI, and PET techniques (1).

In the context of precision oncology, radiomics is emerging as a promising approach to improve patients' stratification and support clinical decisions. Recent works include the use of radiomics to inform about genomic alterations (2, 3), predict gene expressions patterns (4), and genomic heterogeneity (5). Other studies reported on radiomic signatures able to predict outcome in independent cohorts of patients with cancer (6, 7).

The typical workflow of radiomic analysis includes image acquisition, segmentation, feature extraction, and analysis of high-dimensional datasets (8–10). In the past years, a variety of tools to address the first steps of the analysis have been proposed, including software for image processing and visualization (11, 12) and the computation of large sets of texture parameters (radiomic features; refs. 13, 14), even though not providing radiomics capabilities beyond feature extraction. The recent development of public clinical images databases, such as The Cancer Imaging Archive (TCIA; ref. 15), allows researchers to test, implement and validate radiomics-based analyses. However, despite the availability of software and data repositories, the processing and analysis of high-dimensional radiomic datasets remains challenging due to the lack of specific tools.

To address this issue, we developed RadAR (Radiomics Analysis with R), a new package for the statistical language R (https://www.R-project.org/) for the analysis of radiomic datasets. RadAR includes functions allowing both beginners and intermediate/advanced users to perform comprehensive analysis of radiomic datasets, including data import/output, feature data processing, statistics and visualization. To showcase RadAR, we applied our tool to multiple publicly available datasets (7, 16), including CT images of more than 850 patients with non–small cell lung and head and neck cancer and demonstrated that RadAR may represent a valuable tool for the radiomics community.

The Radiomics Analysis with R (RadAR) tool was developed as a package for the statistical language R. RadAR consists of five main modules, including import, preprocessing, statistics, visualization, and output (Fig. 1). The RadAR R package along with a detailed step-by-step tutorial on the processing of radiomic datasets is freely available under MIT license at https://github.com/cgplab/RadAR.

The main functionalities of RadAR are described in the following sections.

The import module includes a set of functions to import features extracted by widely used feature extractor software, including LifeX (11), PyRadiomics (14), and 3DSlicer (12). Features are automatically annotated using a manually curated dictionary with the aim of facilitating the interpretation of radiomic features in downstream analysis.

Also, this module allows the user to import features extracted by other feature extraction tools. In this case, data should be provided in the form of generic table where rows and columns represent radiomic features and region of interests (ROI) of patients, respectively.

When a radiomic dataset is imported, an R object of the class SummarizedExperiment (https://bioconductor.org/packages/release/bioc/html/SummarizedExperiment.html) is automatically generated. The SummarizedExperiment class is used to store rectangular matrices of experimental results and was originally developed to deal with sequencing and microarray experiments. Each object stores observations of multiple patients' ROIs along with additional meta-data describing both the observations (features) and patients (i.e., clinical and molecular data). We termed this object as RadAR object.

This module allows the user to apply different preprocessing techniques to the radiomic features, such as outlier identification, data scaling, and data normalization.

• - Outlier identification implements outlier detection through the interquartile range rule. Patients showing a user selectable fraction of outliers can be therefore identified and filtered out. This procedure is useful to improve the homogeneity of the dataset and therefore improve downstream signature development.

• - Data scaling procedure makes features comparable in terms of their range values. This procedure does not affect the distribution of features. RadAR implements two methods for data scaling, including min-max and median subtraction.

• - Data normalization makes the distribution of radiomic features identical in terms of statistical properties. To do that, we used quantile normalization (17).

The statistics module is the core module of RadAR. It implements a set of statistical techniques to perform comprehensive analysis of radiomic datasets.

To perform feature selection, RadAR implements the following methods.

• 1. Minimum-redundancy-maximum-relevance (mRMR) is a well-established method to perform feature selection (18). This is a supervised feature selection method that therefore requires a response variable (e.g., patient outcome, response to therapy).

• 2. A method based on hierarchical clustering (HC) analysis of the features' correlation matrix. In this strategy, the redundancy between radiomic features is estimated by correlation matrix. Then, a HC analysis of correlation matrix enables to identify groups of highly similar features. User can decide the number of features to retain and the HC method. Finally, a user-selectable criterium is applied to select representative features. Using concordance, the features showing best concordance index within each group are selected. This is a supervised approach as it requires a response variable. Using variability, the features showing the largest dynamic range within each group are selected.

• 3. A method based on k-means clustering of principal component analysis of the radiomic dataset. Dimensionality reduction of principal components is performed by the selection of the most informative components based on cumulative proportion of explained variance criterion (user-selectable parameter). Finally, a k-means clustering is applied to the multidimensional informative components space to identify a user-selectable number of groups of similar features. The selection of representative features within each group is performed using the aforementioned criteria (i.e., concordance or variability).

• 4. Binomial and Cox regression model fitting methods using the GLMNET package (19, 20).

This procedure is useful when patients' response variables are not available. A widely used unsupervised analysis is cluster analysis, which is used to identify groups of patients showing similar radiomic profiles. RadAR implements unsupervised analysis based on HC.

We define differential radiomics as the use of supervised statistical analysis with the goal of identifying statistically different radiomic features between groups of patients. The groups of patients can be defined on the basis of different clinical (e.g., responders vs. nonresponders), molecular (e.g., presence vs. absence of a certain gene mutation) or histopathologic (e.g., presence vs. absence of the expression of a certain cell receptor) conditions. For this statistical procedure, RadAR implements three nonparametric tests: the Wilcoxon-Mann-Whitney (WMW) and the area under the curve (AUC) for two groups comparisons and the Kruskal–Wallis test for comparisons between three or more groups.

A specific function enables to test the performance of a radiomic signature in independent dataset(s) (obtained for instance by upstream feature selection procedure) using Cox regression model.

This module implements functions to facilitate the visualization of radiomic features and entire radiomic datasets. The functions include (i) the generation of box plots and strip charts showing the distribution of radiomic features across different conditions; (ii) the visualization of feature correlation matrix by heatmap or correlation plot; (iii) the visualization of the output of the unsupervised HC analysis through heatmap, dendrogram, and annotation tracks.

This module enables to export and save out of R different kinds of products of the RadAR analysis, including text tables reporting top-ranked features resulting from differential radiomic analyses, figures and the RadAR R object that includes all the information related to the radiomic dataset and its processing with RadAR.

To demonstrate the reliability of RadAR for the analysis of radiomic datasets, we considered three CT cancer datasets (7, 16) available through TCIA (15). The collection includes images and clinical data of 422 patients with non–small cell lung cancer (NSCLC) from the NSCLC-Radiomics dataset (Lung1, https://doi.org/10.7937/K9/TCIA.2015.PF0M9REI), 137 patients with head and neck (HN) from the Head-Neck-Radiomics-HN1 dataset (HN1, https://doi.org/10.7937/tcia.2019.8kap372n) and 606 head and neck patients from the Oropharyngeal Carcinoma dataset (OPC, https://doi.org/10.7937/tcia.2019.8dho2gls). Manual contouring of the gross tumor volume was available for a total of 857 patients, including 388 from Lung1, 118 from HN1 and 351 from OPC. No additional segmentations were performed. DICOM to nrrd conversion was performed using the SlicerRT extension (version 0.18) of 3D Slicer (version 4.10.2). PyRadiomics version 2.2.0 using default parameters was used to perform radiomic feature extraction. For all the datasets, a total of n = 851 features were extracted from original images (n = 107) and 8 wavelet filtered images (n = 93 each).

A typical radiomic workflow is here reported starting from the import of radiomic features, identification of features associated with prognosis, application of unsupervised analysis, radiomic signature development and validation. Similarly to Aerts and colleagues, 2014 (7), we used the Lung1 as training and the HN1 and OPC as validation datasets. Finally, a comparative study of the radiomic signatures generated using the methods implemented in RadAR and the original signature from Aerts and colleagues was carried out. The code along with detailed comments to reproduce the proposed analyses is reported as Supplementary Data.

The .csv files generated by PyRadiomics were imported into RadAR using the import_pyradiomics function. The degree of redundancy of the Lung1 dataset was evaluated by correlation matrix visualization (Fig. 2A) using the plot_correlation_matrix function.

To identify radiomic features associated with patients' outcome, the Wilcoxon-Mann-Whitney test was applied to dead versus alive patients using the calc_differential_radiomics function. A table reporting the top-ranked radiomic features was generated using the select_top_features function (Supplementary Table S1) and four top features were visualized as box plots with overlaid strip charts using the plot_features function (Supplementary Figs. S1 and S2). To facilitate the interpretation of radiomic features across different conditions, quantile normalized values can be used for visualization (Supplementary Fig. S3). The features better discriminating different stages of the disease can be identified using the Kruskal–Wallis test implemented in the calc_differential_radiomics function (Supplementary Table S2; Supplementary Fig. S4).

To identify groups of patients showing similar radiomic profiles, unsupervised analysis of the Lung1 dataset was performed: as first steps, distance matrix computation (using Spearman correlation for patients and Euclidean distance for features) and Ward HC were computed using the do_hierarchical_clustering function. Then, the result of the clustering, along with annotation tracks reporting clinical and histologic data was visualized through the plot_heatmap_hcl function. To investigate whether unsupervised analysis of radiomic features could inform on patients' survival in this dataset, we labeled patients with their cluster membership (assuming n = 2 clusters) by the find_clusters function and estimated statistical significance using a Fisher exact test. Doing that, we found an enrichment of dead events in one of the two clusters [P = 1e-2, OR = 1.8; Supplementary Fig. S5]. Notably, we observed a marked improvement in patients' stratification when the same strategy was applied to quantile normalized data (Fig. 2B and C; P = 3e-3; OR = 2.0).

Finally, we show how to build a radiomic signature and evaluate it in two independent datasets, including HN1 and OPC. A comparative study of the radiomic signatures generated using the 7 methods implemented in RadAR and the original signature from ref. 7 is reported (Fig. 2D). This analysis shows that all the radiomic signatures obtained good performance, with 5 out of 8 methods obtaining Concordance Index (CI) > 0.60 in both the validation datasets.

RadAR is, to our knowledge, the first software specifically designed for the processing and statistical analysis of radiomic datasets. The functions implemented in RadAR have been extensively documented to allow both beginners and intermediate/advanced users to exploit all the functionalities of our software. The analysis of three large datasets demonstrated the reliability of RadAR to perform comprehensive analysis of radiomics datasets. Also, we showed that unsupervised analysis of the Lung1 dataset without the application of upstream feature selection procedures was able to identify clusters of patients with NSCLC reflecting different clinical outcome. Finally, we performed a comparative analysis of different radiomic signatures identified by RadAR, showing its capability for radiomics studies. To conclude, RadAR may represent a valuable resource for the radiomics community, with also potential for facilitating the application of this emerging field in the clinical setting.

M. Benelli reports personal fees from Novartis outside the submitted work. No potential conflicts of interest were disclosed by the other authors.

M. Benelli: Conceptualization, resources, data curation, software, supervision, funding acquisition, methodology, writing-original draft, project administration, writing-review and editing. A. Barucci: Data curation, validation, writing-original draft, writing-review and editing. N. Zoppetti: Resources, data curation, software, writing-review and editing. S. Calusi: resources, data curation, validation, writing-original draft, writing-review and editing. L. Redapi: Data curation, validation, writing-review and editing. G. Della Gala: Data curation, validation, writing-review and editing. S. Piffer: Data curation, validation, writing-review and editing. L. Bernardi: project administration, writing-review and editing. F. Fusi: Project administration, writing-review and editing. S. Pallotta: Supervision, writing-original draft, writing-review and editing.

This study was supported by Fondazione Cassa di Risparmio di Firenze (to M. Benelli and A. Barucci).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.