In cgplab/RadAR: Radiomic Analysis with R
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
-
Aerts HJWL, Velazquez ER, Leijenaar RTH, Parmar C, Grossmann P, Carvalho S, et al. Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach. Nat. Commun. 2014;5:4006.
-
Fedorov A, Beichel R, Kalpathy-Cramer J, Finet J, Fillion-Robin J-C, Pujol S, et al. 3D Slicer as an Image Computing Platform for the Quantitative Imaging Network. Magnetic Resonance Imaging. 2012 Nov;30(9):1323-41.
-
Gillies RJ, Kinahan PE, Hricak H. Radiomics: Images Are More than Pictures, They Are Data. Radiology. 2016;278:563–77.
-
Nioche C, Orlhac F, Boughdad S, Reuzé S, Goya-Outi J, Robert C, et al. LIFEx: A Freeware for Radiomic Feature Calculation in Multimodality Imaging to Accelerate Advances in the Characterization of Tumor Heterogeneity. Cancer Res. 2018;78:4786–9.
-
Shi Z, Zhovannik I, Traverso A, DAnkers JWMF, Deist TM et al. Distributed radiomics as a signature validation study using the Personal Health Train infrastructure. Sci Data 6, 218 (2019).
-
van Griethuysen JJM, Fedorov A, Parmar C, Hosny A, Aucoin N, Narayan V, et al. Computational Radiomics System to Decode the Radiographic Phenotype. Cancer Res. 2017;77:e104–7.
cgplab/RadAR documentation built on Nov. 10, 2021, 1:32 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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
-
Aerts HJWL, Velazquez ER, Leijenaar RTH, Parmar C, Grossmann P, Carvalho S, et al. Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach. Nat. Commun. 2014;5:4006.
-
Fedorov A, Beichel R, Kalpathy-Cramer J, Finet J, Fillion-Robin J-C, Pujol S, et al. 3D Slicer as an Image Computing Platform for the Quantitative Imaging Network. Magnetic Resonance Imaging. 2012 Nov;30(9):1323-41.
-
Gillies RJ, Kinahan PE, Hricak H. Radiomics: Images Are More than Pictures, They Are Data. Radiology. 2016;278:563–77.
-
Nioche C, Orlhac F, Boughdad S, Reuzé S, Goya-Outi J, Robert C, et al. LIFEx: A Freeware for Radiomic Feature Calculation in Multimodality Imaging to Accelerate Advances in the Characterization of Tumor Heterogeneity. Cancer Res. 2018;78:4786–9.
-
Shi Z, Zhovannik I, Traverso A, DAnkers JWMF, Deist TM et al. Distributed radiomics as a signature validation study using the Personal Health Train infrastructure. Sci Data 6, 218 (2019).
-
van Griethuysen JJM, Fedorov A, Parmar C, Hosny A, Aucoin N, Narayan V, et al. Computational Radiomics System to Decode the Radiographic Phenotype. Cancer Res. 2017;77:e104–7.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.