knitr::opts_chunk$set(echo = TRUE)

Introduction

The quantitative analysis of biomedical images, referred to as Radiomics, is emerging as a promising approach to facilitate clinical decisions and improve patients’ stratification. The typical radiomic workflow includes image acquisition, segmentation, feature extraction and analysis of high-dimensional datasets. While procedures for primary radiomic analyses have been established during the last years, the processing of the resulting radiomic datasets remains challenging due to the lack of specific tools.

Here, we present RadAR (Radiomics Analysis with R), a new software to perform comprehensive analysis of radiomic features. RadAR allows the users to carry out the entire processing of radiomic datasets, from data import to feature processing and visualization and implements multiple statistical methods for the analysis of these data. We used RadAR to analyse the radiomic profiles of more than 850 cancer patients from publicly available datasets and showed that it was able to recapitulate expected results, demonstrating its reliability and proving that RadAR may represent a valuable tool for the radiomic community.

Here we show how to perform analysis of radiomic features by RadAR (Radiomics Analysis with R).

Installation

RadAR is freely available under GPL-3 license at https://github.com/cgplab/RadAR. First, install biocViews to facilitate the installation of the package dependecies:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("biocViews")

Then, you can install RadAR from github with:

# install.packages("devtools")
devtools::install_github("cgplab/RadAR")

To include vignettes use:

# install.packages("devtools")
devtools::install_github("cgplab/RadAR", build_vignettes = T)

Importing radiomic features

RadAR is able to import radiomic features extracted using different software, including pyradiomics (van Griethuysen et al., 2017), 3DSlicer (Fedorov et al., 2012) and LifeX (Nioche et al., 2018).

Importing features from Pyradiomics

To work with files generated by pyradiomics, use

rdr <- import_pyradiomics(dir = "PATH/TO/FOLDER/")

where "PATH/TO/FOLDER/" is the path to the folder containing the pyradiomics csv files. This function takes as input one file per patient.

Importing features from 3DSlicer

To work with files generated by 3DSlicer, use

rdr <- import_3dslicer(dir = "PATH/TO/FOLDER/")

where "PATH/TO/FOLDER/" is the path to the folder containing the pyradiomics tsv files. This function takes as input one file per patient.

Importing features from LifeX

To work with files generated by LifeX, there are two options: 1) xls of csv files (one per patient) or 2) session file includind multiple patients:

rdr <- import_lifex(dir = "PATH/TO/FOLDER/")

rdr <- import_lifex_session(dir = "PATH/TO/SESSION_FILE")

where "PATH/TO/FOLDER/" is the path to the folder containing the LifeX xls or csv files. This function takes as input one file per patient. "PATH/TO/SESSION_FILE" is the path to the LifeX session file including radiomic features from multiple patients.

Importing features (generic radiomic table)

RadAR allows users to import generic tables of radiomic features (rows) from multiple patients and/or regions of interest (columns). To do that, use

rdr <- import_radiomic_table(dir = "PATH/TO/RADIOMIC_TABLE")

where "PATH/TO/RADIOMIC_TABLE" is the path to the radiomic table including features from multiple patients. This is a tab delimited file reporting features (rows) for each patient (columns). First row should report the patient ids, first column should report feature ids, first cell should be empty.

Getting started: analysis of the NSCLC-Radiomics dataset

