The expression data, gene mutations, and clinical information were collected from the Cancer Genome Atlas (TCGA) website for 284 patients with LUAD (Tomczak et al., 2015). And 214 LUAD patients’ information were collected through dataset GSE31210 in the GEO database. Single-cell sequencing data (GSE176021) of tumor-infiltrating T lymphocytes from six NSCLC patients were obtained from the GEO database. For data inclusion criteria, we selected patient samples with RNA transcriptome sequencing data and complete clinical data. For data from different datasets de-batching effects were performed and normalized, and we used fragment per kilobase transcript/fragment per million mapping (FPKM) expression values for further analysis. In FPKM, RNA-seq data were normalized to the length of each gene and the total number of aligned reads in the library (Trapnell et al., 2010). FPKM values were transformed using log2 (FPKM + 1). The flowchart of our investigation was displayed in Figure 1.

The molecular signature database provided information on TEX signaling pathways and marker genes (Wherry, 2011; Liberzon et al., 2015). Similar to previous studies, we performed an unsupervised cluster analysis of LUAD patients using IFN-γ, TNF, and IL-2 signaling pathways to represent the TEX pathway (Zhang et al., 2022). The specific method is to use ssGSEA through the “GSVA” R package to estimate the activity score of each patient’s TEX pathway (Hänzelmann et al., 2013). The percentage of patients at different stages in different clusters is also shown.

To further explore the differences among different subgroups, we first drew Kaplan-Meier (K-M) survival curves for different clusters with the mark of 50% survival rate. Then the K-M survival curve of pairwise subgroups was drawn.

To further explore the abundance of immunocytes in different classifications, we used CIBERSORTX algorithm to evaluate 22 immunocytes in samples from different clusters (Steen et al., 2020). After cell infiltration of each sample was obtained, COX regression analysis was conducted to explore the prognostic value of various cells in each TEX cluster. Based on median immune cell content, LUAD was divided into diverse subgroups, and survival rates were compared between groups.

The correlation patterns among genomics can be described using a systems biology method known as a weighted correlation network analysis (WGCNA) (Langfelder and Horvath, 2008). The R package repository had the package WGCNA 3.6.1 that was used for the WGCNA. WGCNA analysis was performed after deunion of the four cluster differential genes. The significance of each gene was taken into account when calculating the association between the gene expression profile and the TEX score, and the relationship between module eigengenes and gene expression profiles was taken into consideration when determining module membership. The soft threshold parameters were set at a power of 4 and a scale-free R2 of 0.9, in order to ensure the topology network was scale-free despite the number of nodes. The analysis consisted of retrieving an initial set of six modules, and the Grey modules that showed the strongest correlations were applied for further investigation.

Using the survival random forest of 1,000 trees by the R package randomSurvivalForest version 3.6.4, it was possible to validate the results and rank the importance of 7 genes obtained from Lasso regression using the R package randomSurvivalForest Version 3.6.4 (Taylor, 2011). The relative importance of gene > 0.2 is considered the ultimate hub gene.

Gene expression tends to show significant collinearity between genes, which means it is necessary to use prognostic models if needed. LASSO regression and other methods reduce the number of variables to further reduce the redundancy of the model and increase the convenience of clinical use. Based on the following formula, we were able to calculate the risk score according to the following (Tibshirani, 1997): R i s k   s c o r e = ∑ i = 1 n Exp i * C o e f i

Here, the TCGA-LUAD data set was employed as the training set to construct the risk model of LUAD patients based on survival random forest screening genes, took OS as the outcome event, and p value less than 0.05 as the limit of statistical significance. We then categorized the patients into diverse subgroups based on the formula generated by the risk model. Meanwhile, the K-M survival curve was drawn for the high and low risk group. Receiver operating characteristic (ROC) curves have a wide range of uses in identifying the diagnostic power of threshold changes. To further analyze the predictive power of prognostic models, we plotted the ROC. An analysis of multivariate cox regression was undertaken in order to determine the independent prognostic significance of risk scores. To ensure the value of external generalization of the prognostic model, we used another LUAD dataset (GSE31210) for validation.

A follow-up analysis looked at genes with differential expression (DEGs) between high and low risk groups. It was established that the cut-off criteria for the study were |FC| > 2 and adj. p.val > 0.05. GSEA software was used to analyze three data sets of HALMARK, KEGG and GO for the biological functional differences among high-risk patients (Shi and Walker, 2007; Subramanian et al., 2007).

