We studied 407 patients with GC whose data (including RNA sequencing data and clinical information) was obtained from The Cancer Genome Atlas (TCGA) (15)—a landmark cancer genomics program that contains a large amount of genomic and clinical datasets. RNA sequencing data were normalized using the fragments per kilobase transcriptome per million reads method. Among the 407 patients with GC, 46 patients with contrast-enhanced computed tomography (CECT) imaging data were obtained from The Cancer Imaging Archive (TCIA) (16)—a service that hosts a large archive of medical images of cancer accessible for public download. Most patients underwent imaging on a multi-detector row CT scanners: 4-, 16-, or 64-slice CT scanners with slice thickness of 2.5~5 mm at 120 kVp and 200~500 mA.

In addition, because the number of the patients in TCGA database was limited (there were only 46 patients who owned imaging data), we obtained external 424 patients with GC with imaging data in Nanfang Hospital to develop a prognostic nomogram model. Thirty-two patients were excluded for the invisible lesion or insufficient image quality, and then 392 patients with GC were included in further research. All of these patients underwent partial or total radical gastrectomy and did not receive other therapy preoperatively. All of the contrast-enhanced abdominal CT images were acquired within 30 days before surgical resection; the detailed inclusion and exclusion criteria were described in our previous study (17). Patient characteristics are detailed in Table 1 .

To avoid artificial and subjective deviation, two professional abdominal radiologists cooperated to delineate each tumor using the Medical Imaging Interaction Toolkit (MITK) (version 3.6.0 Apr. 1, 2017) without referring to any clinical information. Four samples are shown in Supplementary Figure S1 . Forty-four imaging features were extracted, including four first-order statistics and 40 texture features, from the delineated tumor outlines for each person using MATLAB automatically. Four different matrices, including the gray-level cooccurrence matrix (GLCM) (18), the neighborhood gray-tone difference matrix (NGTDM) (19), the gray-level size zone matrix (GLSZM) (20), and the gray-level run-length matrix (GLRLM) (21), were utilized to evaluate the 40 texture features of the regions of interest (ROIs). A detailed description of the four matrices and the features is shown in Supplementary Table S1 .

Intra- and interclass correlation coefficients (ICCs) were used to assess the intra- and interobserver reproducibility of radiomics feature extraction. Thirty images were chosen randomly for ROI segmentation and feature extraction. The ROI segmentation was performed by two radiologists with, respectively, 12 (reader 1) and 8 years (reader 2) of experience in gastric CT interpretation. Reader 1 then repeated the same procedure 1 week later. An ICC greater than 0.75 indicates good agreement of the feature extraction. To ensure the independence of imaging features, we removed the features with high correlation. We calculated the mean absolute correlation using the “findCorrelation” function in the caret package (version 6.0-84) in R software. When the pairwise correlation was greater than 0.8, the feature was identified as the one with redundancy and will be eliminated.

Compared with traditional analysis methods, which focus on a single gene or a few genes, we used weighted gene coexpression network analysis (WGCNA) (22), a systematic biological method to describe the pattern of gene association between different samples and obtain highly synergistic gene sets. The distinct advantage is that WGCNA converts an enormous amount of highly correlated genes into modules defined as clusters of densely interconnected genes to eliminate the problem of multiple hypothesis testing corrections.

WGCNA methodology can be used to construct a gene coexpression network, identify modules using hierarchical clustering and the dynamic tree cut method, study module relationships using eigengene networks, and obtain hub genes—the key regulators of each module. Because of the early advantages, the application of the WGCNA method has therefore been extended to variable types of high-throughput datasets, such as proteomic and metabolomic datasets, in recent years (23–25). To ensure the reliability of the prognostic gene modules, all the 407 patients with gene expression data were included.

The survival analysis of a time-to-event outcome is a widely applied method to predict survival from a set of patient predictors or covariates. Univariate Cox proportional hazards regression (26) was used to identify the gene modules with significant prognostic ability in OS, defined as the time from pathology diagnosis to the date of death or the last follow-up. For the module eigengene of gene modules that were continuous, the contributions to OS were investigated by Kaplan–Meier survival analysis (log-rank test).

We annotated the biological functions of the survival-related modules of GC using Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, allowing the identification of common biologic pathways from a public database. GO and KEGG analyses were performed using the clusterProfiler package (27, 28) in R. In addition, we calculated the false discovery rate (FDR) using the Benjamini–Hochberg method (29) to correct the multiple hypothesis testing of the survival analysis.

