A Four-gene Autophagy-related Prognostic Signature and Its Association With Immune Landscape in Lung Squamous Cell Carcinoma

Background: Recent advances in immune checkpoint inhibitors (ICIs) have dramatically changed the therapeutic strategy against lung squamous cell carcinoma (LUSC). In the era of immunotherapy, effective biomarkers to better predict outcomes and inform treatment decisions for patients diagnosed with LUSC are urgently needed. We hypothesized that immune contexture of LUSC is potentially dictated by tumor intrinsic events, such as autophagy. Thus, we attempted to construct an autophagy-related risk signature and examine its prediction value for immune phenotype in LUSC. Method: The expression prole of LUSC was obtained from the cancer genome atlas (TCGA) database and the prole of autophagy-related genes (ARGs) was extracted. The survival ‐ related ARGs (sARGs) was screened out through survival analyses. Random forest was performed to select the sARGs and construct a prognostic risk signature based on these sARGs. The signature was further validated by receiver operating characteristic (ROC) analysis and Cox regression. GEO dataset was used as an independent testing dataset. Patients were divided into high-risk and low-risk group based on the risk score. Then, gene set enrichment analysis (GSEA) was conducted between the two groups. The Single-Sample GSEA (ssGSEA) was introduced to quantify the relative inltration of immune cells. The correlations between risk score and several main immune checkpoints were examined. And the ESTIMATE algorithm was used to calculate the estimate/immune/stromal scores of the LUSC. Results: Four ARGs (CFLAR, RGS19, PINK1 and CTSD) with the most signicant prognostic values were enrolled to construct the risk signature. Patients in high-risk group had better prognosis than the low-risk group (P < 0.0001 in TCGA; P < 0.01 in GEO) and considered as an independent prognosis factor. We also found that high-risk group indicated an immune-suppression status and had higher levels of inltrating regulatory T cells and macrophages, which are correlated with worse outcome. Besides, risk score showed a signicantly positive correlation with the expression of PD-1 and CTLA4, as well as estimate score and immune score. Conclusion: This risk signature, which could be applied as an independent prognostic indicator for LUSC patients. And the autophagy-related scores are bound up with immune phenotype of LUSC, with higher score indicating an immune-suppression status. This study provides a novel and comprehensive sight to the correlation of autophagy and immune landscape in the tumor microenvironment of LUSC.

tumor immune microenvironment [4]. More effective biomarkers for prognosis and optimal treatment selection are urgently needed. Therefore, it is vital to nd more effective and precise methods for outcome predicting in LUSC.
Autophagy, an intracellular lysosomal degradation pathway that supports nutrient recycling and metabolic adaptation, is supposed to protect cells and tissues from stressors in normal physiological processes [5]. A growing evidence suggested that autophagy also plays an important role in various pathological processes, especially in cancer [5]. In the earliest stages of tumorigenesis, autophagy may limit the development of tumors [6]; however, during the advanced stages of tumors, autophagy is upregulated, and promotes tumor cell proliferation through absorb nutrients and energy driving from degraded proteins and organelles [7]. Autophagy can help cope with intracellular and environmental stresses, such as hypoxia, nutrient shortage, or cancer therapy, thereby favoring tumor progression [5,8,9]. In the tumor microenvironment, autophagy is an important regulator of immune responses, modulating the functions of immune cells and the production of cytokines [10]. Autophagy seems to be a "double-edged sword" in immune cells within the tumor microenvironment. It can promote or suppress tumor development, which depends on the properties of the tumor and cell types.
With the great success of PD-L1 /PD-1 inhibitors and other ICIs in NSCLC, there arise the question of how autophagy would affect cancer treatments in the age of immunotherapy. A comprehensive understanding of autophagy-related genes (ARGs) associated the immune microenvironment would favor the expansion of immunotherapies and explore the predictive biomarkers for the design of patient-tailored combination treatments. At present, few studies have systematically explored the correlation between autophagy and immune microenvironment of LUSC. Moreover, whether ARGs and immune in ltration level could be prognostic in LUSC subtype remains to be fully elucidated.
Therefore, in this study, using the transcriptomic data of LUSC from the TCGA dataset, we developed a prognostic risk signature based on four ARGs in LUSC by random forest algorithm and then explore whether it correlates with the immune landscape. The results shine light on clarifying the association of ARGs and immune escape, and establish a more personalized precision predicting model for immunotherapy.

