Transcriptome profiling and clinicopathologic characteristics of 522 LUAD specimens were retrospectively curated from the Cancer Genome Atlas (TCGA) project utilizing TCGAbiolinks R package (23). Table 1 lists clinicopathological data of 522 LUAD patients. AS data were curated from TCGA SpliceSeq (https://bioinformatics.mdanderson.org/TCGASpliceSeq) (24). Then, Percent Spliced In (PSI) values that ranged from 0 to 1 were determined for AS events across transcripts. AS events were classified into seven forms, containing Alternate Donor site (AD), Alternate Acceptor site (AA), Alternate Terminator (AT), Alternate Promoter (AP), Mutually Exclusive Exons (ME), Exon Skip (ES), and Retained Intron (RI). AS events with PSI value ≥ 75%, and average PSI value ≥ 0.05 were enrolled for subsequent analysis. UpSetR package was employed for visualizing the distribution of AS events in LUAD (25).

OS-relevant AS events were selected across LUAD patients through the survival R package utilizing univariate regression analyses following the criteria of p-value < 0.05. In addition, UpSet and volcano plot were adopted for describing the distribution of OS-relevant AS events. Thereafter, the first 20 AS events in different types of AS were visualized into bubble plots.

SpliceAid project was employed to curate specific splicing factors (26). Furthermore, Pearson correlation test was adopted for analyzing the interactions of splicing factors with OS-relevant AS events. The Cytoscape (version 3.8.0) was utilized for visualizing this interactional network of splicing factors with OS-relevant AS events and correlation coefficient > 0.6 as well as p < 0.05 as the filtering criteria (27).

The glmnet R package was adopted to establish a least absolute shrinkage and selection operator (LASSO) prognostic model based on OS-relevant AS events across LUAD patients (28). The prognostic scoring formula was conducted with this formula: risk score= PSI value of AS event1 × Coef1 + PSI value of AS event2 × Coef2 … + PSI value of AS eventn × Coefn, in which Coefn represented the regression coefficient. Then, we stratified LUAD patients into different risk subpopulations according to median risk score. The receiver operating characteristic (ROC) curve was generated utilizing timeROC R package for showing the specificity and sensitivity of risk score in evaluating prognosis of LUAD. The Kaplan-Meier curves were applied to assess the differences in OS rate with the survival R package. Additionally, Cox regression models were conducted for analyzing the interactions of age, gender, tumor stage, and risk score with OS outcomes.

In order to evaluate OS outcomes, a prognostic nomogram comprised of independently prognostic indicators AS-relevant risk signature as well as stage was conducted for estimating 1‐, 2‐, and 3‐year OS probabilities with the rms R package. Subsequently, calibration curves which showed the survival implications of this nomogram were depicted. The calibration curve close to 45° was considered as an excellent indicator in this nomogram.

The cell type identification by Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) deconvolution algorithm was adopted to estimate the abundances of 22 diverse leukocyte subsets (29). CIBERSORT results for samples with p < 0.05 indicated that the estimated abundances of leukocyte subsets were reliable, which were eligible for subsequent analysis. For each specimen, estimations were standardized to sum up to 1, thereby being interpreted directly as cellular fraction. The ssGSEA from Gene Set Variation Analysis (GSVA) was employed for quantification of the relative abundances of 29 immune cells as well as functions following the special feature gene panels across LUAD specimens (30). The ssGSEA enrichment score was indicative of the relative abundance, which was standardized to range from 0 to 1.

The Estimation of Stromal and Immune Cells in Malignant Tumors using Expression Data (ESTIMATE) R package possesses the significant advantage in estimating the specific features of transcriptome profiles (31). The gene sets of immune checkpoints were downloaded from recent research (32, 33). The mRNA expression of immune checkpoints was quantified across LUAD specimens. Tumor mutational burden (TMB) was employed to predict clinical response to immunotherapy (34). TMB was calculated according to the formula: (entire counts of variants)/(the entire lengths of exons) in line with the variants of LUAD specimens that were extracted from the mutational profiles.

Half-maximal inhibitory concentration (IC50) values for cisplatin, gemcitabine, gefitinib, and erlotinib were estimated with the pRRophetic R package by ridge regression analysis (35, 36). IC50 indicated the treatment response to above chemotherapeutic agents in TCGA cohort.

Spearman’s correlation analysis was conducted to estimate composition differences. Wilcoxon signed rank test was applied for comparisons in two groups. Kaplan-Meier survival curve was implemented for evaluating the survival differences between groups. Cox regression analysis was conducted for verifying the associations of certain indicators with LUAD prognosis. To evaluate the performance of prognosis prediction, time-independent ROC curves were conducted and area under the curve (AUC) was calculated with timeROC R package. Statistical analysis was achieved utilizing R software (version 4.02). P <0.05 was taken into consideration statistically.