We utilized Spearman rank correlations to build a radiogenomic map combining selected imaging features and prognostic gene modules. The FDR was used again to correct for multiple testing. The imaging features that showed high correlation with the prognostic gene modules (P < 0.05, FDR < 0.05) were further analyzed by the univariate Cox proportional hazards model to assess the prognostic significance. Afterward, a bootstrap approach (30) was used to assess whether the prognostic signature had significant power compared with random chance (C-index = 0.5). For 100 times, we calculated the C-index of the imaging signature based on 35 randomly selected samples with correct outcome data and on 35 randomly chosen samples with random outcome data. The Wilcoxon test was applied to assess the differences between the two distributions (imaging signature and random chances).The prognostic signature with significant power (P < 0.05) was identified as the imaging biomarker.

On the basis of the prognostic imaging biomarkers identified from the TCGA database, the radiomics score (RIS) was constructed using a multivariate Cox proportional hazards regression model. The equation of RIS is detailed in the Supplementary Material . Because the median RIS is convenient for clinical application, it was used as a cutoff value to divide the patients into low‐risk and high‐risk groups. The log-rank test (Mantel–Cox test) was used to compare the difference in the survival status between two groups. We calculated RIS both in the TCGA database and Nanfang Hospital database, and log-rank test was used to identify their prognostic value. To discuss the prognostic significance of the imaging biomarkers and RIS, we tried to explain it from the perspective of immunity. We utilized the CIBERSORT algorithm (an algorithm that calculates the ratios of infiltrated immune cells) to investigate the immune infiltration level in patients with GC.

To compensate for the limitation of the number of patients in TCGA database, we utilized the database from Nanfang Hospital to construct a prognostic nomogram model for its large samples. We divided patients from Nanfang Hospital into training and validation cohorts in a ratio of 1:1 using the stratified randomization method. The training cohort was used to construct the nomogram that integrated both the RIS and significant clinicopathological characteristics, and the validation cohort was used to identify the predictive accuracy of the nomogram. The nomogram was constructed using “rms” package (31) of R software. With the application of the bootstrap method (1,000 replicates), a calibration curve was used to visualize the deviation of predicted probabilities from what actually happened. Moreover, decision curve analysis (32) was employed to quantify the clinical utility and compare the nomogram that we created with the TNM staging system. Finally, the receiver operating characteristic curve (ROC) analysis (33) was performed to measure the predictive performance of the nomogram score and the TNM stage. The Delong test was used to assess the differences between the two models.

Only two-sided P-value < 0.05 was considered to have necessary to statistic. The statistical analysis was implemented on R software (version 4.0.4). Several packages including ggplot2, rms, survival, survivalROC, survminer, nomogramFormula, and rmda were used to perform other statistical calculations and graphical work.

We extracted 44 imaging features from each person, and 16 imaging features were removed because their pairwise correlation was greater than 0.8. The correlation matrix is shown in Supplementary Figure S2 . The intraobserver ICCs ranged from 0.812 to 0.978, and the interobserver ICCs ranged from 0.786 to 0.912, indicating a favorable intra- and interobserver feature extraction reproducibility.

We grouped transcripts with correlated expression levels into gene coexpression modules using the WGCNA approach. The cluster dendrogram and gene counts for each module are shown in Supplementary Figure S3A . Analysis of network topology for various soft-thresholding powers is shown in Figures S3B, C . The soft thresholding power β was set at 9 because the scale independence reached 0.8 ( Figure S3B ) and had a relatively high-average connectivity ( Figure S3C ). At last, we identified 33 gene modules among 58,428 genes. Across the 33 modules, four prognostic modules (MEgreen, MEwhite, MEdarkturquoise, and MEroyalblue) were identified using univariate Cox proportional hazards regression ( Table 2 ), and the Kaplan–Meier plots of the four prognostic modules are shown in Figure 1 . Afterward, four modules were tested for the enrichment of specific biological functions in GO terms and KEGG pathways. The results are shown in Supplementary Table S2 . For example, the gene module MEdarkturquoise, enriched in protein kinase activator activity, was strongly correlated with a poor prognosis (P < 0.05). The gene module MEwhite, enriched in NADH dehydrogenase activity, was strongly correlated with positive survival.

We created the radiogenomic map that combined the remaining imaging features and four prognostic gene modules using Spearman rank correlations ( Figure 2 ).