To show how to use RadAR for the analysis of radiomic datasets, we will apply it to the analysis of the NSCLC-Radiomics from Aerts et al. (2014).
CT Images, available RTstruct and clinical data for N=388 patients were downloaded from The Cancer Imaging Archive (TCIA) portal (https://www.cancerimagingarchive.net). Features (N=851) were extracted using pyradiomics v2.2.0.

The pre-processed pyradiomics results of the NSCLC-Radiomics datasets can be downloaded from here. This archive contains 388 .csv files generated by pyradiomics and the clinical data.
In R (we suggest using RStudio interface, see https://rstudio.com), the RadAR R package can be loaded using the following code

## load the RadAR package
library(RadAR)

To import Lung1 dataset into R, the import_pyradiomics function should be used.

## import pyradiomics data
rdr <- import_pyradiomics(dir = "path/to/NSCLC-Radiomics dataset/pyradiomics_out/")
data(package = "RadAR")
rdr <- lung1_rdr
clinical <- lung1_clinical
clinical$deadstatus.event <- as.character(clinical$deadstatus.event)

The function import_pyradiomics imports the csv files generated using pyradiomics and outputs rdr, a R object of the class SummarizedExperiment containing all the information related to the Lung1 radiomic dataset. The function filters out duplicated feature names. After duplicated feature names removal, the dataset includes 388 samples (columns) x 797 features (rows). The following code shows how to access to the information included in the rdr object.

rdr

The available clinical data can be also imported in the rdr object and used for downstream analyses.

## read the clinical file
clinical <- read.csv("path/to/NSCLC-Radiomics dataset/clinical.csv", header = T)
## assign rownames to clinical
rownames(clinical) <- clinical$PatientID
## order patients in clinical as in rdr
clinical <- clinical[colnames(rdr), ]
## print first lines of clinical data
head (clinical)

The clinical data includes the information summarized here.

head (clinical)
## add clinical data to the rdr object
colData(rdr) <- cbind(colData(rdr), clinical)
## print summary of rdr
rdr

Compared to the previous rdr, now colData (meta-data for patients) includes 13 elements (3 original + 10 from clinical).

Datasets included in RadAR

To facilitate the reproducibility of our analyses, the pre-processed rdr objects of Lung1 (lung1_rdr) and HN1 (hn1_rdr) datasets, along with available clinical information (lung1_clinical and hn1_clinical), are included in the RadAR package. These data are automatically loaded with RadAR.

Working with a RadAR object of the class SummarizedExperiment

Data included in rdr

Let's start giving a look at the information included in the Lung1 dataset. First, we can check which image types (i.e., normal or wavelet filtered) features have been extracted from.

## print available image types
print_image_type(rdr = rdr)

To use only a subset of features, the filter_by_image_type function can be used.

## filter by image types
rdr_filt <- filter_by_image_type(rdr = rdr, image_type = c("original"))
rdr_filt

rdr_filt now contains only 101 features. The rdr object can be also filter by feature type.

## print available feature types
print_feature_type(rdr = rdr)

To use only a subset of feature types, the filter_by_feature_type function can be used.

## filter rdr by feature types
rdr_filt2 <- filter_by_feature_type(rdr = rdr, 
                                   feature_type = c("first_order_shape",
                                                    "first_order_statistics"))
rdr_filt2

rdr_filt2 now contains 212 features.

To access to the elements of patients (columns) or features (rows) the following code can be used.

## print available elements for samples (columns)
names(colData(rdr))

## print available elements for features (rows)
names(rowData(rdr))

## how many male and females in my dataset?
table(colData(rdr)$gender)

## how many features per image types?
table(rowData(rdr)$image_type)

A manually curated dictionary of radiomic features

To facilitate the interpretation of the results of radiomic analysis, RadAR includes a manually curated dictionary of radiomic features. The description of each radiomic feature is included in the feature_description element of rowData(rdr). Dictionary is automatically added to the RadAR object when import_3dslicer(), import_lifex() or import_pyradiomics() are used to import the data. To facilitate the access to the dictionary, RadAR implements a dedicated function (features_to_dictionary).

## print dictionary for a selection of radiomic features 
my_features <- c("Elongation", "Flatness", "Maximum2DDiameterRow", "Sphericity")
out <- features_to_dictionary(rdr = rdr, 
                              feature_names = my_features)
out

Visualize the correlation matrix

One of the first steps of the analysis of a radiomic dataset is the study of the correlation matrix to evaluate the degree of redundancy of feature data in the dataset. The following code shows how to visualize the correlation matrix first by heat map and then correlation plot (to facilitate the visualization, only radiomic features extracted from unfiltered images are considered).

## view correlation matrix of features by heatmap
plot_correlation_matrix(rdr = rdr_filt, 
                        method_correlation = "spearman", 
                        view_as = "heatmap", 
                        which_data = "normal")
## view correlation matrix of features by heatmap
plot_correlation_matrix(rdr = rdr_filt, 
                        method_correlation = "spearman", 
                        view_as = "heatmap", 
                        which_data = "normal", 
                        cex_labels = .3)

Here, values are computed using Spearman's correlation and hierarchical clustering generated using Ward's. Different methods for both calculating correlation coefficients and generating dendrograms can be used (digit ?plot_correlation_matrix for more details). The option which_data is used to select which data should be used to generate the plot (original values, scaled or normalized values, see below for more details). The object rdr_filt is the rdr reduced to features related to only original, unfiltered, images and obtained through the upstream application of the filter_by_image_type function.

A different way to visualize the correlation matrix is by using correlation plot, that can be generated using the following code

## view correlation matrix of features by correlation plot
plot_correlation_matrix(rdr = rdr_filt, 
                        method_correlation = "pearson", 
                        view_as = "corrplot", 
                        which_data = "normal")
## view correlation matrix of features by correlation plot
plot_correlation_matrix(rdr = rdr_filt, 
                        method_correlation = "pearson", 
                        view_as = "corrplot", 
                        which_data = "normal", cex_labels = .3)

To save these plots to a file, the save_plot_to option should be used.

## view correlation matrix of features by correlation plot
plot_correlation_matrix(rdr = rdr_filt, 
                        method_correlation = "pearson", 
                        view_as = "corrplot", 
                        which_data = "normal",
                        save_plot_to = "myplot.pdf")

Pre-processing and visualization of radiomic features

Radiomic features are differently distributed and take values in extremely different ranges. To facilitate the processing and visualization of radiomic features, RadAR implements data processing strategies described here, such as scaling and normalization.

Using scale_feature_values() function, features are scaled by min-max or median strategies. This processing facilitates the visualization of features without affecting the distribution of data. This is useful when comparison or visualization of multiple features is desired. In the example below, we plot the distribution of the four features included in the original Aerts et al. radiomic signature across the patients' vital status (dead vs alive) retrieved from available clinical data. As the feature "compactness" had been deprecated in PyRadiomics, we considered "sphericity" (defined as the cube of compactness).

First, normal (unprocessed) feature values are used, as shown in the following box

## select some features to plot
aerts_sign <- c("Energy.original",
                "GrayLevelNonUniformity.original",
                "GrayLevelNonUniformity.wavelet.HLH",
                "Sphericity.original")

## plot selected features between alive and dead patients using normal data
plot_features(rdr = rdr, 
              feature_names = aerts_sign, 
              which_data = "normal", 
              conditions = colData(rdr)$deadstatus.event)

To represent features using comparable scales, plots can be generated using data scaled by min-max strategy.

## scale features by min-max method
rdr <- scale_feature_values(rdr = rdr, method = "minmax")

## plot selected features between alive and dead patients using scaled data
plot_features(rdr = rdr, 
              feature_names = aerts_sign,
              which_data = "scaled", 
              conditions = colData(rdr)$deadstatus.event)

Using quantile normalization, features are transformed to show identical distribution in terms of their statistical properties.

## normalize features by quantile normalization method
rdr <- normalize_feature_values(rdr = rdr, which_data = "scaled")
## plot selected features between alive and dead patients using normalized data
plot_features(rdr = rdr, 
              feature_names = aerts_sign,
              which_data = "normalized", 
              conditions = colData(rdr)$deadstatus.event)

As shown in the previous plot, quantile normalization can be used to improve the visualization of radiomic features.

Differential radiomic analysis by RadAR

RadAR allows the user to carry out differential radiomics analysis with aim of identifying radiomic features better differentiating between different conditions (e.g., prognosis, molecular alterations, pathology, HPV status). RadAR implements three not-parametric tests: the Wilcoxon-Mann-Whitney (WMW or wilcox) and the Area Under the Curve (AUC) for differential analysis between two groups, while the Kruskal-Willis test can be used to perform differential analysis across multiple (n>2) groups.

In the following code, we show the utility of the calc_differential_radiomics to identify the radiomic features associated with patients' outcome (dead vs. alive).

## perform differential radiomic analysis between good vs. poor prognosis groups of patients
rdr <- calc_differential_radiomics(rdr = rdr, 
                                   conditions = rdr$deadstatus.event, 
                                   method = "wilcox",
                                   adjust_pvalues_by = "BH")

The option conditions requires a vector of labels indicating the status of each patient. Top ranked features can be retrieved by the select_top_features() function.

## select top ranked features
top_features <- select_top_features(rdr = rdr, 
                                    which_statistics = "wilcox")
## print (first) top ranked features
rownames(top_features) <- NULL
## print the first top ranked features
top_features
## plot 4 top-ranked features
plot_features(rdr = rdr, 
              feature_names = c("Median.original", 
                                "JointAverage.wavelet.HHH", 
                                "LowGrayLevelEmphasis.wavelet.HHH", 
                                "SmallDependenceLowGrayLevelEmphasis.wavelet.LHH"),
              which_data = "normalized", 
              conditions = colData(rdr)$deadstatus.event)

In the following code, we show the utility of the calc_differential_radiomics to identify the radiomic features associated with cancer stage. Since this involves the comparison across n>2 classes, we will use the Kruskal-Wallis test.

## Print numbers of available stage 
table(rdr$Overall.Stage)
## perform differential radiomic analysis among patients with  different stage of cancer
rdr <- calc_differential_radiomics(rdr = rdr, 
                                   conditions = rdr$Overall.Stage, 
                                   method = "kruskal-wallis",
                                   adjust_pvalues_by = "BH")

## select top ranked features
top_features <- select_top_features(rdr = rdr, 
                                    which_statistics = "kruskal-wallis")

## print the first top ranked features
top_features
## plot two top-ranked features
plot_features(rdr = rdr, 
              feature_names = c("Sphericity.original", 
                                "Kurtosis.wavelet.LLH"),
              which_data = "normalized", 
              conditions = colData(rdr)$Overall.Stage)

Unsupervised clustering of radiomic features

Unsupervised analysis is useful when a response variable is not available. A widely used unsupervised analysis is cluster analysis, which is used to identify group of data showing similar characterististics. RadAR implements unsupervised analysis based on hierarchical clustering (hcl). Different distance measures and hierarchical clustering methods are implemented.

## print available distance methods
print_distance_methods()
## print available hcl methods
print_hcl_methods()

As first steps, distance matrix computation and hierarchical clustering are computed. In the example below Spearman’s correlation as distance measure is used for columns (eg., patients) and Euclidean distance measure is used for rows (features). Ward’s (default) hcl is used for both columns and rows.

## perform Ward's hierarchical clustering on correlation matrix estimated by Spearman correlation  
rdr <- do_hierarchical_clustering(rdr = rdr, 
                                  which_data = "scaled", 
                                  method_dist_col = "correlation.spearman")

Now we can visualize the results of the clustering by the heat map of radiomic features and dendrogram representing the results of the hcl. This can be done through the plot_heatmap_hcl() function. As the dataset includes hundreds of samples, it is useful to add color-coded annotation tracks reporting additional molecular and/or clinical information above the heat map. The annotation_tracks option allows the user to obtain automatic color-coded annotation tracks (exploiting pheatmap R package)

## draw a clusterized heat map of the NSCLC-Radiomics dataset + annotation tracks
plot_heatmap_hcl(rdr = rdr, 
                 annotation_tracks = c("Overall.Stage", 
                                       "Histology", 
                                       "deadstatus.event", 
                                       "gender")
                 )
## draw a clusterized heat map of the NSCLC-Radiomics dataset + annotation tracks
plot_heatmap_hcl(rdr = rdr, 
                 annotation_tracks = c("Overall.Stage", 
                                       "Histology", 
                                       "deadstatus.event", 
                                       "gender"), 
                 annotation_colors = list(Overall.Stage=c(I="#fef0d9", II="#fdcc8a", IIIa="#fc8d59", IIIb="#d7301f"),
                                          Histology=c(adenocarcinoma="#7fc97f", `large cell`="#beaed4", nos="#fdc086", `squamous cell carcinoma`="#ffff99"),
                                          gender=c(male="#66c2a5", female="#fc8d62"),
                                          deadstatus.event=c(`0`="grey90", `1`="black"))
                 )

The result of the clustering seems to indicate an enrichment of dead events in the cluster on the right part. This also suggests that unsupervised analysis of radiomic features could inform on patients’ survival in this dataset. To test if this is true, we can label patients with their cluster membership (assuming N=2 clusters) by the find_clusters() function of RadAR, that attaches cluster labels to the hcl_cols element of the RadAR object. Then, we can estimate statistical significance using a Fisher’s exact test.

## assign clustering membership to each patients based on previous hcl
rdr <- find_clusters(rdr = rdr, n_clusters_cols = 2)
## print contigency table
table(rdr$hcl_cols, rdr$deadstatus.event)
## test if hcl predicts enrichment of dead events / depletion of alive events in the right cluster
fisher.test(table(rdr$hcl_cols, rdr$deadstatus.event))

This analysis shows that unsupervised hiearchical clustering analysis applied to unselected radiomic features is able to infer the survival status of NSCLC patients. Indeed, in the cluster labeled as 1 we observe that 62% of patients are annotated as dead, while in cluster 2 48% are annotated as dead (p-value=1.4e-2, Odds-Ratio(OR)=1.8).

We can now apply the same strategy using quantile normalized feature data.

## perform Ward's hierarchical clustering on correlation matrix estimated by Spearman correlation  
rdr <- do_hierarchical_clustering(rdr = rdr, 
                                  which_data = "normalized", 
                                  method_dist_col = "correlation.spearman")
## draw a clusterized heat map of the NSCLC-Radiomics dataset + annotation tracks
plot_heatmap_hcl(rdr = rdr, 
                 annotation_tracks = c("Overall.Stage", 
                                       "Histology", 
                                       "deadstatus.event", 
                                       "gender")
                 )
## draw a clusterized heat map of the NSCLC-Radiomics dataset + annotation tracks
plot_heatmap_hcl(rdr = rdr, 
                 annotation_tracks = c("Overall.Stage", 
                                       "Histology", 
                                       "deadstatus.event", 
                                       "gender"), 
                 annotation_colors = list(Overall.Stage=c(I="#fef0d9", II="#fdcc8a", IIIa="#fc8d59", IIIb="#d7301f"),
                                          Histology=c(adenocarcinoma="#7fc97f", `large cell`="#beaed4", nos="#fdc086", `squamous cell carcinoma`="#ffff99"),
                                          gender=c(male="#66c2a5", female="#fc8d62"),
                                          deadstatus.event=c(`0`="grey90", `1`="black"))
                 )
## assign clustering membership to each patients based on previous hcl
rdr <- find_clusters(rdr = rdr, n_clusters_cols = 2)
## print contigency table
table(rdr$hcl_cols, rdr$deadstatus.event)
## test if hcl predicts enrichment of dead events / depletion of alive events in the right cluster
fisher.test(table(rdr$hcl_cols, rdr$deadstatus.event))

Using normalized feature data, unsupervised hcl analysis improved patients' stratification. Indeed, in the cluster labeled as 1 we now observe that 63% of patients are annotated as dead, while in cluster 2 46% are annotated as dead (p-value=3.2e-3, Odds-Ratio(OR)=1.97).

We can visualize this result also by Kaplan-Meier plot (the package survminer is required for this step).

## rename clusters based on clinical data
rdr$pts_hcl <- ifelse(rdr$hcl_cols == 1, "poor", "good")
## Create a Surv object with survival R package
surv_obj <- Surv(time = colData(rdr)$Survival.time, 
                 event = as.numeric(colData(rdr)$deadstatus.event))
## calculate Kaplan-Meier curve
kaplan_meier <- survfit(surv_obj ~  pts_hcl, 
                        data = colData(rdr), 
                        conf.type = "log-log")
## test if survminer package is available
if (requireNamespace("survminer", quietly = TRUE)) {
## plot Kaplan-Meier curve 
  survminer::ggsurvplot(kaplan_meier, 
                        risk.table = T, 
                        pval = T, 
                        xlim = c(0, 1500))
}

Feature selection methods

RadAR implements five methods to perform feature selection:

1) Minimum-redundancy-maximum-relevance (mRMR) feature selection using the mRMR.classic function of the mRMRe R package.

2) A method based on the calculus of correlation matrix followed by hierarchical clustering (hcl).

