Development of an Immune-Related Risk Signature for Predicting Prognosis in Lung Squamous Cell Carcinoma

Lung squamous cell carcinoma (LSCC) is the most common subtype of non-small cell lung cancer. Immunotherapy has become an effective treatment in recent years, while patients showed different responses to the current treatment. It is vital to identify the potential immunogenomic signatures to predict patient’ prognosis. The expression profiles of LSCC patients with the clinical information were downloaded from TCGA database. Differentially expressed immune-related genes (IRGs) were extracted using edgeR algorithm, and functional enrichment analysis showed that these IRGs were primarily enriched in inflammatory- and immune-related processes. “Cytokine-cytokine receptor interaction” and “PI3K-AKT signaling pathway” were the most enriched KEGG pathways. 27 differentially expressed IRGs were significantly correlated with the overall survival (OS) of patients using univariate Cox regression analysis. A prognostic risk signature that comprises seven IRGs (GCCR, FGF8, CLEC4M, PTH, SLC10A2, NPPC, and FGF4) was developed with effective predictive performance by multivariable Cox stepwise regression analysis. Most importantly, the signature could be an independent prognostic predictor after adjusting for clinicopathological parameters, and also validated in two independent LSCC cohorts (GSE4573 and GSE17710). Potential molecular mechanisms and tumor immune landscape of these IRGs were investigated through computational biology. Analysis of tumor infiltrating lymphocytes and immune checkpoint molecules revealed distinct immune landscape in high- and low-risk group. The study was the first time to construct IRG-based immune signature in the recognition of disease progression and prognosis of LSCC patients.


INTRODUCTION
Lung cancer is the second most prevalent human malignancy that arises from epithelial cells in both men and women, and is by far the leading cause of cancer death worldwide, accounting for 25% of all cancer deaths. It's estimated that there are about 228,820 new diagnosed cases and 135,720 deaths from this disease in the United States in 2020 (Siegel et al., 2020). Lung cancers are classified into two main types, non-small cell lung carcinoma (NSCLC) and small cell lung carcinoma (SCLC). Lung squamous cell carcinoma (LSCC) is the second most common histologic type of NSCLC followed adenocarcinoma (LADC) (Fernandez-Cuesta and Foll, 2019), causing about 30% of lung cancers. LSCC is strongly correlated with tobacco smoking than other forms of NSCLC (Cipriano et al., 2019). Although the incidence and mortality continues to decrease partly due to less smoking and advances in early diagnosis and treatment, patients still have a poor prognosis and 5-year survival remains at a very low level (Lu et al., 2019;Thomas et al., 2015). The main treatments of lung cancer are surgical resection and chemotherapy, but vary for various factors, including tumor stage, lung function and genomic alterations. Early stage LSCC patients are typically receiving resected surgically and chemotherapy and/or radiation could be as an adjuvant therapy, while advanced LSCC are given the firstline systematic therapy, commonly a platinum-based regimens (Gandara et al., 2015).Compared to lung adenocarcinoma, much less frequent mutations in EGFR and ALK rearrangements for LSCC postponed the development of targeted therapies (Derman et al., 2015;Yang et al., 2019). Therefore, identification of novel and effective biomarkers will contribute to monitor the progression of LSCC and reduce the numbers of patients that not diagnosed before the onset of this aggressive disease. Many online tools have been developed for identification of biomarkers to assist the diangosis of cancer subtypes and provide survival prediction for cancers based on massive genomic data in recent years (Kamps et al., 2017). As an aggressive disease with leading mortality and incidence worldwide, several databases were constructed to find prognosis related biomarkers in lung cancer, while few online survival analysis softwares available for specific subtypes of lung cancer except than Kaplan-Meier plotter (Gyõrffy et al., 2013) and OSluca (Yan et al., 2020). Cancer immunotherapy has become a promising treatment for different types of human cancers in recent years. Some immunotherapies have been utilized to leverage the immune system to fight tumors (Kobold et al., 2018;Popovic et al., 2018). Immune checkpoint inhibitors (ICIs), such as anti-PD-1/PD-L1, have emerged as an effective therapeutic selection for advanced breast cancer, metastatic melanoma and NSCLC Ahern et al., 2019;Liu D. et al., 2019). NSCLC is characterized by several mutations in the immune system, making it possible that these patients may benefit from immunotherapy. Several monoclonal antibodies targeting the immune checkpoint pathways have been approved for treatment of NSCLC. Anti-PD-1 agent (Nivolumb) could improve the clinical outcome for LSCC. Despite the success of immunotherapy, different patients showed diverse responses, and only a small number of patients benefited from treatment with ICIs agents. The prognostic significance of PD-1 expression in patients with tumors are still controversial. This suggests that it is essential to identify additional predictive biomarkers for developing potential immunotherapy. Different studies proposed models regarding the prognostic utility of immune-related genes in various cancers, including papillary thyroid cancer , bladder cancer , head and neck squamous cell carcinoma (Cao et al., 2018) and nonsquamous NSCLC (Li et al., 2017a). However, the clinical relevance and prognostic significance of IRGs in LSCC has not been well illuminated.
This study aimed to investigate the prognostic utility of IRGs on monitoring the prognosis, and identify novel biomarkers in developing potential targeted therapies for LSCC patients. The gene expression profile and clinical information of LSCC cohort were downloaded from TCGA database. Differentially expressed IRGs were identified, and the prognostic landscape of these IRGs were comprehensively assessed through computational biology. Importantly, we proposed a prognostic signature based on the immunogenomic profile for predicting prognosis in LSCC patients. The study may provide new insights into understanding the functional regulatory mechanisms of IRGs, and help to develop immune-related targeted therapy in the treatment of LSCC in further in-depth work.

