Development and validation of an individual alternative splicing prognostic signature in gastric cancer

Overview of AS events in the TCGA GC cohort

We preliminarily detected 119,697 AS events from 14,433 genes, which accounted for approximately 70% of protein-coding genes [18]. A large proportion of AS events could only be detected in a few samples, and specific splicing isoforms were barely detected (PSI value < 0.05). Then, we implemented a series of filters (percentage of samples with PSI values ≥ 75 and average PSI value ≥ 0.05) to obtain a reliable set of AS events in GC. Consequently, a total of 37,139 AS events were obtained from 10,380 genes, which were used for further analysis. Among these AS events, they were comprised of “15,816 ES events from 6,462 genes”, “6,633 AP events from 4,016 genes”, “5,988 AT events from 3,665 genes”, “3,285 AA events from 2,404 genes”, “2,691 AD events from 1,986 genes”, “2,564 RI events from 1,780 genes” and “162 ME events in 156 genes” (Supplementary Figure 1A). Among these seven types of splicing patterns, ES was the most frequent splicing pattern, followed by AP and AT, while ME was the least frequent splicing pattern. Of note, the number of AS events far exceeded the total number of genes, which indicated that one gene might undergo multiple splicing patterns. As shown in Supplementary Figure 1B, most genes contained two or more AS events, and several genes could have up to five different splicing modes.

Identification of GC-associated AS events

To identify GC-associated AS events, we compared the PSI values between 33 paired tumor and normal tissues. According to the defined threshold, a total of 1,287 GCAS events were identified from 837 genes, which included 738 upregulated AS events from 625 genes and 549 downregulated AS events from 509 genes (Figure 1A and Supplementary Table 1). The proportion of AS splicing modes between the GCAS profiles and the entire AS profile was inconsistent (Figure 1B). A larger proportion of ES events (42.6%) were detected among all AS profiles than among GCAS events (32.6%). Moreover, AP events accounted for 17.9% of events among all AS profiles and accounted for 35.1% of GCAS events. Forty percent of the splicing-parent genes possessed more than one GCAS event. However, the proportion of splicing-parent genes derived from more than two GCAS events was very low (Supplementary Figure 1C).

Figure 1. Profiling of AS events identified in the TCGA GC cohort. (A) GCAS events between GC and adjacent normal tissues were visualized in a volcano plot. (B) The number of GCAS events and total AS events were depicted according to the seven AS types. (C) Heat map for the PSI values of 173 OS-GCAS events identified in 354 GC patients. (D) Kaplan-Meier curves for the paired survival-related GCAS events (IL11RA_86208_AP and IL11RA_86209_AP). (E) Violin plots for the PSI values of the paired survival-related GCAS events (IL11RA_86208_AP and IL11RA_86209_AP) between GC and adjacent normal tissues.

Identification of OS-GCAS events

To investigate the relationship between GCAS and OS in GC patients, we performed a univariate Cox regression analysis of 1,287 GCAS events in 354 patients. As shown in Figure 1C, 173 OS-GCAS events were identified from 130 genes (Supplementary Table 2). The proportion of AS splicing modes between the OS-GCAS and GCAS profiles was also inconsistent (Supplementary Figure 1D). The AP pattern (69 cases) contained the most OS-GCAS events, followed by the ES pattern (55 cases) and AT pattern (29 cases), while the AA pattern (2 cases) contained the fewest OS-GCAS events. Although 33.3% of splicing-parent genes had two or more OS-GCAS events, almost all the genes had only one splicing pattern (Supplementary Figure 1E). Of note, among the 130 splicing-parent genes, 38 genes contained paired survival-related GCAS events. For example, the AP of IL11RA in exon site 1 (IL11RA_86209_AP) was a favorable prognostic indicator, whereas the AP of IL11RA in exon site 2 (IL11RA_86208_AP) was a poor prognostic indicator (Figure 1D). Interestingly, the splicing pattern of all the paired survival-related GCAS events was either AP or AT, and the expression of the paired survival-related GCAS events was opposite between tumor and normal tissues (Figure 1E).

