The study was conducted in accordance with the Declaration of Helsinki, and the study protocol was approved by the Ethics Committee of the Istituto Nazionale Tumori “Fondazione G. Pascale (protocol number 1/20). Sixteen patients with diagnosed ESCA were extrapolated by combining the public databases TCIA-TCGA and included in this study[32,33]. All subjects performed CT investigation with iodized contrast medium injection. Of 16 patients, one was excluded due to the presence of artifacts on acquired images. Due to their availability for all 15 patients, the transcriptome profiling and the following clinical variables were extracted and considered as outcomes for this study: Cancer stage, history of significant alcohol consumption and BMI value. We did not perform analyses on smoking measurements (tobacco history, age at starting smoking, pack-year smoked) due to the incompleteness of these data. To perform radiogenomic analyses, we divided patients into two groups according to stages I−II or III−IV, making stage outcome binary. Refer to Table 1 for clinical characteristics and outcomes of included patients.

CT examinations were acquired on GE Medical Systems (9 patients) and Siemens (6 patients) CT scanner, with slice thickness ranging from 1.25 mm to 2.5 mm and pixel size varying from 0.67 to 0.9. three-dimensional (3D) regions of interest (ROIs) encompassing the tumor were manually delineated slice-by-slice by using ITK-SNAP (version 3.6.0, http://www.itksnap.org) on the post-administration of contrast agents CT images. The tumor localization was divided into three regions: 7 patients presented lesions at the middle level, 2 subjects in the middle/distal area and 6 patients in the distal esophagus.

Gray-levels normalization was not performed, since CT gray values reflect absolute world values (Hounsfield units) and should be comparable between scanners. To correct variability from parameters related to voxel size and so unify voxel size across the cohort, radiomics data were extracted from images resampled to isotropic voxels of 1 × 1 × 1 mm³ using B Spline interpolator.

Radiomic features extraction: A total of 1302 radiomics features were extracted from segmented ROIs by using the open source Python package PyRadiomics (https://pyradiomics.readthedocs.io/en/Latest/). The extracted radiomics features were categorized into five groups: (1) Shape features (n = 14); (2) First-order features including 18 intensity statistics; (3) 74 multi-dimensional texture features including 23 gray level co-occurrence matrix (GLCM), 16 gray level size zone matrix (GLSZM), 16 gray level run length matrix (GLRLM), 14 gray level dependence matrix (GLDM) and 5 neighboring gray tone difference matrix (NGTDM) features; 1196 transformed first-order and textural features including; (4) 736 wavelet features in frequency channels LHL, LLH, HHH, HLH, HLL,HHL, LHH and LLL, where L and H are low- and high-pass filters, respectively; and (5) 460 LoG filtered features with sigma ranging from 1.0 and 5.0, with step size = 1. Features of groups (2) and (3) were grouped together and, from now on, this group will be referred to as “original features”. The extracted radiomics features grouped by similarity in four categories are listed in the Supplementary Table 1. The computing algorithms can be found at www.radiomics.io and the image biomarker standardization initiative presented a document to standardize the nomenclature and definition of radiomic features[34].

Radiomic features selection: Features selection was performed separately for shape features, original features, wavelet features and LoG filtered features in two steps, followed by a third step involving the whole set of features passed through the step I and II. In the first step, a correlation filter based on the absolute values of pairwise Spearman’s correlation (ρ) coefficient was used to reduce feature redundancy. Threshold for ρ was set to 0.8. Briefly, if two features had ρ > 0.8, the function looks at the mean absolute correlation of each variable and the variable with the largest mean absolute correlation is removed. The second step varied according to the type of outcome variables. For binary outcomes, a further feature restriction through a univariate analysis was performed by using non-parametric Wilcoxon rank-sum test to investigate the statistical significance with respect to the outcome. Statistical significance was set to P < 0.05. The significantly different features were then selected and further reduced in the third step. For continuous outcomes, the second step consisted of computing ρ between each feature selected through the step I and the reference outcome. Then, all features with a significant ρ > 0.5 and P value < 0.05 were considered. In order to check for redundancies among features belonging to the different four groups, the third step consisted in applying the correlation filter described in step I to the whole feature set passed through step II. All steps were implemented using Matlab R2020a (The MathWorks Inc., Natick, MA, United States).

Predictive models building and analysis for stage assessment: In order to evaluate the predictive power of CT radiomic features taken by them for ESCA staging, a fourth step of feature selection was performed for features that were associated with stage. The latter step consisted in ranking the remaining features based on the mutual information (MI) between the distribution of the values of a certain feature and the membership to a particular class. Features were evaluated independently, and the final feature selection occurred by aggregating the five top ranked ones[35-37]. For the binary stage I-II/stage III-IV classification task, the reduced feature set was used to build logistic regression models of order from 1 to 5 that would best predict ESCA stage by using an imbalanced-adjusted bootstrap resampling (IABR) approach on 1000 bootstrap samples[38]. Specifically, the training set was made up 1000 bootstrap samples randomly drawn with replacement from the available dataset. The testing set consisted of the instances that did not belong to the bootstrap sample. Then, application of the imbalance-adjustment step made the probability of picking a positive and a negative instance in the bootstrap sample the same[39-41].

For each model order, the combination of features maximizing the 0.632+ area under the receiver operating characteristic curve (AUC) within 1000 bootstrap training and testing samples was identified. Finally, IABR on 1000 samples was performed again for all models to assess prediction performances[38,42].

Additional analyses were performed starting from the first two, three and four features surviving after the MI-based feature selection step (which was used to build logistic regression models of order from 1 to 2, 1 to 3 and 1 to 4, respectively) that would best predict ESCA stage. Moreover, given that the overall stage is determined after the cancer is assigned categories describing the tumor (T), node (N) and metastasis (M) categories, we tested the capability of these features for predicting T and N status. Analyses assuming M status as clinical outcome were not performed due to the extremely unbalanced sample. Patients were divided into two groups according to T1-T2 or T3−T4 tumor status, making T stage outcome binary. Similarly, we evaluated if CT radiomic features could assess N status by dividing patients into two groups according to the absence (N0) or presence (N1-N2-N3) lymph node status[43].

RNA-Seq and miRNA-Seq data of esophageal carcinoma of tumor tissues were downloaded from the GDC Data Portal (https://gdc-portal.nci.nih.gov/) considering for TCGA-ESCA project only patients with associated imaging data from TCIA database (see Supplementary Tables 2 and 3 for Clinical and Transcriptomic data, respectively). The mRNAs expression levels were considered as read count based on gene length and the total number of mapped reads (FPKM values). Moreover, miRNA expression quantification was downloaded and normalized counts in reads-per-million-miRNA-mapped were considered.

An integrative study design was defined (see Figure 1 for the flowchart reporting the organization of data and analyses in the study) and reported as radiogenomic workflow in Figure 2 in order to evaluate potential association between significant radiomic features according to selected clinical variables (stage, history of significant alcohol consumption and BMI) with biomarkers and RNA regulators characterizing esophageal cancer. For this purpose, a Spearman’s correlation analysis was investigated between radiomic features selected after the three features selection steps described above and transcriptomic signatures suggested by Zeng et al[22] and Xu et al[25]. Specifically, we calculated ρ between the whole selected feature set and the eight m6A RNA methylation regulators (KIAA1429, HNRNPC, RBM15, METTL3, WTAP, YTHDF1, YTHDC1, YTHDF2), as well as ρ between the whole selected feature set and the five up-regulated miRNAs (miRNA-93, miRNA-21, miRNA-4746, miRNA-196a-1, miRNA-196a-2). The significance level was set to 0.05. All analyses were performed using Matlab R2020a. The statistical methods of this study were reviewed by our bioinformatics and biostatistics group of our research support center.

Open in New Tab   Full Size Figure   Download Figure

Figure 1 Flowchart reporting the organization of data and analyses in the study. BMI: Body mass index; CT: Computed tomography; TCGA: The Cancer Genome Atlas; TCIA: The Cancer Imaging Archive.

Open in New Tab   Full Size Figure   Download Figure

Figure 2 Workflow of radiogenomic analysis implemented in the study. On the first row the radiomic analysis steps. On the second row the radiogenomic analysis for each clinical outcome. BMI: Body mass index; miRNA: MicroRNA; mRNA: Messenger RNA.

Radiomic feature selection: None of the 14 shape features passed the step II of features selection, both considering binary outcomes (namely stage and alcohol history) and BMI. So, radiogenomic analysis was not performed for this feature group. Concerning original, wavelet and Log sigma feature groups, the step I of feature selection reduced the feature sets from 92, 736 and 460 to 23, 127 and 65, respectively. Concerning stage analysis, Wilcoxon rank-sum test used in step II of feature selection revealed significant results for 26 radiomic features, of which 25 belonging to the wavelet feature set and the remaining one was the original first order maximum value. However, the latter feature did not pass the second correlation filter of step III (refer to Table 2). Considering the presence of alcoholic history as clinical outcome, Wilcoxon rank-sum test revealed significant results for five features, of which one belonging to original feature set (Kurtosis) and the remaining ones to wavelet feature set. These five features passed step III, since Kurtosis showed a correlation lower than 0.8 with all four wavelet features (refer to Table 2). Lastly, considering BMI as the clinical outcome, 26 wavelet and three LoG sigma features were selected due to a significant ρ > 0.5 (P < 0.05) with BMI. In total, 29 radiomic features (see Table 3) survived after the second correlation filter and were associated with the three examined clinical outcomes.

Predictive models building and analysis for stage assessment: The top five features selected after the MI-based feature selection step were wavelet LLH GLDM high gray level emphasis, LLH NGTDM complexity, HHH GLCM joint entropy, HLH entropy and HLL GLCM cluster prominence. Prediction performances of multivariable logistic regression models for the stage I-II/stage III-IV classification task were very high for both five model orders. However, by inspecting prediction performances values in Table 4, we determined that the simplest multivariable model with the best prediction performances were reached by the second order model (AUC = 87%, sensitivity = 64%, specificity = 83% and accuracy = 79%), which was based on wavelet LLH NGTDM complexity and HHH GLCM joint entropy. These results were also confirmed by additional analyses (Supplementary Tables 4-6). The top five features were also found to be able to predict T and N staging, with best AUCs (0.79 and 0.80, respectively) reached by second order models (see Supplementary Tables 7 and 8 for analyses involving five features). Results of additional analyses for T and N prediction by using two, three and four features are reported in Supplementary Tables 9-11 (T staging) and Supplementary Tables 12-14 (N staging).

Stage: Overall, radiogenomic analysis revealed that six of the eight mRNA regulators were significantly correlated with 10 of the 25 selected radiomic features, which were all belonging to the wavelet group. In particular, HRNPC and WTAP were positively correlated with wavelet HHL NGTDM strength (ρ = 0.61, ρ = 0.61, P < 0.05, respectively); METTL3 was positively correlated with wavelet LHL GLDM high gray level emphasis, HLL GLRLM gray level variance, HLL GLDM dependence entropy and LLL GLDM small dependence high gray level emphasis (ρ = 0.54, ρ = 0.53, ρ = 0.56, ρ = 0.6, P < 0.05, respectively) and negatively correlated with wavelet HLL GLCM inverse variance (ρ = -0.6, P < 0.05); YTHDF1 reported a positive and significant correlation with the wavelet feature HLL GLCM maximum probability (ρ = 0.6, P < 0.05) and a negative correlation with HLL 90^th percentile (ρ = -0.61, P < 0.05); YTHDF2 was positively correlated with the wavelet feature HHH GLCM contrast and HHH GLSZM zone percentage (ρ = 0.57, ρ = 0.56, P < 0.05, respectively); the latter feature was also positively correlated with YTHDC1 (ρ = 0.6, P < 0.05). Moreover, correlation analysis with the five up-regulated miRNA revealed a significant positive correlation between miRNA-93 and two radiomic features, namely wavelet LHL GLDM high gray level emphasis (ρ = 0.69, P < 0.05) and HHH GLCM joint entropy (ρ = 0.58, P < 0.05). Notably, HHH GLCM joint entropy contributed to building the best predictive models for stage assessment, as well as T and N assessment. Finally, a positive correlation between miRNA-4746 and HHL GLCM cluster shade was found (ρ = 0.53, P < 0.05). The radiogenomic results for stage are shown through a heatmap in the Figure 3.

Open in New Tab   Full Size Figure   Download Figure

Figure 3 Radiogenomic analysis using stage as clinical outcome. Heatmap depicting the correlation matrix between transcriptomic features and radiomic features significantly associated with stage. GLCM: Gray level co-occurrence matrix; GLDM: Gray level dependence matrix; GLRLM: Gray level run length matrix; GLSZM: Gray level size zone matrix; H: High-pass filter; L: Low-pass filter; NGTDM: Neighboring gray tone difference matrix.

Alcohol history: From the integrated analysis, none of the five selected radiomic features was significantly correlated with any of eight RNA regulators (data not shown). Conversely, correlation analysis with the five up-regulated miRNA revealed a significant correlation between the wavelet feature LHH first order Mean and three up-regulated miRNA, namely miRNA-21 (ρ = -0.56, P < 0.05), miRNA-4746 (ρ = 0.64, P < 0.05) and miRNA-93 (ρ = 0.61, P < 0.05), as reported in the heatmap in Figure 4.

Open in New Tab   Full Size Figure   Download Figure

Figure 4 Radiogenomic analysis using alcohol history as clinical outcome. Heatmap depicting the correlation matrix between transcriptomic features and radiomic features significantly associated with alcohol history. H: High-pass filter; L: Low-pass filter; miRNA: MicroRNA.

BMI: Overall, radiogenomic analysis revealed that four of the eight mRNA regulators were significantly correlated with 10 of the 29 selected radiomic features, of which five belonging to wavelet group and three belonging to LoG sigma group as reported in Figure 5. In particular, METTL3 was positively correlated with HLH median, LHL GLDM high gray level emphasis, HLL GLRLM gray level variance and LLL GLDM small dependence high gray level emphasis (ρ = 0.71, ρ = 0.54, ρ = 0.53, ρ = 0.6, P < 0.05, respectively) and negatively correlated with HLL GLCM inverse variance and LoG sigma 10^th percentile calculated with sigma = 5.0 (ρ = -0.6, ρ = -0.66, P < 0.05, respectively); YTHDF1 was positively correlated with wavelet feature HLL GLCM maximum probability (ρ = 0.6, P < 0.05), whereas YTHDC1 was positively correlated with the wavelet HHH GLCM difference average (ρ = 0.67, P < 0.05) and inversely correlated with and LoG sigma root mean squared calculated with sigma = 1.0 (ρ = -0.67, P < 0.05); similarly, YTHDF2 was positively correlated with the wavelet HHH GLCM contrast (ρ = 0.57, P < 0.05) and inversely with the LoG sigma 10^th percentile calculated with sigma = 5.0 (ρ = -0.56, P < 0.05). Moreover, correlation analysis with the five up-regulated miRNA revealed a significant correlation only between miRNA-93 and two radiomic features as reported in Figure 5 and in particular with the wavelet feature LHL GLRLM high gray level run emphasis (ρ = 0.69, P < 0.05) and the LoG sigma 10^th percentile calculated with sigma = 5.0 (ρ = -0.54, P < 0.05).

Open in New Tab   Full Size Figure   Download Figure

Figure 5 Radiogenomic analysis using body mass index as clinical outcome. Heatmap depicting the correlation matrix between transcriptomic features and radiomic features significantly associated with body mass index. GLCM: Gray level co-occurrence matrix; GLDM: Gray level dependence matrix; GLRLM: Gray level run length matrix; GLSZM: Gray level size zone matrix; H: High-pass filter; L: Low-pass filter; miRNA: MicroRNA; mRNA: Messenger RNA; NGTDM: Neighboring gray tone difference matrix.

Esophageal cancer (ESCA) is the sixth most common malignancy in the world, and its incidence is rapidly increasing. Radiogenomics provides clinically useful prognostic predictions by linking molecular characteristics such as gene mutations and gene expression patterns of malignant tumors with medical images and could provide more opportunities in the management of patients with ESCA.

Recently, several microRNAs (miRNAs) and messenger (RNA) targets were evaluated as potential biomarkers and regulators of epigenetic mechanisms involved in ESCA. In addition, the use of computed tomography (CT) radiomics is rapidly increasing in the field of ESCA and plays an important role in different ESCA management steps. Moreover, there are no previous radiogenomic studies on ESCA investigating the association with clinical staging and potential risk factors. This has motivated us to investigate on the relationship between the expression levels of eight N6-methyladenosine (m6A) RNA methylation regulators, as well as the five up-regulated miRNAs, and radiomic features extracted from CT images that were significantly associated to clinical outcomes.

To explore the combination of CT radiomic features and molecular targets associated with clinical outcomes for characterization of ESCA patients.

Fifteen patients with diagnosed ESCA were included in this study, and their CT imaging and transcriptomic data were extracted from The Cancer Imaging Archive and gene expression data from The Cancer Genome Atlas, respectively. Cancer stage, history of significant alcohol consumption and body mass index (BMI) were considered as clinical outcomes. Radiomic analysis was performed on CT images acquired after injection of contrast medium. In total, 1302 radiomics features were extracted from three-dimensional regions of interest by using PyRadiomics. Radiogenomic analysis involved Spearman’s correlation (ρ) analysis between radiomic features associated with clinical outcomes and transcriptomic signatures consisting in eight m6A RNA methylation regulators and five up-regulated miRNA.

Radiogenomic analysis with stage as clinical outcome revealed that six of the eight mRNA regulators and two of the five up-regulated miRNA were significantly correlated with 10 and three of the 25 selected radiomic features, respectively. Assuming alcohol history as clinical outcome, no correlation was found between the five selected radiomic features and mRNA regulators, while a significant correlation was found between one radiomic feature and three up-regulated miRNA. Radiogenomic analysis with BMI as clinical outcome revealed that four mRNA regulators and one up-regulated miRNA were significantly correlated with 10 and two radiomic features, respectively.

Our study revealed interesting relationships between the expression of eight m6a RNA regulators, as well as five up-regulated miRNAs, and CT radiomic features associated with clinical outcomes of ESCA patients.

This preliminary study revealed interesting associations between m6a RNA regulators, as well as miRNAs, and CT radiomic features associated with clinical outcomes of ESCA patients. Further investigations on different ESCA omics data are required. Moreover, multimodal data combined with artificial intelligence techniques characteristics are desirable.