Myc pathway-guided alternative splicing events predict the overall survival of lung squamous cell carcinoma

Abstract

Oncogenic signaling pathways-guided alternative splicing events and their clinical prognostic values were investigated in lung squamous cell carcinoma (LUSC). RNAseq and AS data of LUSC were obtained from TCGA and SpliceSeq database, respectively. Pathway enrichment-guided activity study of alternative splicing (PEGASAS) analysis was performed to evaluate the activity of oncogenic pathways. We identified 4229 AS events in 2311 genes significantly associated with the Hallmark_Myc_Targets_V1 pathway. Twenty-nine prognostic Myc guided-AS events classified LUSC into 3 subtypes, of which Cluster 1 subtype showed the worst prognosis, while cluster 2 subtype showed a favorable prognosis. A prognostic model based on 10 most contributive prognostic AS events identified by LASSO could effectively stratify LUSC patients into high- and low-risk groups. The high-risk group was associated with worse survival and contained more samples in poor prognostic subtype (cluster 1). GSEA revealed that cell adhesion and immune response pathways were significantly enriched in the high-risk group, and AS pathways in the low-risk group. Fifty-nine regulatory SFs were identified to regulate these 10 AS events. Myc-guided AS events could classify LUSC samples into subtypes with significantly different survival. A prognostic model based on 10 Myc-guided AS events could predict the survival of LUSC.

KEYWORDS:

• Lung squamous cell carcinoma
• Myc signaling pathway
• alternative splicing
• prognosis

Introduction

There were 2,206,771 new cases and 1,796,144 deaths globally in 2020 for lung cancer, accounting for 11.4% of the global cancer burden and 18.0 of global cancer-related mortality, respectively (Sung et al. 2021). Non-small cell lung cancer (NSCLC) accounts for approximately 85% of all lung cancers, and lung squamous cell carcinoma (LUSC) is a major subtype of NSCLC (Sabbula and Anjum 2019). It’s the most common type of NSCLC in men, which may be due to its high correlation with smoking. Although around 70% of stage I patients survived more than 5 years, the overall five-year survival rate of LUSC is around 20%, due to the last diagnosis of LUSC at advanced stages (Siegel et al. 2019; Sabbula and Anjum 2019). Traditional chemotherapy is the primary treatment for advanced LUSC, but the response to treatment in individuals varies greatly. In recent years, predictive molecular biomarkers have been considered a powerful addition to the available treatments for LUSC, which can identify the patients who benefit the most and predict the risk of adverse events (Genova et al. 2016; Liao et al. 2020). Therefore, it is necessary to identify novel prognostic biomarkers for LUSC.

Alternative splicing (AS) is a conserved biological process that generates multiple RNA variants from a single transcript (Kornblihtt et al. 2013). AS causes variable incorporation of coding regions, resulting in the production of functionally different proteins and contributing to proteomic diversity, while AS of noncoding regions can modulate the localization, stability and translation efficiency of mRNAs (Kim et al. 2008). Therefore, keeping AS under control is very important, but abnormality of AS is prevalently found in diseases and is associated with disease progression (Garcia-Blanco et al. 2004; Scotti and Swanson 2016). Dysregulation of AS is also identified in cancers, including LUSC (Garcia-Blanco et al. 2004; Bonnal et al. 2020; Wang and Aifantis 2020). AS is emerging as a new hallmark of cancers (Ladomery 2013; Farina et al. 2020), and contributes to tumor development, drug resistance, immune evasion and metastasis (Marzese et al. 2018; Sciarrillo et al. 2020). Lack of the RAS-binding domain of BRAF V600E resulting from alternative splicing in melanoma leads to resistance to RAF inhibitors (Poulikakos et al. 2011). Exon 4 deletion variant of epidermal growth factor receptor associates with metastasis of ovarian cancer (Zhang et al. 2013), while DeltaRon, a constitutively active isoform of tyrosine kinase RON, promotes invasion and motility of cancer cells (Ghigna et al. 2005).