In total, 43,945 AS events were identified across 522 LUAD patients ( Figure 1A ). ES accounted for the most frequent AS signature, followed by AT and AP. Univariate analyses were presented to qualify the impact of each AS event on patients’ OS. Subsequently, 2415 AS events displayed remarked associations with survival outcomes of LUAD patients, in which 1356 were protective factors and 1059 were risk factors ( Figure 1B ). Notably, one gene may possess two or more OS-relevant AS events across LUAD patients, as shown in the UpSet plots ( Figure 1C ). The first 20 significant OS-relevant genes of AS events are separately shown in Figures 1D–J , which indicated that the seven alternative splicing modes exhibit great variability.

Splicing factors act as dominant regulators of AS events, which may affect the splicing of oncogenes as well as tumor suppressors (37). For exploring the underlying interactions of the expressions of splicing factors with AS events, we visualized the splicing-regulatory network, as depicted in Figure 1K . In total, three splicing factors (including SEC31B, CLK1, and DDX39B) displayed prominent associations with 44 OS-relevant AS events. Furthermore, most favorable AS events exhibited positive interactions with the expression of splicing factors and the three splicing factors were in relation to multiple AS events. Thus, splicing factors may act as an indispensable role in modulating AS events during lung carcinogenesis.

For avoiding over-fitting, LASSO Cox analysis was adopted for developing a prognostic model of LUAD on the basis of OS-relevant AS events. Through cross-verification, the optimal parameters were selected ( Figure 2A ) and the coefficients in LASSO regression model were determined ( Figure 2B ). Ultimately, 12 OS-relevant AS events (BEST3|23330|AT, CDKN2A|86004|AP, TTC39C|44852|AP, MEGF6|315|ES, PKIB|77377|AP, CA5B|98313|ES, HNRNPLL|53258|AT, LDB1|12935|AP, C12orf76|24406|AT, AP2B1|40327|AD, LETM2|83398|AT, MRPL33|53046|ES) were identified ( Table 2 ). In line with the regression coefficients and PSI value of 12 OS-relevant AS events, we calculated risk scores of LUAD patients. Thereafter, LUAD patients were classified into different groups with median risk score of 0.8834 ( Figure 2C ). Moreover, we noticed that high-risk subpopulations were often accompanied by high mortality ( Figure 2D ). Heatmap depicted the heterogeneity in PSI values of 12 OS-relevant AS events ( Figure 2E ). Prognostic analyses uncovered that high-risk subpopulations exhibited remarkedly dismal OS outcomes ( Figure 2F ). The validity of the prognostic model in prognosis prediction was verified through ROC analysis. The AUC values at 1-, 3-, and 5-year OS were separately 0.762, 0.770, and 0.725, showing the good effectiveness of this model in prognosis prediction ( Figure 2G ).

Through ROC analysis, we presented the comparisons of AUC values and noticed that risk score displayed the higher AUC values under 1-, 3-, and 5-year survival compared with clinicopathological features (age, gender, and stage; Figures 3A–C ). Additionally, the differences in risk score between distinct clinicopathological features were compared among LUAD patients. No significant differences were observed between age ≤65 and > 65 ( Figure 3D ) as well as between non-metastasis (M0) and metastasis (M1; Figure 3E ). Increased risk score was investigated in male and female patients ( Figure 3F ). As T, N, and stage increased, risk score was gradually elevated ( Figures 3G–I ), indicating that the prognostic model contributed to LUAD progression.

We further verified the prognostic value of clinical characteristics and risk score and found that risk score and stage possessed the potential to independently predict LUAD prognosis ( Figures 4A, B ). Thereafter, a prognostic nomogram containing independent prognostic indicators risk score as well as clinicopathological stage was conducted for forecasting patients’ outcomes ( Figure 4C ). Calibration curves were indicative of the powerful prognostic predictive capacity of this nomogram in 1-, 3-, and 5-year OS ( Figures 4D–F ).

For further applying our findings to clinical practice, this study constructed a nomogram prognostic score system in the prediction of 1-, 3-, and 5-year OS outcomes of LUAD patients ( Figure 4C ). The scoring system included the prognostic model and stage. Thereafter, for verifying the reliability of the prognostic nomogram, calibration plots were conducted and confirmed the practical significance of the model. As depicted in Figures 4D–F , the model possessed the potential in determining survival outcomes with a high predicted accuracy.

