The Gene Expression Omnibus (GEO) database(ncbi.nlm.nih.gov/geo/) contains scRNA-seq data for two purified LUSC samples (GSM3304009 and GSM3304010). The Cancer Genome Atlas (TCGA) (portal.gdc.cancer.gov/) provides bulk RNA-seq data for patients with LUSC, covering 502 tumor cases and 51 normal cases (dataset name: TCGA-LUSC). Following the elimination of cancer patients lacking survival information, a total of 488 cancer patients were incorporated into the study. In addition, the GSE37745, GSE73403 and GSE74777 datasets were sourced from the GEO. These datasets were combined into a single set and adjusted for batch variations using the ‘ComBat’ function in the ‘sva’ package (Version: 3.54.0), with all data transformed to log₂ values (10).

The present study used the R packages ‘Seurat’, ‘Singler’ and ‘magrittr’ to analyze scRNA-seq data (11). Raw arrays for each sample were filtered, excluding genes expressed in fewer than three individual cells, cells with fewer than 50 genes, and cells with more than 5% of genes encoded by mitochondria. First, the scRNA-seq dataset was normalized using the ‘NormalizeData’ property in the Seurat R bundle, and the statistics were transformed into standardized scRNA-seq objects. The feature ‘FindVariableFeatures’ was used to detect the first 1,500 highly variable genes. Second, a principal component analysis (PCA) was performed. Previously, a JackStraw evaluation was used to screen for essential PCs), and 20 PC cells were selected for cluster assessment at P<0.05. The FindNeighbors property was used to compute the closeness of clusters, whilst the FindClusters property was used to analyze mobile clusters. The mobile clusters were subsequently validated with RunTSNE, and the FindAllMarkers property was used for each cluster to compute the differentially expressed genes (DEGs), for which logFoldChange=1 and the expression ratio of the least differential genes=0.25. Subsequently, each cell cluster was annotated with the ‘CreateSinglerObject’ property in the ‘Singler’ package. Finally, using the ‘monocle’ package in R, a quasitemporal analysis was performed (12). Based on the temporal gene expression of each cell, each cell was arranged according to a quasitemporal pattern, and this model was split once into multiple cell groups (states) based on gene expression status to create an intuitive lineage improvement tree.