Many studies have reported the association of AS with the prognosis of multiple cancers (Hu et al. 2019; Chen et al. 2020; Wan et al. 2020; Wu et al. 2021). Epithelial splicing regulatory proteins 1 and 2 (ESPR1/2) are the master regulators for AS, and studies have shown that ESPR1/2 expression promotes AS of CD44 and associates with poor outcomes in breast and ovarian cancer (Yae et al. 2012; Chen et al. 2017). Δ40p53 isoform is elevated in tumor tissues and positively correlated with unfavorable prognosis in triple-negative breast cancer (Avery-Kiejda et al. 2014). AS in LUSC is also related to cancer prognosis (Liu et al. 2020; Yan et al. 2021). Yan et al. detected 1996 AS events in LUSC that correlated with overall survival of LUSC (Yan et al. 2021), while another group, 1991 prognostic AS events (Liu et al. 2020). However, these studies focused mainly on the effect of individual AS events, but did not focus on a set of AS events guided by oncogenic pathways.

Here we aim to investigate the global AS events in LUSC and identify a set of AS events associated with important cancer-related signaling pathways. Its predictive capacity on LUSC prognosis was examined thereafter through the workflow shown in Figure S1. Our study identified oncogenic Myc-related AS events that exhibited a potent prognostic ability for LUSC patients.

Materials and methods

Data collection and processing

The RNAseq data and clinical information of LUSC from The Cancer Genome Atlas (TCGA) database were downloaded. The percent spliced index (PSI) values (range 0-1, an intuitive ratio to quantify different AS events) of AS in TCGA_LUSC samples were downloaded from the SpliceSeq database (Ryan et al. 2012). The data were filtered with the following criterion: ≥75% samples have PSI value and the average PSI value ≥ 0.05.

A total of 540 samples (48 healthy samples and 492 tumor samples with tumor stage information) that had RNA-seq data and PSI value were included in our study. The tumor samples contain 401 early-stage samples and 91 advance-stage samples (Supplementary Table 1). Among them, 485 samples with prognostic information were used for the construction of a prognostic model.

Pathway enrichment-guided activity study of alternative splicing (PEGASAS) analysis

H: hallmark gene sets containing 51 tumor-related pathways were downloaded from the MSigDB database (V7.1) (Liberzon et al. 2011). The activity score of each pathway was calculated by PEGASAS (Phillips et al. 2020) based on the expression profile. Activity comparison between normal samples, early-stage, and advance-stage groups were performed using student’s t-test. The correlation analysis between pathway activity score and PSI value of AS was performed with a cutoff of |correlation coefficient| > 0.2 and p < 0.0002.

KEGG pathway enrichment and gene ontology analysis

The GO biological process (BP) and KEGG pathway enrichment analysis were performed by the online tool DAVID with a cutoff of p < 0.01 and count >2. The redundant GO terms were filtered out by REViGO and displayed in scatterplots.

Screening of prognostic-related AS events

The AS events showed a significant difference between normal, early-stage, and advance-stage samples were screened using t-test, and their prognostic values were evaluated using univariate Cox regression analysis. AS events with p < 0.05 were significantly associated with prognosis.

Unsupervised clustering of AS events

Unsupervised clustering of prognostic-related AS events was performed by heatmap (Version 1.0.8) in R3.6.1 based on centered person correlation to reveal subtypes of tumors. The correlation of each subtype with prognosis was analyzed by Kaplan–Meier (KM) survival curve using survival package (Version 2.41-1) in R3.6.1. The clinical information between different subtypes was compared thereafter.

Construction of the prognostic risk model

Tumor samples were divided into a training set and a validation set randomly and equally. The training set, validation set, and whole set were used for the prognostic model construction and validation. The most contributive prognostic AS events were further screened by the LASSO Cox regression analysis using the glmnet package (Version 4.0-2) in R3.6.1, in which AS events were analyzed with 10-fold cross-validation to obtain an optimized λ value. AS events with a regression coefficient value equal to 0 were the optimized AS events used to construct a prognostic model. The prognostic score (PS) was calculated by the following formula: $\text{Prognostic score }(\text{PS})=\sum \text{Coe}{\text{f}}_{\text{AS}}\times \text{Ex}{\text{p}}_{\text{PSI}}$Coef_AS represents the LASSO prognostic coefficient of AS event. Exp _PSI represents the expression level of AS event in the training set.

Evaluation and comparison of the predictive efficacy of the prognostic model

On the basis of the formula, the PS score for all samples in each data set was calculated, followed by subgrouping samples into High-risk (PS > Median_(PS)) and Low-risk (PS < Median_(PS)) subgroups based on the median value. The KM survival curves were plotted by survival package (Version2.41-1) in R.3.6.1 to evaluate the differences in survival between of High-risk and Low-risk groups.

