Systematic Analysis of Survival-Associated Alternative Splicing Signatures in Gastrointestinal Pan-Adenocarcinomas.

BACKGROUND: Gastrointestinal pan-adenocarcinomas, which mainly include adenocarcinomas of the esophagus, stomach, colon, and rectum, place a heavy burden on society owing to their poor prognoses. Since aberrant alternative splicing (AS) are starting to be considered as efficacious signatures for tumor prognosis predicting and therapeutic targets, systematic analysis of AS events is urgent.
METHODS: Prognosis-related AS events were selected by using univariate COX regression analysis. Gene functional enrichment analysis revealed the pathways enriched by prognosis-related AS. Then, prognostic signatures based on AS events were developed for prognosis prediction. Potential mechanism to regulate splicing events by splicing factors was analyzed via Pearson correlation and regulatory networks were constructed.
FINDINGS: A total of 967, 918, 674, and 406 AS events were identified as prognosis-related AS events in esophagus, stomach, colon, and rectum adenocarcinomas, respectively. Survival-associated AS events were distinguishing in the four subtypes of adenocarcinoma. Furthermore, computational algorithm results indicated that perturbation of ribosome and ubiquitin-mediated proteolysis pathways were the potential molecular mechanisms corresponding to inferior prognoses. Most notably, several prognostic signatures based on AS events displayed moderate performance in prognosis predicting. The area under curve values of the time-dependent receiver operating characteristic were 0.961, 0.871, 0.870, and 0.890 in esophagus, stomach, colon, and rectum adenocarcinomas. Survival-associated splicing factors were submitted to construct the AS regulatory network, which could be an underlying mechanism of AS events.
INTERPRETATION: AS may could be ideal indiactors in the prognosis of gastrointestinal pan-adenocarcinomas. Exploring interesting splicing regulatory networks is conducive to solve the puzzles of AS.

Evidence Before This Study

With the advantage of high-throughput RNA-seq, the TCGA dataset provides multiple sources for the investigation of whole-genome or transcriptome analyses, including genome splicing exploration. TCGA dataset. SpliceSeq is a java program providing a clear view of inclusion level of each exon and splice junction. Recently, Ryan et al. extended the methodology of SpliceSeq and calculated each potential splicing event across 33 types of cancer to establish the TCGA SpliceSeq database. In the present study, we integrate AS events from SpliceSeq and TCGA clinical prognostic parameters together to comprehensively investigate the prognostic value of the alternative splicing events in gastrointestinal pan-adenocarcinomas.

Added Value of This Study

We are the first group to explore the prognostic value of alternative splicing in gastrointestinal pan-adenocarcinomas. More importantly, new trends in cancer research show that alternative splicing has emerging clinical potential in cancer therapy. However, the understanding of prognostic value of alternative splicing is lacking. Hence, findings in the presents study could provide novel insights into the molecular characteristics of gastrointestinal pan-adenocarcinomas.

Implications of All the Available Evidence

Here, we identified prognosis-related alternative splicing events, which could be the targets for cancer therapy. Furthermore, we constructed several prognostic signatures, which could be excellent in predicting the clinical outcome of gastrointestinal pan-adenocarcinomas. The potential splicing-regulated networks we developed help improve our understanding of the mechanisms of alternative splicing in gastrointestinal pan-adenocarcinomas.

Alt-text: Unlabelled Box

Introduction

Alternative splicing (AS), the mechanism by which a single pre-mRNA molecule produces diverse mature mRNAs, provides the potential for remarkable regulatory and functional complexity in cells [1, 2]. AS is the most important mechanism for expanding protein diversity in terms of the number of genes is limited [3]. Experimental studies have shown that AS plays a decisive role in producing receptor diversity and controlling growth and development [[4], [5], [6]]. Besides being decisive in the regulation of cell differentiation and cell-type-specific functions, abnormal variations in AS are also indispensable in multiple pathological processes, including cancers [[7], [8], [9]]. Accumulated evidence highlights the multifaceted AS events in several carcinogenesis hallmarks, including sustaining proliferative signaling, evading growth suppressors, angiogenesis, invasion, and metastasis [10]. More importantly, new trends in cancer research show that AS has emerging clinical potential in cancer therapy [5, [11], [12], [13]]. Therefore, perturbed homeostasis of AS offers a seedbed for tumor cells and could represent a target for therapy.

Gastrointestinal pan-adenocarcinomas, which originate from the columnar epithelium, mainly include adenocarcinomas of the esophagus, stomach, colon, and rectum [14]. Since these adenocarcinomas share similar endodermal developmental origins, the genomic and other molecular characters across these types of cancer can possess many similarities [15]. However, owing to the complex of tumors, they still have many differences [16]. Gastrointestinal pan-adenocarcinomas accounted for an estimated 183,780 new cases and an estimated 77,280 deaths in the United States in 2018 [17]. Clinically, most individuals are diagnosed at an advanced stage, when effective curative therapies are not satisfactory, which renders advanced gastrointestinal pan-adenocarcinomas one of the most lethal cancers globally [[18], [19], [20], [21]]. Therefore, it is imperative to more deeply explore molecular mechanisms in the prognosis of gastrointestinal pan-adenocarcinomas patients. In particular, AS is a new field, which could lead to deeper understanding of the molecular mechanisms of gastrointestinal pan-adenocarcinomas and provide novel insights.

Despite the indispensable role of AS in oncogenesis, there is little understanding of its clinical significance and potential regulatory mechanism in gastrointestinal pan-adenocarcinomas. Interestingly, globally dysregulated AS events in cancer are prone to be orchestrated by several splicing factors (SFs), especially the serine/arginine-rich (SR) and the heterogeneous nuclear ribonucleoproteins (hnRNP) family. SFs assist spliceosome to recognize and bind to specific sequences of pre-mRNA and subsequently result in mature mRNA [22]. Hence, it is imperative to draw a comprehensive regulatory network of SFs [23, 24]. Considering the close connection between AS and SFs and the fact that they are only superficially understood, it is imperative to explore their prognostic value, as well as the regulatory mechanism in gastrointestinal pan-adenocarcinomas.