3) A method based on on Principal Component Analysis (PCA) followed by k-means clustering.

4) Lasso and Elastic-Net Regularized Generalized Linear Model (GLMNET) using cox regression model.

5) GLMNET using binomial regression model.

All methods are available through the function do_feature_selection(). For methods 2) and 3), after that groups of similar features have been identified by hcl or k-means clustering, the option select_by allows the user to choose between different criteria to select the most promising features:

a) Using concordance, the feature showing best concordance index within each cluster is selected. As for mRMR, this is a supervised method as it requires a response variable (a Surv object).

b) Using variability, the feature showing maximum values of (Q3-Q1)/Q2 within each cluster is selected.

c) Using random, a feature per each group is randomly picked up.

Radiomic signature building and evaluation

A widely used strategy to develop and evaluate a radiomic signature is to exploit multiple datasets. The signatures is initially developed on a training dataset and then evaluated for clinical relevance in multiple independent datasets (test). Here we see how to exploit the feature selection methods implemented in RadAR to delevop a radiomic signature using the Lung1 dataset and how to test it on a test dataset (Head & Neck 1 (HN1) dataset, from Aerts et al., 2014). HN1 original images and annotations can be downloaded from TCIA here. To facilitate its usability, we distribute the rdr object of HN1 dataset (only GTV-1 ROIs) within the RadAR package (as R data package). Also, we perform a comparative study of the radiomic signatures generated using the methods implemented in RadAR and original signature from Aerts et al., 2014.