Construction of regulatory network between splicing factor (SF) and prognostic AS events

Seventy-one SFs (corresponding to 67 genes), validated in previous studies, were downloaded from the SpliceAid2 database (Piva et al. 2012). Their expression levels were obtained from the LUSC transcriptomic data. The Spearman correlation coefficient (ρ) between SF expression and prognostic AS events was calculated, followed by Benjamini-Hochberg corrections. The adjusted p < 0.05 and |ρ| > 0.2 were set as cutoff values to screen SF-AS regulatory pairs. The SF-AS regulatory network was constructed by cytoscape (Version 3.4.0) (Shannon et al. 2003).

Gene set enrichment analysis (GSEA) analysis of the high-risk and low-risk group

The differential analysis between the high-risk and low-risk group was performed by limma (Version 3.34.7) based on the expression data to calculate logFC value. The GSEA pathway enrichment analysis was performed by the clusterProfiler (Version 3.16.0) based on reference KEGG pathway gene sets from the MSigDB database. All genes were input and ranked by LogFC descending, followed by Benjamini-Hochberg corrections. The adjusted p.adjust < 0.05 was set as a cutoff value to present significantly enriched pathways.

Validation of key genes in an independent cohort

The expression pattern of the genes involved in key prognostic AS events was analyzed between LUSC tumor and normal samples using a T-test based on the TCGA cohort and two other independent datasets. GSE19188 dataset contained 27 LUSC and 65 normal samples, and GSE30219 contained 61 LUSC and 14 normal samples. Both these datasets were downloaded from the Gene Expression Omnibus database.

Results

Identification of Hallmark gene sets associated AS events

The pathway activity scores of 51 hallmark gene sets in normal, early-stage, and advance-stage were calculated by the PEGASAS analysis (Table S2). The comparison of pathway activity identified 2 differential pathways (Hallmark_Apical_Surface and Hallmark_MYC_Targets_V1) that showed a significant difference between these groups (Figure 1(A,B)). A total of 1868 AS events (Corresponding to 1126 genes) and 4229 AS events (Corresponding to 2311 genes) were significantly associated with the Hallmark_Apical_Surface and Hallmark_MYC_Targets_V1, respectively (Table S3).

Figure 1. Identification of Hallmark gene sets-associated AS events. (A, B) The activities of tumor signaling pathways in normal, early-stage and advance-stage of LUSC were calculated by the PEGASAS analysis. Hallmark_Myc_Targets_V1 (A) and Hallmakr_Apical_Surface (B) were shown. (C, D) GO BP analysis of Hallmark_Myc_Targets_V1 (C) and Hallmakr_Apical_Surface (D)-associated AS events. (E) KEGG pathway enrichment of Hallmakr_Apical_Surface (left) and Hallmark_Myc_Targets_V1 (right)-associated AS events.

The GO BP and KEGG pathway enrichment analyses were performed on the Hallmark_Apical_Surface and Hallmark_MYC_Targets_V1 associated AS genes. Twenty-six GO BP and nine KEGG pathways were obtained from the Hallmark_Apical_Surface-associated AS genes, while 41 GO BP and 10 KEGG pathways were enriched from the Hallmark_MYC_Targets_V1-associated AS genes (Figure 1(C–E)). Notably, RNA splicing-related pathways, including GO BP terms regulation of RNA splicing, mRNA splicing via spliceosome and KEGG pathway spliceosome, were observed in the Hallmark_MYC_Targets_V1-associated AS genes, suggesting this hallmark pathway was closely related to RNA splicing.

Screening of prognosis-related AS events and AS-based LUSC subtype classification

A total of 552 AS events guided by the Hallmark_MYC_Targets_V1 showed significant differences among normal, early-stage, and advance-stage samples (Figure 2(A)). Of which, univariate Cox regression analysis revealed 29 prognosis-related AS events (Table 1).

Figure 2. Screening of prognosis-related AS events and AS-based LUSC subtype classification. (A) Heatmap of PSI value of differential AS events between normal, early-stage and advance-stage LUSC. (B) Unsupervised heatmap of PSI value of prognostic AS events. (C) The Kaplan-Meier survival curve of LUSC clusters. (D–G) distribution of patient tumor stage (D), age (E), gender (F) and pathologic_N (G) in LUSC clusters.