We firstly estimated infiltration levels of immune and stromal cells across LUAD patients. Accordingly, patients with a high-risk score displayed reduced immune score and stromal score ( Figures 5A, B ). Nevertheless, higher tumor purity was investigated in high-risk patients ( Figure 5C ). Then, we determined ESTIMATE and noticed the prominently decreased ESTIMATE score in the high-risk group ( Figure 5D ). Thus, low-risk tumors were accompanied by abundant infiltrations of immune and stromal cells. Then, we systematically investigated the immune cell infiltration landscape across LUAD with CIBERSORT algorithm. We noticed that the low-risk group displayed high infiltration levels of B cells naïve, T cells CD4 memory resting, monocytes, and mast cells resting ( Figure 5E ). Oppositely, the high-risk group exhibited increased infiltration levels in T cells CD4 memory activated, T cells follicular helper, T cells regulatory (Tregs), NK cells resting, macrophages M0, and macrophages M1. Supplementary Figure 1 displays interactions of risk score signature with above immune cell infiltrations. Subsequently, we revealed the activities of immune functions and immune cell infiltrations across LUAD ssGSEA method. In Figures 5F, G , higher abundance levels of aDCs, B cells, HLA, iDCs, mast cells, neutrophils, T helper cells, TIL, and type II IFN response were investigated in the low-risk group while MHC class I and NK cells exhibited higher abundance levels in the high-risk group. We also evaluated the interactions of the prognostic model with immune checkpoints across LUAD. As depicted in Figure 5H , this prognostic model possessed a positive association with CD274. Moreover, we observed that the low-risk group was characterized by increased expression of most immune checkpoint-related genes ( Figure 5I ). Thus, low-risk patients were indicative of higher immune response as well as immune cell infiltration.

The interaction of the prognostic model with TMB was also observed across LUAD. As shown in Figure 6A , the high-risk group exhibited a remarkedly increased TMB score. Moreover, we presented survival analysis among diverse subgroups. We noticed that subpopulations possessing an elevated TMB score as well as a reduced risk score displayed the most favorable survival outcomes while those with a low TMB score and high-risk score exhibited the poorest survival outcomes ( Figure 6B ). Chemotherapy and targeted therapy were gradually applied in treatments for patients with advanced lung adenocarcinoma. It is of great significance to evaluate the responses of certain drugs in different risk subpopulations. Herein, we identified the treatment responses of some drugs that were widely used in the treatment of LUAD. As shown in Figures 6C–F , the high-risk group possessed prominently lowered IC50 values of cisplatin, erlotinib, gefitinib, and gemcitabine, indicating that this subpopulation possessed higher sensitivity to these therapeutic agents. The above findings provide more clues for individualized treatment strategies in LUAD patients.

We found that CDKN2A, PKIB, and TTC39C exhibited a higher expression in LUAD than normal tissues among the 12 AS events-relevant genes in the prognostic models ( Figures 7A–C ). Moreover, survival analysis uncovered that highly expressed CDKN2A and PKIB were in relation to more dismal survival probabilities of LUAD ( Figures 7D, E ). In contrast, high TTC39C expression was indicative of the marked survival advantage ( Figure 7F ).

We further investigated the interactions of the three prognostics AS events-related genes (CDKN2A, PKIB, and TTC39C) with immune response and immune cell infiltration across LUAD. We found that deregulated CDKN2A did not affect estimate, immune, and stromal score as well as tumor purity ( Figures 8A–D ). For Figures 8E–H , high PKIB expression was characterized by increased estimate, immune, and stromal score but reduced tumor purity. Moreover, high TTC39C expression displayed remarkedly decreased estimate, immune, and stromal score but elevated tumor purity ( Figures 8I–L ). In Figure 8M , CDKN2A upregulation was in relation to increased infiltration levels of T cells CD8, T cells CD4 memory activated, and macrophages M1. PKIB deregulation was in relation to infiltrations of B cell naïve, B cells memory, plasma cells, T cells CD8, macrophages M1, macrophages M2, dendritic cells resting and mast cells resting ( Figure 8N ). B cells native, plasma cells, T cells follicular helper, and NK cells activated exhibited the increased infiltration levels in high TTC39C expression group ( Figure 8O ). The ssGSEA results uncovered the increased infiltrations of APC co-inhibition, CD8+ T cells, inflammation-promoting, MHC class I, NK cells, T cell co-stimulation, Tfh, and Th1 cells in high CDKN2A expression group ( Figure 8P ). PKIB upregulation was in relation to most immune functions and immune cell infiltrations ( Figure 8Q ). In Figure 8R , we noticed the prominent interactions of high TTC39C expression with activation of most immune functions and immune cell infiltrations. We also estimated the associations of CDKN2A, PKIB, and TTC39C with immune checkpoint molecules. Most immune checkpoint molecules exhibited positive interactions with CDKN2A, PKIB, and TTC39C ( Figures 8S–U ).