The pattern of alternative splicing in lung adenocarcinoma shows novel events correlated with tumorigenesis and immune microenvironment

Abstract

Lung adenocarcinoma (LUAD) is the leading cause of cancer deaths worldwide due to the lack of early diagnostic markers and specific drugs. Previous studies have shown the association of LUAD growth with aberrant alternative splicing (AS). Herein, clinical data of 535 tumor tissues and 59 normal tissues were extracted from The Cancer Genome Atlas (TCGA) database. Each sample was analyzed using the ESTIMATE algorithm; a comparison between higher and lower score groups (stromal or immune) was made to determine the overall- and progression-free survival-related differentially expressed AS (DEAS) events. We then performed unsupervised clustering of these DEASs, followed by determining their relationship with survival rate, immune cells, and the tumor microenvironment (TME). Next, two prognostic signatures were developed using bioinformatics tools to explore the prognosis of cases with LUAD. Five OS- and six PFS-associated DEAS events were implemented to establish a prognostic risk score model. When compared to the high-risk group (HRG), the PFS and OS of the low-risk group (LRG) were found to be considerable. Additionally, a better prognosis was found considerably associated with the ESTIMATE score of the patients as well as immune cells infiltration. Our analysis of AS events in LUAD not only helps to clarify the tumorigenesis mechanism of AS but also provides ideas for revealing potential prognostic biomarkers and therapeutic targets.

Keywords: Alternative splicing; Immune; LUAD; Prognosis; Tumor microenvironment.

Fig. 1

The relationship between stromal/immune scores and prognosis of LUAD patients. A The K–M survival curves for OS of high and low immune scores groups. B The K–M survival curves for OS of high and low stromal scores groups. C Comparison of immune scores between T stages. D Comparison of immune scores between M stages. E Comparison of stromal scores between M stages

Fig. 2

DEAS events between high and low stromal/immune scores groups. A The UpSet plot of intersections between AS events and the corresponding gene intersections. The heatmaps (B, D) and volcano plots (C, E) of DEAS events between the high and low stromal/immune scores groups. F, G The UpSet plots of DEAS events and the corresponding gene intersections of the immune group and the stromal group

Fig. 3

Intersection DEAS events and GO and KEGG pathway enrichment analysis. The up-regulated (A) or down-regulated (B) DEAS events in both stromal score and immune scores groups by Venn diagram. C GO analysis. D KEGG pathway enrichment analysis

Fig. 4

The immune microenvironment was closely related to the prognosis of LUAD patients. A TCGA LUAD cohort was clustered into three subgroups by unsupervised cluster analysis. B The heatmap of clusters and clinical parameters. C K–M survival curves for OS of three clusters. D, E The comparison of immune scores and stromal scores between three clusters. F The comparison of immune cell infiltration between three clusters. G, H GSVA enrichment analysis showing the activation states of biological pathways in distinct clusters (C1–C2, C1–C3)

Fig. 5

Establishment of a 5-DEAS-based OS signature. A The forest plot of the univariate Cox analyses. B LASSO coefficient profiles of the candidate survival-related DEASs. A coefficient profile plot was produced against the log λ sequence related to OS. C Dotted vertical lines were drawn at the optimal values using the minimum criteria related to OS. D ROC curve analysis to evaluate the sensitivity and specificity of the gene signature. E K–M survival curve to test the predictive effect of the gene signature. F DEAS expression levels of the high-risk and low-risk groups with the OS signatures. G OS scatter plots for LUAD patients. H Risk score distribution of patients with the OS signatures

Fig. 6

Establishment of a 6-DEAS-based PFS signature. A The forest plot of the univariate Cox analyses. B LASSO coefficient profiles of the candidate survival-related DEASs. A coefficient profile plot was produced against the log λ sequence related to PFS. C Dotted vertical lines were drawn at the optimal values using the minimum criteria related to PFS. D ROC curve analysis to evaluate the sensitivity and specificity of the gene signature. E K–M survival curve to test the predictive effect of the gene signature. F DEAS expression levels of the high-risk and low-risk groups with the PFS signatures. G PFS scatter plots for LUAD patients. H Risk score distribution of patients with the PFS signatures

Fig. 7

Development of two nomograms combining DEAS-based signature and independent prognostic clinical variables to predict OS and PFS in LUAD patients. A Nomogram of OS combining the OS signature and three clinical variables of LUAD patients. B–D Calibration curves of the nomogram at 1, 2, and 3 years. E Nomogram of PFS combining the PFS signature and two clinical variables of LUAD patients. F–H Calibration curves of the nomogram at 1, 2, and 3 years

Fig. 8

Regulatory network between SFs and AS events in LUAD. A OS, B PFS

Fig. 9

Analysis of the crucial SOD2-78,301-AT DEAS event. The K–M survival curves for OS (A) and PFS (B) of high and low SOD2-78,301-AT expression level groups. C The relationship between SOD2|78,301|AT and ESRP2. D The correlation between SOD2 expression in LUAD and abundance of immune infiltrates. E The correlation between somatic copy number alterations (SCAN) and abundance of immune infiltrates of SOD2