Table 1. Univariate Cox regression analysis of prognosis-related AS events.

AS events	HR (95% CI)	P value
ID_9625	1.638 (1.242–2.161)	0.000
ID_81187	1.503 (1.138–1.984)	0.004
ID_9624	0.670 (0.509–0.883)	0.004
ID_81185	0.668 (0.506–0.882)	0.004
ID_24704	0.672 (0.510–0.886)	0.005
ID_24705	1.478 (1.121–1.948)	0.006
ID_37798	1.469 (1.115–1.936)	0.006
ID_69089	1.465 (1.112–1.930)	0.007
ID_59424	1.460 (1.106–1.927)	0.008
ID_59422	0.698 (0.529–0.921)	0.011
ID_38130	0.710 (0.540–0.934)	0.014
ID_49932	0.710 (0.539–0.935)	0.015
ID_66486	1.408 (1.069–1.856)	0.015
ID_38131	1.400 (1.064–1.841)	0.016
ID_18268	1.397 (1.058–1.845)	0.018
ID_79157	0.732 (0.557–0.963)	0.026
ID_28797	0.739 (0.561–0.973)	0.031
ID_59273	0.739 (0.561–0.973)	0.031
ID_49685	1.350 (1.026–1.776)	0.032
ID_63079	0.697 (0.500–0.972)	0.033
ID_63415	1.349 (1.023–1.777)	0.034
ID_79158	1.340 (1.019–1.762)	0.036
ID_62294	0.745 (0.565–0.983)	0.037
ID_66487	0.747 (0.567–0.984)	0.038
ID_13995	0.751 (0.570–0.988)	0.041
ID_29320	1.330 (1.011–1.751)	0.042
ID_68258	0.752 (0.572–0.990)	0.042
ID_13994	1.323 (1.005–1.741)	0.046
ID_68257	1.322 (1.005–1.740)	0.046

HR: hazard ratio; CI: confidence interval.

Based on the expression of the 29 prognostic AS events, unsupervised clustering-assigned tumor samples into three clusters (Figure 2(B)). There were 157, 188, and 140 samples in clusters 1, 2, and 3, respectively. The KM survival curve showed a different prognosis of each cluster. Cluster 1 subtype showed the worst prognosis, while cluster 2 subtype showed a relatively good prognosis (Figure 2(C)). Combined with clinical information, we found significant differences in tumor stage, age, gender, and pathologic-N between clusters (Figure 2(D–G)). More early-stage samples and pathologic N0-N1 samples were clustered in the Cluster 3 subtype. Therefore, early-stage and low pathologic grades might be drivers for the better survival of the Cluster 3 subtype.

Construction and evaluation of the predictive efficacy of the prognostic model

Among the 29 prognosis-related AS events, LASSO regression analysis identified 10 most contributive prognostic AS events (Figure 3(A) and Table 2), including AT events (ID_81185 and ID_81187) in DNAJC2 gene, AT event (ID_28797) in EFCAB11 gene, AT event (ID_29320) in PPP2R5C gene, AT event (ID_9625) in C4BPB gene, AP event (ID_38130) in VPS9D1 gene, AP event (ID_63415) in PPARG gene, AP event (ID_66486) in CCDC14 gene, AP event (ID_79157) in AQP1 gene and ES event (ID_59273) in EPB41L1 gene. KM survival curves were plotted to evaluate the survival differences for each AS event between the high PSI (PSI value ≥ median) and the low PSI group (PSI value < median), every AS event showed significant differences between the high PSI and low PSI groups (Figure 3(B–K)). Among them, high levels of several AS events (ID_9625, ID_63415, ID_66486, ID_81187, and ID_29320) were correlated with a poor prognosis of LUSC, while others (ID_28797, ID_38130, ID_59273, ID_79157, ID_81185) were correlated with favorable prognosis.

Figure 3. Construction and evaluation of the predictive efficacy of the prognostic model. (A) λ selection in LASSO model. (B–K) The Kaplan-Meier survival curve of OS-related AS events.

Table 2. Optimized prognostic AS events.

