Systematic profiling of alternative splicing events and splicing factors in left- and right-sided colon cancer

Clinical features of left- and right-sided colon cancer patients

The present study included a total of 434 colon cancer patients, among whom 176 had LC (affecting the splenic flexure of the colon, descending colon and sigmoid colon) and 258 had RC (affecting the ileocecum, ascending colon, hepatic flexure of the colon and transverse colon). The clinical features of the LC and RC patients are listed in Table 1. The proportion of elderly patients was significantly higher in the RC group (50.8% vs. 35.2%, P=0.001). Distant and lymphatic metastases were significantly more common in LC patients (P<0.05). The proportion of patients with microsatellite instability was higher in the RC group than in the LC group (18.9% vs. 2.9%, P=0.064).

DEAS events in left- and right-sided colon cancer

Integrated AS event profiling was performed with data from the 434 colon cancer patients. AS events were quantified based on their percent-spliced-in (PSI) values, which are commonly used for this purpose. We observed extremely low values (PSI <0.05) for certain splicing isoforms. To obtain a reliable a set of AS events in colon cancer, we only included AS events that occurred in at least 75% of the samples with an average PSI value ≥0.05. After the results were thus filtered, a total of 26843 AS events from 8879 genes remained for further analysis. Seven different splicing patterns were identified in colon cancer: 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) (Figure 1A). The ratio of AS events to genes was around 4:1, indicating that each gene underwent four AS events on average. ES was the most frequent splicing pattern, followed by AT and AP, while ME was the least frequent splicing pattern (Figure 1B).

