Establishing a cancer driver gene signature-based risk model for predicting the prognoses of gastric cancer patients

Construction of the CDG signature for GC

Overlapping prognostic CDGs from The Cancer Genome Atlas Stomach Adenocarcinoma (TCGA-STAD) and GSE62254 databases were selected to determine the candidate CDGs. Twenty-four CDGs were identified for final analysis (Supplementary Figure 1A). Next, least absolute shrinkage and selection operator (LASSO)-penalized Cox analysis was performed to narrow down the list of CDGs, and 12 genes were identified for downstream analyses (Supplementary Figure 1B). Multivariate Cox analysis was performed based on the CDGs selected by LASSO analysis (Supplementary Figure 1C). The prognostic risk scores of the CDG signature were determined as follows:

Risk score = (-0.06570) * (expression level of Damage Specific Deoxyribose Nucleic Acid [DNA] Binding Protein 2 [DDB2]) + 0.04589 * (expression level of Aminopeptidase [ENPEP]) + 0.00243* (expression level of Guanine Nucleotide binding protein, Alpha Stimulating activity polypeptide [GNAS]) - 0.10790 * (expression level of Musashi Ribose Nucleic Acid [RNA] Binding Protein 2 [MSI2]) + 0.14947 * (expression level of myosin Va [MYO5A]) + 0.22932 * (expression level of Pleomorphic Adenoma Gene 1 [PLAG1]) - 0.18526 * (expression level of RNA Binding Motif 15 [RBM15]; Supplementary Figure 1D).

CDG expression and its mutations in GC

To study the differences in CDG expression between tumor and normal tissues, we examined CDG messenger RNA (mRNA) levels in samples from TCGA-STAD. The results revealed that tumor tissues had significantly higher expression of DDB2, MSI2, and RBM15 than the normal tissues (Figure 1A). However, no differences were observed in the expression of ENPEP, GNAS, MYO5A, or PLAG1. We also examined the levels of proteins encoded by these CDGs using clinical samples from the Human Protein Atlas (HPA) database and found that the differences in protein levels were consistent with the observed differences in mRNA levels (Figure 1B).

Figure 1. Expression levels of cancer driver genes (CDGs) and their alterations in gastric cancer (GC). (A) CDG mRNA expression levels in GC obtained from The Cancer Genome Atlas Stomach Adenocarcinoma (TCGA-STAD). (B) Expression levels of proteins encoded by CDGs in normal tissues as obtained from the Human Protein Atlas (HPA) database (Data for GLAP1 was not available at HPA database).

The prognostic value of the CDG signature

To evaluate the prognostic value of the CDG signature in the training set, the GC patients were divided into high-risk (n = 167) and low-risk (n= 167) groups according to the median risk score. The corresponding signature risk score survival statuses were ranked and displayed on a dot-plot (Figure 2A, 2B). Individuals exhibited a greater risk of mortality with increasing risk score (Figure 2C). Heatmaps of the seven prognostic CDGs are displayed in Figure 2D.

Figure 2. Prognostic value of the cancer driver gene (CDG) signature in The Cancer Genome Atlas (TCGA) training set. (A) Distribution of risk scores per patient. (B) Relationships between overall survival (OS) status and survival time of gastric cancer (GC) patients ranked on the basis of risk score. (C) Comparison of mortality risk between the two groups in TCGA cohort. (D) Heatmap representing the expression profiles of the seven CDGs. (E) Kaplan-Meier analysis of OS between high- and low-risk groups in TCGA set. (F) Time-dependent receiver operating characteristic (ROC) analysis for OS prediction in TCGA set. (G) Time-dependent ROC analysis for grade prediction in TCGA set of OS. (H) Time-dependent ROC analysis for stage prediction in TCGA set of OS.

Kaplan-Meier analysis revealed that in the training set, patients in the high-risk group had shorter overall survival (OS) than those in the low-risk group (P < 0.001) (Figure 2E). Time-dependent receiver operating characteristic (ROC) analysis demonstrated that the area under the curve (AUC) values for one year, three years, and five years were 0.712, 0.613, and 0.611, respectively (Figure 2F). We also compared the prognostic value to traditional clinicopathological predictors. Time-dependent ROC analyses for grade-based and stage-based prediction of OS (TCGA) are shown in Figure 2G, 2H. Univariate and multivariate Cox regression analyses showed that age and risk score were both independent prognostic factors for GC (Figure 3A, 3B).

Figure 3. Forest plot depicting associations between risk factors and other clinical features and the prognosis of gastric cancer. (A) Univariate Cox regression analysis. (B) Multiple Cox regression analysis.

Validation of the CDG signature

