Identification of potential biomarkers and their clinical significance in gastric cancer using bioinformatics analysis methods

Data acquisition

A total of 407 samples (375 GC samples and 32 normal samples) were enrolled for comprehensive integrated analysis. The data were download from The Cancer Genome Atlas (TCGA) database. In addition, we used the Data Transfer Tool (provided by GDC Apps) to download the level 3 mRNASeq gene expression data and clinical information of those patients. In addition, data of splicing events were obtained from the SpliceSeq database (Ryan et al., 2012). The filter condition is sample percentages with PSI values ≥75 and ΔPSI ≥ 30% (Ryan et al., 2012). Finally, the resulting matrix files and PSI files are used for subsequent analysis.

Gene Set Enrichment Analysis (GSEA)

GSEA abandons the previous method that the analysis software only focuses on a group of up-regulated or down-regulated genes, and focuses on a group of genes with the same or similar biological processes, through a comparative analysis of the overall changes of a group of genes, and then explain the effects of different treatments on the sample or reveal the biological significance (Subramanian et al., 2007). GSEA was used to enrich key Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of high and low ECT2 expression in GC.

Oncomine database analysis

The expression level of the ECT2 gene in various types of cancers was identified in the Oncomine database (Rhodes et al., 2007). The threshold was determined according to the following values: P-value of 0.001, fold change of 1.5, and gene ranking of all.

Quantitative reverse transcription polymerase chain reaction (qRT- PCR) assays

Total RNA from cells or tissues was isolated using TRIzol (Invitrogen, Canada) reagent, the specific operation is carried out with reference to the instructions for the operation of the kit. RNA (1 µg) was converted into cDNA using the RevertAid First Strand cDNA Synthesis Kit (Takara, China). qRT-PCR was performed using SYBR Green Mixture (Takara, China) in the ABI StepOne-Plus System (ABI7500, USA). Target gene expression was normalized against GAPDH. The primer sequences are F: 5′-CAGACTCCGAAGGAAGTTGTATG-3′,R:5′-TCCACTGAGCCGTGGGATGTCA-3′.

Statistical analyses

We used the R packages (“UpSetR”) to get a overview of AS events profiling in GC (Conway, Lex & Gehlenborg, 2017). We combine the survival data with the AS data to obtain survival-related AS data for subsequent analysis, then use the R package (“UpSetR and survival”) to analyze the variable shear events associated with survival. A univariate Cox regression analysis was employed by “survival” package in R to identify survival-associated splicing events. Next, a multivarite Cox regression was analyzed upon meaful genes (p < 0.05) screened from univariate regression. Subsequently, Cytoscape 3.5 was employed to construct the potential regulatory network. R language packages (ggplot2, pheatmap, pROC, and corrgram) are used for other statistical computations and figure drawing.

Result

The landscape of AS events in gastric cancer

60,754 AS events of 22,039 genes were identified. In detail, we detected 31,730 exon skip (ES) in 6,973 genes, 8,393 alternate terminator (AT) in 3,666 genes, 10,005 alternate promoter (AP) in 4,025 genes, 4,006 alternate acceptor site (AA) in 2,799 genes, 3,450 alternate donor site (AD) in 2,401 genes, 2,944 retained intron (RI) in 1,956 genes, and 226 mutually exclusive exons in 219 genes respectively (Figs. 1A and 1B). After strict filtering by survival, 32,129 AS events of 17,913 genes were identified, including 12,894 ESs in 5,598 genes, 5,518 ATs in 3,192 genes, 6,253 APs in 3,620 genes, 2,817 AAs in 2,111 genes, 2,410 ADs in 1,782 genes, 2,182 RIs in 1,469 genes, and 145 MEs (mutually exclusive exons) in 141 genes respectively (Figs. 1C and 1D). These data suggest that a single gene may have multiple types of mRNA splicing events. ES is the primary splicing event and ME is a rare splicing event in GC patients.

