SpliceSeq data of TCGA-LUAD were downloaded from the TCGA database (https://tcga-data.nci.nih.gov/tcga/) (Tomczak et al., 2015). SpliceSeq tool, a java-based application, was usually used to unambiguously quantify the mRNA splicing levels of samples in TCGA. A novel value could be calculated by SpliceSeq based on seven types of AS events about each protein-coding gene provided from the Ensemble gene database (Ryan et al., 2012). For the following seven kinds of AS events, the Percent Spliced In (PSI) value was calculated, quantifying splicing event levels range from 0 to 1: Mutually Exclusive Exons (ME), Exon Skip (ES), Alternate Promoter (AP), Retained Intron (RI), Alternate Acceptor site (AA), Alternate Donor site (AD), and Alternate Terminator (AT). The schematic graph explaining these seven types of AS is shown in Supplementary Figure S1A.

Patients’ clinical parameters of LUAD were downloaded and extracted from the TCGA database. In total, 444 LUAD patients were included in this analysis. Clinical information for these patients and the pathological details obtained from TCGA are provided in Supplementary Table S1. The PSI value of AS events in samples were collected and subjected to univariate Cox analysis. All the AS events screen for a p value <0.05, and these events were considered as candidate prognosis-related events. The “glmnet R″ software package was used to perform the least absolute shrinkage and selection operator (LASSO) analysis to filter out the most valuable and concise AS events in all AS events filtered in univariate Cox analysis (p < 0 .05). Afterward, the prognostic independence of the AS signature was constructed by multivariate Cox analysis. Then, based on the coefficient of each above AS event, each patients’ risk score could be calculated by the signature, respectively. Meanwhile, all patients were divided into distinct subgroups based on the median value of risk scores.

The Kaplan-Meier curve was implemented to evaluate the differential survival status in both groups. The receiver operating characteristic curves (ROC) were conducted to detect both the sensitivity and specificity of prognostic signatures using the “survivalROC” R package (https://www.r-project.org/, v3.5.3).

We developed the Upset intersective plot, a more scalable visualizing diagram than Venn, which was used to explore the interactive sets of AS events. The “UpSet” R package was used to visualize their potential interrelationship. The expression data of the Splicing factors (SFs) was extracted from TCGA-LUAD mRNA-seq data. All SF genes were subjected to univariate Cox analysis when their p < 0.05. These SFs were considered the survival-associated splicing factors. The relationship closeness between SFs expression value and AS’s PSI value were calculated by the Spearman test. At the same time, the interaction network diagram of these SFs and prognosis-related AS events was illustrated using Cytoscape 3.7.0 (https://cytoscape.org/).

22 types of tumor-infiltrating immune cells were estimated using the CIBERSOFT algorithm, which characterized the cellular composition of complex tissues based on normalized gene expression profiles. The comparison of immune cell distribution between high and low-risk groups was made using the Mann-Whitney U-test.

Potential drugs for the treatment of LUAD were selected using the Connectivity Map (CMap) database (https://clue.io/). We uploaded differentially expressed genes between high and low-risk groups into the CMap database for genomic enrichment analysis. We screened the small drug molecules with an enrichment score of <90, and obtained 3D structure through PubChem database (http://www.pubchem.ncbi.nlm.gov), a public repository of small molecules in properties. The IC50 value of these small molecule drugs was provided by the GDSC database (https://www.cancerrxgene.org/).

Comprehensive AS events were examined in a cohort of 444 TCGA-LUAD patients (Supplementary Figure S1A).A total of 30,735 AS events from 9,635 genes were detected, including 11,768 ES events in 5,467 genes, 6,129 AP events in 3,424 genes, 5,782 AT events in 3,372 genes, 2,605 AA events in 1,961 genes, 2,199 AD events in 1,659 genes, 2,103 RI events in 1,461 genes, 149 ME events in 146 genes (Supplementary Figure S1B). In TCGA-LUAD, ES events were the most common AS events, accounting for approximately just over one-third of all events, followed by the number of AP and AT events, while the number of ME events was the least. Remarkably, the number of AS events went far beyond those of their corresponding mRNAs. Furthermore, a subset of overlapping AS events in the seven types of AS in LUAD was displayed by the UpSet plot diagram (Supplementary Figure S1C).

First, we conducted a univariate Cox analysis based on the 30,735 AS events related to 444 patients to appraise the relationship between AS events and overall survival (OS) status in LUAD. Consequently, 1,429 AS events and corresponding 1,125 genes were conspicuously related to the overall survival of LUAD patients (Figure 1A). Figures 1B–H showed the top 20 most important AS events related to OS among these seven types of AS events. Interestingly, some special survival-associated AS genes contain multiple types of AS events. For example, AD, AT, AA, RI, and ES of C1orf159 and AP, AA, RI, AD, and ES of MRPL55 were all related to OS of LUAD patients.

We selected the significant prognostic associated AS events as candidates by univariate Cox analysis, aiming to further screen out the most significant AS events related to patient prognosis by LASSO Cox regression analysis (Supplementary Figure S2). Furthermore, several prediction signatures based on these prognostic associated AS events were constructed by multivariate Cox analyses. Eventually, a combined prognostic model was built, integrated from different types of AS events (Supplementary Table S2). The median value of risk scores was considered as the cutoff criteria for dividing patients into a high-risk group and a low-risk group. The Kaplan–Meier curves shown in Figures 2A–H, showed that LUAD patients in the high-risk group had appreciably shorter OS than patients in the low-risk group, demonstrated that these AS signatures could be powerful biomarkers to distinguish patient prognosis. The combined prognostic signature showed better predictive properties than the single type of AS events (Figure 2H). Then, the ROC curve was performed to appraise the prognostic efficiency of prognostic AS models. The results show that all signatures had a robust predictive property with AUC values from 0.708 to 0.849, except the ME signature (AUC = 0.582, Figure 2I). Conceivably, the combined model contains different types of AS events that had the highest efficiency (AUC = 0.849). The distribution diagram of patients’ risk score, survival status, and expression profiles of related AS model events are shown in Supplementary Figures S3A–H.

To confirm the prognostic value of our combined AS signature, we randomly selected 50% of patients in all LUAD patients as a testing cohort and performed Kaplan-Meier cure and ROC analysis (Supplementary Figure S4). The results were consistent with the above results and demonstrated the combined AS signature had a robust predictive ability (AUC = 0.856).

To explore the predicted performance of AS signatures and other clinical characteristics of survival, univariate and multivariate Cox analyses were performed to find out if these AS signatures were independent prognostic factors for LUAD patients. The univariate Cox analysis showed that almost all risk score signatures (except ME-risk score signature), TNM stage, T, N, and M stage were remarkably related to LUAD patients’ overall survival (Supplementary Figure S5). Furthermore, multivariate Cox analysis results show that most risk score signature, T and N stage still have a predictive ability when all univariate significant factors are considered together, suggesting that risk score, T and N stage were unassisted risk elements (Figures 3A–G). Taken together, all these results demonstrated that AS signatures exhibit powerful predictive performance in LUAD patients. In addition, circos plots were depicted to display the details of AS events and their interacting genes in the chromosome (Figure 3H).

Interestingly, extensive dysregulated AS events in many types of cancers are easily programmed by some specific SFs. Thus, an interesting issue is whether several key SFs could regulate these prognosis-associated AS events in LUAD. To determine those specific SFs which had a close connection with prognosis-associated AS events in LUAD, univariate Cox analysis of all SFs were performed based on the gene expression value of LUAD patients. The results showed that there were 30 SFs related to the OS of LUAD patients (Supplementary Table S3). Furthermore, correlations between SFs and prognostic AS events were tested by Spearman’s test (Figure 4A). In the correlation networks, 27 SFs (purple dots) were related to 248 prognosis-associated AS events, involving 141 favorable AS events (green dots) and 107 adverse AS events (red dots). Intriguingly, there was a positive relationship (red lines) between most of the poor survival-related AS events (red dots) and SFs (purple dots), while there was a negative correlation (green lines) between most of the favorable survival-related AS events (green dots) and SFs. For example, SFs RNF34 and HNRNPK were related to worse survival of LUAD patients (Figures 4B,C). ES events of YPEL5 were an adverse factor, while the AP events of PCNA and AP events of PDHX were related to a favorable prognosis. The relationship between RNF34 and the AP of PCNA or ES of YPEL5 were shown in dot plots, implicating the high expression of RNF34 was negatively related to adverse prognosis (Figures 4D,E). Similarly, the relationship between HNRNPK and the AP of PDHX or ES of YPEL5 were shown in dot plots (Figures 4F,G), illustrating that the high expression of HNRNPK was positively related to poor prognosis.

We finally constructed an mRNA signature for LUAD patients by subjecting differentially expressed mRNAs in LUAD to univariate and multivariate Cox analysis: mRNA risk score = (0.1674*XPR1) + (0.0529*MMP1) + (-0.2108*MAP3K8) + (0.0882*RHOV) + (0.1099*CDH17) + (0.1488*UCA1). Then, Kaplan-Meier and ROC curves were implemented to contrast prognostic ability between AS signature and mRNA signature. Both results from Kaplan-Meier and ROC analyses illustrated that AS signature had significantly better survival and higher ROC than mRNA signature (Figures 5A,B). These data demonstrated that the predictive power of AS signature is greater than the mRNA signature. In general, AS signature could be used as a superior indicator to predict the prognosis of LUAD patients.

We also analyzed tumor microenvironment cell infiltration of signature. We calculated the median absolute score that CIBERSORT gave for 22 cell types in two groups. The results showed that the fraction of T cells follicular helper, NK cells resting, Monocytes, and Macrophages M1 was significantly higher in high-risk than that in the low-risk group. However, B cells naive, Macrophages M0 and Mast cells resting were remarkably higher in the low-risk than in high-risk group (Figure 6A). Furthermore, the results showed that the expression of some immunomodulator agonists was correlated with a risk score. We found that some immune ligands and receptors were significantly higher in the high-risk group, while antigen presentation was mainly higher in the low-risk group (Figure 6B).

To predict small molecule drugs that can inhibit the resistance of LUAD patients, DEGs of high-risk and low-risk groups were assigned into up-regulated and down-regulated groups. Then we matched it to a small molecule drug in the CMap database. Finally, we selected the three most important small molecular compounds, including AZ-628 (score = −95.17), NVP-AUY922 (score = −91.77) and Nomifensine (score = −90.86). The 3D structure of these three small molecular compounds was downloaded from the PubChem database (Figure 7). These small molecules might potentially improve the outcomes of LUAD patients and provide recommendations for the selection of LUAD-targeted drugs, while the specific mechanism and effectiveness need to be further studied. We also explored the IC50 values of these small molecule drugs in LUAD cells through the GDSC database. The results showed that NVP-AUY622 (IC50 = 0.004–10.222 μM) had lower IC50 values than AZ-628 (IC50 = 0.050–46.073 μM), indicating that NVP-AUY622 has a stronger effect in LUAD cells (nomifensine not provided). These results demonstrated that NVP-AUY922 might become a novel drug to improve the survival of LUAD patients.