To validate the CDG signature, patients in the validation set were divided into high- and low-risk groups, according to the median risk score (calculated using the CDG signature). The results were compatible with those obtained in the training set derived from TCGA. Figure 4A shows the heatmap of the seven prognostic CDGs. The corresponding signature risk score survival statuses were ranked and are displayed as a dot-plot (Figure 4A–4C). Individuals were at greater risk of mortality and recurrence with increasing risk scores (Figure 4E, 4F). Kaplan-Meier analysis also showed that the high-risk group had shorter OS and disease-free survival (DFS) than the low-risk group (P < 0.001; Figure 4G, 4H). Time-dependent ROC analysis of OS showed that the 1-, 3-, and 5-year AUC values were 0.662, 0.652, and 0.646 respectively (Figure 4I). The AUC values of one-, three-, and five-year DFS were 0.628, 0.637, and 0.619, respectively (Figure 4J).

Figure 4. External validation of the cancer driver gene (CDG) signature using the gastric cancer (GC) data from the Gene Expression Omnibus (GEO) validation set. (A) Heatmap representing expression profiles of the seven CDGs. (B) Distribution of risk scores per patient. (C) Relationships between overall survival (OS) status and survival time in GC patients ranked by risk score. (D) Relationships between disease-free survival (DFS) status and survival time in GC patients ranked by risk score. (E) Comparison of OS risk between the two groups. (F) Comparison of DFS risk between the two groups. (G) Kaplan-Meier analysis of OS between high- and low-risk groups in GSE62254. (H) Kaplan-Meier analysis of DFS between high- and low-risk groups in GSE62254. (I) Time-dependent receiver operating characteristic (ROC) analysis for OS prediction in the GSE62254 cohort. (J) Time-dependent ROC analysis for DFS prediction in the GSE62254 cohort.

Subgroup analysis of the CDG signature

To further estimate the utility of the CDG signature in predicting survival outcomes, stratification analysis was conducted based on specific clinicopathological characteristics. These subgroups included age (<65 or ≥65 years), gender, grade, stage, and T, N, and M stages (Figure 5A–5P). Stratification of the training and validation datasets revealed that the CDG signature could categorize patients into different survival groups and provide statistically significant prognostic values (Tables 1, 2).

Figure 5. Subgroup analysis of the cancer driver gene (CDG) signature. (A) Age < 65 years, (B) Age ≥ 65 years, (C) Male, (D) Female, (E) G2, (F) G3, (G) Stage I, (H) Stage III, (I) Stage IV, (J) T3, (K) T4, (L) M0, (M) M1, (N) N0, (O) N1, and (P) N3.

Gene set enrichment analysis (GSEA)

To further analyze the functions of the seven CDGs identified, GSEA was conducted for the high- and low-risk patients in the four datasets. GSEA results for the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways indicated that the “calcium signaling pathway,” “cell adhesion molecules,” (CAMs) “extracellular matrix receptor interaction,” “focal adhesion,” and “gap junction” categories were highly enriched in the high-risk group (Figure 6A). GSEA results for Gene Ontology (GO) terms indicated that the “collagen-containing extracellular matrix,” “contractile fiber,” “glycosaminoglycan binding,” “hormone binding,” and “muscle system processes” were highly enriched in the high-risk group (Figure 6B).

Figure 6. Gene set enrichment analysis (GSEA) of high- and low-risk groups. Top five representatives from (A) Gene Ontology (GO) term enrichment analysis and (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis.

Analysis of tumor immunity

TCGA-STAD gene expression matrix was uploaded to the Cell type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) platform to estimate the proportions of the 22 immune cell types. The high-risk group had high proportions of activated natural killer (NK) cells, monocytes, M2 macrophages, resting dendritic cells, and resting mast cells (Figure 7A). The low-risk group had high proportions of CD8 T-cells, CD4 memory-activated T-cells, follicular helper T-cells, resting NK cells, and M1 macrophages (Figure 7B).

Figure 7. Tumor immunity analysis based on the CDG signature. (A) Relative proportion of immune cells between high- and low-risk groups. (B) Violin plot depicting differences in the abundances of 22 types of immune cells between the high- and low-risk groups.

Establishment and evaluation of the nomogram

To predict the survival of GC patients, we constructed a nomogram based on the training set, which included the CDG signature risk score, age, gender, and pathological stage (Figure 8A). Time-dependent ROC analysis was performed to assess the predictive accuracy of the nomogram. Plotting the one-, three-, and five-year ROC values of OS in the training set of the nomogram revealed AUC values of 0.696, 0.639, and 0.632, respectively (Figure 8B). Time-dependent ROC analyses from the nomogram for one-, three- and five-year OS probabilities in the validation set returned AUC values of 0.825, 0.784, and 0.767, respectively (Figure 8C). In addition, ROC analysis of the DFS predictions in the validation set revealed that the nomogram was highly discriminatory, with AUC values of 0.825, 0.784, and 0.767 for one-, three-, and five-year DFS levels, respectively (Figure 8D).