Gene Expression Datasets and IRGs
Level 3 RNA-seq raw count data from 502 lung squamous cell carcinoma patients (LSCC) and 49 non-tumor samples were downloaded from the TCGA database. Clinical information of these patients was also downloaded and extracted. Patients with overall survival (OS) time less than 30 days were removed, and finally 475 patients were employed in further analysis. Immunology Database and Analysis Portal (ImmPort) is a database that updates and shares immunology data accurately and timely (Bhattacharya et al., 2014). More importantly, the database provides a list of IRGs, which were identified to actively participated in the process of immune activity, and were classified into functional categories. A total of 2498 IRGs were derived from the ImmPort for this study (Supplementary Table S1).

Identification of the Differentially Expressed Genes (DEGs) and IRGs
The differentially expressed genes were identified between the LSCC and adjacent non-tumor samples by edgeR R package (Robinson et al., 2010). Genes with | log FC| > 2 and False discovery rate (FDR) < 0.01 were selected as the DEGs. The differentially expressed IRGs were extracted from the DEGs list.

Heatmap and Clustering Analysis
Heatmap and clustering were performed by the pheatmap R package.

Gene Functional Enrichment Analysis
Functional enrichment analysis of the DEGs and the differentially expressed IRGs were conducted using clusterProfiler R package (Yu et al., 2012) to identify significantly enriched GO terms, including biological process (BP), molecular function (MF), and cellular components (CC). The pathway analysis with reference from KyotoEncyclopedia of Genes and Genomes pathways (KEGG) was also performed. The p-value was adjusted by Benjamini and Hochberg method and less than 0.05 was considered as statistically significant.

Identification of Overall Survival (OS)-Related IRGs
A log2 (normalized value + 1) expression matrix was used to identify OS-related IRGs by univariate Cox regression analysis using survival R package. IRGs that significantly associated with the OS of patients were delivered to further functional analysis and construct the prognostic risk signature.