Enrichment and interaction analysis of OS-GCAS events

All 130 splicing-parent genes derived from 173 OS-GCAS events were assessed by functional and pathway enrichment analysis to explore the underlying mechanisms of the OS-GCAS events. The results revealed that genes were enriched in GO categories closely related to GC development, such as “microtubule polymerization”, “centrosome localization” and “cytokinesis” (Figure 2 and Supplementary Table 3). In addition, some canonical pathways involved in GC metastasis and recurrence were also enriched, such as the “Notch signaling pathway”, “adherens junction” and “bacterial invasion of epithelial cells” (Figure 2 and Supplementary Table 3).

Figure 2. Potential biological functions of 173 OS-GCAS events in GC.

Since AS events could inevitably affect protein characteristics, it is necessary to investigate these OS-GCAS events at the protein level [19]. A protein interaction network of OS-GCAS events was constructed not only to provide an overview of the interactions under normal conditions but also to reveal the potential influence of OS-GCAS events on the entire network. After removing the isolated nodes, 79 genes were mapped in the protein interaction network, and these parental genes were significantly correlated with each other (Supplementary Figure 2A). The top three hub genes were SORBS1, BPTF and SEPT2, which were determined by the maximal clique centrality method. The MCODE algorithm identified three individual modules from the whole protein interaction network (Supplementary Figure 2B–2D).

The network of OS-GCAS events and SFs

SFs are critical regulators of AS events that bind to pre-mRNAs and influence exon selection and splicing site choice [20]. It has been determined that dysregulated AS events within the tumor microenvironment may be regulated by a limited number of SFs [21]. Thus, correlation analysis was performed to explore the potential interaction networks of OS-GCAS events and SFs. As shown in Figure 3A, 103 of 173 OS-GCAS events were significantly associated with 41 experimentally validated SFs in the regulatory network (|R| ≥ 0.4 and P-value < 0.05). Among the networks, most SFs were significantly correlated with more than one OS-GCAS event, while a single OS-GCAS could be regulated by up to 21 different SFs (Figure 3A). Moreover, we found that some SFs played opposite roles in the regulation of different AS events. For example, the expression of SRSF9 was negatively correlated with ATP2B4_9450_ES but positively correlated with APOLD1_20517_AT (Figure 3B).

Figure 3. The regulatory splicing correlation network in GC. (A) The correlation of OS-GCAS events with SFs is shown in network plots. The circular node represents the OS-GCAS event. The diamond node represents the SF. (B) Representative positive correlations between OS-GCAS events and SFs are shown in scatter plots. (C) Representative negative correlations between OS-GCAS events and SFs are shown in scatter plots.

Construction and evaluation of the AS prognostic signature

To facilitate the application of AS events in the clinical monitoring of the prognosis of GC patients, we applied a LASSO penalty with a multivariate Cox regression analysis to the 173 OS-GCAS events in the training set of 261 patients. A total of 23 features were identified with nonzero coefficients (Figure 4A, 4B). These 23 LASSO-selected features were used to build the AS prognostic signature. The coefficients and PSI values of these 23 OS-GCAS events were used together to calculate the risk scores for both the training and validation datasets (Figure 4C).

Figure 4. Construction and evaluation of the 23-AS event prognostic signature for GC patients. (A, B) LASSO regression analysis of OS- GCAS events. (C) The 23 OS-GCAS events included in the signature. Corresponding coefficients from multivariate Cox regression using LASSO and log10(HRs) are depicted by horizontal bars and dots, respectively. (D) Time-dependent ROC curves for the 23-AS event signature in the training and validation datasets. (E, F) Kaplan-Meier curves with difference detection by log-rank test for patients from the training and validation datasets.

This 23-AS event prognostic signature showed excellent performance in both datasets (Figure 4D). The C-index values of the signature were 0.725 and 0.762 in the training and validation datasets, respectively. The signature significantly stratified patients into low- and high-risk groups in the training set and validation set (Figure 4E, 4F). Significant RMS time ratios were also observed in the two datasets (Table 1). Moreover, the 23-AS event prognostic signature remained an independent prognostic factor in multivariate analyses after adjusting for clinical and pathologic factors, such as age, TNM stage and microsatellite instability (MSI) status (Table 2).