Figure 8. Construction of the nomogram. (A) A nomogram for predicting one-, three-, and five-year overall survival (OS) generated by integrating the risk score, age, gender, and stage. (B) Time-dependent receiver operating characteristic ROC curves of the nomogram for OS prediction from the training set. (C) Time-dependent ROC curves of the nomogram for OS prediction from the validation set. (D) Time-dependent ROC curves of the nomogram for DFS predictions from the validation set.

Calibration and decision curve analysis (DCA) revealed the reliability of the nomogram for predicting prognoses in the training and validation sets. The calibration plot revealed that the predictions made using the nomogram were consistent with the true observations (Supplementary Figure 2A–2C). The DCA curves for the predictive nomogram revealed that it had high net benefit (Supplementary Figure 2D–2F). The web-based calculator (https://prognosis.shinyapps.io/STAD/) could predict the OS of GC patients based on the established nomogram (Supplementary Figure 3A, 3B) and is convenient in terms of its usage and visualization of the prognostic nomogram.

Data collection and processing

Publicly available transcriptomic and clinical data associated with GC samples were obtained from TCGA (https://tcga-data.nci.nih.gov/tcga/) and the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/gds/) and analyzed retrospectively. We used the RNA-Seq fragments per kilobase of transcript per million mapped reads (FPKM) data from TCGA. After excluding cases with follow-up times of <30 days, 634 patients were enrolled in the study, including 334 patients from TCGA-STAD project and 300 from the GSE62254 cohort [35]. Data from the GSE62254 cohort were obtained from the Asian Cancer Research Group (ACRG) Gastric cohort, which included 199 male and 91 female GC patients. The median age was 64 years and the range was 24-86 years.

Data for 568 CDGs (Supplementary Table 1) were downloaded from the Integrative OncoGenomics (IntOGen) pipeline (https://www.intogen.org/search) [7]. The immunohistochemical data associated with proteins encoded by each CDG in GC and normal tissues were obtained from the HPA (https://www.proteinatlas.org/).

Construction of the CDG signature

To narrow down the screening range, overlapping prognostic CDGs were selected from TCGA-STAD and GSE62254 cohorts via Cox univariate analysis. TCGA-STAD and GSE62254 were then used for model training (n = 334) and validation (n = 300), respectively. Previously selected CDGs were further screened and confirmed by LASSO Cox regression analysis (with the penalty parameter estimated by 10-fold cross-validation) using the “glmnet” package. A formula was developed using the CDG signature constructed above, where β corresponds to the correlation coefficient:

Risk score = β1 × (expression of RNA1) + β2 × (expression of RNA2) + ··· + βn × (expression of RNAn)

The patients in each dataset were assigned to a high- or a low-risk group using the median risk score as a cutoff. The ROC curves were created using the “survivalROC” package, and the AUC values were calculated to evaluate the predictive potential of the CDG signature.

Validation of the CDG signature

To validate the CDG signature, the patients in the validation set were separated into high- or low-risk groups according to the median risk score, which was calculated according to the CDG signature. Kaplan-Meier curve and time-dependent ROC analyses were conducted to assess CDG signature categorization.

Subgroup analysis of the CDG signature

To validate the effectiveness of the prognostic CDG signature, stratification analysis was performed on the training and validation sets using different demographic and clinical characteristics. The GC cases were divided into two risk groups according to their characteristics and risk scores, and Cox regression analysis was performed to analyze differences between the subgroups.

Estimation of immune cell infiltration

To analyze the relationship between the CDG signature and immune cell characteristics, CIBERSORT was used to estimate the fractions of immune cell types between the high- and low-risk groups [36]. Statistical analysis of the proportions of 22 immune cell types in each of the 334 GC samples was performed using the Wilcoxon rank-sum test.

GSEA

GSEA was performed to explore the GO terms and KEGG pathways that were significantly enriched in high-risk GC samples (http://www.broadinstitute.org/gsea). Gene sets were considered significantly enriched when FDR < 0.05 and |NES| > 1.

Construction and evaluation of the nomogram

We designed a novel nomogram model containing the CDG signature and clinicopathological predictors to establish a quantitative clinical tool to monitor and predict outcomes of GC patients. Subsequently, we developed a web-based calculator based on this model for clinical applications. Time-ROC curves and calibration plots were generated, and DCA was performed to evaluate the clinical utility of the novel nomogram.

Statistical analysis

All statistical analyses were performed using R software (version 4.0.5, R Development Core Team, 2021) and GraphPad Prism (version 8.3.0, GraphPad software, Inc., San Diego, CA, USA). All statistical tests (two-tailed) with P < 0.05 were considered statistically significant.