By using the GSVA analysis, it was possible to explore the differences between subgroups in signaling pathways for disease development (Hänzelmann et al., 2013). Moreover, a correlation analysis was conducted between the partial signal pathway score and the risk score.

A boxplot was made to illustrate the expression of 11 immune checkpoints between different subgroups. An analysis of 22 different types of immune cell infiltration was performed using the CIBERSORTX algorithm.

The detailed gene mutation statuses of the subgroups were displayed using the R package “maftools” to make comparisons between the two subgroups (Mayakonda et al., 2018). With the IMvigor210 package, we were able to determine gene expression and immunotherapeutic effectiveness in the IMvigor210 cohort (Mariathasan et al., 2018). IMvigor210 cohort is widely used to analyze the efficacy of immunotherapy.

It is a standard processing procedure that is used to do downstream processing on scRNA-seq data which is carried out using Seurat R software package, version 3.0.2, and a standard downstream processing package for this analysis (Stuart et al., 2019). In addition, genes detectable in fewer than 3 cells and genes detected in fewer than 200 cells were excluded, and the percentage of mitochondria detected was limited to no more than 20% of the total number of genes. Then, t Data was normalized using LogNormalize. A non-linear method used for reducing the dimensions of a sample is t-distributed stochastic neighborhood embedding (t-SNE) that is used for unsupervised clustering and unbiased visualization of cell populations on a two-dimensional map after principal component analysis (PCA) (Van Der Maaten and Hinton, 2008). A minimum fraction of 0.25 cell population fraction was used in both populations in order to identify marker genes in each cluster using the “FindAllMarkers” function. The filtering criterion was filter value of absolute log2 fold change (FC) ≥1. To visualize each marker gene’s expression patterns within the cluster, the “DotPlot” function in Seurat was used. Then, the SingleR package (version 1.0.0) was utilized for annotating cell types based on marker-based information (Aran et al., 2019).

The A549 and BEAS-2B cell lines were obtained from Dr Liu. A549 cell lines were cultured in RPMI‐1640 (Invitrogen) and BEAS-2B cells were cultured in DMEM medium, The medium was supplemented with 10% FBS (Gibco).

The expression of CCL20 in tumor and normal tissues of LUAD patients was compared through GEPIA2 online website, and we analyzed the overall survival rate of high expression group and low expression group (Tang et al., 2017).

The total RNA was extracted using the Trizol reagent. RT was performed with DNA-free total RNA in Revert Aid First Strand cDNA Synthesis Kit (Thermo). For PCR amplification, specific primers were used to amplify the transcribed cDNA. CCL20 Forward: ATG​TGC​TGT​ACC​AAG​AGT​TTG​C; CCL20 Reverse: CCA​ATT​CCA​TTC​CAG​AAA​AGC​C.

Integrated DNA technologies (Coralville, IA, United States) provided us the synthetic siRNA and the scrambled negative control siRNA. This experiment consisted of transfecting cells with LipofectamineTM RNAiMAX (Thermo Fisher Scientific, Massachusetts, United States) in opti-MEM according to the procedure given by the manufacturer.

The manufacturer’s instructions were followed when performing flow cytometry. Apoptosis was detected with the Annexin V-PE/7-ADD Apoptosis Detection Kit (Vazyme, A213-01). The B525 nm wavelength was selected for the Fluorescein (FITC) signal channel, and the B610 nm signal channel was selected for the ECDPE-TR (ECD) signal channel.

It was determined that two groups with normally distributed variables and those with variables that were not normally distributed were statistically significant using independent t-tests and Mann-Whitney U tests. In order to make a comparison between the two groups on the basis of differences between the groups, we conducted an analysis of variance (ANOVA) and a Kruskal–Wallis test (Hazra and Gogtay, 2016). We performed Spearman correlation and distance correlation analyses using the R package Hmisc 4.4.1. To analyze the correlation between the objects. Those objects whose coefficient was greater than 0.5 were considered highly correlated (Faul et al., 2007). For the purpose of identifying the prognostic factors, Cox regression analyses were conducted. A survival curve with all survivorship curves generated by the R package survminer was also used to determine the overall survival (OS) and TEX riskScore values before generating any survival curves with the R package survminer. As a means of plotting the heatmaps, the R package Complex Heatmap 2.4.3 was used. R package ggplot2 was used for visualizing data comparisons. There were two-sided statistical analyses conducted using R software, which was used for all statistical analyses.