Furthermore, we performed sensitivity analyses according to age, sex and TNM stage to evaluate the robustness of the 23-AS event prognostic signature in different clinical subgroups. The prognostic signature could stratify patients into significantly different prognostic groups in all the subgroups, which indicated that this prognostic signature might predict OS independently of clinical characteristics in GC (Supplementary Figure 3).

Clinical, molecular and immune features underlying the 23-AS event prognostic signature

To explore the relationship between patient characteristics and the 23-AS event prognostic signatures, patients in the entire TCGA dataset were stratified into high- and low-risk groups according to the optimal cut-off of 0.046 (Figure 5A). Among the 327 patients, 137 patients were assigned to the high-risk group, and 190 patients were assigned to the low-risk group. The low-risk group had a significantly more favorable prognosis of OS than the high-risk group (Figure 5B). We found that almost all 23 LASSO-selected features were significantly differentially expressed between the two risk groups (Figure 5C, 5D). Moreover, the distribution of MSI status, TNM stage and molecular subtype were significantly different between the high- and low-risk groups (Figure 5C and Table 3).

Figure 5. Clinical and molecular features underlying the 23-AS event signature. (A) The optimal cut-off of the 23-AS event signature in the entire TCGA dataset. (B) Kaplan-Meier curve with difference detection by log-rank test for patients in the entire TCGA dataset. (C) Heat map for the expression patterns of the 23 OS-GCAS events for the entire TCGA sample set sorted by the risk score in ascending order. (D) The differential expression of the 23 OS-GCAS events between the high- and low-risk groups. (E, F) GSEA of the 50 hallmark gene sets between the high- and low-risk groups.

We then performed GSEA to verify the difference in biological function and pathway between the high- and low-risk groups. The results revealed that some well-known pathways related to GC development, such as “apoptosis”, “hypoxia” and “epithelial-mesenchymal transition”, were significantly enriched in the high-risk group (Figure 5E). Intriguingly, immune-related pathways, such as “inflammatory response”, “interferon-alpha response” and “interferon-gamma response”, were also significantly enriched in the high-risk group (Figure 5F).

Thus, we compared the tumor immune micro-environment between the high- and low-risk groups to comprehensively characterize their different immune features. Both the immune and stromal scores were significantly higher in the high-risk group (Figure 6A). Correspondingly, a lower tumor purity was also found in the high-risk group (Figure 6A). In detail, our studies of immune cell infiltration revealed that many types of immune cells were not randomly distributed across the two groups (Figure 6B, 6C). We found a significantly higher proportion of T cells, monocytes, myeloid dendritic cells, endothelial cells and fibroblasts in the high-risk group (Figure 6B, 6C). However, the proportions of the other five types of immune cells were comparable between the two groups (Figure 6B, 6C). These results indicated the activation of stromal and immune components in the tumor immune microenvironment of high-risk patients together with the activation of oncogenic pathways based on the proposed signatures, which likely contributed at least partially to the poor prognosis of these patients.

Figure 6. Immune microenvironment features underlying the 23-AS event signature. (A) Violin plots for the immune scores, stromal scores, and tumor purity between the high- and low-risk groups. (B) Heat map for immune cell infiltration between the high- and low-risk groups. (C) The differential expression of immune and stromal cells between the high- and low-risk groups.

Integrating the 23-AS event prognostic signature with clinical characteristics

In addition to the 23-AS event prognostic signature, age and TNM stage were also determined as independent prognostic factors in the training dataset, which suggested their complementary value (Table 2). To further improve the prognostic accuracy, we integrated the 23-AS event prognostic signature with the other two independent prognostic factors, age and TNM stage, using the coefficients generated from the multivariate Cox regression model in the training dataset and derived an AS-clinical prognostic model as follows: Risk score = (1.948 × RSASS) + (0.348 × stage) + (0.022 × age). This integrated model was further validated in the validation dataset. Significant improvement in the assessment of survival was achieved with the AS-clinical model relative to the 23-AS event signature only (Figure 7A, 7B).