AS ID	Gene symbol	AS type	LASSO coefficient	HR
ID 28797	EFCAB11	AT	−0.122	0.739
ID 29320	PPP2R5C	AT	0.506	1.33
ID 38130	VPS9D1	AP	−0.785	0.71
ID 59273	EPB41L1	ES	−0.542	0.739
ID 63415	PPARG	AP	0.096	1.349
ID 66486	CCDC14	AP	0.916	1.408
ID 79157	AQP1	AP	−0.392	0.732
ID 81185	DNAJC2	AT	−2.373	0.668
ID 81187	DNAJC2	AT	0.004	1.503
ID 9625	C4BPB	AT	0.129	1.638

The prognostic model was established based on these 10 AS events, and the PS scores of all samples in the training set, validation set, and whole set were calculated. According to the median PS score, samples were obviously divided into high-risk and low-risk groups. The distribution of PS scores and PSI values in the training set was shown (Figure 4(A)). The PSI values of AS events (ID_66486, ID_63415, ID_9625, ID_81187, ID_29320) were gradually elevated along with increase of PS score, while others (ID_79157, ID_81185, ID_59273, ID_28797, ID_38130) were decreased. Samples in each set were divided into high-risk group and low-risk group based on the median of PS scores, and KM survival curves were plotted. Patients in the high-risk group showed worse survival than those in the low-risk group (Figure 4(D)). These observations were further confirmed in the validation set (Figure 4(B,E)) and whole set (Figure 4(C,F)), demonstrating our prognostic model can effectively predict the prognosis of LUSC.

Figure 4. Evaluation of the predictive efficacy of the prognostic model. (A–C) PS scores distribution (Up), patient survival information (Middle), and heatmap for PSI of 10 prognostic AS events (Low) in the training set (A), validation set (B), and whole set (C). (D–F) The Kaplan-Meier survival curve of high- and low-risk groups in the training set (D), validation set (E), and whole set (F). (G) distribution of LUSC clusters in high- and low-risk groups.

To further interrogate the associations between risk groups and LUSC subtypes, we examined the samples’ distribution of each cluster in high-risk and low-risk groups (Figure 4(G)). The high-risk group contains more cluster 1 samples (46%), while the low-risk group was dominated by cluster 2 patients (53%). This result was consistent with the survival analysis that cluster 1 subtype had worse survival, while cluster 2 subtype had favorable survival than other subtypes. Altogether, these data demonstrated that the higher proportion of cluster 1 subtype patients might explain the worse survival of the high-risk group.

Expression pattern of genes involved in the prognostic AS events

Nine genes were involved in the 10 most prognostic AS events in the prognostic model, including DNAJC2, EFCAB11, PPP2R5C, C4BPB, VPS9D1, PPARG, CCDC14, AQP1 and EPB41L1. The expression pattern of these 9 genes was analyzed. In the TCGA cohort, all 9 genes were differentially expressed in LUSC tumor samples compared with normal samples (Figure S2(A)). Additionally, their expression was also validated in two independent cohorts. In both datasets, in addition to CCDC14, EPB41L1, PPP2R5C and VPS9D1, the other five genes were all differentially expressed in LUSC tumor samples compared with normal samples, with similar expression pattern in two datasets and TCGA cohorts (Figure S2(B,C)). The results suggested that these prognostic AS events might play important roles in LUSC by affecting the expression pattern of genes.

Construction of regulatory network between SF and prognosis-related AS events

SFs are major contributors of AS and play important roles in oncogenesis. A total of 149 SF-AS interaction pairs between 10 AS events and 59 SFs were identified (Figure 5). Seven of them were regulated by many SFs. Notably, SF CELF2 were significantly correlated with the 10 AS events, while HNRNPC, NOVA2, SRSF3 were associated with 9 of 10 AS events. This close correlation between AS events and SFs suggests the involvement of SFs in the regulation of AS in LUSC, and contributes to the regulation of LUSC progression.

Figure 5. Regulatory network of AS and splicing factors. The diamond nodes represent splicing factors (SF); the ellipse nodes represent favorable alternative splicing (AS) events, and the circle nodes represent unfavorable AS events; the lines represent the correlation between AS events and SFs, respectively.

GSEA pathway enrichment analysis