The outcome of unsupervised cluster analysis were displayed in Figures 2A, B, where the best result was to classify the 284 LUAD patients into four TEX cluster (consensus matrix k = 4) according to the GSVA scores of the three IFNG/TNFA/IL-2 pathways. Figure 2C displayed the GSVA scores of the three pathways within the four clusters. Cluster B has the highest IFNG/TNF/IL-1 pathway score, and cluster D has the lowest score. The number of patients with stage1 in cluster b was larger than in cluster A, B and C. The number of patients with stage1 in cluster B was larger than that in cluster A, B and C (Figure 2D), suggesting that the overall survival of patients in cluster b may be better than that in other cluster populations. We plotted K-M survival curves for the overall survival of the four cluster populations, and there was no statistically significant difference (p = 0.116) in survival between the four clusters (Figure 2E). We plotted the K-M survival curves separately for cluster B and cluster A C D, and the results suggested that the overall survival rate of cluster B was higher than that of cluster A C D (Figures 2F-H; Supplementary Figure S1) (Cluster C vs. D, p = 0.042; cluster A vs. D, p = 0.025; cluster B vs. D, p = 0.073).

To further analyze the abundance of immune cells in different cluster, we applied the CIBERSORTX algorithm to evaluate the 22 immune cells in the samples in different clusters. After cell infiltration score was obtained for each sample, COX regression analysis was performed to explore the prognostic value of various cells in each TEX cluster. Mast cell resting was a protective factor in both TEXa and TEXc (Figures 2I, J).

For further analysis of the transcriptome differences between different cluster, we will contrast between different cluster differences in gene analysis, analysis of the standard is greater than or equal to | logFC | = 0.5, rectify the p value is less than 0.05, and mapped the volcano map is used to display the results of the analysis (Pearson correlation coefficient = 0.2, p value < 0.001) (Figures 2K–P).

WGCNA analysis results suggested that the grey module was most relevant to survivals related information, and the grey module was selected for subsequent analysis (Figures 3A–D). There were 36 genes chosen as hub genes in the Grey module since they had absolute values of module membership [MM] that were greater than 0.5 and absolute values of gene significance [GS] that were greater than 0.5 within the module (Supplementary Table S1). Variable selection based on minimum depth values above the threshold (0.001) and importance values above the threshold (0.2) yielded seven tentative (SOX9, CD109, CCL20, DUSP5, DKK1, TNS4, and LCAL1) candidate prognostic markers for LUAD. This suggests that these seven genes are most relevant to the prognosis of LUAD (Figures 3E–G).

A TEX risk model that includes five genes was constructed using lasso regression analysis. The formula for the risk score is as follows: risk socre = (0.2628*CD109 + 0.0464*CCL20 + 0.0163*DKK1 + 0.0359*TNS4 + 0.0348*TRIM29). The TCGA-LUAD patients were divided into high- and low-risk groups based on their risk scores. The K-M survival curve between high and low risk groups suggested that the high-risk group had worse overall survival (p < 0.001) (Figure 4A). The AUC values of the TEX risk model were 0.823 in the first year, 0.688 in the third year, and 0.619 in the fifth year (Figure 4B). Multivariate COX analysis showed that TEX Score was an independent prognostic factor (p < 0.05, Hazard Ratio :1.625 [1.329−1.986]) (Figure 4C). In the validation set GSE31210, we also found that high-risk LUAD patients had worse OS (Figure 4D). The AUC values of TEX risk model in the validation set were 0.643 in the first year, 0.655 in the third year, and 0.700 in the fifth year. These results TEX risk model have good predictive power (Figure 4E).