To evaluate the potential of AS events in the prognosis prediction of gastrointestinal pan-adenocarcinomas and the AS prognostic network regulated by SFs, we calculated SFs and AS events in Gastrointestinal Adenocarcinomas (GIAC) tissues by analyzing data provided from The Cancer Genome Atlas (TCGA) database. We also identified several moderate prognostic marker panels enriched in gastrointestinal pan-adenocarcinomas. Prognosis-related SF-AS networks in gastrointestinal pan-adenocarcinomas were also constructed to reveal the underlying mechanisms corresponding to this phenotype. These findings have, for the first time, demonstrated novel molecular characteristics by which gastrointestinal pan-adenocarcinomas possess similarities and differences among the four subtypes, according to their AS events.

Materials and Methods

Data Acquisition and Pre-processing

TCGA SpliceSeq includes the mRNA splicing patterns of 33 types of cancers and adjacent normal samples, when available, across a dataset of >10,000 samples; using this tool, investigators can conveniently analyze alternative mRNA splicing patterns (http://bioinformatics.mdanderson.org/TCGASpliceSeq/) [25]. SpliceSeq, a java program, can also be used to calculate the Percent-Spliced-In (PSI) value, which can provide a clear view of the splice junction and the proportion of exons included in different samples. PSI values were used to quantify seven types of AS events [26]: Exon Skip (ES), Retained Intron (RI), Mutually Exclusive Exons (ME), Alternate Donor site (AD), Alternate Acceptor site (AA), Alternate Promoter (AP), and Alternate Terminator (AT); the PSI values of these seven types of AS in patients with esophageal carcinoma (ESCA), stomach adenocarcinoma (STAD), colon adenocarcinoma (COAD), and rectum adenocarcinoma (READ) were download from TCGA SpliceSeq. Only AS events with a PSI value >75% and a standard deviation >0.1 were included in the present analysis. For ESCA, only those patients who were diagnosed with esophageal adenocarcinoma (ESAD), proven by pathology results, were submitted for further analysis. Simultaneously, patient clinical information was also downloaded and extracted from TCGA database.

Survival Analysis

A total of 83 ESAD, 357 STAD, 410 COAD, and 146 READ patients were included in this study. Patients who died within 90 days postoperatively were omitted. Univariate Cox regression analysis was used to calculate the relationships between the PSI values and the overall survival (OS) of patients. AS events were listed as candidate prognosis-related events when the p-value < 0.05. Next, multivariate Cox analysis was employed to identify the possibility of events as independent prognostic factors, and prognostic indexes (PI) were constructed. According to the median PI value, samples were separated into two groups to observe whether they suffered diametrically distinct prognoses. The Kaplan-Meier (K-M) survival analysis was performed to analyze the difference between the two groups. The prognostic value of the PIs was assessed through a commonly used binary response model, the survivalROC package in R software, which is useful for characterizing the predictive accuracy of prognostic markers [27]. Furthermore, concordance index (C-index) was calculated to validate the performance of prognostic predictors we proposed.

UpSet Plot and Gene Interactions Network

The UpSet package of R software was used to visually reveal the interactive events between the seven types of AS. Genes of candidate prognosis-related AS events were submitted to the String 10.5 online database (https://string-db.org) for protein-protein interaction (PPI) analysis [28]. For the sake of more reliable interaction results, we set the threshold to 0.9. Subsequently, the information on the interacting effects of these genes was downloaded and visualized via Cytoscape software. The hub genes in the PPI network were selected based on the number of connections. Then, genes included in the PPI networks were submitted for the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis by using the Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 online database (https://david.ncifcrf.gov/).

Splicing Factor Genes and the Potential Regulatory Network

The global prognosis-related AS events were regulated by a limited number of SFs. A list of SFs was extracted from the SpliceAid 2 (www.introni.it/spliceaid.html) database [29]. Then, the level 3 mRNA-seq expression profiles of the SFs were downloaded from TCGA data portal. The primitive count values were converted into transcripts per million (TPM), which is considered as a more reasonable data format for RNA-seq [30]. univariate Cox regression analysis was conducted to selected survival-associated SFs.

To find potential regulatory relationships between the prognosis-related SFs and the prognosis-related AS events, Pearson's correlation test was performed. Then, the regulatory network between SFs and AS events was generated using Cytoscape.

Results

Survival-Related AS Events

The design mainly includes three parts in the present study: identification of prognosis-related alternative splicing events, development of prognostic signatures and construction of splicing regulatory network (Fig. 1). First, survival analysis was performed in four types of gastrointestinal pan-adenocarcinomas, including 11,470 AS events of ESAD, 12,336 AS events of STAD, 6883 AS events of COAD, and 6985 AS events of READ (Supplementary Fig. 1). For ESAD, a total of 967 AS events in 693 genes were identified as candidate prognosis-related genes, including 62 AAs in 61 genes, 48 ADs in 47 genes, 293 APs in 179 genes, 128 ATs in 73 genes, 353 ESs in 297 genes, 4 MEs in 4 genes, and 79 RIs in 76 genes (Fig. 2A). For STAD, a total of 918 AS events in 649 genes were identified as candidate prognosis-related genes, including 36 AAs in 36 genes, 66 ADs in 62 genes, 302 APs in 184 genes, 120 ATs in 65 genes, 337 ESs in 282 genes, 11 MEs in 11 genes, and 46 RIs in 42 genes (Fig. 2B). For COAD, a total of 674 AS events in 484 genes were identified as candidate prognosis-related genes, including 49 AAs in 49 genes, 51 ADs in 49 genes, 129 APs in 76 genes, 205 ATs in 119 genes, 174 ESs in 154 genes, 3 MEs in 3 genes, and 63 RIs in 60 genes (Fig. 2C). For READ, a total of 406 AS events in 310 genes were identified as candidate prognosis-related genes, including 22 AAs in 22 genes, 27 ADs in 27 genes, 75 APs in 48 genes, 110 ATs in 64 genes, 131 ESs in 122 genes, 2 MEs in 2 genes, and 38 RIs in 35 genes (Fig. 2D). The prognosis related alternative splicing events account for 8.43, 7.44, 9.79 and 5.81% of the total in ESAD, STAD, COAD and READ. And genes with prognosis related alternative splicing events account for 14.06, 12.51, 14.28 and 9.05% of the total in ESAD, STAD, COAD and READ.

Fig. 1

Flowchart of the present study.

Supplementary Fig. 1

The proportion of prognosis-related alternative splicing events. Red represents the number of risk factors while green represents the number of protective factors. And blue was alternative splicing events that have not prognostic value.

(A) Esophageal adenocarcinoma; (B) stomach adenocarcinoma; (C) colon adenocarcinoma; and (D) rectum adenocarcinoma.

Fig. 2

The number of prognosis-related alternative splicing events and involved genes. Red indicates the number of prognosis-related alternative splicing events and green indicates the number of genes with prognosis-related alternative splicing events. Green columns are equally high or higher than the red columns owing to one gene may have up to one or more types of prognosis-related alternative splicing events.

(A) Esophageal adenocarcinoma; (B) stomach adenocarcinoma; (C) colon adenocarcinoma; and (D) rectum adenocarcinoma.

Remarkably, one gene could possess more than one AS event, which was closely related with survival. The number of prognosis-related AS events are visually depicted in Fig. 3.

Fig. 3

UpSet plot of alternative splicing events. One gene may have several types of alternative splicing to be associated with patient survival. UpSet plot of interactions between the seven types of prognosis-related alternative splicing events. One gene may have up to four types of prognosis-related alternative splicing events.

(A) Esophageal adenocarcinoma; (B) stomach adenocarcinoma; (C) colon adenocarcinoma; and (D) rectum adenocarcinoma.

Gene Network Construction and Functional Enrichment Analysis

To reveal the interacting effects among genes of prognosis-related AS, genes with prognosis-related AS events were used to PPI networks construction for ESAD (Fig. 4A), STAD (Fig. 4B), COAD (Fig. 4C), and READ (Fig. 4D), respectively. In the PPI networks, RNA Polymerase II Subunit L (POLR2L), Anaphase Promoting Complex Subunit 11 (ANAPC11), and Transcription Elongation Factor B (SIII), Polypeptide 1 (TCEB1) were the top three hub genes in the ESAD network (Supplementary Fig. 2A). POLR2L, Phosphoinositide-3-Kinase Regulatory Subunit 1 (PIK3R1), and FYN Proto-Oncogene (FYN) were the top three in the STAD network (Supplementary Fig. 2B). Nuclear Cap Binding Protein Subunit 2 (NCBP2), and TCEB1 were the top two in the COAD network (Supplementary Fig. 2C), and TCEB1, Guanine Nucleotide Binding Protein Beta Polypeptide 2-Like 1 (GNB2L1), and Tyrosine 3-Monooxygenase/Tryptophan 5-Monooxygenase Activation Protein Epsilon (YWHAE) were the top three hub genes in the READ network (Supplementary Fig. 2D).

Fig. 4

Protein-protein interaction networks of genes of alternative splicing events. Information of interactions were extracted from online database STRING (http://string-db.org/). The larger, the brighter circles are more important in the network. The thicker lines between two nodules represents the higher combined score.

(A) Esophageal adenocarcinoma; (B) stomach adenocarcinoma; (C) colon adenocarcinoma; and (D) rectum adenocarcinoma.

Supplementary Fig. 2

Hub genes of the protein-protein interaction network. The number represents the connected nodes of the correspondence genes.

(A) Esophageal adenocarcinoma; (B) stomach adenocarcinoma; (C) colon adenocarcinoma; and (D) rectum adenocarcinoma.

KEGG pathways with a p value < 0.05 were considered significant. For ESAD, KEGG analysis revealed that a class of pathways, including “Endocytosis,” “Spliceosome,” and “Ubiquitin mediated proteolysis,” were most significantly enriched by these genes (Fig. 5A). In STAD, there were several essential oncogenic pathways, such as the “AMPK signaling pathway,” the “Wnt signaling pathway” and the “ErbB signaling pathway” (Fig. 5B). “Ribosome,” “Spliceosome,” and “Ubiquitin mediated proteolysis” were the three significant pathways in COAD (Fig. 5C). Similarly, “Ubiquitin mediated proteolysis,” “mRNA surveillance pathway,” “Ribosome,” and “Purine metabolism” were significant pathways in READ (Fig. 5D). These details have been summarized in Table 1.

Fig. 5

Kyoto Encyclopedia of Genes and Genomes analysis. Red circles represent the enriched pathways. The size of the circles represents the number of the gene enriched in the pathway. A greater size indicates a larger number. The colour depth displays P value of pathways. A darker colour indicates a smaller P value. Green circles represent genes.

(A) Esophageal adenocarcinoma; (B) stomach adenocarcinoma; (C) colon adenocarcinoma; and (D) rectum adenocarcinoma.

Table 1

Kyoto Encyclopedia of Genes and Genomes pathways of genes with prognosis associated alterative splicing events (P < 0.05).

Subtype	Term	Count	PValue	Genes
ESCA	hsa04144:Endocytosis	14	3.59E-04	RAB7A, CHMP3, FGFR3, USP8, ERBB3, CHMP7, PML, ARRB2, GRK6, NEDD4L, BIN1, ARAP3, CLTCL1, DNM2
ESCA	hsa04120:Ubiquitin mediated proteolysis	8	0.008014712	FBXW7, HUWE1, PML, SIAH1, TCEB1, NEDD4L, ANAPC11, MID1
ESCA	hsa05410:Hypertrophic cardiomyopathy (HCM)	6	0.009752513	PRKAG2, DAG1, ITGB4, SGCD, TPM2, TPM4
ESCA	hsa00920:Sulfur metabolism	3	0.011513671	SUOX, ETHE1, MPST
	hsa04512:ECM-receptor interaction	6	0.015163369	LAMA4, LAMA3, COL6A3, DAG1, ITGB4, FN1
ESCA	hsa03010:Ribosome	7	0.026488508	RPL18A, RPL32, RPLP0, RPS15, RPL37A, RPL28, RPS7
ESCA	hsa04062:Chemokine signaling pathway	8	0.03651303	AKT1, PLCB3, GNAI2, ARRB2, PTK2B, BCAR1, GRK6, ELMO1
ESCA	hsa05100:Bacterial invasion of epithelial cells	5	0.04163531	BCAR1, CLTCL1, ELMO1, FN1, DNM2
ESCA	hsa05145:Toxoplasmosis	6	0.047911017	IRAK4, AKT1, IRAK1, LAMA4, LAMA3, GNAI2
ESCA	hsa04921:Oxytocin signaling pathway	7	0.049384014	MYL6, PLCB3, GNAI2, PRKAG2, PPP3CB, CAMKK1, NFATC1
STAD	hsa04520:Adherens junction	9	1.13E-04	MAP3K7, ACTG1, SORBS1, FYN, LMO7, CTNND1, CDH1, INSR, TCF7L2
STAD	hsa04152:AMPK signaling pathway	10	0.001052524	MAP3K7, PRKAG2, TSC2, PPARG, STRADA, ADRA1A, TBC1D1, INSR, PIK3R1, CAMKK2
STAD	hsa05410:Hypertrophic cardiomyopathy (HCM)	8	0.001228388	ACTG1, DMD, PRKAG2, ITGB4, CACNB3, TPM2, TPM1, SGCA
STAD	hsa04810:Regulation of actin cytoskeleton	12	0.004917302	FGFR2, ARHGEF4, ACTG1, ENAH, FGFR4, ARHGEF7, BCAR1, FGF11, ITGB4, PIP5K1A, BDKRB2, PIK3R1
STAD	hsa05100:Bacterial invasion of epithelial cells	7	0.005913666	ACTG1, BCAR1, GAB1, ILK, CDH1, SHC1, PIK3R1
STAD	hsa04660:T cell receptor signaling pathway	8	0.005952153	MAP3K7, PTPRC, FYN, NFATC2, GRAP2, MAP2K7, PIK3R1, NFATC1
STAD	hsa04921:Oxytocin signaling pathway	10	0.006097389	MYL6, ACTG1, CAMK2G, PRKAG2, CACNB3, NFATC2, PLA2G4B, PIK3R1, CAMKK2, NFATC1
STAD	hsa04380:Osteoclast differentiation	9	0.006365485	MAP3K7, IL1R1, FYN, PPARG, NFATC2, MAP2K7, IFNGR2, PIK3R1, NFATC1
STAD	hsa05416:Viral myocarditis	6	0.006812455	ACTG1, FYN, DMD, CASP8, ABL2, SGCA
STAD	hsa05414:Dilated cardiomyopathy	7	0.008444102	ACTG1, DMD, ITGB4, CACNB3, TPM2, TPM1, SGCA
STAD	hsa05412:Arrhythmogenic right ventricular cardiomyopathy (ARVC)	6	0.016761239	ACTG1, DMD, ITGB4, CACNB3, TCF7L2, SGCA
STAD	hsa04310:Wnt signaling pathway	8	0.026545061	MAP3K7, CTBP2, CCND3, CAMK2G, NFATC2, TCF7L2, NFATC1, DVL1
STAD	hsa04015:Rap1 signaling pathway	10	0.033404171	FGFR2, ACTG1, FGFR4, BCAR1, EFNA3, FGF11, CTNND1, CDH1, INSR, PIK3R1
STAD	hsa04012:ErbB signaling pathway	6	0.036477432	CAMK2G, GAB1, SHC1, MAP2K7, ABL2, PIK3R1
STAD	hsa05166:HTLV-I infection	11	0.044097891	IL1R1, CCND3, RAN, CREM, RANBP3, ANAPC11, KAT5, NFATC2, PIK3R1, NFATC1, DVL1
STAD	hsa04014:Ras signaling pathway	10	0.049482855	FGFR2, FGFR4, EFNA3, GAB1, FGF11, SHC1, ABL2, INSR, PLA2G4B, PIK3R1
COAD	hsa03010:Ribosome	8	0.001348748	RPS25, RPL30, RPL34, RPL21, RPL35, RPS21, MRPL35, RPS24
COAD	hsa03040:Spliceosome	6	0.023574399	NCBP2, SRSF5, PRPF3, U2AF1L4, PRPF40A, TXNL4A
COAD	hsa04120:Ubiquitin mediated proteolysis	6	0.02638476	UBE2D3, UBE2Z, MGRN1, TCEB1, UBE2C, MID1
READ	hsa04120:Ubiquitin mediated proteolysis	7	6.43E-04	UBE3B, UBA1, TCEB1, NEDD4L, UBOX5, ITCH, CUL4B
READ	hsa03015:mRNA surveillance pathway	5	0.00538823	PAPOLA, SMG7, PPP2R5C, CPSF4, WDR33
READ	hsa03010:Ribosome	5	0.0212527	RPS25, FAU, RPS9, RPS20, RPS3
READ	hsa00230:Purine metabolism	5	0.047997843	POLR3G, ENTPD8, AK2, PDE9A, PDE4D

Splicing Characteristics of Gastrointestinal Pan-Adenocarcinomas Subtypes

Similarities and differences in the four subtypes of gastrointestinal pan-adenocarcinomas are worth exploring. Surprisingly, no commonality was found in the AS events of these four subtypes of gastrointestinal adenocarcinomas (Fig. 6A). No AS events was significant in all the four subtypes. The AS events of two genes, Thioredoxin Related Transmembrane Protein 2 (TMX2) and METTL23 (Methyltransferase Like 23), were present in all four subtypes concurrently (Fig. 6B). In the KEGG pathway analysis, the “Ubiquitin mediated proteolysis” and “Ribosome” pathways were closely related to the prognosis of ESAD, COAD and READ, but not STAD.

Fig. 6

Venn plot of alternative splicing events. Red stands for colon adenocarcinoma (COAD), light blue means rectum adenocarcinoma (READ), purple means stomach adenocarcinoma (STAD), and green means studious and green represents esophageal carcinoma (ESCA).

(A) A Venn diagram shows the overlap of prognosis-related alternative splicing events in COAD, READ, STAD and ESCA; and (B) A Venn diagram shows the overlap of genes with prognosis-related alternative splicing events in COAD, READ, STAD and ESCA.

Construction of AS-Based Prognostic Signatures

To facilitate the application of AS when clinically monitoring the prognosis of GIAC patients, we developed risk prediction formulas based on the top ten most significant prognosis-related events among the seven types of AS. First, the signatures based on the seven types of AS events in ESAD, STAD, COAD and READ were constructed, respectively. The formulas are summarized in Table 2. The Kaplan-Meier analysis showed that these prognostic signatures could significantly separate patients with distinct prognosis (Fig. 7).

Table 2

General characteristics of prognostic predictors for GIACs.

Tumor type	Alternative splicing	Formula	Hazard ratios (95% CI)	ROC	C-index
Esophageal adenocarcinoma	AA	“U2AF1L4-49,280-AA” * 0.02105 + “TICRR-32428-AA” * (−0.01414) + “RSRC2-24968-AA” * (−0.02971) + “PREPL-53439-AA” * 0.02192 + “PPIL2-61247-AA” * 0.03579 + “FAM135A-76,637-AA” * 0.03004 + “CDV3-66839-AA” * 0.04154 + “ABCB7-89517-AA” * (−0.02951)	4.492 (2.333–8.647)	0.897	0.798
	AD	“ZNF384-19927-AD” * (−0.05812) + “RPP14-65434-AD” * (−0.03077) + “PQBP1-89028-AD” * 0.02678 + “MFSD11-43690-AD” * 0.02587 + “COX6C-84,682-AD” * 0.04491	4.490 (2.362–8.534)	0.745	0.800
	AP	“ZNF623-85469-AP” * (−0.04535) + “KIAA0513-37876-AP” * 0.01779 + “FAM19A5-62,732-AP” * 0.04833 + “ALDH6A1-28,367-AP” * 0.03834	4.366 (2.276–8.376)	0.817	0.794
	AT	“TRIM4-80864-AT” * (−2.67e-02) + “RNASEH2B-25,927-AT” * (−2.61e + 02) + “RNASEH2B-25,926-AT” * -2.61e + 02 + “MCPH1-82574-AT” * 1.72e-02 + “ARL6-65732-AT” * 4.29e-02 + “AHI1-77886-AT” * (−2.30e-02)	4.143 (2.219–7.734)	0.831	0.778
	ES	“TNC-87345-ES” * 0.03180 + “PML-31651-ES” * -0.01033 + “NBPF15-91080-ES” * (0.07165 + “MYL6-22384-ES” * (−0.04943) + “MRPL43-12857-ES” * -0.03976 + “IRF9-117161-ES” * (−0.01943)	4.235 (2.242–8.001)	0.904	0.843
	ME	“SDR39U1-27,012-ME” * (−0.01054) + “KLHL2-71038-ME” * (−0.02045) + “CMC2-37707-ME” * 0.02191	3.522 (1.858–6.674)	0.864	0.690
	RI	“ZNF131-71926-RI” * (−0.03666) + “SLC52A3-58,464-RI” * 0.03577 + “PPARGC1B-74,051-RI” 0.01544 + “PCGF3-68404-RI” * 0.02876 + “MDK-15570-RI” * 0.05833 + “MAF-37687-RI” * 0.05197 + “FAM9C-88,504-RI” * (−0.02045)	5.643 (2.97–10.72)	0.827	0.795
Stomach adenocarcinoma	AA	“RPLP0-24727-AA” * (−0.01196) + “NAT6-64990-AA” * 0.01181 + “MRVI1-14373-AA” * 0.00731 + “LMO7-26065-AA” * 0.01256 + “BDKRB2-29192-AA” * (−0.00776)	2.299 (1.649–3.205)	0.718	0.639
	AD	“YIPF2-47605-AD” * (−0.01345) + “SPHK2-50793-AD” * 0.01230 + “SENP1-21411-AD” * 0.00842 + “PGAP2-14004-AD” * (−0.01311) + “NFATC1-46241-AD” * 0.00822 + “CCDC51-64653-AD” * (−0.00669)	2.048 (1.468–2.858)	0.773	0.630
	AP	“RCAN1-60494-AP” * 0.01222 + “PLCD1-64009-AP” * 0.00891 + “LTBP1-53179-AP” * (−0.00835) + “FAM65B-75,537-AP” * 0.02503 + “ABL2-9101-AP” * (−0.00857)	2.493 (1.790–3.472)	0.646	0.665
	AT	“ZNF846-47399-AT” *7.44e-03 + “ZFYVE28-68559-AT” (−2.01e-02) + “STEAP4-80362-AT” * (−5.58e + 01) + “STEAP4-80361-AT” * (−5.58e + 01) + “KIF1B-602-AT” *5.58e + 01 + “KIF1B-601-AT” *5.58e + 01 + “CXCL12-11344-AT” * (−1.50e-02) + “CLDN11-67617-AT” * (−1.41e-03) + “ABCB5-78909-AT” *1.67e-02 +	2.647 (1.893–3.700)	0.773	0.658
	ES	“UBXN11-1263-ES” * (−0.02472) + “TMEM230-58637-ES” * 0.01826 + “SRSF3-75985-ES” * 0.01602 + “SORBS1-12641-ES” * 0.00880 + “P4HA2-73,263-ES” * (−0.01086) + “CREM-11245-ES” * (−0.01342)	2.461 (1.765–3.432)	0.801	0.679
	ME	“N4BP2L1-25,590-ME” * (−0.00781) + “KDM6A-98,323-ME” * (−0.00884) + “FYN-77273-ME” * 0.00995 + “CCDC53-106010-ME” * 0.00973	1.491 (1.070–2.077)	0.631	0.590
	RI	“TREX1-64682-RI” * 0.01573 + “SRSF7-53276-RI” * (−0.02285) + “RPS15-46490-RI” * (−0.01897) + “LDHA-14642-RI” * 0.00548 + “BICD2-86883-RI” * (−0.00777) + “ALS2CL-64,462-RI” * 0.00983	2.842 (2.036–3.968)	0.800	0.656
Colon adenocarcinoma	AA	“RASSF7-13691-AA” * 0.02897 + “PTGR1-87219-AA” * 0.02644 + “FAM173A-32,964-AA” * 0.02221 + “DPP3-17040-AA” * (−0.01817) + “CDV3-66842-AA” * 0.02860	3.983 (2.623–6.048)	0.737	0.723
	AD	“RNF14-73855-AD” * 0.01562 + “IP6K2-64,759-AD” * 0.01261 + “HPS4-61506-AD” * 0.02158 + “HDGF-8323-AD” * 0.01968 + “ANKRD46-84712-AD” * 0.02147 + “ADPGK-31594-AD” * 0.01645	3.224 (2.126–4.888)	0.680	0.676
	AP	“TUBB3-38167-AP” * (−0.01030) + “RAB3IP-23,345-AP” * (−0.03643) + “MAZ-35938-AP” * 0.02384 + “FADS2-16289-AP” * 0.01878 + “ENO2-20011-AP” * 0.00736	3.436 (2.268–5.206)	0.710	0.681
	AT	“ZNF765–51718-AT” * (−4.08e-02) + “UPK3B-80,182-AT” * 2.05e-02 + “RASEF-86677-AT” * 1.37e + 02 + “RASEF-86676-AT” * 1.37e + 02 + “NRG4-31911-AT” * 3.04e-02 + “AIG1-77972-AT” * (−1.62e-02)	3.262 (2.147–4.958)	0.752	0.720
	ES	“VTI1B-28,083-ES” * (−0.01045) + “STRN3-27098-ES” * (−0.02262) + “RHOC-4236-ES” * 0.01902 + “PRMT1-51042-ES” * 0.01818 + “PLEKHM2-767-ES” * 0.01972 + “DMWD-50528-ES” * 0.02082 + “D2HGDH-58,423-ES” * 0.02800	3.013 (1.981–4.582)	0.769	0.730
	ME	“CNOT10-63822-ME” * 0.0236	1.479 (0.976–2.242)	0.585	0.561
	RI	“ZNF226-50290-RI” * 0.01131 + “NPIPA5-34148-RI” * 0.02178 + “ELP5-38889-RI” * (−0.03011) + “ALS2CL-64,463-RI” * 0.01968	2.711 (1.788–4.110)	0.696	0.679
Rectum adenocarcinoma	AA	“ZNF467-82205-AA” * −0.01766 + “RNPC3-3907-AA” * −0.03513 + “GGT1-61440-AA” * −0.01812 + “BTN3A1-75,660-AA” * −0.01810	4.321 (1.902–9.819)	0.728	0.766
	AD	“OSBPL9-2975-AD” * 0.01167 + “METTL23-43637-AD” * 0.03522 + “BCS1L-57,522-AD” * (−0.02063)	3.177 (1.401–7.205)	0.794	0.718
	AP	“TADA2B-68,732-AP” * (−0.01496) + “PTCH1-86955-AP” * 0.02493 + “DAB2IP-87,442-AP” * 0.02235	4.897 (2.162–11.09)	0.860	0.694
	AT	“PUS10-53676-AT” * 0.0635 + “NOTCH2NL-4437-AT” 0.0541	4.592 (2.015–10.46)	0.912	0.738
	ES	“SPAG9-42496-ES” * −0.0325 + “SERPINA1-29134-ES” * 0.0455 + “PHB2-20048-ES” * (−0.0223) + “FGFR1OP2-20,856-ES” * 0.0529	8.372 (3.577–19.6)	0.832	0.797
	ME	“RBMS2-22465-ME” * (−0.01912)	1.742 (0.769–3.947)	0.574	0.591
	RI	“ZNF692-10557-RI” * −0.05799 + “WDR33-55246-RI” *0.03283 + “TMEM91-50046-RI” * 0.00968 + “SIDT2-18886-RI” * 0.03883 + “EXOSC9-70501-RI” * (−0.02814) + “ADARB1-60863-RI” * (−0.03719)	8.766 (3.841–20)	0.937	0.767

Hazard ratios were estimated between high- and low-risk groups.

Fig. 7

Kaplan-Meier survival analysis of prognostic signatures based on seven types of alternative splicing events in gastrointestinal tract adenocarcinomas.

Prognostic signatures based on Alternate Acceptor site for esophageal adenocarcinoma (A), stomach adenocarcinoma (B), colon adenocarcinoma (C), and rectum adenocarcinoma (D).

Prognostic signatures based on Alternate Donor site for esophageal adenocarcinoma (E), stomach adenocarcinoma (F), colon adenocarcinoma (G), and rectum adenocarcinoma (H).

Prognostic signatures based on Alternate Promoter for esophageal adenocarcinoma (I), stomach adenocarcinoma (J), colon adenocarcinoma (K), and rectum adenocarcinoma (L).

Prognostic signatures based on Alternate Terminator for esophageal adenocarcinoma (M), stomach adenocarcinoma (N), colon adenocarcinoma (O), and rectum adenocarcinoma (P).

Prognostic signatures based on Exon Skip for esophageal adenocarcinoma (Q), stomach adenocarcinoma (R), colon adenocarcinoma (S), and rectum adenocarcinoma (T).

Prognostic signatures based on Mutually Exclusive Exons for esophageal adenocarcinoma (U), stomach adenocarcinoma (V), colon adenocarcinoma (W), and rectum adenocarcinoma (X).

Prognostic signatures based on Retained Intron for esophageal adenocarcinoma (Y), stomach adenocarcinoma (Z), colon adenocarcinoma (AA), and rectum adenocarcinoma (AB).

Then, we integrated the seven kinds of AS events and developed comprehensive prognostic signatures for ESAD, STAD, COAD, and READ. The final survival risk signatures integrating all types of AS events displayed moderate performance in predicting prognosis. With respect to ESAD, patients in the high-risk group presented a low OS (hazard ratio (HR) = 8.206, 95% confidence interval (CI): 4.264–15.79, p < 0.0001; Fig. 8A). Regarding STAD, patients in different groups showed a statistically significant association with OS (HR = 4.449, 95% CI: 3.179–6.228, p < 0.0001; Fig. 8B). AS-based signatures were well practiced in separating patients with COAD into two groups with distinct outcomes (HR = 8.624, 95% CI: 5.666–13.13, p < 0.0001; Fig. 8C) and READ (HR = 7.717, 95% CI: 3.4–17.51, p < 0.0001; Fig. 8D). The time-dependent receiver operating characteristic curve (ROC) analysis validated that the integrated prognostic signatures displayed a more ideal performance than other models based on a single specific type of AS event. The area under curve (AUC) value of prognostic signatures with all types were 0.961, 0.871, 0.870, and 0.890 in ESAD, STAD, COAD, and READ, respectively (Supplementary Fig. 3). C-indexes were 0.871, 0.781, 0.840 and 0.936, in ESAD, STAD, COAD and READ respectively. These findings suggested that the four prognostic signatures could be ideal prognostic predictors. The four computational prognostic models displayed moderate ability in the classification of patient survival (Fig. 9).

Fig. 8

Prognostic signatures based on all types of alternative splicing events in gastrointestinal tract adenocarcinomas. Patients were divided into high- and low-risk groups according to prognostic signatures. Each panel contains three parts: [1] survival differences were estimated by Kaplan-Meier survival curve; [2] number of patients in different groups; and [3] number of censoring at different times.

(A) Esophageal adenocarcinoma; (B) stomach adenocarcinoma; (C) colon adenocarcinoma; and (D) rectum adenocarcinoma.

Supplementary Fig. 3

Time-dependent receiver operating characteristic curves of four prognostic signatures in gastrointestinal tract adenocarcinomas.

(A) Esophageal adenocarcinoma; (B) stomach adenocarcinoma; (C) colon adenocarcinoma; and (D) rectum adenocarcinoma.

Fig. 9

The four prognosis-relevant splicing event groups influencing the survival of gastrointestinal tract adenocarcinoma patients. X-axis represents the orders of patients based on prognostic signatures and Y-axis represents the survival days of patients. Dotted lines were used to distinguish patients in high-risk group and low-risk group.

(A) Esophageal adenocarcinoma; (B) stomach adenocarcinoma; (C) colon adenocarcinoma; and (D) rectum adenocarcinoma.

The prognostic signatures stratified patients into high- and low-risk groups. We observed that the four prognostic signatures could divided patients with early-stage (I and II) cancers and advanced stage (III and IV) into significantly different prognostic groups (Fig. 10).

Fig. 10

Kaplan-Meier curves of overall survival (OS) among patients with low stage and advanced tumor stage.

Esophageal adenocarcinoma (A) early stage, (B) advanced stage. Stomach adenocarcinoma (C) early stage, (D) advanced stage. Colon adenocarcinoma (E) early stage, (F) advanced stage. Rectum adenocarcinoma (G) early stage, (H) advanced stage. Patients in high-risk group suffer poor OS whether with early or advanced tumor stage.

Survival-Associated SFs

A total of 66 SFs collected from the SpliceAid 2 database were submitted for survival analysis. Respectively, 16, 2, 1, and 4 survival-associated SFs were obtained in ESAD, STAD, COAD and READ (Fig. 11). Then, a Pearson correlation coefficient was used to estimate the prospective relationships between prognosis-related SFs and AS events. In ESAD and STAD, most SFs were negatively correlated with protective AS events, while positively correlated with risk AS events. However, in READ, the results were the exact opposite (Fig. 12).

Fig. 11

The prognostic analysis of splicing factors in four subtypes of gastrointestinal pan-adenocarcinoma. Red circles represent splicing factors with hazard ratio >1 and blue circles represent splicing factors with hazard ratio <1. Circles with black edge indicate significant (P < 0.05).

(A) Esophageal adenocarcinoma; (B) stomach adenocarcinoma; (C) colon adenocarcinoma; and (D) rectum adenocarcinoma.

Fig. 12

Survival-associated splicing factors and the splicing correlation network in gastrointestinal tract adenocarcinomas. Expression of prognosis-related splicing factors (green dots) were positively (red line) or negatively (blue line) correlated with the PSI values of prognosis-related alternative splicing events.

(A) Esophageal adenocarcinoma; (B) stomach adenocarcinoma; (C) colon adenocarcinoma; and (D) rectum adenocarcinoma.

Discussion

Preliminary investigations revealed that AS perturbation was involved in the initiation and progression of several cancers [[31], [32], [33]]. In the present study, we identified prognosis-related AS events in ESAD, STAD, COAD, and READ patients. We also found that the splicing characteristics in four subtypes of gastrointestinal pan-adenocarcinoma have significant individual variations. Moreover, we observed that genes of survival-associated AS events were mainly involved in several metabolism pathways and interacted closely with each other. Most notably, we performed subdividing according to unique individual AS events and found that AS events could be potential prognostic factors for GIAC patients. SFs-AS event networks provided a potential regulatory mechanism for the abnormal changes in gastrointestinal pan-adenocarcinomas. Our comprehensive and integrated computational investigation first focused on the AS event characteristics of gastrointestinal pan-adenocarcinomas and then broadened to the novel visual field of prognostic and molecule-targeted implications.

Normal tissues can precisely control genomic stability, thereby maintaining the usually low spontaneous mutation and keeping normal physiological function. However, the development of genomic instability is a typical hallmark of cancerous cells [34]. Particularly, a high proportion of human genetic disorders result from AS events [35]. Hence, abnormal splicing variants actively participate in the development of cancers [36]. For example, VEGFB is an antiangiogenic isoform of VEGF, and it subverts the understanding of the angiogenic role of VEGF. The balance of pro- and antiangiogenic isoforms of VEGF is disturbed between tumor and non-tumor tissues [37]. Further, ZAKα and ZAKβ are two isoforms of ZAK, and they toward contrary way in the proliferation of cancers. ZAKα exerts an anti-neoplastic effect, while ZAKβ has essential antiproliferation properties [38]. However, the onset and progression of gastrointestinal pan-adenocarcinomas is supposed to suffer substantial aberrant AS events [12].

The recent advance in high-throughput technology is conductive to providing a valuable overview of aberrant AS events on a genome-wide scale. In the present study, a series of survival-associated AS events were identified, which provide several prognosis monitoring indexes for gastrointestinal pan-adenocarcinomas. This study is comparatively novel in its topic and methodology. We are the first group to explore the relationships between AS events and prognoses of gastrointestinal pan-adenocarcinomas patients.

We attempted to reveal the similarities and differences of gastrointestinal pan-adenocarcinomas according to their splicing characteristics. Surprisingly, no commonality was achieved among the AS events of the four subtypes of gastrointestinal adenocarcinomas studies. Notably, for those genes that generate prognosis-related AS events, gene functional enrichment was performed, which indicated that “ribosome” and “ubiquitin mediated proteolysis” may be the most significant interfered pathways related to AS. The protein degradation process occurs in two main ways: the lysosomal degradation pathway and ubiquitin mediated proteolysis [39, 40]. Abnormal proteolytic activity is associated with many diseases, especially in cancers [41]. Weatheritt et al. [42] proposed that principal components of human exon AS events, which have been detected in transcripts with medium-to-high abundance, are engaged by ribosomes and therefore likely translated. Ubiquitin mediated proteolysis is a major process to degrade protein in cells. It plays a crucial role in maintaining cellular homeostasis and metabolism [43, 44]. More importantly, the stability of ubiquitin mediated proteolysis is crucial for cell cycle and apoptosis. Therefore, we have speculated that the clinical outcome of patients results from AS events may be disturbed by confused protein degradation.

In view of the significant morbidity and mortality of gastrointestinal pan-adenocarcinomas, deeper mining and the development of prognostic signatures is urgent. Several researchers have proposed prognostic signatures for gastrointestinal pan-adenocarcinomas based on several types of molecules, such as lncRNAs [45], mRNAs [46], and miRNAs [47, 48]. With the advances in high-throughput RNA-seq, TCGA dataset provides multiple resources for the investigation of AS events at the genome-wide level. To illustrate the prognostic value of splicing events in cancers, several researchers have identified different prognostic subtypes for non-small cell lung cancer [49] and ovarian cancer [50] based on AS events. They have proposed that splicing events could be preeminent biomarkers for predicting cancer prognoses. Hence, we integrated AS events and clinical outcome data into the comprehensive mining of prognosis-related AS events in gastrointestinal pan-adenocarcinomas. The signatures we proposed are ideal for predicting prognoses.

We also constructed an SF-AS regulatory network. SFs are the regulatory elements of AS events [51]. We constructed the regulatory network and proposed several prognostic SFs. Indeed, the entire regulatory network of AS events is quite complex, and AS events have far more regulators than do SFs. Altered SFs play a crucial driving role in pathological splicing events [52]. However, the potential molecular mechanisms are unclear and lack the comprehensiveness of regulatory relationships between SFs and AS events. This phenomenon indicates that multiple splicing factors can affect the survival of gastrointestinal pan-adenocarcinomas patients by synergistically regulating the alternative splicing events of genes.

However, several limitations should be considered. First, the number of patients included in the ESAD and READ cohorts were limited. Second, no another independent cohort of gastrointestinal pan-adenocarcinomas patients has been used to show that the prognostic models being proposed here are reproducible. Third, the present in-silico analysis should be verified in the future.

In summary, these results highlight the prognostic value of AS events and SFs and explore potential regulatory mechanisms. An extensive amount data at the genome-wide level has uncovered the general implications of AS events in several aspects of gastrointestinal pan-adenocarcinomas, particularly in prognosis prediction, which may reveal new opportunities for targeted therapies for these cancers.

The following are the supplementary data related to this article.