Molecular Characteristics of the Differentially Expressed IRGs
To explore the interplay between these OS-related IRGs, the protein-protein interaction (PPI) network was assessed from the STRING database (Franceschini et al., 2013). The PPI network was reconstructed by the Cytoscape software (Shannon et al., 2003).
In addition, transcription factors (TFs) have been known to directly mediate the expression levels of the respective genes. Cistrome Cancer database is an online resource that integrates cancer genomics data from TCGA with over 23,000 profiles of ChIP-seq and chromatin to provide the regulatory interactions between TFs and transcriptomes (Mei et al., 2017). To investigate the potential ability of TFs in regulating these clinically relevant IRGs, a total of 318 TFs were downloaded from Cistrome. The correlation between these IRGs with TFs was calculated, and the Person correlation coefficient greater than 0.3 was set as the cutoff value to construct the regulatory network of the IRGs and potential TFs.
The database of Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST) was used to identify the key regulated factors of OS-related IRGs. TRRUST is a reliable curated portal for human, and mouse transcriptional regulatory networks, which contains 8,444 TFs-target regulatory relationships of 800 human TFs ( Han et al., 2018).
Construction and Validation of the Immune Related Prognostic Risk Signature for LSCC OS-related IRGs were employed to construct the prognostic risk signature using multivariate Cox stepwise regression analysis. The minimum number of genes that comprised of the final signature was determined by the Akaike information criterion (AIC) criterion (Vrieze, 2012). The model discrimination performance was assessed by the receiver operating characteristic (ROC) curve using survivalROC R package. The patient's prognostic risk score was calculated based on the corresponding gene expression data multiplied by the Cox regression coefficient. Patients were divided into high-and low-risk group according to the median risk score. The predictive utility of the prognostic signature was evaluated by Kaplan-Meier curve and log rank test.
Subset analysis was conducted to see the utility of the risk score for OS prediction of patients in different clinical parameters set, including age, gender, tumor stage, TNM stage.
To validate the prognostic capability of the immune related risk signature, two independent LSCC cohorts with clinical information, including GSE4573 (n = 130) and GSE17710 (n = 56), were downloaded from the GEO database. The risk scores and clinical pathological characteristics for LSCC patients from TCGA and these two validation cohorts were summarized in Supplementary Table S2. Association of Risk Score With Tumor Immune Landscape Using CIBERSORT and TIMER Database The CIBERSORT was developed to accurately quantify the abundance of distinct cell types within a complex mixture of gene expression data using deconvolution algorithm (Newman et al., 2015). We performed CIBERSORT analysis to calculate the proportions of 22 immune cell subtypes, including seven T cell types, naïve and memory B cells, plasma cells and NK cells, of LSCC patients in high-and low-risk groups, the samples with p< 0.05 were selected for further analysis.
Tumor Immune Estimation Resource (TIMER) is an online database to estimate the abundances of 6 subtypes of tumor infiltrating immune cells in 32 cancer types from TCGA database, including B cells, CD4 T cells, CD8 T cells, macrophages, neutrophils, and dendritic cells via gene expression data. Immune infiltration levels of LSCC patients were calculated, and the Pearson correlation between the risk score and immune cell infiltration was analyzed. Additionally, the correlation of sevenmodel genes expression with the immune cell infiltration levels in LSCC patients was assessed (GraphPad Prism 8.3.0).

Identification of the Differentially Expressed IRGs
To delineate gene expression profiles between normal and lung squamous cell carcinoma, 6,678 DEGs were identified using the edgeR algorithm (Robinson et al., 2010). Among these DEGs, 4,905 genes were up-regulated and 1,703 genes were downregulated ( Figure 1A). A distinct gene expression pattern was observed in the normal and tumor cases ( Figure 1C). 250 IRGs were referenced from ImmPort, including 128 downregulated and 122 up-regulated ( Figure 1B). A similar gene expression difference that defines by IRGs was also apparent in normal and tumor groups ( Figure 1D). GO terms analysis showed that these DEGs were significantly enriched in epidermal cell differentiation, keratinization, extracellular matrix and cellcell injections (Supplementary Figures S1A-C). Cytokinecytokine receptor interaction and alcoholism ranked the top pathways (Supplementary Figure S1D). As to the IRGs, inflammatory processes were the most frequently implicated through functional enrichment analysis, such as "leukocyte migration, " "cell chemotaxis, " "receptor ligand activity, " "cytokine activity, " and "cytoplasmic vesicle lumen" (Figures 2A-C). "Cytokine-cytokine receptor interaction, " "PI3K-AKT signaling pathway, " and "MAPK signaling pathway" were most enriched pathways by the differentially expressed IRGs ( Figure 2D).

Correlation of Individual Differentially Expressed IRGs With OS
To determine the potential prognostic utility of individual differentially expressed IRGs for patients, the clinical information of LSCC patients was downloaded. 27 differentially expressed IRGs were found to be significantly associated with the OS of LSCC patients (P < 0.05, Table 1) using univariate Cox regression analysis. As expected, we found that these OSassociated IRGs were significantly involved in the similar GO terms and pathways that seen in the enrichment analysis using the differentially expressed IRGs, such as "epithelial cell proliferation, " "mesenchyme development, " and "regulation of ERK1 and ERK2 cascade" (Figure 3A). "Hormone activity, " "growth factor activity, " "cytokine receptor binding" (Figure 3B) were the most significant molecular functions. The "Rap1, Ras signaling pathway" and "MAPK signaling pathway" were the top enriched pathways ( Figure 3C).
PPI network analysis of these IRGs identified three modules that defined by the number of nodes greater than 10, and named as SIPR1, EDNRB, and FGFR4, respectively. These three module genes were most significantly increased expression in LSCC cases (Figure 4), and have been implicated in cancer cell proliferation and migration (Tanaka et al., 2014; Huang et al., 2015). In addition, a comprehensive exploration of the molecular characteristics of these OS-related IRGs found that amplification, deep deletion, and mRNA high expression were the three most commonly types of mutations ( Figure 5). FGF4 was the gene with the highest mutation frequency, and there were 12 genes with a mutation rate greater or equal to 5%.

