Prognostic Significance of Alternative Splicing Genes in Cervical Squamous Cell Carcinoma and Endocervical Adenocarcinoma

Introduction

In the population of developing countries, cervical cancer incidence ranks third among all tumors and first among gynecological tumors.¹ The rising incidence of cervical cancer is related to the following factors: early sexual activities, group sex, infrequent use of condoms, multiple pregnancies with Chlamydia infection, and immunosuppression with HIV.² Effective diagnosis in the early stage can significantly reduce the incidence of cervical cancer.³ It has been shown that timely treatment of dysplasia can effectively reduce the incidence and mortality of cervical cancer.⁴ However, the molecular mechanism of cervical cancer has not been fully understood. If lesions can be detected early on the micro level, the prevention and treatment of cervical cancer will become more efficient.

Alternative splicing (AS) allows selective removal of introns and ligation of exons, two processes during which multiple isoforms of mRNA are produced and post-transcriptional proteome diversified.^5,6 AS is currently categorized into 7 transition types: Alternate Acceptor site (AA); Alternate Donor site (AD); Alternate Promoter (AP); Alternate Terminator (AT); Exon Skip (ES); Retained Intron (RI); and Mutually Exclusive Exons (ME). A growing number of studies have shown that tumorigenesis was associated with the abnormal expression of genes and proteins which were induced by AS.^7–12 Referring to different locations and effect on the usage of a splice site, cis-regulatory sequences could be classified into exonic splicing enhancers, exonic splicing silencers, intronic splicing enhancers and intronic splicing silencers, which determined their affinities to cognate splicing factors (SFs).¹³ Differences in gene expression levels of splicing regulatory factors have been observed in many cancers, and these proteins often affect the splicing patterns of many genes that function in certain cancer-specific biological pathways, including cell cycle progression, cellular proliferation, migration, and RNA processing.^14,15

In the current study, integrated bioinformatics analysis of the AS and gene expression profiles was conducted to identify the splicing events that have prognostic value in CESC. The Univariate Cox proportional hazards regression, LASSO Cox regression and multivariate Cox proportional hazards regression models were performed to screen for the potential prognostic splicing events and genes, uncovering a number of survival‐associated AS events. Regulatory network of between AS events and SFs was constructed to elucidate the underlying mechanisms. However, the molecular mechanism of AS in cervical cancer has not been elucidated. This study aimed to screen AS-generated genes related to CESC and provide new candidate molecules for targeted therapy.

Materials and Methods

Process of AS Data Curation