Figure 7. Construction and evaluation of the AS-clinical nomogram for GC patients. (A, B) RMS curves for the 23-AS event prognostic signature and the AS-clinical signature in the training and validation datasets. P-values represent the difference between the two signatures in terms of the C-index. (C) Nomogram prediction of 1-year, 3-year, and 5-year OS. For stage, 0 means TNM stage I, 1 means TNM stage II, 2 means TNM stage III, and 3 means TNM stage IV. (D) Time-dependent ROC curves for the nomogram at different time points. (E) Calibration curves of observed and predicted probabilities for the nomogram.

Moreover, a nomogram was established for model visualization and clinical application (Figure 7C). The ROC curve confirmed that the nomogram showed great performance in predicting the prognosis of GC (Figure 7D). The calibration curve also presented an optimal prediction for 1-year, 3-year and 5-year OS compared with the actual observations (Figure 7E). All these findings indicated that the nomogram built in our study might contribute significantly to the prediction of prognosis for patients with GC.

Data acquisition and curation process

Patients who met the following criteria were included in The Cancer Genome Atlas (TCGA) GC cohort: (1) histologically confirmed primary GC; (2) alternative RNA splicing data available; (3) detailed clinicopathological and follow-up information available, including sex, age, TNM stage, neoplasm grade, microsatellite status, molecular subtype, and OS status; and (4) an OS time of over 90 days; this latter criterion was applied to, avoid immortal time bias. The alternative RNA splicing data of the TCGA GC cohort were downloaded from the TCGASpliceSeq dataset [37]. Corresponding clinical information and RNA-seq data were downloaded from the Genomic Data Commons data portal using the "GDCRNATools" package [38].

Splicing events in the TCGASpliceSeq dataset were divided into seven categories: exon skip (ES), retained intron (RI), alternate promoter (AP), alternate terminator (AT), alternate donor site (AD), alternate acceptor site (AA), and mutually exclusive exons (ME). Each splicing event was quantified by the percent spliced in (PSI) value, which ranges from 0 to 1 and represents the ratio of normalized read counts to indicate the inclusion of a transcript element over the total normalized reads for that event [39]. To generate a reliable set of AS events, we implemented a series of stringent filters, which included "percentage of samples with a PSI value ≥ 75%" and "average PSI value ≥ 0.05" [28, 29]. The missing PSI values were further filled in using the k nearest neighbors (KNN) algorithm [40].

Identification of GCAS events and OS-GCAS events

To identify GCAS events, we applied the Wilcoxon matched-pairs signed-rank test to compare the PSI values between tumor tissues and matched adjacent normal tissues. The P-value was adjusted by the Benjamini-Hochberg (BH) method. GCAS events were defined as a median PSI value that varied more than 0.1 between tumor tissues and matched adjacent normal tissues with a BH-adjusted P-value < 0.05 [30, 41].

To determine the survival-associated GCAS events, we performed univariate Cox proportional hazards regression analysis to estimate the PSI values of GCAS events with OS. GCAS events with a P-value < 0.05 in the univariate Cox regression analysis were selected as OS-GCAS events.

Functional enrichment analysis and protein interaction network

To investigate the potential functions of OS-GCAS events, we subjected their parent genes to enrichment analyses of Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using the ClueGO plug-in of Cytoscape [42]. For the functional enrichment analysis, a BH-adjusted P-value < 0.05 was considered statistically significant.

To observe the PPIs among the OS-GCAS events, we mapped their corresponding genes to the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database [43]. A minimum required interaction score of 0.4 was used to identify the protein interaction results in STRING, which were further visualized by Cytoscape [44]. In addition, hub genes and specific modules of the protein interaction network were identified by the CytoHubba and Molecular Complex Detection (MCODE) plug-ins of Cytoscape [45, 46]. The maximal clique centrality (MCC) algorithm of CytoHubba was used to predict the top three hub genes in the PPI network. The MCODE options were set as degree cutoff = 2, K-core = 2, and node score cutoff = 0.