Construction of Transcription Factors (TFs) Regulatory Network
To explore the potential molecular mechanisms corresponding to the clinical significance of the OS-related IRGs, the regulatory network of these genes with TFs was investigated. We examined the expression profile of 318 TFs, and found 50 TFs were differentially expressed between LSCC and normal samples (Figures 6A,B). Among these 50 TFs, 3 genes (TCF21, HNF1B, and SOX2) were significantly correlated with the OS of LSCC patients (Supplementary Table S3). A regulatory network based on the Pearson correlation between 27 OS-related IRGs and 50 differentially expressed TFs was constructed using Cytoscape software. A correlation score more than 0.3 was set as the cut-off value. The TFs based regulatory network illustrated the regulatory relationships among these IRGs ( Figure 6C). We identified the key regulated factors of OS-related IRGs using the TRRUST database. Seven key transcription factors (RELA, NFKB1, SP1, VDR, FOS, PPARG and STAT3) were found to be associated with the regulation of these IRGs ( Table 2).

Development and Validation of the Immune Related Prognostic Signature
To establish an optimal prognostic immune related signature to define the patients' risk, multivariate Cox stepwise regression analysis was performed that employed the 27 OS-related IRGs, Based on the patients' risk score calculated by the signature, the patients were divided into high-and low-risk groups according to the median value of risk score (Figures 7B,C), and three model genes (GCGR, FGF8 and NPPC) were significantly up-regulated in low-risk patients ( Figure 7D and Supplementary Figure S2), while two genes (SLC10A2 and CLEC4M) were markedly down-regulated (Supplementary Figures S2C,E). Additionally, five model genes were up-regulated in LSCC patients (Supplementary Figure S3), while the remaining two genes (SLC10A2 and CLEC4M) were observed to be increased expression in normal cases (Supplementary Figure S3C,E). Kaplan-Meier curve showed that patients in high-risk group have worse OS than that of patients in low-risk group (P < 0.0001, Figure 8A). The area under curve (AUC) value of the receiver operating characteristics (ROC) curve was 0.7 ( Figure 8B) for 5 years, and 3-year AUC was 0.67 (Figure 8C), suggesting the prognostic signature based on IRGs has moderate capacity for monitoring prognosis. Furthermore, in order to minimize potential over-fitting, we used the least absolute shrinkage and selection operator (Lasso) regression model to select the model genes from 27 IRGs. A signature that comprises of seven IRGs (BMP2, GCGR, PTH, SLC10A2, PPBP, FGF9, and AMH) was constructed using multivariate Cox stepwise regression analysis. The patients were divided into low-and high-risk groups according to the median risk score. The patients in high-risk group have significant shorter OS than that of patients in lowrisk group (P = 4.2e-03), while the AUC of ROC curve of prognostic utility of this risk signature for 3 and 5 years were 0.62 and 0.61, respectively. This indicated that our signature show better predictive performance for LSCC patients. In addition, we made an attempt to develop several immune related signatures by increasing or decreasing the number of OS-related IRGs. The predictive performance of our signature is superior to these signatures (data not shown).
To evaluate predictive capability of the signature for patients' prognosis independently, univariate Cox regression analysis found that age, pathological M and T stages, Tobacco smoking history, tumor stage, and patient's risk score were significantly associated with the OS of LSCC patients ( Table 3). Multivariate Cox regression analysis showed that the risk signature could serve as an independent predictor after adjusting for other clinical parameters, including age, gender, tumor stage, pathological stage, tobacco smoking history and cigarette exposure per day ( Table 3).
The proposed signature was tested using external LSCC cohorts from GSE4573 (n = 130) and GSE17710 (n = 56). The risk scores of patients were calculated. Patients were divided into high-and low-risk group according to the median risk score in both validation datasets. Patients in high-risk group have significant shorter OS than that of patients in low-risk group in GSE4573 (Figure 8D, P = 0.0405). Similar observation was seen in GSE17710 LSCC cohort (Figure 8E, P = 0.0168) both validation cohorts. This demonstrated that our signature has good ability of risk stratification for LSCC patients.
In addition, the tumor infiltration levels among the somatic copy number alternations category (deep deletion, arm-level deletion, diploid/normal, arm-level gain, and high amplification) for the model genes varies in each immune subset compared with the normal using two sided Wilcoxon rank sum test (Supplementary Figures S4A-G).