As the feature "compactness" included in the original radiomic signature in Aerts et al. (2014) has been deprecated in pyradiomics, we derived the mathematical equivalent of compactness by taking the cube of the shape feature "sphericity" (Shi et al., 2019). For a question of convenience, we generate a new rdr replacing "sphericity" with "compactness".

## generate new rdr objects where sphericity is replaced by compactness 
## in training and validation datasets

rdr_new <- rdr 
assays(rdr_new)$values["Sphericity.original", ] <- assays(rdr_new)$values["Sphericity.original", ]^3
rownames(rdr_new)[which(rownames(rdr) == "Sphericity.original")] <- "Compactness.original" 
rownames(assays(rdr_new)$values)[which(rownames(rdr) == "Sphericity.original")] <- "Compactness.original" 

hn1_rdr_new <- hn1_rdr
assays(hn1_rdr_new)$values["Sphericity.original", ] <- assays(hn1_rdr_new)$values["Sphericity.original", ]^3
rownames(hn1_rdr_new)[which(rownames(hn1_rdr) == "Sphericity.original")] <- "Compactness.original" 
rownames(assays(hn1_rdr_new)$values)[which(rownames(hn1_rdr) == "Sphericity.original")] <- "Compactness.original" 

## define the Aerts et al., 2014 signature

