Identification of an Innate Immune-Related Prognostic Signature in Early-Stage Lung Squamous Cell Carcinoma

Introduction

Lung cancer is the second most commonly diagnosed cancer (11.4%) and the leading cause of cancer-related deaths (18%) worldwide.¹ Lung cancer can be classified generally into two types: non-small cell lung cancer (NSCLC), which accounts for approximately 85% of tumors, and small cell lung cancer (SCLC), which accounts for approximately 15%.² About 30% of patients with NSCLC are diagnosed with lung squamous-cell carcinoma (LUSC),³ with these patients having a median overall survival (OS) approximately 30% shorter than in patients with other NSCLC subtypes.⁴ Poor prognosis in patients with LUSC may be attributed to older age, more advanced disease at diagnosis and a more central location predisposing to the invasion of larger blood vessels.^5,6 Additionally, many novel small-molecule inhibitors shown to provide clinical benefit to patients with lung adenocarcinoma are ineffective in patients with LUSC.⁷ Early-stage diagnosis and more effective treatment of LUSC can therefore enhance patient outcomes.

The innate immune system is the first line of defense against invading microbes. Pathogen recognition receptors on immune system cells activate innate immune signaling pathways, leading to the production of various proinflammatory cytokines and chemokines, subsequently inducing adaptive immune responses.^8,9 Innate immune responses also indirectly influence the tumor microenvironment (TME) by controlling T cell fate and critically sculpting the TME.¹⁰ Distinct signal regulation modes in the TME can alter the biological functions of certain innate immune system cells, which may lead to tumor proliferation and evasion of immune surveillance. Additionally, specific immunophenotypes of partial innate immune cell compositions depend on the type, pathological type, and even stage of malignancy. Therefore, understanding the characteristics of the innate immune landscape in specific tumor types and stages may provide a theoretical basis for subsequent monitoring and treatment.

Immune-checkpoint inhibitors (ICIs) have shown promise in the treatment of NSCLC, with some of these agents approved for first-line treatment. These include pembrolizumab monotherapy for advanced NSCLC patients with PD-L1 TPS ≥ 1% and without EGFR/ALK gene alterations, atezolizumab monotherapy for NSCLC patients with high PD-L1 expression and without EGFR/ALK genetic alterations, and combinations of pembrolizumab plus chemotherapy and nivolumab plus ipilimumab.^11–14 Some patients, however, do not benefit clinically from ICI treatment. For example, ICI treatment has not been shown to prolong survival in patients with early-stage LUSC, and some patients develop resistance to ICIs after an initial or develop immune-related adverse events. Owing to the complexity of the tumor immune microenvironment (TIME), no single indicator, including immunohistochemical expression of PD-L1, tumor mutational burden, number of CD8+ tumor-infiltrating lymphocytes or IFN-γ-related gene signature, is optimal for screening responders.^15,16 These drawbacks may be overcome by elucidating the detailed mechanisms of action of ICIs. Although all currently available ICIs target the adaptive immune system to activate T cells,¹⁷ immune sensing and escape in the TIME depend on interactions between the innate and adaptive immune systems. Metabolic changes in innate immune system cells can alter immune activity and antitumor effects.¹⁸ Alternatively, newly developed ICIs targeting the innate immune system could overcome resistance to currently approved ICIs.¹⁷ In addition, screening for specific gene signatures related to innate immunity may also enhance the effects of immunotherapy. Several signature models related to innate immunity, as determined by RNA-sequencing and microarray assays, have successfully predicted OS in patients with hepatoblastoma,¹⁹ head and neck squamous cell carcinoma,²⁰ ovarian cancer²¹ and colorectal cancer.²² Less is known, however, about the association between LUSC and innate immunity.

The present study systematically investigated the relationships between innate immune-related genes (IRGs) and prognosis in patients with early-stage LUSC. Our findings showed that innate IRGs can be used to classify patients with early-stage LUSC based on their clinical and molecular features. Furthermore, this study assessed differences in the immune microenvironment in high- and low-risk patients, providing a potential theoretical basis for future immunotherapy.

Materials and Methods

Mining from GEO and TCGA Databases