Immune-related genes were retrieved from two comprehensive databases(Immport gene list and InnateDB Innate Immunity Genes), Immport (https://www.immport.org/home) and InnateDB (https://www.innatedb.ca/). By evaluating the expression of epithelial marker genes filtered in the preceding step, 278 EIGs were identified. To further assess the clinical relevance of these identified EIGs, univariate Cox regression analysis was performed. Genes exhibiting P<0.05 were deemed prognostic indicators. The protein interaction community was then evaluated using the protein-protein interaction network and the hyper-linkages between the DEGs was consequently determined. To improve the precision of the predictive gene screening, an extensive evaluation using both least absolute shrinkage and selection operator (LASSO) and Cox was performed. This analysis was executed using the ‘glmnet’ software package, which is tailored for fitting generalized linear models via penalized maximum likelihood. The LASSO model was used to evaluate survival according to Cox, and the ‘cv. glmnet’ capability was subsequently used to verify the first-class model (13). LASSO and Cox regression analysis were applied to elucidate genes with significant prognostic value. The differentially expressed EIGs coefficients are presented in Table SI. This technique enabled the systematic analysis of gene expression and risk factors, ultimately leading to the creation of predictive models. Patients were stratified into distinct risk categories using a risk score threshold of 0.91. The basic information of patients in the test group and the training group are presented in Table I. To assess the predictive accuracy of the EIGs, the area under the curve (AUC) was calculated using the ‘survival ROC’ package. A nomogram for 1-, 3- and 5-year outcomes was developed based on age, sex, tumor stage and risk model with the ‘RMS’ package. Decision analysis was performed using the ‘ggDCA’ package. The ‘survminer’ package and Kaplan-Meier methods were used to perform survival analysis. To mitigate the impact of confounding factors on survival curves, the two-stage (TS) test was used to address potential issues of curve crossover in survival analysis (14). Finally, survival analysis using an independent dataset sourced from the GEO (accession nos. GSE37745, GSE73403 and GSE74777; ncbi.nlm.nih.gov/geo) corroborated the predictive capability of the model regarding survival outcomes.

Detailed and in-depth gene set enrichment analysis (GSEA) was performed on the gene sets ‘c5.go.Hs.symbols.gmt’ and ‘c2.cp.kegg.Hs.symbols.gmt’ obtained from the MSigDB database (docs.gsea-msigdb.org/#MSigDB/MSigDB_FAQ) (15). Furthermore, using the ‘Cluster Profiler’ R package, Gene Ontology (GO) enrichment analysis was performed to assess the cellular and molecular characteristics of EIGs, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was used to evaluate their roles in several metabolic and signaling pathways (16). The GO annotation was based on the Bioconductor Project's complete genome annotation package (org.Hs.eg.db). KEGG elucidated the advanced functions, networks and applications of biological systems, providing a foundation for further research. The ‘cluster profiler’ tool was utilized, which uses the Web API to query the most up-to-date KEGG database, obtain path data and perform functional analysis. The threshold of P<0.05 was applied to determine enrichment. Distinct GO and KEGG plots were generated, illustrating variations in gene numbers and selection criteria.

Using the ‘biolinks’ software tool, LUSC mutation information from TCGA database was analyzed to assess the disparities in TMB between risk cohorts. A waterfall plot was subsequently generated using the ‘maftools’ tool to display the top 20 genes with notable mutation differences. This visualization enabled the comparison of the distinct variations in mutation burden among different risk cohorts. Subsequently, the expression data were processed using the ESTIMATE algorithm (1.0.13, bioinformatics.mdanderson.org/), leading to the assessment of the stromal score, immune score and overall ESTIMATE score in LUSC (17). CIBERSORT(0.1.0, http://cibersortx.stanford.edu/) was subsequently used to identify 22 invasive immune cells. Subsequently, immune typing analysis was performed using the R package ‘RColorBrewer’ (18). single-sample GSEA (ssGSEA) was then used to measure immune function in different risk groups. Furthermore, TIDE score files were retrieved from the TIDE site (tide.dfci.harvard.edu/) to compare the different scores among different risk groups.

An analysis of the expression of 79 immune checkpoint genes was performed across different groups to assess tumor immunity, and the ‘ggpubr’ software suite was used to scrutinize disparities in immunosuppressive molecules associated with immune checkpoints (19). Furthermore, gene transcription data from 298 patients with uroepithelial carcinoma was retrieved from the IMvigor210 cohort and their responsiveness to immunotherapy was forecasted to appraise the immunotherapeutic efficacy of the risk model in the present study and to prognosticate the outcome of the IMvigor210 immunotherapy regimen (20).

The drug susceptibility of LUSC was predicted using the R-editing language ‘parallel’ and the ‘oncoPredict’ package. Moreover, using the ‘ggplot2’ package, the susceptibility of the high-risk group to 198 drugs was plotted, considering P<0.05 (21).

To assess the proposed model, RT-qPCR analysis was performed. RNA was extracted from human LUSC H2195, H711 and H1522 cells (Binhui Biotechnology) and from normal bronchial epithelium 16HBE cells (Binhui Biotechnology) using TRIzol reagent (Takara Bio, Inc.). RT-qPCR was performed using the Probe One Step RT-qPCR Kit (NM_004048, Beijing Quality Biotechnology Co., Ltd.). For RT-PCR, begin with a 70°C incubation for 10 min, then synthesize cDNA at 37°C for 65 min, hold at 15.8°C for 150 min, and finally maintain at 4°C. For qPCR, use SYBR Green I (Vazyme Medical Co., Ltd.) as a fluorescent marker. Thermocycling conditions were as follows: Initial denaturation at 95°C denaturation for 2–10 min, followed by 40 cycles of 95°C for 10 sec and 60°C for 30 sec. Finally, the 2^−ΔΔCq method (22) was used to normalize sample gene expression for the final result calculations. Normalization involved using ΔCT and the efficiency measured for the reference gene to calculate the ETΔCT/ERΔCT ratio (23). To compare the expression levels between LUSC cells and normal cells, the unpaired t-test was used for statistical analysis. The aforementioned method enabled precise measurement and comparison of significant gene expression variations across several cell types. The primers used in the RT-qPCR experiment are listed in Table SII.

To assess the protein expression differences in patients with LUSC and healthy individuals, a thorough analysis was performed using data obtained from the Human Protein Atlas (HPA) database (AREG, MUC1 and FABP6; proteinatlas.org/). This comprehensive dataset was used to assess the protein expression profiles of the chosen EIGs in both LUSC and normal tissues and the analytical method enabled the comprehensive understanding of the expression patterns and variations of these genes across different samples.

Statistical analyses were carried out using R version 4.2.1, and the corresponding R packages were also based on this version. A P-value threshold of less than 0.05 was used to determine significance, except where explicitly stated otherwise.

As illustrated in Fig. 1, which presents the study flowchart, scRNA-seq was performed on a sequence of 2,470 LUSC samples. Furthermore, Fig. 2A demonstrates the range, sequencing depth and % of mitochondria for each sample. After a quality control filtration process was used to remove poor-quality cells, the remaining cells were subjected to subsequent analysis. After sample normalization, 7,950 low-mutation genes were excluded and 1,500 highly variable genes were selected (Fig. 2B). Subsequently, the PCA method was applied to simplify the variables within the dataset (Fig. 2C and D). Moreover, to perform a more detailed and in-depth analysis, 20 PCs) were selected, each exhibiting P<0.05, ensuring their statistical significance. Subsequently, the t-SNE algorithm (2.0, geeksforgeeks.org) was used, which allowed the visual mapping and distinguishing of the 12 clusters present within the sample, thus providing a clearer representation of the data distribution and relationships (Fig. 2E). Using the ‘single-R’ package in R (bioconductor.org/packages/release/bioc/html/SingleR.html), it was demonstrated that Clusters 0, 2, 3, 7 and 10 were classified as epithelial cell subpopulations. The expression patterns of distinct genes in individual cell clusters are shown in Fig. 2F. Finally, the developmental trajectory of epithelial cells was assessed using pseudotime analysis (Fig. 2G).

After the overlap of the immune gene set with genes related to the epithelium was confirmed, 278 immune-related genes were identified (Fig. S1A and Table SIII). A single-variable Cox regression analysis was performed using the TCGA LUSC cohort as a training set, and 24 epithelium-related genes (Table SIV) were revealed to be significantly related to overall survival (OS; Fig. S1B). Furthermore, a protein-protein interaction network was constructed using Cytoscape (3.9 http://cn.string-db.org/), and the relationships among the cox genes were demonstrated (Fig. S1C). A total of 8 EIGs were then identified using LASSO analysis, and 4 genes [amphiregulin (AREG), mucin (MUC)1, fatty acid-binding protein 6 (FABP6) and thymic stromal lymphopoietin (TSLP)] were selected for model construction after comparison with Cox regression analysis (Fig. S1D and E). The risk score was calculated as follows: (0.134 × AREG expression) + (0.129 × MUC1 expression) + (0.249 × FABP6 expression) + (−0.199 × TSLP expression). This was followed by an evaluation of the prognostic value of the model. According to the risk score, clinical information of different risk groups was recorded (Table II). Survival rates notably differed between the two risk groups, as indicated by the survival analysis curves, with high-risk patients exhibiting a significantly reduced survival rate, compared with that of the low-risk patients (P=0.007; Fig. 3A). However, the generated curves intersected at the tail end S(t)=0.2. The TS test demonstrated that there was still a significant difference in survival between different risk groups (P=0.00723), indicating that the intersection of the curves did not affect the results. ROC curves were also used to evaluate the prognostic ability of the model. Fig. 3B presents the AUC values at 1, 3 and 5 years, with OS rates for patients with LUAD of 0.698, 0.707 and 0.629, respectively. Similar results were demonstrated for both the experimental group and the entire cohort (Fig. 3C-F), confirming the feasibility of the model. Specifically, for the test cohort, it was revealed that individuals in the high-risk category had a significantly decreased OS rate, and this trend was consistent across all the experimental results. Additionally, through analysis of the expression levels of EIGs, it was revealed that AREG, MUC1 and FABP6 were expressed at higher levels in high-risk patients, indicating increased susceptibility to disease, whilst TSLP was more highly expressed in low-risk patients, suggesting greater resistance (Fig. S2A-C). Moreover, scatter plot analysis indicated that OS was markedly reduced and that the mortality rate was increased (Fig. S2D-I). Kaplan-Meier curves for survival analysis and differential expression of patients with LUSC patients with four model-associated EIGs are shown in Fig. S3. AREG, FABP6, and MUC1 showed decreased expression in the high-risk group compared to the low-risk group, while their expression was elevated in the latter. Conversely, TSLP exhibited lower expression in the low-risk group but was relatively higher in the high-risk group. The survival curves for the AREG and MUC1 genes intersected; therefore, the TS test was performed. The resulting P-values were 0.0314459 and 0.02161114 for AREG and MUC1, respectively. This indicated that despite the intersection, the genes still demonstrated a significant difference in their prognostic implications for survival outcomes.

Univariate and multivariate Cox regression analyses revealed that the risk score held independent prognostic significance in both types of analyses (Fig. S4A and B). Moreover, the clinical circle and heatmap indicated that age and T stage may be independent prognostic factors (Fig. S4C and D). Furthermore, the risk models value was notably greater than that of clinical traits in predicting patient prognosis, according to the nomogram (Fig. 4A). The ROC survival curve also demonstrated that the risk was more predictive than the majority of clinical data (Fig. 4B). Additionally, the decision curve revealed that the risk model had a markedly greater net benefit than the other clinical factors (Fig. 4C). The calibration curve also demonstrated that having a model had a notably greater C-index than having no risk model (Fig. 4D and E). Subsequently, a prospective estimate of patients with LUSC was performed according to age, risk, sex and disease stage (Fig. S5A-C) and survival rates were compared among patients of different ages, sexes, T stages and N stages (Fig. S5D-H) to determine the predictive effect. Notably, there were pronounced differences in risk scores among those aged ≥65 and above, as well as between patients with stage I and stage II tumors. The survival curves plotted by combining clinical features such as age, sex, and tumor stage with the model-predicted risk scores show significant differences between the groups. GEO external data was validated and it was demonstrated that the model is feasible (Fig. S5I).

The GSEA enrichment results indicated that the model gene had a positive association with the ‘granulocytes migration’ and ‘leucocytes migration’, while it showed a negative association with ‘chromosome segregation’ and ‘DNA replication’. Pathway enrichment revealed a high association of the model genes with three pathways related to ‘cytokine-cytokine receptor interaction’, ‘type I diabetes mellitus’ and ‘complement and coagulation cascades’, whilst the association with the ‘DNA replication’ pathway was low (Fig. 5A-D). GO analysis revealed the involvement of immune system genes, T cells, cytokines and white blood cells (Fig. 5E and F). Moreover, KEGG analysis confirmed the close relationship between these genes and the way in which antigen are processed and presented (Fig. 5G and H). According to KEGG analysis, antigen processing and presentation pathways were enriched, and the high-risk group revealed contamination with human T-cell leukemia virus 1. The detailed findings are presented in Table SV. Finally, the potential mechanisms initiated by relevant model genes during the progression of tumors were illustrated (Fig. 6).

The tumor mutation load was assessed in different risk populations. It was demonstrated that the most common mutations in the population were titin and MUC16 (Fig. 7A and B). When comparing TMB across different risk groups, no significant difference was observed (Fig. 7C). Consequently, the outcomes of patients with high and low TMB levels were analyzed. A more favorable prognosis was revealed in patients with a lower TMB than in those with a higher TMB (Fig. 7D). Subsequently, after incorporating the TMB-based risk model, high-risk and low TMB subgroups demonstrated a significantly improved prognosis compared with that of the TMB and low-risk subgroups (Fig. 7E). Although TMB did not show significant differences when comparing different risk groups individually, and even though survival was better in the H-TMB group than in the L-TMB group, after integrating TMB with the high- and low-risk groups defined by our model, survival differences still persisted between different risk and TMB groups (P<0.001). This shows the feasibility of the model in the present study. Subsequently, differences between TIDE and EIGs were assessed. A significant difference in TIDE scores was noted between the two risk groups (Fig. 7F). Due to the significant involvement of epithelial cells in the immune response and migration of tumors, the infiltration of these cells along with immune cells was assessed using data from patients with LUSC. The high-risk cohort exhibited significant variances in matrix scores, immune scores and estimated scores (Fig. 7G). Furthermore, the two risk populations demonstrated markedly different immune types (P=0.033; Fig. 7H). CIBERSORT analysis revealed that in high-risk cases, the proportion of activated dendritic cells, neutrophils and M2 macrophages was markedly higher when contrasted with the levels found in other immune cell populations. (Fig. 7I). Finally, the ssGSEA algorithm revealed a markedly high level of resting dendritic cells and neutrophils in the high-risk population, in comparison with the low-risk population (Fig. 7J).

To further assess the role of epithelial cells as genetic predictors of immunotherapy efficacy, 210 patients in the IMvigor210 cohort were analyzed. Survival analysis revealed significantly lower survival rates after immunotherapy in high-risk patients, compared with low-risk patients (Fig. S6A). Subsequently, the IMvigor210 immunotherapy analysis demonstrated that both high- and low-risk groups exhibited poor sensitivity to immunotherapy (Fig. S6B and C); therefore, immune checkpoints were searched for again. Immune checkpoint genes, which are critical for bypassing autoreactivity, represent new targets for cancer therapy (18). Immune checkpoint proteins poliovirus replication cell adhesion molecule and tumor necrosis factor receptor Superfamily Member 14 (TNFRSF14) exhibit significant differences across various risk groups, and their expression levels correlate with risk scores. (Fig. S6D-G). Significant disparities were observed in the correlations between risk scores and a selection of immune checkpoint proteins, such as cytotoxic T-lymphocyte-associated protein 4, Inducible T Cell Costimulator (ICOS), and CD20, etc. (Fig. S7).

The present study also evaluated differences in susceptibility between patients at low risk and those at high risk. Reduced susceptibility to anticancer drugs, including gefitinib, savolitinib, ipatasertib, erlotinib, uprosertib, afuresertib, alpelisib and nilotinib (Fig. 8A-H), was demonstrated in patients at low risk of cancer. The pharmaceuticals exhibiting heightened sensitivity within the low-risk cohort are listed in Table SVI. High-risk patients, on the other hand, had lower susceptibility to anticancer agents such as dasatinib, fludarabine, leflunomide, and selumetinib (Fig. 8I-L). Table SVII lists the drugs that show greater sensitivity in the low-risk group. These results indicate that EIGs may be good predictors of the efficacy of anticancer drugs (Table SVIII).

Images of immunohistochemical tissue sections of LUSC and normal tissues obtained from the HPA database demonstrated the differential expression of the selected EIGs genes in both tissue types (Fig. S8A). Analyzed using immunohistochemistry techniques, these sections revealed marked expression differences between cancerous and normal tissues. Furthermore, the present study used RT-qPCR to determine the levels of proteins in LUSC epithelial cells (AREG, FABP6, MUC1 and TSLP). The results revealed a significantly higher expression of AREG and MUC1 in normal tissues compared with that in tumor tissues, whereas there was a significantly higher expression of TSLP and FABP6 in tumor tissues compared with that in normal tissues (Fig. S8B). Overall, the RT-qPCR results were consistent with the aforementioned bioinformatics results.