aerts_signature <- c("Energy.original",
                     "GrayLevelNonUniformity.original",
                     "GrayLevelNonUniformity.wavelet.HLH",
                     "Compactness.original")

We generate different radiomic signature using the methods implemented in RadAR.

## set number of features 
## (n=8, roughly estimated from the heat map of correlation matrix)
nfeatures <- 8 

## generate a radiomic signature using mRMR
rdr_mRMR <- do_feature_selection(rdr = rdr,
                                 n_features = nfeatures, 
                                 method = "mRMR",
                                 response = as.numeric(rdr$deadstatus.event),
                                 which_data = "normal")$signature

## generate a radiomic signature using HC+variability
rdr_HC_var <- do_feature_selection(rdr = rdr,
                                   n_features = nfeatures,
                                   method = "hcl",
                                   select_by = "variability",
                                   which_data = "normal")$signature

## generate a radiomic signature using HC+concordance index
rdr_HC_conc <- do_feature_selection(rdr = rdr,
                                    n_features = nfeatures,
                                    method = "hcl",
                                    select_by = "concordance",
                                    surv_obj = surv_obj,
                                    which_data = "normal")$signature

## generate a radiomic signature using PCA+variability
rdr_PCA_var <- do_feature_selection(rdr = rdr,
                                    n_features = nfeatures,
                                    method = "pca",
                                    select_by = "variability",
                                    which_data = "scaled")$signature