Cohorts I (GSE157009) and II (GSE157010) contain large numbers of specimens and datasets from multiple institutions, based on assays by independent laboratories and locked test parameters, thus ensuring objective quality standards.²³ Microarray Suite version 5.0 (MAS 5.0), containing data on normalized mRNA expression and relevant clinical information (including age, gender, T stage, OS and disease-free survival [DFS]), were acquired from the Gene Expression Omnibus (GEOhttps://www.ncbi.nlm.nih.gov/geo/). Samples were included if (1) gene expression data were available; (2) clinical information, including age, gender and T stage, were available; and (3) follow-up time was greater than 30 days. The training cohort consisted of 245 patients in GSE157009, and the validation cohort consisted of 232 patients in GSE157010. Transcriptome profiles in FPKM format were derived from The Cancer Genome Atlas (TCGA) database. Based on the above-described criteria, 438 patients with early stage (ranging from IA-IIIA) LUSC were selected as a second validation cohort.

Establishment of a Prognostic Signature of Innate IRGs

Univariate Cox regression analysis was performed to identify possible prognostic genes associated with OS in the training cohort, with a P-value < 0.001 considered the filtering criterion. A comparison of these genes with innate IRGs obtained from InnateDB was used to select overlapping genes, which were considered candidate genes for the construction of a prognostic signature. The least absolute shrinkage and selection operator (LASSO) method was then applied using the glmnet package (https://CRAN.R-project.org/package=glmnet) in R²⁴ for feature selection to narrow the number of IRGs with non-zero coefficients. Genes with non-zero coefficients were included in multivariate Cox regression analysis, and genes with P-values < 0.05 included in a prognostic signature model. The prognostic risk score of each sample was calculated by the “predict” function in “stats” R package using the formula:

Risk score=h0 (t)*exp (β1X1+β2X2+β3X3+β4X4+β5X5+β6X6),

where h0(t) represents the baseline risk value (unknown constants), β represents the regression coefficient calculated by multivariate Cox regression analysis, and X represents the gene expression level.

Evaluation and Validation of Innate IRG Signatures for Predicting Survival

Patients in the training cohort were divided into low- and high-risk groups based on the median risk score. Kaplan–Meier (K-M) survival curves and time-dependent receptor operating characteristic (time-ROC) curves were used to evaluate the accuracy of this signature. Moreover, the above-described formula and statistical methods were applied to the independent validation cohorts to validate the prognostic capacity of the gene signature.

Construction and Evaluation of a Predictive Nomogram

Risk scores were dichotomized into high- and low-risk score groups based on median risk score and treated as an independent clinicopathological parameter. Univariate and multivariate Cox regression analyses were performed to determine the clinical characteristics (age, gender, T stage and risk score) associated with OS. Factors differing significantly were visualized using a nomogram and their predictive accuracy and clinical efficacy determined. Decision curve analysis (DCA) was performed to estimate their possible clinical applications,²⁵ and ROC curves were used to measure the accuracy of the model.

Functional Enrichment Analysis

Log fold-changes (logFCs) in the levels of gene expression in the high- and low-risk groups were calculated using the “limma” package in R, with these genes sorted from the highest to lowest logFCs. Gene set enrichment analysis (GSEA) was performed using the “clusterProfiler” package in R,²⁶ with results considered significant if |NES| was >1, NOM p-value was <0.05, and FDR q-value was <0.25. Function and pathway enrichment were visualized using Ridgeplot.

Immune Profile Analysis

TIME differences between the high- and low-risk groups in the training and validation cohorts were evaluated using the “estimate” package in R. The relative proportions of 22 tumor-infiltrating immune cells (TICs) in the training cohort were calculated using the CIBERSORT algorithm.²⁷

Statistical Analysis

Univariate and multivariate Cox regression analyses were performed using the “survminer” package in R, and LASSO analysis was performed using the “glmnet” package in R. Kaplan–Meier analyses were performed using the “survival” and “survminer” packages in R, and ROC analysis was performed using the “timeROC” package in R. The “rms”, “ggplot2”, “ggforest”, “cowplot”, “gseaplot2”, “pheatmap” and “VennDiagram” packages in R (version 4.0.3) were used for visualization. A p-value <0.05 was considered statistically significant.

Results

Cohort Characteristics

The flow chart of the present study is shown in Figure 1A. After screening for eligibility criteria, the training cohort consisted of 245 patients in GSE157009, and the independent validation cohort consisted of 232 patients in GSE157010, with a second validation cohort from TCGA consisting of 438 patients. The 245 patients in the training cohort included 88 women and 157 men, of median age 70 years (range: 43–92 years). The GEO validation cohort of 232 patients included 81 women and 151 men, of median age 68 years (range: 46–89 years), and the TCGA validation cohort of 438 patients included 116 women and 322 men, of median age 68 years (range: 39–90 years).

Figure 1 Process used to construct the prognostic signature model based on innate immune-related genes. (A) Work flow. (B) Overlap of 147 genes associated with patient prognosis and 1376 genes associated with innate immunity, followed by selection of 29 genes by LASSO Cox regression analysis. (C and D) Construction of a LASSO Cox regression model from the 29 innate immune-related prognostic genes, followed by calculation of the tuning parameter (λ) based on partial likelihood deviance with 10-fold cross-validation. The vertical black line indicates the optimal log λ value.

Establishment of Innate IRG Signatures Model in Training Cohort

Univariate Cox regression analyses identified 147 significantly prognostic genes (P < 0.001) in the training cohort. Of the 1376 innate IRGs downloaded from the innateDB database, 29 overlapped with genes in the training cohort (Figure 1B). LASSO-penalized Cox analysis identified 18 genes based on the optimal value of λ (Figure 1C and D). Finally, multivariate Cox regression analysis identified six genes (SREBF2, GP2, BMX, NR1H4, DDX41, GOPC) that were prognostically significant (P < 0.05). These genes were used to build a prognostic signature and calculate corresponding regression coefficients. Evaluation of risk scores in each sample yielded the formula:

The median risk score served as the boundary that separated patients in the training cohort into high- and low-risk groups (Figure 2A), with a scatter diagram showing poorer survival outcomes in the high-risk group (Figure 2B). Kaplan–Meier analysis confirmed that OS was significantly shorter in the high- than in the low-risk group (Figure 2C, P < 0.0001). A heatmap showed model genes differentially expressed in the high- and low-risk groups of the training cohort (Figure 2D). Evaluation of the predictive performance of the model using time-dependent ROC curves showed that the areas under the curve (AUCs) at 1, 3, and 5 years were 0.713, 0.722, and 0.708, respectively (Figure 2E).

Figure 2 Risk score distribution, Kaplan-Meier survival analysis and time-dependent ROC curves in the training cohort (GSE157009). (A and B) Distribution of risk scores and survival status. (C) Kaplan-Meier analysis, showing that overall survival (OS) rates were significantly higher in the low- than in the high-risk group (P <0.0001 by Log rank test). (D) Distribution of risk scores as a function of gene expression characteristics. (E) AUC of time-dependent ROC curves measuring the ability of the model to predict OS.

Validation of the Prognostic Model in an Independent Cohort

Two independent validation cohorts were utilized to further investigate the prognostic value of the innate IRGs signatures model. The 232 samples in the GEO cohort GSE157010 were categorized into high- and low-risk groups based on the median risk score, calculated as described above. The risk curve (Figure 3A) and the scatter diagram for survival (Figure 3B) were similar to those observed in the training cohort. Kaplan–Meier analysis showed that OS was significantly shorter in the high-risk group (P = 0.00021; Figure 3C). A heatmap showed the expression of the six genes in the high- and low-risk groups (Figure 3D), and time-dependent ROC curves showed that the average AUCs at 1, 3, and 5 years were 0.661, 0.643, and 0.640, respectively (Figure 3E).

Figure 3 Risk score distribution, Kaplan-Meier survival analysis and time-dependent ROC curves in the validation cohort (GSE157010). (A and B) Distribution of risk scores and survival status. (C) Kaplan-Meier analysis, showing that overall survival (OS) rates were significantly higher in the low- than in the high-risk group (P= 0.00021 by Log rank test). (D) Distribution of risk scores as a function of gene expression characteristics. (E) AUC of time-dependent ROC curves measuring the ability of the model to predict OS.

The same classification methods were applied to the TCGA cohort. Kaplan–Meier analysis showed that OS was markedly shorter in the high-risk group (Supplementary Figure 1A), and time-dependent ROC analysis showed that the average AUCs in 3, 5, and 7 years, were 0.482, 0.580 and 0.634, respectively (Supplementary Figure 1B).

Independent Prognostic Value of the IRG Signature Model

The independence and integration effects of the IRG signature model among multiple clinical parameters were assessed by drawing Forest plots and a nomogram, followed by univariate and multivariate Cox regression analysis. Univariate Cox regression analyses showed that age (hazard ratio [HR]=1.480, p = 0.019), T-stage (HR = 1.740, P < 0.001), and risk score (HR = 2.250, p < 0.001) were significantly associated with OS in the training cohort. Subsequent multivariate Cox regression analysis showed that age (HR = 1.560, p = 0.008), T-stage (HR = 1.600, P < 0.001) and risk score (HR = 2.060, p < 0.001) remained significantly associated with OS (Figure 4A and B). In the validation cohort, univariate Cox regression analysis showed that T-stage (HR = 1.450, p = 0.011) and risk score (HR = 2.010, p < 0.001), but not age, were significantly associated with OS. Multivariate Cox regression analysis also showed that T-stage (HR = 1.440, p = 0.014) and risk score (HR = 1.930, p < 0.001) were significantly associated with OS (Figure 4C and D). Taken together, these findings showed that both risk score and T-stage were independently predictive of OS in the training cohort. These factors were incorporated into a nomogram to quantify the predicted individual survival probability at 1, 3 and 5 years (Figure 5A). The possible application of this model to future clinical decision-making was further estimated by DCA, which quantifies the clinical value of a nomogram by analyzing the clinical results at different periods. This analysis showed that the combined model was more predictive of outcomes than any single model at 1, 3 and 5 years (Figure 5B), and ROC analyses showed that the average AUCs at 1, 3, and 5 years were 0.752, 0.727, and 0.715, respectively (Figure 5C).

Figure 4 Univariate and multivariate Cox regression analysis of the risk scores and other clinical characteristics. (A) Univariate Cox regression analysis in the training cohort. (B) Multivariate Cox regression analysis in the training cohort. (C) Univariate Cox regression analysis in the validation cohort. (D) Multivariate Cox regression analysis in the validation cohort.

Figure 5 Construction and validation of a predictive nomogram. (A) Nomogram predicting the 1-, 3-, and 5-year OS of patients with early stage LUSC. (B) DCA curves showing that the nomogram is superior to a single model in providing optimal clinical decision-making benefits. (C) Time–ROC curves evaluating the predictive efficiency of the combined model at 1-, 3-, and 5-year.

Functional Analysis in the Training and Validation Cohorts

To identify the biological functions and pathways associated with risk scores, samples in the high- and low-risk groups were subjected to differential analyses, and total logFC values were sorted from highest to lowest and included in GSEA. The top 20 biological functions and pathways enriched in both cohorts were visualized by Ridgeplot (adjusted P < 0.05, Figure 6A–D). After integrating, we found that carcinoma-related pathways, such as cell cycle, cellular senescence and microRNAs, were up-regulated in cancers in the high-risk group of the training cohort, while immune and cell death-related pathways, such as cytokine–cytokine receptor interactions, chemokine signaling pathways and phagosome pathways, were down-regulated (Figure 6E). Similar results were observed in the validation cohort, with cell cycle, thermogenesis and chemical carcinogenesis-reactive oxygen species pathways up-regulated in the high-risk group, while apoptosis, cytokine–cytokine receptor interactions and chemokine and NOD-like receptor signaling pathways being down-regulated (Figure 6F).

Figure 6 Enrichment analyses of biological functions and pathways in the high- and low-risk group. (A and B) Ridgeplots showing biological function enrichment in the training and validation cohorts, with the horizontal coordinate corresponding to the peak being the enrichment fraction. (C and D) Pathway enrichment in the training and validation cohorts. (E and F) Gene sets with |NES| > 1, NOM p-value <0.05, and FDR q-value < 0.25 were considered significantly enriched.

TIME Score and Enrichment Analysis of 22 TICs

GSEA results showed that the differential functions and pathways between the high- and low-risk groups mainly included cytokine–cytokine receptor interactions and the chemokine signaling and carcinoma-related pathways. Evaluation of the contents of the immune-stromal components in TME of the training cohort showed that the high-risk group had a higher TumorPurity score but lower StromaScore and ImmuneScore than the low-risk group (Figure 7A), with similar results in the validation cohort (Figure 7B). Evaluation of the proportions of tumor-infiltrating immune subpopulations using the CIBERSORT algorithm showed the relative distributions of 22 TICs in the training cohort (Figure 7C), as well as correlations among these 22 TICs (Figure 7D). A violin plot (Figure 7E) showed that the numbers of naïve B cells, memory activated CD 4 T cells, and gamma delta T cells were markedly higher in the low- than in the high-risk group, whereas the numbers of memory B cells, resting NK cells and activated mast cells were significantly lower in the low-risk than in the high-risk group.

Figure 7 Component analysis of the tumor immune microenvironment. (A and B) Evaluation of tumor purity, Stroma score and Immune score by the “estimate” package in R and visualization by boxplots in the (A) training and (B) validation cohorts. (C) Relative content distribution of 22 TICs. (D) Analysis of correlations among the 22 TICs. (E) Violin plot showing the ratio of each of the 22 TICs between the high- and low-risk groups. (p < 0.05 was significant difference by Wilcoxon rank sum tests).

Validation of the Six-Gene Signature by Evaluating Expression and Mutations

The expression patterns of the six gene signature in LUSC and normal lung tissue were further confirmed using the Human Protein Atlas (HPA; https://www.proteinatlas.org/). Immunohistochemical results from the HPA database showed that levels of expression of SREBF2, DDX41, GOPC and GP2 were higher to varying degrees in LUSC than in normal lung tissue samples, whereas the expression of BMX did not differ significantly (Figure 8A) and NR1H4 was not included in the HPA. The levels of mRNAs encoded by these six signature genes were downloaded from the Gene Expression Profiling Interactive Analysis (GEPIA) database. Except for BMX mRNA, the level of which was markedly lower in LUSC samples, there were no significant differences in expression of the other five genes (Figure 8B). Evaluation of mutations in these six signature genes, including missense, splice and truncating mutations, in 1176 samples from three independent TCGA cohorts using the cBioPortal instrument (http://www.cbioportal.org/) showed that GP2 and NR1H4 were mutated in 5% of these samples (Figure 8C).

Figure 8 Differential expression of the six-genes in LUSC and normal lung tissue and gene mutations in LUSC. (A) Representative immunohistochemical results of LUSC and normal tissue samples from the Human Protein Atlas (HPA) assessing levels of expression of signature genes. (B) Levels of expression of mRNAs encoded by the six signature genes from the GEPIA database, with expression from the TCGA-LUSC and standardized by log2 (TPM+1). (C) Type and frequency of mutations of signature genes in three TCGA cohorts from the cBioPortal website.

Discussion

The present study identified the prognostic innate immune-related signature associated with OS in early-stage LUSC using one cohort (GSE157009) and validated this signature in a second cohort (GSE157010). LASSO and subsequent multivariate Cox regression analyses identified a signature composed of six innate IRGs (SREBF2, GP2, BMX, NR1H4, DDX41, GOPC) as being significantly associated with OS in early-stage LUSC. SREBF2,^28,29 BMX,^30,31 GP2,³² DDX41^33–35 and NR1H4^36–39 have been reported to correlate positively and GOPC negatively with tumor-related phenotypes and adverse events in several types of cancer.⁴⁰ Bone marrow X kinase (BMX) was associated with chemoresistance in SCLC through its modulation of autophagy.³⁰ Moreover, BMX signaling was found to mediate radiation resistance in the vascular endothelium of lung tumors in mice.³¹ A mutation in the gene encoding nuclear receptor subfamily 1 group H member 4 (NR1H4) was found to correlate with distant metastasis of LUSC.³⁶ An SNP, rs12296850 at 12q23.1, may be associated with the expression of NR1H4 and high risk in LUSC.³⁷ Indeed, this study found that the NR1H4 mutation frequency in LUSCs was 5%. NR1H4 has been reported to participate in colorectal cancer tumorigenesis, cell proliferation and cell survival via Wnt/β-catenin and the regulation of c-MYC stability.^38,39 Additional studies are needed to assess the potential role of NR1H4 in LUSC development.

After integrating these six innate IRG genes into a prognostic signature, we calculated the risk score of each sample. Univariate and multivariate Cox regression analyses showed that both risk score and T stage were independently prognostic of OS. DCA and time-ROC curve showed that risk alone had a better capacity than T-stage alone in estimating 3- and 5-year OS in early-stage LUSC.

GSEA of both the training and validation cohorts showed that the high-risk group was enriched in genes associated with tumor-related pathways, such as the cell cycle, microRNAs in cancer and thermogenesis, in agreement with the evaluation of tumor purity in the TME. Moreover, genes associated with cytokine–cytokine receptor interactions and chemokine signal pathways were significantly negatively correlated with risk scores in both the training and validation cohorts. Specific cytokines, such as IFN, IL-2, IL-4, IL13, and IL15,⁴¹ can modulate both innate and adaptive immunity by interacting with their receptors on various immune cells, with chemokines playing crucial roles in inflammation and immunity. Owing to their expression in both tumor cells and leukocytes in the TME, chemokines represent an ideal target for immunotherapy.⁴² These functional analysis results were in accordance with immune microenvironment scores, with the low-risk group having a higher proportion of immune cells (quantified as ImmuneScore) and stromal cells (quantified as StromaScore). The TIME has been confirmed important during tumor development and progression. LUSCs have the highest leukocyte fractions among all solid tumors analyzed,^16,43 making further exploration of the TIME immunity.

Although the adaptive immune system is primarily responsible for cytotoxicity and tumor suppression, various cytokines secreted by innate immune components within the TME can activate or suppress the immune system. For example, IL-4 in the TME initiates STAT6 signaling, polarizing macrophages to immune suppressive M2 macrophages;¹⁰ Neutrophils promote angiogenesis, tumor progression, and invasion by secreting neutrophil elastase, oncostatin-M and PGE2^44,45 factors that can modulate the curative effects of ICIs. Innate immune components and immunophenotypes in the TIME can change dynamically during tumor development, therefore, suggesting that efforts to develop therapeutic approaches should consider the specific conditions of innate immune components in different types and stages of cancer.

Genes representative of specific immune-cell types can predict TME components. The evolutionary trajectory of the immune response, from precancerous lesions to pre-invasive stages of LUSC⁴⁶ has been found to include increases in the numbers of activated T cells, total neutrophils, M1 macrophages, memory B cells and the myeloid signature, accompanied by decreases in the abundance of naive B cells. In addition, resting mast cells were found to be more abundant in early than in later developmental stages, whereas activated mast cells increased during development, with both resting and activated NK cells remaining at low levels. The present study found that some innate immune components differed significantly in the high- and low-risk groups, as estimated by the CIBERSORT algorithm-based TIC analysis. Specifically, the numbers of γδ T cells were higher, while the numbers of resting NK cells and activated mast cells were lower. γδ T cells are considered a rapid lymphoid stress-surveillance system, rapidly responding to stress-induced tissue perturbation, and a bridge between innate and adaptive immunity.^47,48 Moreover, increased numbers of γδ T cells have been found to correlated with prolonged OS in patients with LUSC.¹¹ NK cells have potent antitumor and antimetastatic activity. Although NK cell infiltration did not correlate with clinical outcomes in NSCLC,⁴⁹ activated NK cells in the TIME were found to correlate with OS in patients with LUSC.¹¹ These findings may be due to NK cell dysfunction cells in the TIME,⁵⁰ a hypothesis based on the higher ratio of resting NK cells in high-risk patients. Mast cells in the TME have both pro- and anti-tumorigenic activities, depending on tumor type and developmental stage.⁵¹ Mast cells infiltrating into tumor islets enhance the expression of TNFα, prolonging survival time in patients with NSCLC.⁵² In contrast, our results demonstrated that activated mast cells were more abundant in the high-risk group, suggesting a need for further studies to determine whether activated mast cells have pro- or anti-tumorigenic activities in early-stage LUSC. Despite the innate IRG signature identified in this study, components of adaptive immunity also differed in the high- and low-risk groups. For example, levels of activated memory CD4+ T cells and naïve B cells were higher, while levels of memory B cells were lower. These phenomena may be due to crosstalk between innate and adaptive immunity, which could determine the immune states in the TME and affect patient prognosis.

The present study had several limitations. First, clinical information about patients in these two cohorts in the GEO database was not complete. Although all samples were from patients with stage IA-IIIA LUSC (AJCC 8th),²³ accurate clinical staging for each patient was unavailable. Thus, T stage was included for independent prognostic analysis. Second, further studies in larger numbers of samples are needed to validate the present results. Third, the functions of IRGs in early-stage LUSC have not been fully determined.

Conclusion

Based on bioinformatics analysis, a prognostic model that included six innate immune-related genes was constructed to predict OS in patients LUSC. The efficacy of the model was further validated in two independent cohorts. Risk scores, as determined by this model, were negatively associated with immune cell infiltration and functional performance. These findings may improve understanding of the roles of innate immune-related genes in early stages of LUSC and to further explore dynamic changes in the TIME.