To better understand the prognostic role of AS event in GC, we used a multivariate Cox regression analysis. Multivariate Cox regression showed that 2,042 AS events of 1,615 genes, including 805 ESs in 660 genes, 297 ATs in 203 genes, 461 APs in 304 genes, 157 AAs in 153 genes, 174 ADs in 164 genes, 130 RIs in 113 genes, and 18 MEs in 18 genes respectively, were significantly associated with OS (P < 0.05) (Figs. 1E and 1F). To better understand the biological processes of the OS-associated genes, we annotated their function using gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. The results showed that OS-associated genes are enriched in GTPase activity and PI3K-Akt signaling pathway (Figs. 1G and 1H). Obviously, most of the top 20 significant OS-associated AS events were better prognostic factors (Z < 0) (Figs. 2A–2G). For instance, AA of ECT2, LMO7, STAT3, CBX7, TRAPPC2L, TSC2,TROAP, ZNF410 and HNRNPR were adverse prognostic factors in GC patients, however others were better prognostic factors (Fig. 2A). Table 1 shows the top 15 most significant AS events for up- and down-regulation.

Construction of the prognostic predictor in GC patients

Next, risk score constructed using the top 20 significant OS-related AS events of the eight types, identified by multivariate Cox proportional hazards regression. As the patient’s risk score increases, the number of dead patients increases, indicating that the risk score is related to survival. At the same time, as the patient’s risk score increases, the PSI value of AS increases (Figs. 3A–3X). For instance, the PSI value of AP RCAN1-60494 increases as the risk value increases (Fig. 3C).

Next, the area under the curve (AUC) of ROC was generated, the result showed that risk score exhibited the AUC of 0.841 in AA, followed by AD, ES, All, ME, AP, AT and RI model with AUC of 0.827, 0.818, 0.801, 0.765, 0.759, 0.756 and 0.716, respectively (Figs. 4A–4P). The K-M curve was used to analyze the survival time of patients in the low-risk and high-risk groups. The results show that in each Cox regression model constructed from eight types of AS events, high risk score had a poor survival (Figs. 4A–4P). In addition, Both the univariate (HR: 1.71; 95% CI [1.51–1.95]) and multivariate Cox regression analyses (HR: 1.64; 95% CI [1.44–1.86]) results indicated that the risk score and age were all correlated with the OS (Figs. 5A and 5B).

Next, to determine which SF is associated with AS events associated with survival in the GC, we performed a survival analysis of SF. The results showed that 26 SF was significantly associated with overall survival. In addition, the correlation between the PSI value of significant AS events and the expression of survival-related SF was investigated using the Spearman test (Fig. 5C).

Construct survival-associated alternative splicing events

The counterpart genes of risk score (AA) were ST5, PLK4, BDKRB2, NAT6, APOBEC3B, ECT2, MORF4L2, STAT3, PARPBP, CBX7, TRAPPC2L, C19orf60, DHPS, TROAP, ZNF410, and HNRNPR (Table 2). Compared with normal tissues, GC patients generally contain a lower proportion of ST2 and CBX7 (P < 0.05), while the mRNA expression of PLK4, NAT6, APOBEC3B, ECT2, MORF4L2, STAT3, PARPBP, TRAPPC2L, TROAP, and HNRNPR were opposite (P < 0.05). Compared with normal tissues, GC patients have no considerable difference of BDKRB2, ZNF410, and DHPS (Figs. 6A and 6B). K-M analysis was used to evaluate the association between mRNA expression of 15 genes and OS, and the result indicate high expression of ECT2 predicts poor prognosis in GC patients (Figs. 6C–6G).

High expression of ECT2 predicts poor prognosis in GC patients

Next, both the univariate (HR: 1.32; CI [1.11–1.56]) and multivariate Cox regression analysis (HR: 1.26; CI [1.06–1.51]) results indicated that high ECT2 expression correlated significantly with a poor overall survival (Fig. 7A). To identify signaling pathways that are differentially activated in GC, we conducted GSEA between low and high ECT2 expression data sets. We selected the most significantly enriched signaling pathways based on their normalized enrichment score. The GSEA shows that cancer pathway, prostate cancer pathway, and wnt signaling pathway are differentially enriched in ECT2 high expression phenotype (Figs. 7B–7F), and parkinson’s disease, ribosome, oxidative, and Huntington disease are differentially enriched in ECT2 low expression phenotype (Fig. 8B). Next, we further validated ECT2 using data from Oncomine, TIMER, and GEO database. The results showed that mRNA levels of ECT2 were significantly upregulated in GC patients compared with normal samples, and high expression of ECT2 predicts poor prognosis in GC patients (Figs. 8A–8C). To further validate ECT2 in GC, RT-qPCR was used to detect the ECT2 mRNA expression in GC, and paired adjacent normal tissue (PANT). Compared with the PANT group, the ECT2 mRNA level was significantly higher in the GC group (Fig. 8D).