## generate a radiomic signature using PCA+concordance index
rdr_PCA_conc <- do_feature_selection(rdr = rdr,
                                     n_features = nfeatures,
                                     method = "pca",
                                     select_by = "concordance",
                                     surv_obj = surv_obj,
                                     which_data = "scaled")$signature

## generate a radiomic signature using glmnet-cox
rdr_glmnet_cox <- do_feature_selection(rdr = lung1_rdr,
                                       method = "glmnet-cox",
                                       surv_obj = surv_obj,
                                       which_data = "normal")$signature

## generate a radiomic signature using glmnet-binomial
rdr_glmnet_binomial <- do_feature_selection(rdr = lung1_rdr,
                                            method = "glmnet-binomial",
                                            response = rdr$deadstatus.event,
                                            which_data = "normal")$signature

To evaluate the performance of a radiomic signature in independent radiomic datasets, RadAR implements the function test_radiomic_signature().

## generate the Surv object for the validation (HN1) dataset
surv_obj_hn1 <- Surv(time = hn1_clinical$overall_survival_in_days, 
                     event = hn1_clinical$event_overall_survival)

## evaluate the radiomic signature from Aerts et al., 2014 in training and validation datasets
training_Aerts <- test_radiomic_signature(rdr = rdr_new, 
                                         signature = aerts_signature, 
                                         surv_obj = surv_obj, 
                                         which_data = "normal")$concordance