Figure 1. Overview of AS event profiling in left- and right-sided colon cancer. (A) Illustration of the 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). (B) The number of AS events and involved genes for each AS type in the 434 colon cancer patients. (C) The number of (DEAS events upregulated in left- or right-sided colon cancer. (D) UpSet plot of overlapping genes among the seven patterns of DEAS events that were upregulated in right-sided colon cancer. One gene may have up to three splicing patterns. (E) UpSet plot of overlapping genes among the seven patterns of DEAS events that were upregulated in left-sided colon cancer.

To investigate the distinctions between LC and RC at the level of AS, we performed a differential expression analysis. Given the small range of PSI values (from zero to one), we filtered the results based on an adjusted P-value <0.05. Ultimately, 1248 DEAS from 836 genes were identified. Among the DEAS events, 646 AS events from 557 genes were upregulated in RC, while 602 AS events from 550 genes were upregulated in LC (Supplementary Figure 1). The proportion of different types of AS events differed significantly between LC and RC (Figure 1C, P=4.69×10⁻¹²); for instance, the ES pattern was significantly more common in RC, while the RI pattern was significantly more common in LC (P<0.05). Considering one gene had more than one AS pattern, we used UpSet plot to visualize the intersecting sets of different AS pattern. We found that, among the DEAS that were upregulated in RC, genes with the ES pattern occupied the largest number (210 cases), followed by those with the AT pattern and the AP pattern (Figure 1D). However, among the DEAS that were upregulated in LC, genes with the AT pattern occupied the largest number (176 cases), followed by those with the AP pattern and the ES pattern (Figure 1E). All these findings suggested that AS events contribute to the heterogeneity of LC and RC.

Enrichment and interaction analysis of DEAS events

AS can directly alter protein function. To explore the potential functions and pathways of the DEAS events, we performed an enrichment analysis of the differentially spliced genes (DSGs). As shown in Figure 2A, the DEAS that were upregulated in RC were significantly enriched for genes involved in axon guidance (P=0.02), adenosine monophosphate kinase signaling (P=0.009), FoxO signaling (P=0.04), VEGF signaling (P=0.03) and colorectal cancer (P=0.04). The DEAS that were upregulated in LC were significantly enriched for genes involved in B cell receptor signaling (P=0.004), natural killer cell-mediated cytotoxicity (P=0.005), ribosomes (P=0.01) and colorectal cancer (P=0.04) (Figure 2D). Gene ontology (GO) molecular function enrichment analysis indicated that cadherin binding (P=5.4×10⁻³) and cell adhesion molecule binding (P=1.3×10⁻³) were enriched in LC and RC, respectively (Figure 2B and 2E). GO biological process enrichment analysis indicated that coenzyme metabolic processes (P=5.1×10⁻⁶), cofactor metabolic processes (P=7.7×10⁻⁶) and spliceosomal small nuclear ribonucleoprotein assembly (P=1.2×10⁻⁵) were enriched in RC, while nucleobase-containing compound catabolic processes (P=1.3×10⁻⁵) and DNA catabolic processes (P=2.4×10⁻⁵) were enriched in LC (Figure 2C and 2F).

Figure 2. Functional GO analysis and KEGG analysis of DSGs between left- and right-sided colon cancer. The vertical axis represents GO or KEGG pathway annotations. The horizontal axis represents the number of genes assigned to the corresponding annotation. (A–C) right-sided colon cancer; (D–F) left-sided colon cancer. (A and D) KEGG pathways; (B and E) GO molecular functions; (C and F) GO biological processes.

Since AS inevitably alters the translation and features of proteins, we performed a protein-protein interaction network analysis of the proteins encoded by the DSGs in LC and RC. Nodes with over 14 degrees were identified as hub genes with potentially vital regulatory functions in the network. The protein-protein interaction network of DSGs in LC is displayed in Supplementary Figure 2. There were 211 nodes and 405 edges in the network. UBB, RNPS1, RPS29 and PPP2R2A were hub genes in the network. Supplementary Figure 3 depicts the protein-protein interaction network of DSGs in RC. There were 213 nodes and 389 edges in the network. Interestingly, most of the hub genes in RC encoded ribosomal proteins such as RPL37A, RPLP0, RPL24, RPL30 and RPL15, suggesting that the AS of genes encoding ribosomal proteins promotes the development of RC.

Survival-associated DEAS events in left- and right-sided colon cancer

To investigate the relationship between DEAS and overall survival in colon cancer patients, we performed a univariate Cox regression analysis of the 1248 DEAS events in the 434 patients. As shown in Figure 3A, 114 survival-associated DEAS were identified (P<0.05). Six of the seven splicing patterns (all but the ME pattern) contained survival-associated DEAS. The AT pattern contained the most survival-associated DEAS (48 cases), followed by the ES pattern (31 cases). The AD pattern contained the fewest survival-associated DEAS (2 cases). For each splicing pattern, the hazard ratios of the five AS events with the smallest P-values are visualized in Figure 3B–3G.

Figure 3. Forest plots for subgroup analyses of survival-associated DEAS. (A) Volcano plot depicting the P-values from univariate Cox regression analysis of the 1248 DEAS. Log FC: the log-transformed fold-change in the PSI value of a DEAS. (B–G) Forest plots of hazard ratios for the five AS events with the smallest P-values in the AA, AD, AP, AT, ES and RI splicing patterns, respectively. P-values are indicated by the color scale on the side. Horizontal bars represent the 95% confidence intervals.

Next, we sought to identify independent prognostic DEAS in colon cancer. Since the univariate Cox regression was only a preliminary screening, we used a relatively loose filter (P-value <0.15) to select variables for multivariate Cox regression analysis. We performed separate multivariate Cox regression analyses for the six splicing patterns. The multivariate Cox regression analysis results for each pattern of AS events are shown in Figure 4A–4F. Colon cancer patients were divided into high-risk and low-risk groups according to the median risk scores predicted by the prognostic models. Five of the six prognostic models exhibited significant power to distinguish good from poor outcomes in colon cancer patients. The prognostic model based on the ES pattern was the most powerful, with a P-value <0.0001. To further assess the discriminatory abilities of these prognostic models, we generated receiver operating characteristic (ROC) curves and calculated the area under the curve (AUC) (Figure 4G). The prognostic model based on the ES pattern displayed the greatest discriminatory ability, with an AUC of 0.755.

Figure 4. Kaplan-Meier plots and ROC curves of prognostic models for different AS patterns. (A–F) Kaplan-Meier curves of prognostic models built with the AA, AD, AP, AT, ES and RI patterns of AS, respectively. The red line indicates the high-risk group, whereas the blue line indicates the low-risk group. (G) The ROC curves of the predictive models for the different AS patterns.

To obtain the final prognostic model, we selected independent prognostic DEAS events from the multivariate Cox regression analysis of each splicing pattern, and further assessed them by multivariate Cox regression analysis. Ten independent prognostic DEAS were selected, and their hazard ratios and P-values are summarized in Figure 5A and 5B. Colon cancer patients were divided into high-risk and low-risk groups according to the median risk score predicted by the final prognostic model. Survival analysis demonstrated that the final prognostic model had significant power to distinguish good from poor outcomes in colon cancer patients (P<0.001) (Figure 5C and 5G). Subgroup analysis indicated that the final prognostic model could efficiently distinguish good from poor outcomes in patients with either LC or RC (RC: Figure 5D and 5H; LC: Figure 5E and 5I). ROC curve analysis revealed that the final prognostic model was more efficient than any of the individual splicing-pattern-based prognostic models in distinguishing good from poor outcomes in colon cancer patients (AUC for the final prognostic model: 0.84). The final prognostic model exhibited greater accuracy in RC than in LC (AUC = 0.90 for RC; Figure 5F). Figure 6 displays the expression of the 10 independent prognostic DEAS in LC and RC. Detailed information on the 10 AS events in the prognostic model is provided in Supplementary Table 1.

Figure 5. The prognostic value of the DEAS signature. (A) Univariate analysis of the 10 predictive factors for overall survival. P-values are indicated by the color scale on the side. Horizontal bars represent 95% confidence intervals. (B) Multivariate analysis of the 10 predictive factors for overall survival. P-values are indicated by the color scale on the side. Horizontal bars represent 95% confidence intervals. (C) Kaplan-Meier curves of the final prognostic model. Patients were divided into the high-risk and low-risk groups according to the median risk score. (D) Kaplan-Meier curves of the final prognostic model in right-sided colon cancer patients. (E) Kaplan-Meier curves of the final prognostic model in left-sided colon cancer patients. (F) The ROC curves of the final prognostic model in all, right-sided and left-sided colon cancer patients. (G–I) Construction and analysis of risk scores. The top panels indicate the risk scores of the patients. The middle panels depict the survival statuses and survival times of the patients distributed by risk score. The bottom panels display the heatmap of the PSI values for the 10 predictive factors distributed by risk score. (G) all patients; (H) right-sided colon cancer patients; (I) left-sided colon cancer patients.

Figure 6. The expression of the 10 independent prognostic DEAS events in left- and right-sided colon cancer.

Survival-associated DEAS events expression in colon cancer and protein structure prediction

To explore survival-associated DEAS events expression in colon cancer, we next used reverse-transcription quantitative PCR (RT-qPCR) to assess the expression of three independent prognostic AS events in clinical specimens. Among the three selected AS events, two (FIP1L1-ES and SATB2-AP) were members of the final prognostic model, while SMAGP-AP were independent prognostic AS events from AP pattern models. We designed two pairs of primers to quantify each AS event. One pair of primers specifically amplified the included fragment in the isoform of interest, and was used to quantify the expression of a specific AS product. The other pair of primers amplified the common fragment among the different isoforms, and was used to quantify the total expression of the various isoforms. The expression of SMAGP-AP did not differ significantly between cancer samples and adjacent tissues (Supplementary Figure 4). On the other hand, as shown in Figure 7A and 7B, FIP1L1-ES and SATB2-AP were significantly downregulated in colon cancer samples compared with adjacent tissues (P<0.05). Notably, the expression of the common fragment did not differ significantly between the cancer samples and adjacent tissues, indicating that the differences in AS events between cancer samples and adjacent tissues were not caused by expression changes at the whole-gene level.

Figure 7. Survival-associated DEAS events expression in colon cancer and protein structure prediction. (A) The schematic diagram (top panel) depicts the ES of FIP1L1, where exon sequences are denoted by boxes and intron sequences are denoted by a horizontal line. The excluded exon is marked with a white stripe. The pair of red arrows indicates the primers amplifying the excluded exon, while the pair of cyan arrows indicates the primers amplifying the common exon among the different isoforms. The left panel displays the expression of FIP1L1-ES in cancer and adjacent tissues. The right panel displays the expression of FIP1L1 in cancer and adjacent tissues. (B) The schematic diagram (top panel) depicts the AP of SATB2, where exon sequences are denoted by boxes and intron sequences are denoted by a horizontal line. The excluded exon is marked with a white stripe. The pair of red arrows indicates the primers amplifying the excluded exon, while the pair of cyan arrows indicates the primers amplifying the common exon among the different isoforms. The left panel displays the expression of SATB2-AP in cancer and adjacent tissues. The right panel displays the expression of SATB2 in cancer and adjacent tissues. (C) Predicted structures of FIP1L1. The shorter variant (left panel) and longer variant (middle panel) were predicted by I-TASSER. The black circle in the right panel indicates the structure that could not be aligned.

Different isoforms of mRNA can be translated into proteins with different structures, thus increasing the diversity of the proteome and the complexity of cells. Therefore, we explored the effects of AS on protein structures. The AP of SATB2 is located in the 5′ untranslated region, where it can influence the translation process, but not the protein structure. Thus, we only used the ES of FIP1L1 as an example. The ES of FIP1L1 occurred in the 11^th exon and resulted in two isoforms of FIP1L1. The longer variant could be translated into a sequence of 588 amino acids, while the shorter variant could be translated into a sequence of 552 amino acids. The predicated three-dimensional structures of the two variants are shown in Figure 7C. The shorter variant lacked a structure that could have altered the protein functional domain and protein-protein interactions.

Differentially expressed SFs in left- and right-sided colon cancer

SFs are vital regulators of AS events. SFs bind to pre-mRNAs at specific positions and subsequently process them into mature RNAs. A single SF can induce numerous AS events, and it has been suggested that a limited number of SFs orchestrate the dysregulated AS events in the tumor microenvironment [32]. Given the marked differences in AS events between LC and RC (as illustrated above), we considered it important to explore the relationship between AS events and SFs in LC and RC.

First, we identified differentially expressed SFs (DESFs) in LC and RC. The mRNA levels of 71 experimentally validated SFs were obtained and used for differential expression analysis. We identified 10 DESFs in LC and RC. As shown in Figure 8, five SFs were upregulated in RC and five SFs were upregulated in LC. ELAVL2 exhibited the most significant difference in expression between RC and LC.

Figure 8. The expression of the 10 DESFs in left- and right-sided colon cancer.

Correlation network of DESFs and DEAS evens

Next, we analyzed the correlations between the DESFs and DEAS events, and constructed a correlation network from the significantly related pairs (|R| > 0.4 and adjusted P-value <0.05). As shown in Figure 9A, 218 DEAS events were significantly associated with 7 DESFs. Among the 218 DEAS events, 119 were upregulated in RC (orange dots) and 99 were upregulated in LC (green dots). RBM25 was a hub SF in the correlation network, and was significantly associated with 121 of the 218 DEAS, indicating that it was a key determinant of the distinct AS events in LC and RC. Representative correlations between DESFs and DEAS events are presented as dot plots (Figure 9B). For example, SRSF6 expression correlated positively with RBM39 expression, suggesting that RBM39 could be a potential target of SRSF6.

Figure 9. Correlation network of DESFs and DEAS events. (A) Correlation network. The correlations between the mRNA levels of the 10 DESFs and the PSI values of each DEAS events were analyzed, and a splicing regulatory network was built among the significant correlations. Quadrate nodes indicate SFs that were upregulated in right-sided (red nodes) or left-sided colon cancer (green nodes). Circular nodes indicate AS that were upregulated in right-sided (red nodes) or left-sided colon cancer (green nodes). Orange lines indicate positive correlations, while deep pink lines indicate negative correlations. (B) Representative dot plots of correlations between the mRNA levels of DESFs and the PSI values of DEAS (P<0.05).

We then analyzed the clinical significance of RBM25 in colon cancer. RBM25 expression was significantly greater in colon cancer tissues than in adjacent tissues (P=2.4×10^-8). To assess the diagnostic value of RBM25 in colon cancer, we used an ROC curve to analyze its sensitivity and specificity. The AUC of RBM25 was 0.672, and its sensitivity and specificity values at a cut-off value of 6.28 variance stabilizing transformation (VST) value were 0.66 and 0.71, respectively. However, a survival analysis indicated that the overall survival did not differ significantly between patients with high and low RBM25 expression (Supplementary Figure 5).

Clinical specimens

In total, 14 colon cancer patients who underwent colectomies at the Guangxi Medical University Cancer between June and July of 2019 were included in this study. Colon cancer and adjacent tissue specimens were collected. All patients had primary colon cancer and had not undergone chemotherapy or radiotherapy before the collection of their tissues. The patients included 10 men and 4 women with a mean age of 55.5 years (range: 34–78 years). Detailed clinical information on the 14 patients is shown in Supplementary Table 2.

Written informed consents were obtained from all patients. The study was approved by the Ethics and Human Subject Committee of Guangxi Medical University Cancer Hospital. All experiments and methods were performed according to relevant guidelines and regulations.

Data acquisition

We downloaded AS data on colon cancer from The Cancer Genome Atlas (TCGA) SpliceSeq, a web-based resource for exploring the AS patterns of 33 different tumor types [50]. PSI values, which range from zero to one, were used to quantify AS events. Given that the PSI values of many AS events were relatively small, we filtered the results (based on ≥75% of samples having an AS event, with an average PSI value ≥0.05) to generate a set of AS events [31]. The level 3 RNA-Seq data and corresponding clinical information from colon cancer patients were downloaded from the Genomic Data Commons data portal (https://portal.gdc.cancer.gov/). The barcodes from TCGA were used to match the AS data, RNA-Seq data and clinical data with each other. Patients who met the following criteria were included in the study: 1. Patients with complete clinical parameters, including sex, age, and information on the cancer location, local invasion, lymph node metastasis, distal metastasis, pathologic stage and survival, and 2. Patients with corresponding RNA-Seq data and AS data. The list of 71 SFs was obtained from SpliceAid-F, a database of experimentally validated SFs [51].

Identification of DEAS and enrichment analysis

A t-test was performed to identify DEAS events between LC and RC, and P-values were adjusted by the Benjamini-Hochberg method. Given that PSI values are small, we used an adjusted P-value of <0.05 to identify statistically significant DEAS events. The parent genes of these DEAS events were then subjected to GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses in clusterProfiler [52]. Terms with P-values <0.05 were selected for further analysis. The interactions between the parent genes of these DEAS events were downloaded from the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) 9.1 database [53]. A required interaction score of 0.9 was used for the protein interaction networks generated in STRING, and the default parameters were used for other settings. Cytoscape v3.4.0 was used to depict the gene interaction network [54].

Survival analysis

Survival was initially assessed in a univariate Cox regression analysis based on the PSI value of each DEAS product. DEAS events with P-values <0.15 in the univariate Cox regression analysis were entered into the multivariate Cox regression analysis. We first performed multivariate Cox regression analyses based on different AS patterns and constructed corresponding predictive models. Independent DEAS events from the multivariate Cox regression analyses of the different AS patterns were entered into the final multivariate Cox regression analysis. Then, independent DEAS events from the final multivariate Cox regression analysis were used to construct the final prognostic model. Kaplan-Meier curves were used to determine whether the prognostic models could distinguish good from poor outcomes in colon cancer patients. The discriminatory ability of each prognostic model at five years (1825 days) was further assessed by ROC curve analysis in the survival ROC package.

Identification of DESFs

The count values of the SFs were derived from RNA-seq data. DESFs in LC and RC were identified through the DESeq2 package [55], and the normalized mRNA levels were calculated with the variance-stabilizing transformation function of this program. P-values were adjusted by the Benjamini-Hochberg method. The threshold for DESFs was set at an adjusted P-value <0.05.

Correlation network construction

The correlations between the normalized mRNA levels of DESFs and the PSI values of DEAS events were calculated with the cor.test function in R. P-values were adjusted by the Benjamini-Hochberg method. Adjusted P-values <0.05 and correlation coefficients with absolute values >0.4 were considered significant. The correlation plots were generated with Cytoscape. The topology structure of the network was analyzed with the NetworkAnalyzer in Cytoscape. A node with over 14 degrees was regarded as a hub node in the network.

RT-qPCR validation of AS events

RT-qPCR was performed to validate the differential expression of selected AS events. Total RNA was extracted with Trizol reagent (Invitrogen, USA) according to the manufacturer’s instructions. Then, M-MLV Reverse Transcriptase (Promega, USA) was used to reverse-transcribe 2-6 μg of the total RNA into cDNA. RT-qPCR was performed on a qTOWER3 G Real-Time PCR system (Analytik Jena, Germany) in a 20-μL reaction mixture including 0.1 μM primers, 10 μL of GoTaq® qPCR Master Mix (Promega, USA) and 20-100 ng of the cDNA sample. The PCR conditions included denaturing at 95 °C for 10 min, followed by 40 cycles of denaturing at 95 °C for 15 s and annealing and extension at 60 °C for 1 min.

We quantified the expression of specific AS events in a method similar to PSI value calculation, which was the percentage of include exon. We designed two pairs of primers for each gene. One pair specifically amplified the included fragment, and was used to quantify the expression of a specific AS product. The other pair amplified the common exon among different isoforms, and was used to quantify the total expression of the various isoforms. The primers used in the current study are listed in Supplementary Table 3. The relative expression of each gene was calculated by the 2^-ΔΔCT method and normalized to that of the reference gene GAPDH [56].

Three-dimensional structural modeling

We employed I-TASSER to predict the three-dimensional structures of different protein isoforms. I-TASSER is a fully automated three-dimensional structure prediction tool that employs a hierarchical approach [57]. The amino acid sequences of different isoforms were entered into I-TASSER, and the three-dimensional structures of the proteins were predicted with the default parameters. PyMol was used for structure visualization, and the ‘align’ function of PyMol was used to compare protein structures.