The Utility of the Prognostic Signature in OS Prediction for Patients With Different Clinicopathological Factors
The subset analysis was performed to determine the utility of our signature in predicting patient's OS in different clinicopathological parameters. According to the Kaplan-Meier analysis, the risk score has potential prognostic values for the different subsets of LSCC patients (Figures 9A-F), such as patients with pathological T2 and T3 stages, N0-1 stages, M0-1, early tumor clinical stages (stage 1/2), age greater than 60 years old, and male patients. Patients with high-risk patients did have significantly worse outcome than those of low-risk groups. Furthermore, the risk score in high-and low-risk groups with different clinical parameters show significant difference (Supplementary Figures S5A-D).

The Relevance of the Prognostic Signature and Tumor Immune Landscape
Much attention has been paid to develop anti-tumor immune therapies for lung cancers recently, and major advances have been made, especially for immune checkpoint blockade, such as anti-PD1 antibodies Nivolumab and Pembrolizumab (Tanaka et al., 2014;Sul et al., 2016), and anti-PD-L1 Atezolizumab (No Authors Listed, 2016). The expression of immune checkpoint molecules, involving cytotoxic T-lymphocyte-associated protein 4 (CTLA4), programmed cell death 1 ligand (PD-L1), lymphocyte activation gene-3 (LAG-3), and T cell immunoglobulin-3 (TIM-3) of LSCC patients stratified by the prognostic signature, showed that CTLA4 and TIM-3 were significantly increased expression in high-risk patients (Figures 10A-D), suggesting patients in high-risk group might have poor response to the targeted molecular immunotherapy.
Tumor infiltrating lymphocytes (TILs) have been proposed to play a vital role in regulating tumor immune microenvironment (TME), treatment response and clinical outcome. CIBERSORT was applied to calculate the proportions of 22 immune cells types in high-and low-risk LSCC patients. Low-risk patients had remarkably higher fraction of plasma cells, memory activated CD4 T cells, and follicular helper T cells, whereas CD4 memory resting T cells, monocytes cells, macrophages M2 and neutrophils were in high levels (Figure 11).
To see if the immunogenome accurately reflected the status of TME, we found that the infiltration level of five immune cell types (Figures 12A-F), including CD4 T cell, CD8 T cell, neutrophil, macrophage and dendritic cells, are significantly associated with patient' risk score. In addition, infiltration levels of macrophage and dendritic cells ranked the top among the immune cells of the seven model genes (Supplementary Figures S6A-G).