validation_Aerts <- test_radiomic_signature(rdr = hn1_rdr_new, 
                                           signature = aerts_signature, 
                                           surv_obj = surv_obj_hn1, 
                                           which_data = "normal")$concordance


## evaluate the radiomic signature generated using mRMR in training and validation datasets
training_mRMR <- test_radiomic_signature(rdr = rdr, 
                                         signature = rdr_mRMR, 
                                         surv_obj = surv_obj, 
                                         which_data = "normal")$concordance

validation_mRMR <- test_radiomic_signature(rdr = hn1_rdr, 
                                           signature = rdr_mRMR, 
                                           surv_obj = surv_obj_hn1, 
                                           which_data = "normal")$concordance

## evaluate the radiomic signature generated using HC+variability in training and validation datasets
training_HC_var <- test_radiomic_signature(rdr = rdr, 
                                           signature = rdr_HC_var,
                                           surv_obj = surv_obj, 
                                           which_data = "normal")$concordance

validation_HC_var <- test_radiomic_signature(rdr = hn1_rdr,
                                             signature = rdr_HC_var,
                                             surv_obj = surv_obj_hn1, 
                                             which_data = "normal")$concordance

## evaluate the radiomic signature generated using HC+concordance in training and validation datasets
training_HC_conc <- test_radiomic_signature(rdr = rdr, 
                                            signature = rdr_HC_conc, 
                                            surv_obj = surv_obj, 
                                            which_data = "normal")$concordance