Splicing factors (SFs) and the splicing correlation network

A total of 78 genes that participated in the process of alternative RNA splicing (GO:0000380) were obtained from the Molecular Signatures Database (MSigDB) [47]. The read count values of these SFs were extracted from the RNA-seq data, normalized by the trimmed means of M (TMM) method and further transformed by "voom" [48, 49]. Differentially expressed SFs were identified through the "limma" package, and a BH-adjusted P-value < 0.05 was used as the significance threshold [50]. Spearman's correlation analysis was conducted to explore the underlying correlation between the expression of the SFs and the PSI values of the OS-GCAS events. The absolute value of the correlation coefficient ≥ 0.4 with a P-value < 0.05 was considered statistically significant. The splicing correlation network of SFs and the OS-GCAS events was visualized by Cytoscape [44].

Identification of the AS prognostic signature

To enhance the robustness of the AS prognostic signature, the entire TCGA cohort of 327 samples was first randomly distributed into two datasets (7:3), namely, the training and validation datasets. Next, to avoid overfitting in the multivariate model, the least absolute shrinkage and selection operator (LASSO) penalty was applied in the training dataset to build an optimal prognostic signature with the minimum number of OS-GCAS events. Ten-fold cross-validation was conducted to determine the optimal value of the penalty parameter λ, which gives the minimum partial likelihood deviance. Finally, a set of OS-GCAS events and their nonzero coefficients were identified, which was used to build the AS prognostic signature. The time-dependent receiver operator characteristic (ROC) curve and concordance index (C-index) were applied to assess the performance of the prognosis model.

The risk score for the AS signature was calculated for each sample via a linear combination of the selected features, weighted by the corresponding coefficients. Patients were divided into high- and low-risk groups using the cohort-specific median risk score as the cut-off for each dataset. Kaplan-Meier survival analysis was performed to compare patient prognosis between the high- and low-risk groups. Furthermore, the restricted mean survival (RMS) time was also used to compare the prognostic differences between the two groups [51]. The higher the RMS value was, the greater the prognostic difference. A P-value < 0.05 was considered statistically significant.

Association between different risk groups and clinical, molecular, and immune features

Patients in the entire TCGA cohort were divided into high-risk and low-risk groups according to the optimal cut-off value, which was determined by the maximally selected rank statistics test. For the difference between the two groups, clinical features (age, sex, TNM stage, neoplasm grade, and survival status), molecular features (microsatellite status and molecular subtype), and immune features (immune score, stromal score, tumor purity, and immune cell infiltration) were analyzed.

Tumor purity was measured by “ABSOLUTE” [52]; immune and stromal scores were calculated by “ESTIMATE” [53]; and immune cell infiltration was estimated with “MCPcounter” [54]. Kaplan-Meier curve with log-rank test was used for survival data; Mann-Whitney test was used for continuous data; chi-square test was used for categorical data. A P-value < 0.05 was considered statistically significant. In addition, to further explore the difference in biological functions and pathways between high-risk and low-risk groups, gene set enrichment analysis (GSEA) was performed using the “hallmark gene sets” downloaded from MSigDB [47, 55]. A BH-adjusted P-value < 0.05 was considered statistically significant.

Development and verification of an AS-clinical nomogram

Based on the results derived from multivariate analyses, we integrated age, TNM stage, and AS prognostic signature to generate an AS-clinical prognostic model by applying Cox proportional hazard regression in the training dataset. The AS-clinical prognostic model was then applied to the validation dataset for further validation. Next, the prognostic value of the AS-clinical model was compared with that of the AS signature in terms of C-index, and the results were revealed by the RMS curve [18]. RMS represents the life expectancy at five years for patients with different risk scores. Finally, based on the AS-clinical prognostic model, a nomogram was developed to estimate the individual survival probability of patients. The ROC and calibration curves were used to assess the discrimination and calibration ability of the AS-clinical nomogram.