Data Collection
We collected all transcriptome pro les of LUSC available in the database of TCGA (https://portal.gdc.cancer.gov/).Corresponding clinical information of these patients including gender, age, pathological stage, TNM stage, follow-up time, survival status was also obtained from TCGA (Samples missing any clinical characteristics were excluded and samples of which OS ≤ 30 days were also excluded because these patients probably died of unpredictable factors). Our study included the expression pro le of 32 normal samples and 326 LUSC samples. As the status of distant metastasis is missing in lots of samples, we didn't take the status of M stage into consideration in the current study.
From the GEO dataset (https://www.ncbi.nlm.nih.gov/geo/), transcriptome data of GSE41271 were obtained. A total of 78 LUSC subtype samples from GSE41271 with available clinical and survival data was utilized as an independent validation cohort.

Acquisition of ARGs
The Human Autophagy Database (HADb, http://www.autophagy.lu/) provides a complete and real-time updated list of human genes related with the biological processes of autophagy reported in PubMed or other common databases [11]. We got 232 ARGs form HADb, among which 224 genes were available in the expression pro le from TCGA (The gene list is showed in Table S1).

Differentially Expressed Analysis of ARGs and Functional Enrichment Analysis
To compare the expression level of ARGs between tumor and normal sample, differentially expressed analysis on all ARGs was conducted based on Wilcoxon-test with edgeR package under R environment (version 3.6.3). The cut-off criterion for differentially expressed genes (DEGs) was set as p < 0.05 and log2 |fold change| > 1. The results are displayed with pheatmap package. Then, We performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis to explore the main roles of disordered differentially expressed ARGs (DE-ARGs) in LUSC. The cut-off criterion was p < 0.05 and Benjamin-Hochberg adjusted p < 0.05. The analyses and the visualization of results were conducted under R environment with ClusterPro ler package.

Survival-related ARGs (sARGs)
The ARGs associated with clinical outcomes in LUSC patients were identi ed as sARGs. sARGs were selected by Kaplan-Meier method and veri ed using a log-rank test conducted by the R packages survival and survminer (The cut-off value of the expression level of each gene was 50%, P < 0.05). Final 54 survival-related ARGs were found. These sARGs were speci ed for subsequent research.

Construction and validation of Autophagy-Related Prognostic Signature
We used random survival forest (RSF) method for developing a prognosis model based on sARGs. The random forest algorithm is a machine learning strategy, which is based on the construction of many classi cation (decision) trees that are used to classify the input data vector [12]. Using "RandomForestSRC" package for R, four genes were considered important survival relevant variables for subsequent analysis. The developed RSF prognosis model based on the optimal parameters was then evaluated on the independent dataset where the RSF-based score was derived for each sample. Then, to verify the validity and robustness of the RSF prognosis model, standard Kaplan-Meier survival curves were generated for different risk patient groups on the basis of the RSF-based scores. Based on the optimal cutoff values obtained by the "survminer" R package, LUSC patients were classi ed as low-risk and high-risk according to their risk score. To appraise the prognostic performance of the model, Kaplan-Meier analysis and the log-rank test were employed. Time-dependent receiver operating characteristic (ROC) curves were depicted to evaluate the sensitivity and speci city using the "timeROC" R package [13].
Area under the curve (AUC) values were calculated from the ROC curves. We then performed univariate and multivariate Cox regression analyses to verify the prognostic value of the risk score. We took age, gender, TNM stage (M stage excluded) as candidate risk factors for regression analyses. We evaluated if all these factors are risk factors for poor prognosis by univariate Cox regression analysis, and then by multivariate Cox regression, further determined if the risk score could be utilized for predicting the prognosis of LUSC patients independently.

Gene Set Enrichment Analysis (GSEA)
GSEA was conducted to explore signi cant immune phenotypes between the two groups (high-and lowrisk groups) [14]. The gene sets involved in the negative-regulation of immune response was imported from MSIGDB gmt le from Broad institute. We performed GSEA using "Clusterpro ler" R package and "GSEABase" R package. The cutoff criterion for statistically signi cant terms was set as p < 0.05.

Immune Cell In ltration
The Single-Sample Gene Set Enrichment Analysis (ssGSEA) was introduced to quantify the relative in ltration of 28 immune cell types in the tumor microenvironment. "GSVA'' R package was used in the analysis. Feature gene panels for each immune cell type were obtained from previous studies [15,16].
The relative abundance of each immune cell type was represented by an enrichment score from ssGSEA analysis.

ESTIMATE Algorithm
The ESTIMATE((Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data) ) algorithm was applied to the normalized expression matrix for estimating the stromal and immune scores using "estimate" package (http://r-forge.r-project.org; repos = rforge, dependencies = TRUE) for each LUSC sample [17].

Statistical Analysis
Statistical analyses of this study were conducted using the R software (version 3.6.3), and P < 0.05 was regarded as signi cant.

Identi cation of Differentially Expressed Autophagy-
Related Genes(ARGs) 20 up-regulated and 28 down-regulated genes, respectively (p < 0.05 and log2 |fold change| > 1). The results are shown in a heatmap ( Figure 1A). These differentially expressed ARGs (DE-ARGs) were then selected for KEGG pathway enrichment analysis with a threshold of p < 0.05. As shown in Figures 1B, the differentially expressed ARGs were mainly enriched in regulation of autophagy, apoptosis and mitophagy. It is also worth noting that "PD-L1 expression and PD-1 checkpoint pathway in cancer" were enriched, implicating the correlation between these ARGs and immunosuppression pathway.

Construction of the Autophagy-related Prognostic Signature
The Kaplan-Meier method and log-rank test were applied to examine the relationship between overall survival (OS) and all ARGs (rather than DE-ARGs) in 326 LUSC samples from TCGA. As shown in Table   S2, 54 genes were signi cantly related with OS of LUSC patients (p < 0.05) and these genes were de ned as survival-related ARGs (sARGs). Meanwhile, we noticed that not all sARGs are DE-ARGs. Next, Random Forest was performed to select sARGs with the best prognostic value and to build an autophagy-related risk score model in the TCGA cohort. Four sARGs including CFLAR, RGS19, PINK1 and CTSD were selected as most important survival relevant variables. The Kaplan-Meier overall survival (OS) curves of these four genes were displayed in Figure 2 (The cutoff value of gene expression level was 50%). The random survival forest -based score was derived for each sample and de ned as risk score. Patients in the TCGA cohort (training dataset) then were assigned to a high-risk or low-risk group using the optimal cut-off value obtained with the "survminer" R package. Risk score distribution and corresponding fourgene expression patterns were shown in Figure 3A. With the increase of risk score, the expression levels of the four genes were elevated. The overall survival time of the high-risk group was shorter than that of the low-risk group.

Validation of the Autophagy-related Prognostic Signature
For further validation, 78 samples of LUSC subtype from GSE41271 (GEO dataset) were derived as an independent testing dataset. Then, Kaplan-Meier analysis were applied to the two cohorts and demonstrated that patients with a high-risk score were correlated with worse outcomes in the two cohorts ( Figure 3B-C). Besides, the time-dependent ROC curve analysis of the risk score model in the TCGA cohort indicated a promising prognostic ability for OS (1-year AUC = 0.948, 3-year AUC = 0.965, 5-year AUC = 0.97, Figure 3D). Meanwhile, the ROC curve of OS prediction was drawn in GEO cohort (1-year AUC = 0.649, 3-year AUC = 0.559, 5-year AUC = 0.655, Figure 3E). The autophagy-related model showed a less promising prognostic ability in GEO cohort than the TCGA cohort, possibly due to the small sample size of GEO dataset (78 in GEO vs. 326 in TCGA), which deserved further large samples to validate. Furthermore, in TCGA cohort, we took age, gender, TNM stage (M stage excluded) as candidate risk factors for Cox regression analyses. As shown in Table 1, univariate Cox regression analysis indicated that N stage and risk score are risk factors for poor prognosis, and multivariate Cox regression demonstrated the independence of risk score of this signature in prognosis prediction from other clinical factors.

High Risk Group Indicated an Immune-suppression status
To explore the immune phenotype and the potential immune-related mechanisms underlying our constructed prognostic gene signature, GSEA was performed to nd enriched immune-related gene sets annotated by the gene ontology (GO) term. Results unveiled that seven signi cantly altered immunerelated pathways were enriched in high-risk group, and no immune-related pathways were enriched in the low-risk group (Figure 4). High risk score was signi cantly associated with the macrophage activation involved in immune response (P < 0.05), negative regulation of adaptive immune response (P < 0.01), negative regulation of immune effector process (P < 0.01), negative regulation of immune system process (P < 0.01), negative regulation of innate immune response (P < 0.01), negative regulation of adaptive immune response (P < 0.01), and negative response of immune response (P < 0.01). Consistently, high-risk group exhibited an immunosuppressive phenotype.

High Risk Group showed an increased In ltration of suppressive immune cell populations
To further investigate whether there is an association between autophagy-related score and the regulation of suppressive immune cell populations, normalized ssGSEA scores of 28 in ltrating immune cell populations in each sample were calculated and represented their in ltration level. We drew a heatmap to visualize the relative abundance of 28 in ltrating immune cell populations ( Figure 5). And the boxplot of these immune cells labeled with the P value of Wilcoxon rank test between low-risk and high-risk score group was shown in Supplementary Figure 1. We observed a positive correlation between risk score and the in ltration levels of most immune cell subtypes, including cell types executing anti-tumor immunity (Activated CD4 T cell, Activated CD8 T cell, Central memory CD4 T cell, Central memory CD8 T cell, Effector memory CD4 T cell, Effector memory CD8 T cell, Type 1 T helper cell, Type 17 T helper cell, CD56bright natural killer cell) and cell types executing pro-tumor and immune suppressive functions (Regulatory T cell, Type 2 T helper cell, CD56dim natural killer cell, Immature dendritic cell, Macrophage, MDSC, Neutrophil, Plasmacytoid dendritic cell). For further investigation, Pearson's correlation analysis was applied to the abundances of these two categories of cells. Results showed that the abundances of anti-tumor immune cells were positively associated with the abundances of immunosuppression cells in tumor microenvironment ( Figure 6A). And It seemed that immune suppression is stronger in high-risk group than low-risk group, which is consistent with GSEA results. 3.6 Macrophage and Regulatory T cell (Treg) correlated with worse outcome Then, univariate Cox regression and Kaplan-Meier analysis plus the log-rank test were employed to select the prognostic immune cells among these 28 in ltrating immune cell subtypes. Results revealed that only macrophage and regulatory T cell have a signi cant correlation with overall survival. As shown in Figure   6D-E, higher in ltration of macrophage and regulatory T cell correlated with shorter overall survival time (P < 0.05). As above, high-risk group was associated with higher in ltration of macrophage and regulatory T cell (P < 0.01) ( Figure 6B-C). This observation suggested macrophage and regulatory T cell may accounted for the poor prognosis induced by immunosuppression in high-risk group.

Correlation between genes expression and immunocyte in ltration
The correlations of the expression of four genes in this autophagy-related signature with the in ltration level of macrophage and regulatory T cell as well as risk score were investigated. Figure7A demonstrated overall correlations among these arguments. The four genes were strongly interrelated and exhibited signi cantly positive association with the in ltration of two immunocyte. Besides, the genes expression and immunocyte in ltration were increasing along with the growth of risk score.

Immune Checkpoints Analysis
Immune checkpoints inhibitors are becoming the promising strategies in the treatments of lung cancer. Therefore, we investigated the association of the risk score and the main immune checkpoints, including PD-1/PD-L1 and CTLA4. As Figure 7B-D shown, risk score showed a signi cantly positive correlation with PD-1 and CTLA4 (P < 0.01). While, there is no signi cant correlation between risk score and PD-L1 expression. In summary, high-risk group had higher levels of PD-1 and CTLA4. Then, we also examined the correlations between the expression of the four ARGs (CFLAR, RGS19, PINK1, CTSD) and the expression of immune checkpoints (PD-1, PD-L1, CTLA4). Results showed that the expressions of four ARGs had signi cant associations with the immune checkpoints consistently (only the correlation between PINK1 and PD-L1 was not signi cant) (Supplementary Figure 2).
3.9 Estimate score, immune score and stromal score ESTIMATE (Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data) is a newly developed algorithm that uses the transcriptional pro les of cancer tissues to infer the level of in ltrating stromal and immune cells based on speci c gene expression signatures of stromal and immune cells in the speci c cancer type. In the present study, the stromal score, immune score, and estimate score of each included sample of was calculated by applying Estimate R package. Next, the association of autophagy-related risk scores with the estimate/immune/stromal scores was examined. As shown in Figure 8A-C, estimate score and immune score were found to have a positive correlation with the risk scores of our established model (P < 0.05). Stromal score was also rising with the increase of risk score, but not signi cantly (P=0.072). The prognostic value of the three types of score were also investigated. Survival analysis showed that the patients with higher estimate/immune/stromal score had a poorer overall survival, although the statistical signi cance were only observed in stromal score (P < 0.05).

Discussion
Over the past decade, new treatments for patients with LUSC have evolved dramatically, including immune checkpoint inhibitors and combination strategies. Thus, the more re ned risk-strati cation to guide treatment for LUSC patients beyond conventional TNM staging are needed. With the rapid development of high-throughput next generation sequencing, machine leaning, and bioinformatics analyses, prognostic signatures composed of a set of genes have been designed to meet this need.
Compared to the single molecular predictors, signatures integrated the prognostic value of several genes and seemed to have a better outcome prediction ability. In previous studies, several risk signatures have shown satisfactory effects in outcome prediction and treatment guiding in lung cancer, such as miRNAbased signature [18], lncRNA-related signature [19], immune-associated signature [20]. These risk signatures could be used as the tailored algorithms to individualize therapy.
Autophagy, which induce degradation of proteins and organelles or cell death upon cellular stress, are crucial in the pathophysiology of malignant tumor [9,21]. The apoptotic effect of autophagy is controversial as both inhibitory and stimulatory effects have been reported in NSCLC [21]. Consequently, an integrated study of multi-molecules genes model is needed to gure out the exact effects of autophagy in lung cancer. A previous study has reported an autophagy-associated multiple-gene signature that correlated with survival in NSCLC [22]. However, in the era of immunotherapy, the immune contexture of tumor microenvironment is essential for the treatment selection and survival prediction [23]. The correlations between autophagy and immune pathway should not be ignored [10]. Thus, in this study, we conducted a complete and meaningful analysis of autophagy-related genes involved in outcome prediction and associated immune landscape, which may also enable the identi cation of patients who are more likely to respond to immunotherapy. Moreover, different with previous studies, random forest instead of Lasso regression was applied in the current analyses. Random forest is a process of machine learning with accurate algorithms and high e ciency in classi cation prediction.
In the present study, four ARGs including CFLAR, RGS19, PINK1 and CTSD were screened out as most important survival relevant variables and were selected for the construction of an autophagy-related prognostic risk signature. CFLAR, also known as c-FLIPL, is a critical anti-apoptotic protein that inhibits cell death mediated by the death receptors Fas, DR4, DR5, and TNF-R1 [24,25]. CFLAR was found as an independent adverse prognostic biomarker in different cancer types [26,27]. Besides, elevated expression of CFLAR is associated with tumor cells escaping from immune surveillance in vivo, correlates with a more aggressive tumor, and is also considered to be the main cause of immune escape [25]. A recent study reported that down regulation of CFLAR could enhance the antitumor response of T cells and enhances the PD-1 blockade e cacy in melanoma tumor model [28]. And knockdown of CFLAR could also decrease the expression of PD-L1 [28]. Consistent with these results, our analyses found that high expression of CFLAR was signi cantly associated with some key immune checkpoint including PD-1, PD-L1 and CTLA4 (Supplementary Fig. 2). These ndings may provide a new combined therapeutic target for further improving the e cacy of ICIs. RGS19 is a member of the regulators of G protein signaling (RGS proteins) [29]. Rapid termination of G protein signals by RGS proteins can potentially modulate growth signals and hence promote tumorigenesis, although no clear mechanism of RGS19 in lung cancer was reported [29]. Besides, physiological functions of most RGS proteins in immune response are largely unknown [30]. Our study rst demonstrated RGS19 could be a risk factor for poor prognosis, and also correlated with the expression of several immune checkpoint and the in ltration of Tregs and Macrophages ( Supplementary Fig. 2, Fig. 7A). This discovery provides a new perspective that RGS19 may participate in immune process regulation, which worth further studies. PINK1 (PTEN induced kinase 1) mediates the recruitment of Parkin to mitochondria, which facilitates the elimination of the injured mitochondria by autophagy. It also has been proved to play crucial roles in the genesis and development of tumor [31]. In lung carcinoma tissue, the expression of PINK1 was raised, which was associated with a poor prognosis [32]. In recent years, PINK1 was also found as a repressor of the immune system and played a key role in adaptive immunity by repressing presentation of mitochondrial antigens [33]. Similar with previous studies, our study reported PINK1 was a risk factor for worse outcome. Moreover, PINK1 was found to be associated with the expression of PD-1 and CTLA4 as well as the in ltration of Tregs and macrophages. CTSD (Cathepsin D) is a key protein for lysosomal function that is necessary for autophagy in cancer cells [34]. Although CTSD is expressed at high levels in many cells of the immune system, but its role in immune function is still not well understood. In the present investigation, we observed that, similar with the other three genes, CTSD was correlated with poor survival and immunosuppressive status as well.
Taken together, it is reasonable to believe that the combination of these four genes has robust prediction value in LUSC, and the four-gene signature could stratify patients with different immune contexture, which is critical in the context of current immunotherapy. Moreover, GSEA analysis was employed to make the functional annotation, and we found the more abundant negative-regulation of immune responses and processes in the high autophagy-related score group. Later, survival analyses were employed to examine whether the in ltration of 28 in ltrating immune cell subtypes were prognostic and it turned out that only macrophage and Treg cells have a signi cantly negative correlation with overall survival. While, these two types of cells delivering pro-tumor suppression were consistently enriched in high-risk group. Recent studies have shown that autophagy signi cantly controls immune responses by modulating the functions of immune cells and the relationship between autophagy and immunity are complicated [35]. On one hand, autophagy has been shown to be important for priming of tumor-speci c CD8 + T cells, and inhibition of autophagy would impair systemic immunity [36]. However, on the other hand, the induction of autophagy may also bene t tumor cells escape from immune surveillance [35]. A recent report indicated that immunosuppressive Treg cells are critically dependent on autophagy [37]. Results showed that autophagy is active in Treg cells and supports their lineage stability and survival tness [37]. While, for macrophage, autophagy regulates cellular development of monocytes, resulting in the disturbance of macrophage differentiation. Thus, our ndings that high-risk group with higher expression of the four autophagy-related genes had an increased in ltration level of Treg cells and macrophages should be noticed. Moreover, based on estimate/immune/stromal scores from the ESTIMATE algorithm, we found that both the estimate score and immune score of high-risk group are higher and patients with higher scores tended to have a worse outcome.
The high-risk group presented a comparatively suppressed immune status, the phenotype of adaptive immune evasion. Actually, the in ltration of almost all immune cells were upregulated in high-risk group, including both the immune-suppressive subtypes and immune-stimulatory subtypes. This observation possibly suggested the presence of a feedback mechanism such that the recruitment or differentiation of cells specialized for immune suppression may be facilitated by anti-tumor in ammation. This could also be explained by our results that the immunosuppressive molecules points such as PD-1, CTLA4 was also higher in high-risk group, which means the presence of a suppressed pre-existing antitumor immunity that could be re-invigorated by anti PD-1/PD-L1 immunotherapy [38]. Besides, a previous review has been proposed that four different types of tumor microenvironment exist based on the presence or absence of tumor-in ltrating lymphocytes (TILs) and PD-L1 expression [39]. Type I cancers (with higher PD-L1 + and TILs) are more likely to bene t to anti-PD-1/L1 therapy [39]. Thus, we speculated that the high-risk group are more likely to respond to ICIs.

Conclusion
In conclusion, this study developed a new autophagy-related four-gene prognostic risk signature, which could be applied as an independent prognostic indicator for LUSC patients. And the autophagy-related scores are bound up with immune phenotype of LUSC, with higher score indicating an immunesuppression status. This study provides a novel and comprehensive sight to the correlation of autophagy

Consent for publication
Not applicable.

Availability of data and materials
Publicly available datasets were analyzed in this study. These can be found here: TCGA database (https://portal.gdc.cancer.gov/) and the NCBI Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/).

Competing interests
The authors declare that the research was conducted in the absence of any commercial or nancial relationships that could be construed as a potential con ict of interest.    The results of survival analyses.  The results of survival analyses.      Heatmap of immune cell in ltration level of LUSC tumor samples in TCGA cohort. A Single-sample gene set enrichment analysis identifying the relative in ltration of 28 immune cell populations for LUSC tumor samples. Samples in the heatmap were ranked by risk score of each patient. The ssGSEA score which represents the relative in ltration of each cell type was normalized to unity distribution, for which zero is the minimal and one is the maximal score for each immune cell type. (red represents high and blue represents low in ltration). The three parts of the heatmap exhibited the three types of immune cells (antitumor immunity, pro-tumor immune-suppression, and other unclassi ed immune cells).

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.