validation_HC_conc <- test_radiomic_signature(rdr = hn1_rdr, 
                                              signature = rdr_HC_conc, 
                                              surv_obj = surv_obj_hn1, 
                                              which_data = "normal")$concordance

## evaluate the radiomic signature generated using PCA+variability in training and validation datasets
training_PCA_var <- test_radiomic_signature(rdr = rdr, 
                                            signature = rdr_PCA_var, 
                                            surv_obj = surv_obj, 
                                            which_data = "scaled")$concordance

validation_PCA_var <- test_radiomic_signature(rdr = hn1_rdr, 
                                              signature = rdr_PCA_var, 
                                              surv_obj = surv_obj_hn1, 
                                              which_data = "scaled")$concordance

## evaluate the radiomic signature generated using PCA+concordance in training and validation datasets
training_PCA_conc <- test_radiomic_signature(rdr = rdr, 
                                             signature = rdr_PCA_conc,
                                             surv_obj = surv_obj, 
                                             which_data = "scaled")$concordance

validation_PCA_conc <- test_radiomic_signature(rdr = hn1_rdr,
                                               signature = rdr_PCA_conc,
                                               surv_obj = surv_obj_hn1,
                                               which_data = "scaled")$concordance

## evaluate the radiomic signature generated using glmnet-cox in training and validation datasets
training_glmnet_cox <- test_radiomic_signature(rdr = rdr,
                                               signature = rdr_glmnet_cox,
                                               surv_obj = surv_obj,
                                               which_data = "normal")$concordance

validation_glmnet_cox <- test_radiomic_signature(rdr = hn1_rdr,
                                                 signature = rdr_glmnet_cox,
                                                 surv_obj = surv_obj_hn1,
                                                 which_data = "normal")$concordance

## evaluate the radiomic signature generated using glmnet-cox in training and validation datasets
training_glmnet_binomial <- test_radiomic_signature(rdr = rdr,
                                                    signature = rdr_glmnet_binomial,
                                                    surv_obj = surv_obj,
                                                    which_data = "normal")$concordance

validation_glmnet_binomial <- test_radiomic_signature(rdr = hn1_rdr,
                                                      signature = rdr_glmnet_binomial,
                                                      surv_obj = surv_obj_hn1,
                                                      which_data = "normal")$concordance

To compare the performace of the radiomic signature reported above, we report them in a dot chart, ordering them from the best (top) to the worst (bottom) performing signatures in the in HN1 (validation dataset)

data_signatures <- rbind(c(training_Aerts[1], validation_Aerts[1]),
                         c(training_mRMR[1], validation_mRMR[1]),
                         c(training_HC_conc[1], validation_HC_conc[1]),
                         c(training_HC_var[1], validation_HC_var[1]),
                         c(training_PCA_conc[1], validation_PCA_conc[1]),
                         c(training_PCA_var[1], validation_PCA_var[1]),
                         c(training_glmnet_cox[1], validation_glmnet_cox[1]),
                         c(training_glmnet_binomial[1], validation_glmnet_binomial[1])
                         )

rownames(data_signatures) <- c("Aerts et al. (2014)", 
                               "mRMR", 
                               "HC+concordance", 
                               "HC+variability", 
                               "PCA+concordance", 
                               "PCA+variability",
                               "glmnet-cox",
                               "glmnet-binomial")

colnames(data_signatures) <- c("Training (Lung1)", 
                               "Validation (HN1)")

data_signatures <- data_signatures[order(data_signatures[, 2], decreasing = T), ]

dotchart(t(data_signatures),  
         col =  c("#E41A1C", "#377EB8"), 
         pch = 19, xlab = "Concordance Index", 
         xlim = c(.53, .67))

data_signatures

Save the rdr object

To save the rdr, use the following code.

## save a rdr object
save (rdr, file = "radar_object.rda")

References



cgplab/RadAR documentation built on Nov. 10, 2021, 1:32 a.m.