TCGA data portal (https://portal.gdc.cancer.gov/) provided the RNA sequencing data and clinical information of CESC cohorts. mRNA splicing profiles in CESC was analyzed with the aid of SpliceSeq, a java program that explicitly quantifies RNA-Seq reads and identifies its possible functional changes as a consequence of AS in the context of transcript splice graphs. We downloaded the percent spliced in (PSI) value of seven types of AS events: Exon Skip (ES), Mutually Exclusive Exons (ME), Retained Intron (RI), Alternate Promoter (AP), Alternate Terminator (AT), Alternate Donor site (AD) and Alternate Acceptor site (AA). PSI (range: 0–1) is a common value for scoring AS events.

Survival-Associated AS Events in CESC

A total of 306 CESC tissues and were analyzed in this study and the overall survival (OS) of the 306 CESC patients were at least 90 days. We used UpSetR (version 1.3.3)¹⁶ to present the intersections between seven types of AS. Interactive sets among the seven types of AS events were visualized by UpSet plot. Patients were divided into the low-risk group and the high-risk group based on the median levels of risk score.

Univariate Cox regression analysis was conducted to identify survival-associated AS events (P < 0.05). Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways involving the two lists of genes were determined using an online tool Enrichr database (http://amp.pharm.mssm.edu/Enrichr/).¹⁷ The biological functions of the overlap genes were interpreted. P < 0.05 was considered statistically significant. Since PPI (protein–protein interaction) analysis help identify the hub genes with core functions, we utilized STRING (Search Tool for the Retrieval of Interacting Genes Database) (http://www.string-db.org/) to evaluate PPI of hub modules.¹⁸ To visualize the findings, we exported the results from STRING database to Cytoscape software.¹⁹ The topological analysis method Degree and the centrality analysis methods Closeness and Betweenness were used separately to identify the hub nodes in the PPI network. Among the top 10 hub nodes identified by each method, only genes with consistent high Degree, Closeness, and Betweenness values (larger than the median value) were selected as hub genes. MCODE plug-in was employed for searching clustered sub-networks. Default parameters were set: node score cut-off = 0.2, k-core = 2, degree cut-off = 2, and max. depth = 100. Enrichr was conducted for functional analysis of hub clusters

Prognostic Models Featured by AS Events for CESC

The result of univariate Cox regression analysis in identifying survival-associated AS events was also shown by Volcano plot using “ggplot2” in R. Seven types of AS events were revealed by Bubble chart respectively. Each chart contained the top 20 significant survival-associated AS events of corresponding types. Following univariate Cox, we used the least absolute shrinkage and selection operator (LASSO) Cox analysis. LASSO method is used to calculate the regression of high-dimensional predictors.²⁰ In Cox regression survival analysis, LASSO is ideal for high-dimensional data. Then, we used multivariate Cox regression to construct prognostic feature for each AS type. The prognostic model was constructed based on a linear integration of the expression level multiplied regression model (β) using the following formula: Risk score = β_gene1 × expr_gene1 + β_gene2 × expr_gene2 + ··· + β_genen × expr_genen. Patients were then divided into two groups based on the median levels of risk score. This prognostic model and patient survival information were merged. Kaplan–Meier survival curves were used to compare the prognostic ability of prediction models. Area under the curve (AUC) value for the ROC curves of each prognostic model was calculated by survival ROC package in R. Besides, univariate and multivariate Cox regression analysis were performed to compare the hazard ratios (HRs) of prognostic models and important clinical features for CESC. Finally, Chi-square tests were used to compare the tumor status, age, grade, clinical stage and histological type of the two risk groups. Risk analysis related to tumor-infiltrating immune cells was also performed.

Correlation Network of SFs and Survival-Associated AS Events

The data of SFs was obtained from SpliceAid 2.²¹ The expression profiles of SF genes were curated from the TCGA dataset and further coped with univariate Cox analysis. We selected axes between the expression value of SFs and PSI values of prognosis-related AS events to construct the SF-AS regulatory network according to the following conditions: P < 0.05, Pearson correlation coefficient >0.5. Then, we built the correlation plots via Cytoscape version 3.6.1. Univariate and multivariate Cox regression analyses and LASSO Cox regression analysis were further performed for SF genes.

Results

Association Between AS Events and Cervical Cancer Survival

In total, 41,426 AS events were detected from 19,696 genes of 256 CESC patients and classified into seven AS types (Table 1). The events were mainly enriched in ES, followed by AT and AP. In detail, we detected 15,942 ESs in 6271 genes, 8395 ATs in 3663 genes, 8066 APs in 3256 genes, 3424 AAs in 2397 genes, 3017 ADs in2106 genes, 2373 RIs in 1801 genes, and 209 MEs in 202 genes (Figure 1A). Overall, the results reveal that one gene might possess up an average of 2.5 AS events. To explore the relationship between AS events and OS of CESC patients, all AS events mentioned above were subjected to the univariate Cox regression analysis, and 3708 survival-associated AS events in seven AS types in 3355 genes (Table 2) which we identified as significant AS events (Figure 1B) were screened out.

Table 1 AS Events Were Detected and Classified into Seven as Types

Table 2 3708 Survival-Associated as Events in Seven as Types in 3355 Genes

Figure 1 (A) UpSet plots of AS events in CESC (B) UpSet plots of top significant AS events in CESC. The horizontal axis and vertical axes represent number of genes in the interacting sets of one or multiple AS types and number of genes in each AS type, respectively.

Enrichment Analysis of AS-Related Genes in Cervical Cancer

Since we identified 3708 survival-associated AS events in seven AS types in 3355 genes, to explore possible pathways for these AS events, the parent genes of survival‐associated AS events were subjected to functional enrichment analyses. GO analysis was performed to investigate the biological functions and pathways enriched in these genes as well as the interaction network among them. In the aspect of biological processes, the genes were mostly enriched in Golgi to lysosome transport, negative regulation of mitotic nuclear division and negative regulation of nuclear division (Figure S1A). In the aspect of cellular component, the genes were mostly enriched in lysosome, mitochondrion and lysosomal lumen (Figure S1B). In the aspect of molecular functions, the genes were mostly enriched in death domain binding, nucleoside-diphosphatase activity and GTPase activator activity (Figure S1C). KEGG analysis showed that the genes were mostly enriched in glycosaminoglycan degradation, glyoxylate and dicarboxylate metabolism and NF-kappa B signaling pathway (Figure S1D). Molecular Complex Detection (MCODE) was used to screen out the modules in the protein-to-protein network. PPI network propped by these genes revealed that the AS events were engaged in nine MCODE components (Figure S2A). Among these components, the top related cluster was composed of 12 nodes and 36 edges (Figure S3A), the second of 9 nodes and 17 edges and the third of 4 nodes and 6 edges (Figure S3B and C). The top related cluster was transferred to functional analysis. In the aspect of biological process, the genes were mostly enriched in regulation of oxidative stress-induced neuron intrinsic apoptotic signaling pathway, regulation of osteoclast development and regulation of DNA endoreduplication (Figure S4A). In the aspect of cellular component, the genes were mostly enriched in AP-2 adaptor complex clathrin coat of endocytic vesicle and endolysosome membrane (Figure S4B). In the aspect of molecular function, the genes were mostly enriched in endocytic adaptor activity, clathrin adaptor activity and nitric-oxide synthase binding (Figure S4C). KEGG analysis showed that the genes were mostly enriched in ubiquitin mediated proteolysis, endocytosis and endocrine and other factor-regulated calcium reabsorption (Figure S4D). The top 10 hub genes by PPI network were screened according to their degree values, which were significantly associated with AS events (Figure S2B). Among them, UBA52 connecting 38 nodes showed the most prominence. The distributions of survival-associated AS events are displayed in Figure 2A. Among them, the MEs have no statistical significance. The 20 most significant prognosis-related AS events were also shown (Figure 2B–G).

Figure 2 Top 20 most significant AS events in CESC. (A) The red dots represent AS events that are significantly correlated with patient survival. The blue dots represent AS events without correlation. The top 20 AS events correlated with clinical outcome based on acceptor sites (B) AA (C) AD (D) AP (E) AT (F) ES (G) RI.

Prognostic Models Constructed with Seven Types of AS Events for CESC

To investigate the relationship between each AS event and the prognosis of CESC patients, we constructed proportional hazards regression analysis on genes related to AS events. With a cutoff of p < 0.05, 402 genes with possible prognostic significance were screened out, including 33 of AD, 155 of AP, 103 of ES, 71 of AT, 15 of RI, and 25 of AA (Table S1).

These prognosis-related AS-event genes were further analyzed with the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm (Figure S5). Then, multivariate Cox proportional hazards regression analysis was further performed to build the risk signature. Using the seven types of AS events, we constructed seven prognostic models. The risk scores were calculated (Table 3). Based on the risk score, CESC patients were divided into low-risk group and high-risk group based on the median levels of risk score. In each type of model, Kaplan–Meier survival analysis indicated that low-risk patients had significantly better overall survival than high-risk patients (Figure 3). We also performed the receiver operating characteristic curve (ROC) analysis. As shown in Figure 4, PI-AA demonstrated the highest capacity of estimating the prognosis of CESC patients (AUC 0.92), followed by PI-RI (AUC 0.871). The risk score, survival status and the expression of 7 types of AS-event genes in prognostic model for each patient was displayed in Figures S6 and 7, Figure 5. To assess whether these models were the independent predictor of CESC, univariate and multivariate Cox regression analyses were performed, including clinical factors and risk scores. The results showed that this prognostic model showed moderate and independent prognostic power for all seven AS events (Figures 6 and 7).

Table 3 Seven Prognostic Models Based on 7 Types of as Events

Figure 3 Survival analysis based on seven types of AS-event genes in TCGA cohort. (A) ALL (B) AA (C) AD (D) AP (E) AT (F) ES (G) RI.

Figure 4 ROC analysis based on seven types of AS-event genes. (A) ALL (B) AA (C) AD (D) AP (E) AT (F) ES (G) RI.

Figure 5 Heatmap of the 7 types of AS-event genes expression profiles. (A) ALL (B) AA (C) AD (D) AP (E) AT (F) ES (G) RI.

Figure 6 Univariate Cox regression analysis was performed including clinical factors and risk score based on 7 types of AS-event genes. (A) ALL (B) AA (C) AD (D) AP (E) AT (F) ES (G) RI.

Figure 7 Multivariate Cox regression analysis was performed including clinical factors and risk score based on 7 types of AS-event genes. (A) ALL (B) AA (C) AD (D) AP (E) AT (F) ES (G) RI.

Association Between AS Events’ Risk Scores and Clinical Features

We analyzed the association between AS event risk score and the clinical features of CESC (Figure S8). Interestingly, we found that all the AS genes had significant correlation with tumor status in both the high- and low-risk groups. Among them, AT-related genes were associated with tumor status and grade, RI-related genes with tumor status and stage.

Association Between Splicing Factors and AS Events

SFs play an important part in the occurrence and development of AS events via changing exon selection and splicing site. Therefore, it is necessary to uncover the correlation between SFs and AS events. We got 385 SFs from SpliceAid 2, and downloaded the expression profiles of these 385 genes from the TCGA database. Then, we used the Pearson’s correlation matrices to calculate the relationship between the SFs and AS-events. With a cutoff of p < 0.05 and |Cor| > 0.5, we constructed the SF-AS regulatory network with 43 significant SFs (Figure 8A) (Table S3). The relationship between 43 SFs and CESC prognosis was determined with univariate Cox proportional hazard regression analyses. With a cutoff of p < 0.05, seven genes (SNRPA, PQBP1, CCDC12, SNRPF, CELF2, HSPA8 and TXNL4A) were screened out (Table S2). Then, seven genes were further submitted to LASSO Cox regression algorithm. The optimal values of the penalty parameter lambda were determined through 10-times cross-validations. The simplest (smallest parameter) model of prognostic genes signature within one standard error of the best lambda value was screened out (Figure 8B and C). Finally, multivariate Cox proportional hazards regression analysis screened out three genes (SNRPA, CELF2 and HSPA8) (Figure 8D). Risk score=−0.033* SNRPA-0.34*CELF2+0.003*HSPA8. Based on the median levels of risk score, CESC patients were divided into low-risk and high-risk groups. Kaplan–Meier survival analysis indicated that low-risk group had better overall survival than high-risk group (Figure 9A). In ROC curve analysis, risk score of 0.685 reaped the highest specificity and sensitivity in predicting the 5-year survival (Figure 9B).

Figure 8 (A) Correlation network between significant survival‐associated SFs and significant survival‐associated AS events. The connecting line represents the association between the SF and the AS event. Triangles represent SF, dots represent AS events, red dots and red connecting lines represent AS events positively correlated with SF, green dots and green connecting lines represent AS events negatively correlated with SF. (B and C) Construction of prognostic signatures based on LASSO COX analysis target significant survival‐associated SFs. (D) Multivariate Cox proportional hazards regression analysis based on the hub significant survival‐associated SFs.

Figure 9 (A) Survival analysis of hub significant survival‐associated SFs. (B) ROC analysis of hub significant survival‐associated SFs. (C) Univariate Cox regression analysis was performed including clinical factors and risk score based on hub significant survival‐associated SFs. (D) Multivariate Cox regression analysis was performed including clinical factors and risk score based on hub significant survival‐associated SFs.

Meanwhile, we analyzed the association between the clinical parameters and the risk score of three genes. The univariate and multivariate Cox proportional hazards regression showed that only the risk score calculated based on three genes was an independent prognostic indicator of CESC (Figure 9C and D). Of the three genes, SNRPA and CELF2 were protective genes to CESC. Interestingly, the AS most closely related to CELF2 was TC2N, type AP. Our analysis showed that TC2N was negatively correlated with CELF2, and univariate Cox proportional hazard regression analysis showed that TC2N was a causative factor, and the two results were consistent. Similarly, CALCOCO2, PLS3 and FAM76A developed from ES, NGLY1 from AP, COX15 and PPOX from AT; all these AS-developed genes were positively correlated with SNRPA. The results of univariate Cox proportional hazard regression analysis showed that these AS-developed genes were protective factors. PLEKHB2, PPP1CC and GALNT4 developed from AP, LYRM5 and NPIPA7 from AT, and these genes were negatively correlated with SNRPA. Univariate Cox proportional hazard regression analysis proved that they were all pathogenic factors.

Association Between AS Event and Tumor Immune Microenvironment

In order to explore the role of AS-related genes in the tumor immune microenvironment, We divide the samples into high- and low-risk groups, combined with comprehensive analysis of immune cells. The correlation of risk score and tumor immune microenvironment was shown in (Figure 10A). Figure 10B–G indicated that dendritic cells activated, dendritic cells resting, Macrophages M0, Mast cells resting, T cells CD4 memory activated and T cells CD8 were most correlative with risk core.

Figure 10 (A) The correlation of risk score and tumor infiltrating immune cells. (B) The correlation of risk score and dendritic cells activated. (C) The correlation of risk score and dendritic cells resting. (D) The correlation of risk score and Macrophages M0.(E) The correlation of risk score and Mast cells resting. (F) The correlation of risk score and T cells CD4 memory activated. (G) The correlation of risk score and T cells CD8.

Discussion

AS during mRNA editing diversifies the transcription of proteins. Studies have verified the regulatory role of variable scission in the development of tumors.^12,22,23 In the present study, we shifted our focus to cervical cancer.

In this paper, we firstly investigate the performance of AS events in CESC based on the TCGA data. We divided AS events into seven types, established CESC prognostic models using AS-related genes, and screened those with the strongest correlation with CESC. GO analysis showed that those genes were mostly enriched in Golgi to lysosome transport, lysosome and death domain binding. KEGG analysis showed that those genes were mostly enriched in glycosaminoglycan degradation and NF-kappa B signaling pathway. Perera et al have found that lysosome drives pancreatic cancer metabolism.²⁴ Li et al have found that glycosaminoglycan can mark lung cancer.²⁵ Studies have not yet probed into mechanism of Golgi to lysosome transport and death domain binding in tumorigenesis. The NF-kappa B signaling pathway engages tumorigenesis with inflammation.²⁶ Tang et al found TNF-Alpha promoted the invasion and metastasis via NF-Kappa B pathway in oral squamous cell carcinoma.²⁷ Pacifico et al found that f NF-kappa B pathway also acted in human thyroid carcinomas.²⁸ AS events in our model are mainly active in proliferation, migration, apoptosis and tumor immune-related pathways, which may indicate that abnormal AS mainly affects tumor biological processes through these pathways.

We next analyzed the core modules in the PPI network. GO analysis showed that the core AS-related genes were mostly enriched in regulation of oxidative stress-induced neuron intrinsic apoptotic signaling pathway, AP-2 adaptor complex, and endocytic adaptor activity. KEGG analysis showed that these core genes were mostly enriched in ubiquitin mediated proteolysis, endocytosis and endocrine and other factor-regulated calcium reabsorption. Among these AS-related genes, UBA52 showed the strongest intimacy with CESC. Instructive is that the biological terms we teased out were all related to tumor diseases. For example, ubiquitin-mediated proteolysis plays in the regulation of clear cell renal cell carcinoma.²⁹ Zhou et al have found that UBA52 regulates colorectal cancer tumorigenesis.³⁰ Our results demonstrate a range of genes and pathways pertaining to CESC.

We also found that each type of AS event could specify certain features of CESC patients. These AS events, either combined or single, performed well in predicting overall survival time of CESC patients (AUC > 0.8). Kaplan–Meier survival curves also proved that all the AS-events-based models could efficiently predict the risks and survival outcome. Meanwhile, each prognostic model could serve as an independent prognostic factor. To be more specific, Prognostic models were established with nine genes screened out from the model that incorporated all AS events, including FCF1|28425|AD, C1QTNF1|43985|AP, TBC1D22A|62728|ES, ECE1|963|AP, SLC39A6|45210|ES, GEMIN7|50397|ES, ACSS1|58860|AT, ENTPD6|58863|AA and DYX1C1|30732|AT. The prognostic efficiency of these prognostic models verified that the nine genes could be taken as independent prognostic factors for CESC. Han et al have found that C1QTNF1 affects the progression of human hepatocellular carcinoma.³¹ Xu et al have found that ECE-1 overexpression in worsens tumor differentiation and clinical outcome of head and neck cancer and colon cancer.³² SLC39A6 functions in hepatocellular carcinoma and esophageal carcinoma.^33,34 Yun et al have found that ACSS1 decides the survival of hepatocellular carcinoma.³⁵ Tada et al have demonstrated that ENTPD6 monitors cisplatin resistance in testis and testicular cancer.³⁶ Rosin et al have found that DYX1C1 can indicate the poor survival in breast cancer.³⁷ How FCF1 and TBC1D22A act in CESC is waiting to be answered by future research. However, due to the limited database, we cannot use another cohort to verify the efficiency of the model.

SFs are indispensable for AS events. We constructed an SF-AS regulatory network based on SFs. Cox regression analysis identified that AS was most closely related to CELF2 and TC2N. Ramalingam et al found RNA binding protein CELF2 might be a putative tumor suppressor gene in colon cancer.³⁸ CELF2 can inhibit the development of glioma and non-small cell lung carcinoma.^39,40 TC2N can accelerate lung cancer progression by suppressing p53 signaling pathway.⁴¹ Our research validated these opinions. As our new finding, the role of CELF2 in CESC needs new evidence.

We found that SNRPA interacted with CALCOCO2, PLS3, FAM76A, NGLY1, COX15, PPOX, PLEKHB2, PPP1CC, GALNT4, LYRM5 and NPIPA7. Among them, CALCOCO2, PLS3, FAM76A, NGLY1, COX15, PPOX and PLEKHB2 were all protective factors for their positive association with SNRPA. In contrast, PPP1CC, GALNT4, LYRM5 and NPIPA7 were pathogenic factors for their negative association with SNRPA. SNRPA functions in the progression of gastric cancer.⁴² Kurumi et al have found that PPOX can change the fluorescence intensity in 5-aminolevulinic acid-mediated laser-based photodynamic endoscopic diagnosis for early gastric cancer.⁴³ GALNT4 is implicated in clear cell renal cell carcinoma and colon cancer.^44,45 However, how these AS genes cooperate with SFs to regulate tumor development has not been fully studied, which is a task lying in our future research. In addition, we found that HSPA8 as SF negatively regulates the GEMIN4-38259-AP in CESC. HSPA8 protein integrates compensatory mechanisms to drive cell growth and is dysregulated in a variety of chronic stress diseases, including cancer.⁴⁶ In addition, HSPA8 could be used as candidate biomarkers and potential target for treatment in multiple tumors.^47,48 The results of this study showed that SNRPA, CELF2, HSPA8 affected the prognosis by affecting AS events. The relevant mechanisms and relationship involved in SNRPA, CELF2, HSPA8 have not been found yet, which can be further explore in the further.

Research on prognostic alternative splicing signature by Wang et al reveals the landscape of immune infiltration in pancreatic cancer.⁴⁹ In fact, Immunotherapy is currently widely used in tumor treatment options.⁵⁰ In our study, dendritic cells activated, dendritic cells resting, macrophages M0, Mast cells resting, T cells CD4 memory activated and T cells CD8 were found to be the most correlative with risk score. Dendritic cells are a special type of antigen presenting cells, which play a key role in the process of innate and adaptive immune response. The literature has proved that it is closely related to cancer immunology and immunotherapy.⁵¹ Dendritic cells have also been shown to be closely linked to cervical cancer treatment.^52,53 Macrophages M0 was shown to be associated with the risk of cervical cancer in our previous study.⁵⁴ T cells CD4 and CD8 are the most common immune cells. There are not a few studies about their application in the treatment of cervical cancer.^55–58 All of the above support our analysis. The role of Mast cells has not been fully explored in cervical cancer. This will become a new research direction.

The study also has the following shortcomings. First, no other cohorts can be used to re-validate the prognostic model. In addition, more in vivo and in vitro functional experiments are needed to further explore the impact of dysregulated AS events and SFs on cancer development.

Conclusion

This study screened out the AS events related to CESC and their most enriched functions and pathways. SNRPA and CELF2 were two CESC-engaged SFs cooperating with a list of AS events. This scientific proposition may provide direct guidance for future biological experiments and contribute to the management and individualized treatment of CESC patients.