One first-order statistic (volume.eccentricity) and 14 texture features (GLCM.Energy, GLCM.Contrast, GLCM.Entropy, GLCM.Variance, GLRLM.GLN, GLRLM.RP, GLRLM.SRLGE, GLRLM.LRLGE, GLSZM.GLN, GLSZM.ZSN, GLSZM.ZP, GLSZM.HGZE, GLSZM.LZLGE, and NGTDM.Complexity) were identified to have a strong correlation with the gene modules (P < 0.05, FDR < 0.05). There were 19 pairs of significant correlations among the features and gene modules. For example, the gene module MEdarkturquoise was negatively correlated with GLCM.Energy (P = 0.0379), GLRLM.GLN (P = 0.0161), GLSZM.ZSN (P = 0.0181), and GLSZM.ZP (P=0.0224) and positively correlated with GLCM.Contrast (P = 0.0079) and GLRLM.SRLGE (P = 0.0115).

Univariate Cox proportional hazards regression was employed to assess the survival relationship of the 15 imaging features with survival outcomes. Three texture features, namely, NGTDM.Complexity, GLRLM.LRLGE, and GLSZM.ZP (shown in Table 3 ), were significantly correlated with OS. By using the bootstrap method, the Wilcoxon test showed significant differences between the imaging features and random chances. Thus, the three imaging features were identified as the imaging biomarkers of GC. In the TCGA database, the RIS was constructed using a multivariate Cox proportional hazards regression model based on three imaging features and their corresponding coefficients ( Table S3 ). Then, patients were divided into a low-RIS group (n = 22) and a high-RIS group (n = 23) according to the median RIS. Compared with the low-RIS group, high-RIS group had a significantly lower OS rate [Hazard Ratio (HR) = 3.05, 95% confidence interval (CI): 1.35–6.92] ( Figure 3A ). To verify the ability of the RIS in predicting the OS of patients with GC, we further validated our findings in the database from Nanfang hospital, which yielded the similar results as above (HR = 3.41,95% CI: 2.41–4.82) ( Figure 3B ). Furthermore, multivariate Cox regression analysis was performed to determine whether the RIS was an independent prognostic factor for patients’ OS in both databases. In TCGA database, the P-value of RIS was 0.003 (HR = 1.99, 95% CI: 1.27–3.12), and in the database from Nanfang hospital, the P-value of RIS was less than 0.001(HR = 2.04, 95% CI: 1.52–2.75). Thus, RIS could be identified as an independent prognostic factor. Correlation between the imaging biomarkers/RIS and immune infiltration level is shown in the Supplementary Material .

In the training cohort of Nanfang Hospital database, the multivariate Cox analysis was employed to identify the independent prognostic factor in the RIS and clinicopathological characteristics. RIS, tumor metastasis, pathological lymph node positive detection rate, primary tumor, and age were significantly associated with OS ( Table 4 ), and they were integrated to construct a nomogram ( Figure 4A ). After adding the sum of the points of the five variables, we can obtain the patients’ OS probabilities. Patients were divided into high– and low–nomogram point groups based on the median point. The Kaplan–Meier curve showed that the OS of the high–nomogram point group was significantly poorer [ hazarad ratio (HR) = 8.58, 95% CI: 5.17–14.23, P < 0.001] ( Figure 4B ). The calibration curve of our nomogram is shown in Figures 4C, D , which demonstrated the accurate predictive ability of the nomogram. Moreover, the nomogram model showed a better net benefit and broader threshold probability than the AJCC TNM staging system in the decision curves ( Figure 4E ), and if the threshold probability is between 0 and 0.53, then using the radiogenomics nomogram to predict OS adds more benefit than AJCC stage. ROC analysis revealed that the nomogram model (AUC = 0.838) had a better prognostic value compared to the TNM stage (AUC = 0.765) ( Figure 4F ). Moreover, the AUC value of 1-, 3-, and 5-year survival were 0.801 (95% CI: 0.707–0.895), 0.829 (95% CI: 0.755–0.904), and 0.809 (95% CI: 0.688–0.931), respectively ( Figure 4G ). It indicated that the nomogram model was a powerful tool for predicting patients’ survival.

To further identify the predictive ability of the nomogram, we performed similar analysis in the validation cohort. The nomogram score for each patient was calculated in the same way (detailed in Table S4 ). Consisted with the findings mentioned above, Kaplan–Meier analysis showed that the patients in high–nomogram point group had a poor OS (HR = 4.93, 95% CI: 3.05–7.95, P < 0.001) ( Figure S4A ). Moreover, ROC analysis demonstrated that the nomogram (AUC = 0.791) model improved prognostic value in the validation compared with the TNM stage (AUC = 0.686) ( Figure S4B ). In addition, the time-dependent ROC showed that the nomgram model had a good accuracy with 0.843 (95% CI: 0.772–0.913) in 1 year, 0.772 (95% CI: 0.686–0.857) in 3 years and 0.767 (95% CI: 0.664–0.869) in 5 years ( Figure S4C ).