{'fraction_non_alphanumeric': 0.031900027947458776, 'fraction_numerical': 0.014532678564299118, 'mean_word_length': 5.2231055900621115, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 4, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Alternative Splicing Events in Immune Infiltration of Lung Adenocarcinoma

Background The exact mechanisms of lung adenocarcinoma (LUAD) etiology and pathogenesis remain unclear. Material/Methods In this study, we evaluated alternative splicing (AS) events in LUAD. We separately analyzed LUAD data from The Cancer Genome Atlas (TCGA) database and AS-related feature lists from SpliceSeq to develop risk models for AS events and further validated risk models for AS events. The association between immune-related features and the risk model of AS events was evaluated. Results We found that exon skip (ES) and mutually exclusive exons (ME) exhibited the most and least AS events, respectively. The risk score and stage of LUAD patients were strongly associated with prognosis and were independent influences on the prognosis of LUAD. The pathological stage (stage, T, N) and risk model were incorporated to construct a column line graph with better predictive ability. Risk models were associated with tumor microenvironment, and higher risk score values were associated with higher M2 macrophage content and lower M0 macrophage and helper T lymphocyte content. The correlation between core genes and immunity was further assessed by analyzing differentially-expressed genes between risk models. These results suggest an association between the level of risk for AS events and the density of immune cell infiltration. Conclusions Our findings suggest a clear association between AS risk model and patient prognosis, and was performed to confirm the biological relationship between AS and immunity.


Background
Lung cancer (LC) has the highest incidence and mortality rates among malignant tumors worldwide, with lung adenocarcinoma (LUAD), the predominant pathological type of LC, causing a great burden of disease [1]. Due to the lack of sensitivity of diagnostic tools and the limitations of diagnostic conditions, 75% of patients are at advanced stages at the time of diagnosis [2]. Prognostic assessment of LUAD and individualized treatment are very important to reduce metastasis and recurrence and to prolong survival [2]. It is gradually being accepted that LUAD recognition is not a sudden change, but the result of a combination of genetic and epigenetic changes playing a role [3]. Therefore, there is an urgent need to understand the molecular mechanisms of LUAD. With the advent of next-generation sequencing technologies, genes involved in the onset, progression, and metastasis of LUAD, such as KRAS, EGFR, ALK, and ROS1, are currently an important topic of research that can improve early screening and more precise individualized treatment of LUAD [1,3].
Recent studies have gradually established that alternative splicing (AS) is also responsible for carcinogenesis [4]. LUAD tissues exhibit more aberrant activation patterns of AS; however, the exact role of AS in lung cancer is unknown [5]. In recent years, epigenetic and transcriptomic studies have identified many valuable tumor markers that may provide new ideas for future diagnosis and treatment [6][7][8]. Immune response affects tumor development and sensitivity to drugs, while AS influences the type and function of immune cells [9,10]. In this study, we established an AS-related risk model by use of the bioinformatics method to further analyze the correlation between the risk model and immunity, and provide new ideas for future in-depth research. The workflow of this study is shown in Figure 1.

Data Acquisition
The expression data and corresponding clinical information of LUAD were downloaded from the TCGA database (https:// portal.gdc.cancer.gov/). Patients who met the following criteria were included in the study: (1) LUAD patients with RNAseq data; (2) patients with AS data; (3) patients with detailed clinicopathological and follow-up information, including sex, age, clinical stage, overall survival, and survival status. Preprocessing such as gene ID conversion and batch correction proofreading was performed in R software.
AS events can be divided into 7 basic types: alternate acceptor site AA), alternate donor site (AD), alternate promoter (AP), alternate terminator (AT), exon skipping (ES), exon mutual exclusion (ME), and retained intron (RI) [11]. Seven types of AS events of LUAD were downloaded from the TCGA SpliceSeq database (https://bioinformatics.mdanderson.org/TCGASpliceSeq/). Each AS event was quantified by the percentage splice (PSI) value, which ranged from 0 to 1 [11]. AS events that met the following criteria were included for the next step of analysis, with PSIvalue ³75%. The visualization of interaction of AS events was performed through UpSet plots using the UpSet package.
Further, the AS dataset was processed with survival information (survival status and overall survival) by removing blank samples and merging. To identify AS events associated with survival in LUAD, univariate Cox proportional risk regression analysis was performed. P<0.05 was selected as the significance level for candidate prognostic AS events. We created volcano plots using the ggplot package to show details of AS associated with survival. The UpSetR package creates UpSet diagrams to show the interaction of the 7 types of AS events. Bubble plots were created using the ggplot package to show prognosis-related AS events in each of the 7 types of AS, and P<0. 05 and P values ranked in the top 20 from smallest to largest are presented as conditions.

AS-Related Prognostic Risk Model
The prognosis-related AS genes identified above were further analyzed.
Minimum absolute shrinkage and selection operator (LASSO) and multivariate Cox proportional risk regression analysis were performed using the glmnet package and survival package to narrow down the candidate genes, calculate the coefficient, and build a prognostic risk model. The penalty parameter (l) was determined by the minimum parameter. The risk score was calculated using the following equation: Risk score=eSi(Coefi*PSIi) Coefi: correlation coefficient of the identified AS events, PSIi: PSI of the AS identified events.
The risk score in each sample was calculated and divided into a high-risk score subgroup (>median value) and a low-risk score subgroup (<median value). The results were constructed by visualizing distribution curves for risk score subgroup, survival distributions, and model gene expression heatmap. Kaplan-Meier analysis was used to construct survival curves using the survival package and SurvMiner package between 2 subgroups, and the log-rank test was performed. P<0.05 was considered statistically significant. The accuracy of the prognostic prediction model was assessed using the timeROC package to construct receiver operating characteristic (ROC) e934491-2 curves. An area under the curve (AUC) >0.06 was considered as acceptable prediction.

Association of the Risk Score with Clinical Characteristics
To determine the relationship with prognosis in LUAD, univariate Cox proportional risk regression analysis was performed for several factors, including age, sex, clinical stage (TMN), and risk score. To exclude the interference of confounding factors, multivariate Cox proportional risk regression analysis was performed for the above-mentioned multiple influencing factors. The above process was performed in the survival package, and P<0.05 was considered statistically significant.
The Limma package was used to analyze the differences between the risk score subgroups and several factors, including age, sex, and clinical stage (TMN). The Wilcoxon rank sum test was used to compare the 2 subgroups, and the Kruskal-Wallis test was used to compare 3 or more subgroups. P<0.05 was considered statistically significant.
A nomogram can be used as a multiple-indicator combination to diagnose or predict disease onset or progression [12]. The rms package was used to plot the column line graph constructs, and then the mean of the predicted survival probability and the corresponding actual survival probability (calculated by the Kaplan-Meier method) were calculated separately for each group of study subjects, and the 2 were combined to plot the 4 calibration points, and finally the 4 calibration points were connected to obtain the prediction calibration curve. If the prediction calibration curve is closer to the standard curve, the better predictive ability of the column line plot is indicated.

Correlation of Risk Score Subgroup with Immunity
The ESTIMATE package was used to assess tumor microenvironment (TME) while generating 4 measures (stromal score, immune score, estimate score, and tumor purity). The higher the score of immune cells or stromal cells, the more immune or stromal components in TME. The Wilcoxon rank sum test was used to compare the differences in stromal score, immune score, estimate score, and tumor purity between the high-risk score group and low-risk score group using the Limma package, and P<0.05 was considered statistically significant.
The CIBERSORT package was used to obtain the content of immune infiltrating cells in each LUAD sample, including 22 immune cells such as CD8T cells, activated dendritic cells, macrophages, natural killer T cells, and regulatory T cells. The Wilcoxon rank test was used to compare the difference between risk score subgroup and immune infiltrating cells. The correlation between risk score and significantly different immune cells was further compared by Spearman's method. P<0.05 was considered statistically significant.
Single-gene enrichment analysis (ssGSEA) was used to quantify the relative abundance of immune infiltrating cells and immune-related functions in each sample using the immune score. The samples were first analyzed by ssGSEA using the GSVA package to obtain an immune score in the LUAD dataset. Then, the Wilcoxon rank test was used to compare the differences in immune score between the risk score subgroups using the Limma package. P<0.05 was considered statistically significant.

Differential Genes Between Risk Score Subgroup
The Limma package was used to apply differential analysis to identify differentially-expressed genes (DEGs) between highand low-risk score groups, and DEGs with statistically significant differences were screened by logFC >1 and P<0.05. The Limma package was further used to compare the differential expression of DEGs in different clinical traits. The Wilcoxon rank sum test was used to compare 2 groups and the Kruskal-Wallis test was used to compare 3 or more groups. Samples were divided into a high-expression group (>median value) and a lowexpression group (<median value) according to the expression of DEGs. Survival curves were generated by Kaplan-Meier method for predictive analysis between groups, and the log-rank test was used to determine the significance of differences. P<0.05 was considered statistically significant.

Correlation Between Differentially-Expressed Genes And Immunity
The Wilcoxon test was applied using the Limma package to compare the association between DEGs and immunity, including immune checkpoint, tumor microenvironment, 22 immune infiltrating cells, and immune score. The preparation materials for the above process were the same as the previous steps, and P<0.05 was regarded as statistically significant.

Alternative Splicing-Alternative Splicing Factors Network
The expression of alternative slicing factors in transcriptional sampling data was extracted and processed uniformly. The correlation of survival-related AS and alternative splicing factors were analyzed using the Spearman method with P<0.001 and a coefficient of 0.6 as screening conditions. Cytoscape software was used for alternative splicing-alternative splicing factors network graph construction.

Statistical Analysis
All data were processed using R software (version 4.0.3). Data extraction was performed using PERL (version 5.8.3). Network plots were created using Cytoscape (version 3.7.2). P<0.05 was considered statistically significant.

Identification of AS-Related Events of LUAD
Seven AS events types were included in this study: AA, AD, AP, AT, ES, ME, and RI. A total of 513 cohorts of LUAD tissues and 59 cohorts of paraneoplastic tissues from the TCGA database were included in this study (Supplementary Table 1 (Figure 2A).  A significant proportion of these genes contained more than 2 AS events, and up to 7 different AS events occurred in a single gene. To identify the relationship between AS and survival, further univariate Cox proportional risk regression analysis was performed (Supplementary Table 2). We found that 19 45 AS events and the corresponding 1408 genes were significantly associated with the overall survival of LUAD patients ( Figure 2B). According to the priority ranking of AS events, ES events and AP events were the prognosis-related AS events that occupied a larger proportion. Most genes have only 1 AS event, and a few genes have more than 2 different splicing patterns. Volcano plots were used to identify AS prognosisrelated genes in LUAD (Figure 2C), and the top 20 AS events of these 7 types of events are demonstrated by bubble plots in Supplementary Figure 1.

Construction of Risk Model
The prognosis-related AS events were reduced by LASSO regression analysis (Supplementary Figure 2). Further, multivariate Cox regression analysis was performed to obtain a list of prognostic characteristics of AS events, and risk scores were calculated according to the weight of each event (Supplementary Table 3). As shown in Figure 3A, the number of patient deaths increased with higher risk score values. The heatmap in Figure 3B demonstrates the differential expression of identified AS events in the risk score subgroups. The identified AS events were unevenly distributed between high-and low-risk subgroups. Kaplan-Meier survival curve analysis showed that the high-risk score group was associated with poorer OS (P<0.001, Figure 3C). The ROC curve analyzed the relationship between risk model and survival at 1, 2, and 3 years, which were 0.771, 0.770, and 0.784, respectively, indicating that the risk prognostic model had high predictive ability ( Figure 3D).

Association of Risk Model with Clinical Characteristics
To better assess the correlation between the risk score and clinical characteristics, the independent variables were included in the Cox regression analysis, and the continuous variable was risk score, and the categorical variables were sex, age, and clinical stage (Stage). Then, univariate and multivariate Cox regression analyses were performed for the above variables. The results of the univariate Cox proportional risk regression analysis showed that risk score and clinical stage (Stage) were the risk factors affecting overall cancer survival (P<0.001, Figure 4A). The results of multivariate Cox regression model analysis showed that risk score and clinical stage (Stage) were independent factors influencing the prognosis of LUAD patients (P<0.05, Figure 4B). In conclusion, the risk score showed strong potential for clinical application.
Further analysis of differences in risk scores and clinical characteristics (age, sex, and clinical stage) was performed, and the risk score was significantly associated with sex, clinical stage (T), clinical stage (N), and clinical stage (Stage) (P<0.05, Figure 4C). The variables clinical stage (T), clinical stage (N), clinical stage (Stage), and risk score were included to construct a risk prediction model using column line plots ( Figure 4D). A calibration chart shows that the closer the prediction calibration curve is to the standard curve, the better the predictive ability of the column line plot is.

Association of Risk Model with Immunity
The ESTIMATE package was used to calculate the immune score, stromal score, estimate score, and tumor purity. Compared to the high-risk score group, the low-risk score group had more stromal and immune components in TME and lower tumor purity (P <0.05, Figure 5A).
The CIBERSPOT package was used to assess the contents of 22 immune infiltrating cells and further compare the relationship between the risk model and immune infiltrating cells, showing significant differences between the risk score subgroup in terms of M0 macrophages, M2 macrophages, and helper T cells (P<0.05, Figure 5B). Further correlation analysis showed that higher risk score values were associated with higher M2 macrophage content and with lower M0 macrophage and helper T lymphocyte content (P<0.05, Supplementary Figure 3).
Calculating immune scoring showed the risk score subgroups were significantly different in dendritic cells, B lymphocytes, major histocompatibility complex class I, type II interferon response, tumor-infiltrating lymphocyte, helper T lymphocytes, regulatory T cells, human leukocyte antigen, mast cells, neutrophils, and natural killer cells (P<0.05, Figure 5C).

Differential Gene Analysis
The expression of TTC39C, CDKN2A, and PKIB in paraneoplastic tissues was significantly lower than that in LUAD tissues, and the difference was statistically significant (P<0.001, Figure 6A). Patients with highly expressed PKIB and CDKN2A had shorter survival, while those with high expression of TTC39C had longer survival (P<0.001, Figure 6A). Further, PKIB was significantly different between N0 and N1 (P<0.05, Figure 6B), CDKN2A did not show significant differences in sex, age or clinical stage (Figure 6C), and TTC29C was significantly different in sex, T2, and T3 (P<0.05, Figure 6D).  T T1 T2 T3 T4   T  T1  T2  T3  T4   15  First, the relationship between DEGs and tumor microenvironment was assessed. Compared with the low-expression group, the immune cell score and stromal cell score were higher in the high-expression PKIB group (P<0.001, Figure 7A). CDKN2A was not significantly correlated with immune cell score and stromal score (Figure 7B). Compared with the high-expression group, the immune cell score and stromal score were higher in the low-expression TTC29C group (P<0.001, Figure 7C).

Construction of Alternative Splicing-Alternative Splicing Factors Network
A number of alternative splicing factors were included, such as ACIN1, AGGF1, ALYREF, AQR, ARGLU1, BAG2, BCAS1, and BCAS2. By performing correlation analysis between AS and alternative splicing factors, we found that alternative splicing factors DDX39B, SEC31B, and CLK1, and a series of AS exist in a network of interactions that together play a role in tumor development (Figure 9).

Discussion
In this study, an AS-related risk model was constructed using a series of prognosis-related AS signatures. First, survival-related AS signatures were screened by the LASSO and Cox regression model constructs, which were more stable and reliable compared with a single method. The AUCs of 0.771, 0.770, and 0.784 at 1, 2, and 3 years, respectively, were confirmed by internal tests, demonstrating the strong predictive ability of the AS-related risk model. Further, risk score and clinical staging (Stage) were found to be significantly associated as independent influences on the prognosis of LUAD patients. Clinical staging (TNM) is widely used to assess the prognosis of LUAD patients, but lacks individuality. Therefore, the construction of columnar maps to predict the prognosis of LUAD patients based on the patient's profile can be helpful for clinical prediction. Unfortunately, columnar maps have been shown to have good predictive power for LUAD, but due to their limitations, they cannot accurately project survival in LUAD and still need further refinement.
Immune infiltrating cells are associated with tumor growth, infiltration, and invasion [20]. Using the CIBERSPOT package to demonstrate the content of 22 immune infiltrating cells, it was found that the content of immune infiltrating cells was different in the AS risk score subgroups, with the low-risk score group having more M0 macrophages and helper T lymphocytes and fewer M2 macrophages. It was further found that the low-risk score subgroup was more associated with prognosis than was the high-risk score subgroup. The relevance of macrophage infiltration in the tumor microenvironment has long been known. According to the "macrophage homeostasis e934491-20 theory", there are be 2 types of macrophages with different functions in the tumor microenvironment: M1 type, the classically activated macrophages, and M2 type, the alternative activated macrophages. The M1 type is involved in antigen delivery and tumor cell killing by secreting various antigenic molecules and inflammatory factors and expressing MHC II and B7 molecules, while the M2 type can have the effect of promoting the pathological process of cancer cells. Helper T lymphocytes can fight against tumors by secreting cytokines (eg, IL-2, IFN-r, and TBF-B), killing tumor cells independently of CTL, and expressing CD40 ligands to stimulate CTL by binding to antigen-presenting cells CD40 [18]. Differences in immune infiltrating cells in high-and low-risk score groups are associated with the prognosis of LUAD patients, further suggesting that AS may influence the composition of the tumor microenvironment.
Immunotherapy is one of the more successful approaches in cancer treatment in recent years. In this study, the low-risk score subgroup had more immune checkpoints and a better response to immunotherapy compared to the high-risk score subgroup. The risk model possessed a good predictive power for immunotherapy. In addition, a noteworthy point is that CD276 was more highly expressed in the high-risk score subgroup. CD276 is overexpressed in tumors, including melanoma, leukemia, breast, prostate, and colorectal cancers, which is associated with poor prognosis and worse clinical outcome, and is a promising target for tumor immunotherapy and a new hope for immunotherapy-insensitive LUAD [21].
The core genes in the risk model were further screened by differential analysis. The tetratricopeptide repeat domain 39C (TTC39C) is associated with diabetic maculopathy with vision loss and skeletal muscle atrophy [22]. The mechanism by which TTC39C plays a role in cancer is unclear, but there is an association. Overexpression of TTC39C can lead to impaired ERK1/2 MAP kinase and MAP transduction signaling [23]. The relationship between the MAPK signaling pathway and a variety of tumors is becoming known with the discovery of multiple proteins in the MAP transduction signaling pathway. cAMP-dependent protein kinase inhibitor-B (PKIB), a protein kinase inhibitor, is a class of proteins that inhibit cAMP-dependent protein kinase activity [24][25][26]. One investigator found that PKIB expression is upregulated in NSCLC and acts as an important effector of NSCLC cell proliferation and metastasis through the PI3K/Akt pathway [27]. The cyclin-dependent kinase inhibitor 2A (CDKN2A), located on chromosome 9p21, encodes the p16 protein. p16 acts as a cell cycle regulator that inhibits cell cycle proteins, thereby preventing premature transition from G1 to S phase [28]. CDKN2A is associated with familial atypical melanoma multiforme syndrome and familial pancreatic cancer [28]. Further internal datasets validated that TTC39C, CDKN2A, and PKIB are associated with survival and may be prognostically relevant therapeutic targets for LUAD. In addition, the core genes and immune presence associations among risk models were assessed, showing that AS plays an important role in the tumor microenvironment.
Song et al used LASSO-Cox regression analysis to construct a multiple alternative splicing prognostic risk model of LUAD and a single alternative splicing prognostic risk model of LUAD, respectively, and found that the multiple alternative splicing prognostic risk models prognostic risk models possessed better predictive ability [29]. The risk model constructed for multiple alternative splicing events alone in this study also yielded satisfactory results. In addition, multiple quantities including risk scores were incorporated to construct and validate the reliability of column line graphs. Li et al constructed the alternative splicing risk prognostic model of LUAD by Cox risk regression analysis and also obtained more satisfactory results [30]. Different methods exist for constructing risk prognostic models, and the reliability of the risk prognostic models of alternative splicing-related LUAD constructed under different methods can be compared in the future. There are still shortcomings in this study, and the constructed risk model needs a larger sample to confirm the predictive ability in clinical practice. There is little research on the mechanism of action of the identified core genes in LUAD, and further studies are needed to confirm the role of core genes.

Conclusions
In this study, a dataset from public databases was used to build and evaluate AS risk models and further elucidate the relationship with prognosis of LUAD patients. In addition, an in-depth mechanistic analysis was performed to confirm the biological relationship between AS and immunity.