Sixty-seven KEGG pathways showed significant differences between high-risk group and low-risk groups. The top 5 pathways with the largest normalized enrichment score in high-risk group and low-risk groups were plotted (Figure 6(A,B)). Pathways involved in cell interaction (cell adhesion molecules) and immune response (complement and coagulation cascade, Leishimania infection, Lysosome, and Viral myocarditis) were observed in the high=risk group, while RNA splicing-related pathways (spliceosome and RNA degradation) were found in the low-risk group.

Figure 6. GSEA pathway enrichment analysis. The top 5 differential KEGG pathways enriched in the high-risk group (A) and low-risk group (B).

Discussion

LUSC is a common type of lung cancer with a poor prognosis, and many new cases are diagnosed every year. Evaluation and prediction of the clinical outcome of LUSC are critical for patient care and disease management. Abnormality in AS regulation is prevalent in various diseases and contributes to disease progression. AS is emerging as a new hallmark of cancer and contributes to tumor development, drug resistance, immune evasion and metastasis. Dysregulation of AS is also identified in cancers, including LUSC. Our current study demonstrated that the oncogenic MYC-related AS events exhibited potent prognostic power for the overall survival of LUSC. Furthermore, three subtypes of LUSC were classified based on the expression levels of these AS events. Each LUSC subtype showed a different clinical outcome. Further study is required to interrogate the relationship of Myc with AS in LUSC and understand the roles of these AS events in the pathogenesis of LUSC.

AS is conservative and tightly regulated and it showed considerable flexibility to adapt various surrounding environments. AS is found in multiple diseases and it participates in the pathogenesis of many diseases (Kelemen et al. 2013; Scotti and Swanson 2016). While consistent results were also observed in our study, we further interrogated the association of these AS with cancer-related hallmark gene sets, and our data showed that these differential AS events in LUSC were significantly associated with Hallmark_Apical_Surface and Hallmark_Myc_Targets_V1 pathways. Myc signaling is under tight control in normal cells but abnormally activated in tumors. Myc signaling is implicated in the initiation and maintenance of various cancers. Studies have shown that Myc contributes to the malignancy of LUSC (Chanvorachote et al. 2020). Myc can be regulated by many upstream signaling including JAK/STAT pathway, Wnt/β catenin signaling, Notch signaling, and Ras/PI3 K/AKT/GSK-3 pathway to promote cancer cell proliferation, survival, metastasis and drug resistance through direct or indirect modulation of target gene expression (Chanvorachote et al. 2020). Intriguingly, our results demonstrated that pathways involved in alternative splicing regulation were enriched in Hallmark_Myc_Targets_V1 pathway-associated AS genes. This agrees with previous study that demonstrated that oncogenic Myc signaling regulates alternative splicing in different cancers (Koh et al. 2015; Phillips et al. 2020). These results suggest that AS is an important target of Myc signaling in LUSC and insinuates the contribution of AS in the carcinogenesis of LUSC under the control of oncogenic Myc signaling. Furthermore, Myc-guided AS events stratified LUSC patients into 3 clusters with different prognoses. Among them, cluster 1 contained more advanced patients and showed a poorer prognosis than other clusters. These results further strengthened the importance of Myc-associated AS events in the malignancy of LUSC and highlighted their prognostic values.

Among those MYC-guided prognostic AS events, we obtained 10 most valuable AS events after LASSO regression analysis, involving 9 genes, including DNAJC2, EFCAB11, PPP2R5C, C4BPB, VPS9D1, PPARG, CCDC14, AQP1 and EPB41L1. PP2A B subunit isoform B56-gamma (PPP2R5C) is a regulatory subunit of protein phosphatase 2A (PP2A) and a tumor suppressor by activating p53 or suppressing ERK signaling (Letourneux et al. 2006; Li et al. 2007; Shouse et al. 2008). AS of PPP2R5C was reported in gastrointestinal adenocarcinomas and ovarian cancers (Lin et al. 2018; Zhu et al. 2018), but its functional consequences remain unknown. cancer cells may employ the AS machinery to inactivate the tumor suppressor activity of PPP2R5C. Peroxisome proliferator-activated receptor gamma (PPARG) is a transcription factor that regulates fatty acid β-oxidation and plays an important role in cancer development (Tachibana et al. 2008), AS of PPARG is also identified in tumors with unknown effects (Sebestyen et al. 2015). AS of other tumor-relevant genes (VPS9D1, EPB41L1, AQP1) also showed prognostic ability in various cancers (Wang et al. 2017; Song et al. 2019; Liang et al. 2020; Du et al. 2021). Besides, our study identified several novel AS events in cancers. AS of EF-hand calcium-binding domain 11 (EFCAB11) was reported in macrophages upon Mycobacterium tuberculosis infection but not in tumors (Kalam et al. 2017). AS of Coiled-coil domain-containing protein 14 (CCDC14) is not reported yet. DnaJ heat shock protein family member C2 (DNAJC2, also known as ZRF1) is a chromatin regulator and participates in epigenetic modulation of gene transcription in cancers (Aloia et al. 2015), but AS of DNJC2 is not found in other cancers or diseases. Complement component 4 binding protein beta (C4BPB) regulates the activation of the complement cascade and circulating C4BPB may serve as indicators for cancer prognosis and metastasis (Battistelli et al. 2008; Liu et al. 2012). AS of C4BPB is not reported before. These AS events in cancer-associated genes insinuate their functions in carcinogenesis, and the novel AS events may also contribute to tumor progression rather than act as prognosis biomarkers.