DISCUSSION
Squamous cell carcinoma is one of the most frequently diagnosed histologic subtype of NSCLC, representing 20-30% of all the NSCLC cases (Travis, 2020). The majority of patients with this disease were closely correlated with the history of cigarette smoking, although the condition is improved alongside cancer incidence. The significance of TME in cancer development, progression and treatment response has been attracted great attention (Chen and Mellman, 2017). Patients with LSCC have benefited from an expanding immunotherapies, such as programmed cell death-1 (PD-1)/programmed cell death ligand 1 (PD-L1) and cytotoxic T-lymphocyte-associated protein 4 (CTLA-4) inhibition, VEGFR inhibition, and targeted therapies matched to fibroblast growth factor receptor (FGFR) and PI3K-AKT (Derman et al., 2015). However, only a small number of patients showed effective remission to immune checkpoint blockade (Wei et al., 2018) or chimeric antigen receptor (CAR) T cell therapy (Mohanty et al., 2019). Identification of strong gene signature to trace the immune status of cancer patients would be vital to establish reliable immune prognostic biomarkers, and to enable stratification of patients into highand low-risk groups that might be beneficial or responsive to immunotherapy. Several gene signatures that could be representative of tumor immune status have been proposed, and demonstrated to have potential prognostic utility in some cancers, such as non-squamous NSCLC (Li et al., 2017a), hepatocellular carcinoma (Long et al., 2019), renal papillary cell carcinoma , and renal cancer (Geissler et al., 2015). Although insights into the roles of IRGs in tumor progression and immunotherapeutic have been seen in many types of cancers, a comprehensive and transcriptomewide profile investigation for their clinical significance and molecular mechanisms in LSCC are yet to be described. In the present study, we developed a prognostic risk signature that comprised of seven IRGs, which demonstrated to be a reliable predictor in identifying LSCC patients with an unfavorable prognosis. In addition, the clinical significance and potential tumor immune landscapes of these IRGs in LSCC patients were also well illuminated.
GO terms enrichment analysis indicated that these differentially expressed IRGs are mainly involved in the process of inflammatory responses, such as "leukocyte migration, " "cell chemotaxis, " "cytokine activity." "PI3K-AKT signaling pathway" and "cytokine-cytokine receptor interaction" were significantly enriched in the LSCC progression, which is line with previous studies reported that these biological processes and pathways play a crucial role in the proliferation, angiogenesis, immune responses and progression in various cancers (Pons-Tostivint et al., 2017;Zhang et al., 2017;. Cancer progression is associated with a pro-inflammatory environment (Lin and Karin, 2007;Lippitz, 2013). Dysregulation of cytokine interactions was involved in the pathogenesis of lung cancer (Zhang et al., 2016;Ozawa et al., 2019). High cytoplasmic RAP1 can increase cisplatin resistance of NSCLC combined with increased NF-κB activity (Xiao et al., 2017), and RAS-RAF-MEK-ERK signaling was a vital pathway that mediates ALK positive tumor cell survival in lung cancer (Hrustanovic and Bivona, 2016). ALK inhibitors, such as crizotinib and ceritinib, have been applied to treat the ALK positive subset of patients (Hrustanovic et al., 2015). These differentially expressed IRGs in our study provided the clues of densely infiltrated inflammatory microenvironment that occurs often in the initiation of cancer cells, and the correlations between activation of immune-related pathways and disease progression and treatment response. Univariate Cox regression analysis showed that 27 IRGs were significantly correlated with patient's OS, suggesting that these genes could be prognostic biomarkers for their potential predictive utility. Further exploration of these OS-related IRGs observed that amplification, deep deletion, and mRNA high expression were the dominantly molecular traits. For example, Bone morphogenetic protein-2 (BMP2) is overexpressed in the majority of human lung carcinomas independent of cell types (Langenfeld et al., 2005). Previous study indicated that the expression pattern induced by BMP2 in lung fibroblasts may predicts patients' prognosis in lung adenocarcinoma (Rajski et al., 2015). Consistent with previous evidence, interleukin-33 (IL33) and interleukin 1 receptor-like 1 (IL1RL1) were found to be increased in lung cancer and associated with disease clinical stage . IL-33/ST2 signaling pathway has been implicated in tumor-associated immune response and inflammatory disease of the lung (Casciaro et al., 2019), and IL-33 could significantly promote the migration and invasion of lung cancer cells through alpha serine/threonine-protein kinase (AKT) pathway activation . As an attractant, IL-33 can promote the recruitment of Th2-associated cytokines, which act as key mediators for the recruitment of neutrophils, monocytes, NK cells, dendritic cells or T lymphocytes in inflammatory conditions at the site of tumor (Brabcova et al., 2014). Some chemokines were increased in LSCC samples, such as CCL2, CCL15, CCL16, and these chemokines were implicated in the angiogenesis or angiostasis balance and promoted tumor infiltrating hematopoietic cells in the pathophysiology of NSCLC (Rivas-Fuentes et al., 2015). The interaction network among these IRGs and the impact of the mutations might shape the LSCC immune microenvironment. In turn, the TME mediated disease progression and the response to treatment.
To explore the potential molecular mechanisms, a TFmediated network was used to identify vital TFs that might regulate these IRGs. S1PR1, EDNRB and FGFR4 were predominantly signatures that represented the modules. They are tightly correlated with cancer development and progression (Tanaka et al., 2014;Lang and Teng, 2019;Liu Y. et al., 2019), and involved in the egress of T cells from lung tissue in tumor infiltrating lymphocytes conditions (Mackay et al., 2015). Seven TFs were found to be key regulators of these OS-related IRGs. Of these TFs, as a subunit of NF-κB complex, the phosphorylation of RELA was associated with disease progression, inflammatory regulation and various cancers through NF-κB signaling pathway (Lu and Yarbrough, 2015). NFKB1 exerted an inhibitory function in the tumorigenesis and progression of different types of cancers through alleviating the abnormal activation of the NF-κB signaling (Concetti and Wilson, 2018). NF-κB activity in lung cancer was significantly associated with T cell infiltration, suggesting this pathway may mediate immune surveillance and promote antitumor T cell response (Hopewell et al., 2013). Our results may provide some clues that these IRGs (FLT4, AMH, ICAM1, BMP2, and FGF8) regulated by NFKB1 could be potentially implicated in the immune regulation of LSCC. Growing evidence has indicated that STAT3 signaling was involved in carcinogenesis, and immunotherapy response through regulating cancer cell differentiation and proliferation in human squamous cell carcinoma (Zhou et al., 2015;Zhao et al., 2018). Growing evidence showed that ICAM-1, a specific ligand for leukocyte-function associated antigen-1, interleukin-1, tumor necrosis factor-1 and interferon-γ, participates in various inflammatory immune processes through stabilizing T cell receptor-mediated binding between antigen-presenting cells and T-lymphocytes, including cell-cell interaction and leukocyte transmigration (Kotteas et al., 2014). ICAM-1 blockade inhibited lung cancer cell invasion in vitro and tumor metastasis in vivo. Tumor necrosis factor-α (TNF-α) induced ICAM-1 can be suppressed by thalidomide administration through inhibition of NF-κB binding to the ICAM-1 promoter (Lin et al., 2006). This suggested targeting ICAM-1 could be a potentially effective therapy for LSCC patients. Previous studies provided limited direct information about these TF-mediated differentially expressed IRGs in LSCC development. Hence, the immunomodulatory role of these genes in monitoring LSCC progression and prognosis prediction remain to be fully investigated.    Recent years have witnessed a boom of gene signatures in clinical practice. Oncotype DX provided a breast cancer recurrence score based on the expression of 21 genes (Siow et al., 2018), and Colorant also represented a recurrence score tool for colon cancer patients constructed by the expression of 18 genes (Tan and Tan, 2011). These signatures developed based on gene expression profiles suggested that screening new prognostic cancer biomarkers is a promising approach to identify high-risk patients for subsequent clinical decisionmaking, such as therapy selection, outcome prediction and disease progression tracking. As regard to lung cancer, Wan et al. (2012), identified a smoking-associated 7-gene signature for predicting patient's diagnosis and prognosis in smokers. Similarly, the authors also developed a smoking-associated 6gene signature that predicts patient's risk and survival by modeling crosstalk with major lung cancer signaling pathways (Guo and Wan, 2012). Zhang et al. (2019), constructed a 5-microRNA signature based on serum microRNA profiling to predict survival for patients with advanced-stage NSCLC. However, there are few rigorously classifiers that use to predict prognosis for LSCC patients. Although 14-gene signature was developed for LSCC prognosis, the moderate large number of genes make it difficult for clinical use (Li et al., 2017b). Here, we developed an immune related prognostic risk signature with seven OS-related IRGs for LSCC patients.
Association analysis combining with the model risk score and clinicopathological characteristics indicated that a higher risk score was significantly correlated with advanced stage, older age and poor prognosis. This is in keeping with higher risk score representing an immunosuppressive tumor microenvironment. This, in turn, promotes tumor progression and worse response to therapies. The seven-gene signature remained an independent prognostic predictor after adjusting for clinicopathological variables. Considering that a robust prognostic signature could be risk-stratify indicator in other independent cohorts, we used two LSCC datasets to validate the predictive efficacy, and found that the signature performed well in classifying highand low-risk prognostic groups (P = 0.0405 for GSE4573, and P = 0.0168 for GSE17710). Additionally, our prognostic signature has a moderate high AUC using seven genes when compared to other signatures developed by greater than or less than seven IRGs, which make it more implementable in clinical practice.
Immune escape of tumor cells is a big bugbear to antitumor immune response during cancer progression (Dunn et al., 2002). This process was mediated by several immunosuppressive mechanisms, including increased immunosuppressive cells, and overexpression of immune checkpoint molecules (e.g., PD-1, CTLA-4, LAG3, and TIM-3) in the TME. In this study, the immune cell infiltration landscape in high-and low-risk FIGURE 11 | Relative infiltrating proportion of immune cells in high-and low-risk groups. Blue violin reflects high-risk groups, and red violin represents low-risk groups.
LSCC patients and the correlation between risk score with the expression of immune checkpoint genes, including PD-L1, CTLA-4, LAG3, and TIM-3 (Woo et al., 2012), were assessed. Consistent with TIM-3 was up-regulated in LSCC and positively correlated with malignant parameters . Our analysis indicated that CTLA-4 and TIM-3 expression were significantly higher in high-risk LSCC groups. Patients in the low-risk group possessed higher plasma cells, CD4 memory activated T cells, T follicular helper cells, macrophage M0 and M1 as well as other immune cells infiltration. Clinical studies demonstrated that tumor-infiltrating lymphocytes (TILs) have a major impact on the disease progression in various cancers (Al-Shibli et al., 2008;Granetto et al., 2017), and increased infiltration of TILs, such as cytotoxic T cells, memory T cells, and T helper cells 1, was associated with favorable prognosis (Bremnes et al., 2016). Previous reports revealed that different subsets of cells differentiating from CD4 + T-cells can promote the inactivation of CD8 T cells, the killing effect of NK cells, enhance the cytotoxicity of CD8 T cells in immune response (Pinto et al., 2018). Significant positive association of six types of immune cell infiltration with risk score indicated that highrisk LSCC patients tended to have more CD4, CD8, neutrophil, macrophages and dendritic cells than patients in low-risk group. This may mean the anti-tumor response in high T cell infiltration is neutralizing through immunosuppressive TME shaped by increased expression of immune checkpoint molecules. These results might help to explain poor prognosis of high-risk patients, and the seven-gene risk signature may provide the potential immunotherapeutic insight for immune checkpoint inhibitor therapies, while the underlying interactions of the IRGs and immune mediators in LSCC need further identification.
Among the seven genes, the expression level of four genes (GCGR, FGF8, CLEC4M, and SLC10A2) were associated with patients' OS using the median expression value as cutoff in OSluca database (Supplementary Figure S7). CLEC4M, SLC10A2 and FGF4 have been found to be involved in lung cancer progression and the regulation of treatment resistance. CLEC4M is associated with worse prognosis and inhibition of CLEC4M showed potential clinical relevance in counterbalancing cisplatin resistance in NSCLC patients (Tan et al., 2019). Fibroblast growth factor (FGF) isoforms, such as FGF4, FGF19 and FGF7, promoted epithelial-mesenchymal transition (EMT), cell proliferation and migration during cancer progression by inducing store-operated calcium entry in lung carcinoma (Qi et al., 2016). In addition, combination therapy of tyrosine kinase inhibitors, such as FGFR and CFS1R inhibitors, with anti-PD-1 or anti-CTLA-4 antibodies showed promising benefit for cancer patients (Katoh, 2016). As a protective factor, increased expression of SLC10A2 is closely related to suppression of NSCLC cell proliferation and migration, and promote apoptosis under bexarotene treatment (Ai et al., 2018). Consistent with evidence that reduction of glucagon receptor, GCGR, in papillary thyroid carcinoma resulted in the inactivation of EMT and P38/ERK pathways (Jiang et al., 2020), GCGR was downregulated in LSCC and may serve as a potential prognostic biomarker and therapeutic target for LSCC. NPPC and PTH have not been previously reported being related to lung cancer, PTH related protein mediated energy wasting in fat tissues, and neutralization of this protein can ameliorate cancer cachexia, which improves patient's survival (Kir et al., 2014). Our study is the first time to suggest their predictive potential as prognostic markers for LSCC patients.

CONCLUSION
In the present study, we developed a seven-gene risk signature for LSCC patients based on differentially expressed IRGs, which could serve as an independent prognostic predictor, and an indicator of tumor immune landscape. The patients divided by the risk score have markedly different prognoses. Correlation analysis of the risk score with clinicopathological factors showed that the signature has potential utility of estimating LSCC patients' prognosis and guiding clinical use. The data presented may provide insight into developing novel therapeutic strategies. Further work should concentrate on uncovering the molecular mechanisms of these IRGs in regulating LSCC progression and outcome, as well as validating our signature in clinical setting.

AUTHOR CONTRIBUTIONS
DF: conceptualization, design, and writing -review and editing. DF, BZ, and SH: data acquisition. DF, BZ, LY, SH, and WX: methodology. DF and BZ: data analysis, interpretation, and writing -original draft. DF and WX: project administration. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported by the Natural Science Foundation of Jiangxi Province (20192BAB215001, 20192BAB215011, and 2020BAB206067) and the National Natural Science Foundation of China (81360447).