Based on the HALLMARK, KEGG, and GO datasets, we performed an enrichment analysis of biological functions in high-risk patients using GSEA software. The results showed that the five HALLMARK pathways with the highest enrichment were bile acid metabolism, heme metabolism, MYC target v1, peroxisome, and protein secretion (Figure 4F). The five most enriched pathways in the GO database were DNA conformational changes, negative regulation of cellular macromolecular biosynthetic processes, ribonucleoprotein complex biogenesis, mitochondrial matrix, and vacuolar membrane (Figure 4G). The top five enriched KEGG pathways were insulin signaling pathway, melanoma, peroxisome, T cell receptor signaling pathway, and vascular smooth muscle contraction (Figure 4H). The results showed that high-risk patients were highly associated with many tumor proliferation and metabolism-related pathways, suggesting that targeted TEX affects the prognosis of LUAD patients mainly through tumor proliferation and metabolic pathways.

We selected several gene sets for GSVA analysis based on the above GSEA results and found that TEX score was positively correlated with glycosaminoglycan degradation, linoleic acid metabolism, o glycan biosynthesis, leukocyte transendothelial migration, focal adhesion, ECM receptor interaction and p53 signaling pathway (Figures 5A, B). This suggests a potential pathway through which TEX exerts its effects.

Considering the great potential of TEX for immunotherapy, In both high and low risk groups, we plotted the expression levels of 11 immune checkpoints. In the high-risk group, CD44 expression was higher, which may be a therapeutic target in the future (Figure 5C). Based on 22 immune cell infiltrations, the high-risk group had a higher percentage of resting Dendritic cells and a lower percentage of activated Dendritic cells (Figure 5D). Further correlations of TEX risk scores and 22 immunocytes are shown in Figure 5E. The radar chart further showed the contents of 22 immunocytes in the high-risk group (Figure 5F), suggesting that TEX may affect the prognosis of LUAD by regulating the state of Dendritic cells.

We used a map to determine the landscape of gene mutations in high and low-risk subgroups of patients (Figures 6A, B). There was no statistically significant difference between the gene mutation frequencies between the groups that were analyzed, but it was noted that TP53 and TTN had the highest mutation frequencies.

Based on the TCGA-LUAD data set and the IMvigor210 data set, subgraph analyses were conducted to evaluate immunotherapy and chemotherapy in high-risk and low-risk groups. The high-risk group demonstrated a lower percentage of responders to immunotherapy (Figure 7A), and TEXscore was lower in those with low response (Figure 7B). These results suggest that patients with lower TEX risk scores are able to achieve a better immunotherapy response. Specific immunotherapy responses fall into four types: CR: complete response; PR, partial response; SD: stable disease; PD: progressive disease. There was no statistically significant difference in TEX risk scores among the four types of response, suggesting that the specific immunotherapy response was not related to the risk score (Figure 7C).

It was found that the high-risk group had also a poorer overall survival rate than the low-risk group when they received immunotherapy, regardless of the median risk score (Figure 7D). These results suggest that the TEX risk score has a role in predicting the efficacy of immunotherapy and the prognosis of patients receiving immunotherapy.

To further search for the potential therapeutic target-cell interactions of TEX, a total of 8 cell subtypes were identified in the single-cell sequencing dataset of T cells (B, CD4 Tconv, CD8T cell, CD8Tex, DC, NK, T prolif, Treg) (Figure 8A). The key gene in the risk model, CCL20, was most highly expressed on TEX cells (Figures 8B,C), CCL20, a key gene in the risk model, was most highly expressed on TEX cells, suggesting that CCL20 plays an important role in the TEX process in LUAD patients and is a potential therapeutic target. GSEA analysis showed that TEX cells were mainly enriched in cell adhesion (Figure 8D). The results of cell communication showed that TEX mainly interacted with CD8T cells (Figures 8E,F). These results provide new explanatory theories and therapeutic targets for TEX depletion in LUAD.

The GEPIA2 website contained 483LUAD patients and 347 normal lung tissues, and we found the expression level of CCL20 was higher in tumors tissues. Subsequently, the LUAD patients were classified into diverse subgroups based on CCL20 expression value (Figure 9A), and the results also showed that the high expression group had a shorter overall survival (p = 0.022) (Figure 9B). In addition, we detected CCL20 mRNA values in both normal and tumor cell lines. The results displayed that the mRNA expression level of CCL20 in A549 was more than twice that in BEAS-2B (Figure 9C). Subsequently, we knocked down CCL20 in A549 cells by siRNA (Figure 9D), and the CCL20 knockdown cells had more apoptosis than the control cells (Figure 9E).