The 10 AS events-based prognostic model stratified LUSC patients into high- and low-risk groups, and patients in the high-risk group showed significantly poorer prognosis than those in the low-risk group. Intriguingly, the distribution of LUSC clusters in the high-risk and low-risk groups varied greatly. Samples in Cluster 2 subtype (favorable survival) were dominant in the low-risk group, while samples in Cluster 1 subtype (worse survival) were highly enriched in the high-risk group, which was consistent with the prognosis of LUSC clusters. These data demonstrated that our prognosis model could effectively predict the risk of LUSC and may be helpful for the clinical management of LUSC. A previous study also investigated AS in LUSC and identified a survival-associated AS signature (Liu et al. 2020). While consistently discovering the prognostic value of AS events for LUSC, our study interrogated oncogenic Myc pathway-guided AS events rather than global AS in LUSC. Furthermore, the MYC-guided AS events classified 3 different subtypes of LUSC, which revealed the internal diversity of LUSC and deepened our understanding and further management of LUSC. Given the regulatory roles of oncogenic MYC signaling in AS, the MYC-associated AS events may not only act as prognosis markers but also contribute to the malignancy of LUSC. Further study is warranted to explore the oncogenic roles of MYC-associated AS events.

Our current study identified a set of oncogenic Myc signaling-guided AS events in LUSC and classified 3 subtypes of LUSC with different clinical outcomes based on these AS events. We established a powerful Myc-guided prognostic AS-based prognosis model stratifying LUSC patients with different risks. Our findings deepened our understanding of the prognostic power of Myc-associated AS events and suggested their essential roles in the pathogenesis of LUSC. Further studies are required to elucidate the underlying mechanism(s) of Myc-related AS in the progression of LUSC.

Acknowledgements

YL contributed to the conception and design of the research and drafted the manuscript ; KW contributed to the acquisition of data: FL and WL contributed to the analysis and interpretation of data ; YS contributed to statistical analysis; ; JX revised the manuscript for important intellectual content. All authors read and approved the final manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The (log2(fpkm+1) RNAseq data and phenotype data of Lung Squamous Cell Carcinoma (LUSC) used in this study is available publicly from The Cancer Genome Atlas (TCGA) database at https://xenabrowser.net/. The percent spliced index (PSI) values of alternative splicing events of TCGA-LUSC samples were available from the SpliceSeq database http://projects.insilico.us.com/TCGASpliceSeq/. GSE19188 and GSE30219 datasets were available from the Gene Expression Omnibus database at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc.

Additional information

Funding

This work was supported by Kunming Municipal Health Committee Health Research Project [grant number 2020-04-02-115]; Science Research Fund Project of Yunnan Provincial Department of Education [grant number 2018JS226]; Study of genetic risk of the inflammatory body associated genes in Xuanwei lung cancer and its pathogenic mechanism [grant number D-2017013]; The research about PDGFRB functions on lung squamous cell carcinoma progression and its potential usage as a clinical lung squamous cell carcinoma maker [grant number 2017BS029]; Yunnan Provincial Health and Family Planning Commission Reserve Talent Project [grant number H-2017013]; and 2018 CSCO-Qilu Tumor Project [grant number Y-Q201802-011].