Identification of Autophagy- and Ferroptosis-Related lncRNAs Functioned through Immune-Related Pathways in Head and Neck Squamous Carcinoma

The interplay between autophagy and ferroptosis has been highlighted as an important event to decide cancer cell fate. However, the underlying mechanisms remain largely unclear. In this study, we systematically explored the expression, prognostic value and functional roles of lncRNA in autophagy and ferroptosis. By a set of bioinformatics analyses, we identified 363 autophagy- and ferroptosis-related lncRNAs (AF-lncRNAs) and found 17 of them are dramatically related to the prognosis of head and neck squamous cell carcinoma (HNSC) patients, named as prognosis-related AF-lncRNAs (PAF-lncRNAs). Based on six key PAF-lncRNAs, a risk score model was developed and used to categorize the TCGA-retrieved HNSC patients into two groups (high-risk vs. low-risk). Functional analysis showed the differentially expressed genes (DEGs) between the two groups were mainly enriched in immune-related pathways and regulated by a PAF-lncRNA-directed ceRNA (competitive endogenous RNA) network. Combined with a variety of immune infiltration analyses, we also found a decreased landscape of immune cell infiltration in high-risk groups. Together, by revealing PAF-lncRNAs with tumor prognostic features functioned through immune-related pathways, our work would contribute to show the pathogenesis of a lncRNA-directed interplay among autophagy, ferroptosis and tumor immunity in HNSC and to develop potential prognostic biomarkers and targets for tumor immunotherapy.


Introduction
Autophagy-dependent cell death was described as a form of regulated cell death (RCD) that mechanistically depends on the autophagic machinery or components [1]. An increasing number of discoveries have built strong links between autophagy and various types of RCD, including ferroptosis [2][3][4]. Ferroptosis is an iron-dependent form of regulated cell death, which mainly depends on the accumulation of iron and reactive oxygen species (ROS) [5]. Studies have found that autophagic machinery could contribute to ferroptosis by mediating the degradation of ferritin and some genes that are involved in the crosstalk of ferroptosis and autophagy, and thus contribute to ferroptotic cancer cell death [4,6]. Thus far, emerging studies have implicated that the interplay between Life 2021, 11, 835 3 of 17 2020) [28]. We selected 493 tumors with complete follow-up information and survival time longer than 30 days. Phenotype and survival information of these patients are shown in Table S1. A total of 376 autophagic genes were gathered from the overlap of the Autophagy database (http://tp-apg.genes.nig.ac.jp/autophagy, accessed on 22 March 2021) [29] and Human Autophagy Modulator database (HAMdb, http://hamdb.scbdd.com, accessed on 22 March 2021) (Table S2) [30]. A total of 184 human ferroptotic genes were obtained from the FerrDb database (www.zhounan.org/ferrdb, accessed on 19 March 2021) [31] and literature studies (Table S2) [32,33].
In order to construct an AF lncRNA prognostic model, 493 patients were divided into a training set (60%, 294 samples) and testing set (40%, 199 samples). Univariate Cox proportional hazards analysis was used to identify prognosis-related AF-lncRNAs in the training set. Subsequently, multivariate Cox analysis was used to construct the prognostic risk model by employing the survival R package.

Development and Assessment of the Prognostic Risk Model
According to multivariate Cox analysis, the risk score was calculated by the following formula: risk score = coef (lncRNA1) × exp (lncRNA1) + coef (lncRNA2) × exp (lncRNA2) + . . . + coef (lncRNAn) × exp (lncRNAn) [35][36][37]. The coefficient values (coef) were calculated following previously reported methods [37]. The coef represents the coefficient of the corresponding lncRNA, which was calculated by using the survival coxph function of the R package. The "exp" represents the expression of corresponding lncRNA. Based on the median risk-score, HNSC patients in the training set and testing set were divided into highand low-risk groups, respectively. The univariate and multivariate Cox regression analyses were used to calculate the prognostic value of risk-score and clinical features by survival R package. The landscape of survival status was described between the high-and low-risk patients by using the ggplot2 and pheatmap R packages. The principal component analysis (PCA) was applied to verify the classification between high-and low-risk groups. Then Kaplan-Meier survival analysis was performed to estimate the survival difference between these two groups by using the survival and survminer R packages. In order to estimate the sensitivity and specificity of PAF-lncRNAs signatures, we employed time-dependent receiver operating characteristic (ROC) curves and multi-factor ROC curves by using the timeROC (Version:0.4) and survivalROC R package (Version:1.0.3).

Functional Analysis of the DEGs between High-and Low-Risk Groups
The differentially expressed genes between high-and low-risk groups were screened out by the limma R package with the criteria of |log2 (fold change)| > 1 and p < 0.05. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were performed by limma and clusterProfiler R packages. We used GOplot and ggplot2 R packages to display the result of the functional enrichment analysis.

Estimation of Immune Cell Infiltration Proportion
An estimation algorithm was used to calculate the stromal score, immune score, and estimate score by estimate and utils R packages [38]. ImmuCellAI and CIBERSORT algorithms were used to estimate the proportion of immune cell infiltration, respectively. For ImmuCellAI, the estimated immune cell proportion of 24 types were calculated by the Immune Cell Abundance Identifier website (ImmuCellAI) [39]. For the CIBERSORT algorithm, the 22 human hematopoietic cell phenotypes file (LM22.txt) was obtained and the immune cell proportion of 22 types were estimated by performing "source ("CIBERSORT.R")" [40].

Statistical Analysis
All statistical analyses were performed by using R software (version 4.0.3). For each analysis, statistical significance was set at p < 0.05 without special description. Nonnormally distributed variables were analyzed using the Wilcoxon test, and the Benjamin Hochberg method was used to calculate FDR. The survival curves between high-and lowrisk groups were assessed by using the Kaplan-Meier survival analysis with the log-rank test.

Identification of Autophagy-and Ferroptosis-Related lncRNAs (AF-lncRNAs) in HNSC
To identify the autophagy and ferroptosis correlated lncRNAs in HNSC, we used autophagic and ferroptotic genes to construct co-expression networks in HNSC patients ac-  Table S2). risk groups were assessed by using the Kaplan-Meier survival analysis with the log-rank test.

Identification of Autophagy-and Ferroptosis-Related lncRNAs (AF-lncRNAs) in HNSC
To identify the autophagy and ferroptosis correlated lncRNAs in HNSC, we used autophagic and ferroptotic genes to construct co-expression networks in HNSC patients according to the criteria of |Correlation Coefficient| > 0.4 and p < 0.05 by Pearson correlation analysis. Our results revealed 363 autophagy-and ferroptosis-related lncRNAs, namely, AF-lncRNAs ( Figure 2a, Table S2).

Construction of a Prognostic Risk Model in the TCGA Cohort Based on the Identified Prognostic AF-lncRNAs in HNSC
To investigate the prognostic value of AF-lncRNAs in HNSC patients, we performed univariate Cox analysis to estimate the prognostic relationship between AF-lncRNAs and overall survival (OS) in 493 HNSC samples from TCGA. Our result showed that 17 lncRNAs, referred to as prognostic AF-lncRNAs (PAF-lncRNAs), were significantly associated with the survival of HNSC patients (p < 0.01, Figure 2b, Table S3). We subsequently performed the multiple stepwise cox regression analysis to investigate the impact of these 17 prognostic-associated AF-lncRNAs on patient survival time and clinical outcomes. PCED1B-AS1, AL512274.1, AL354836.1, MIR9-3HG, MIR4435-2HG and LINC02541 were found to be independent factors in HNSC ( Figure 3a, Table S3). Then, we identified a network between the six key PAF-lncRNAs and autophagy and/or ferroptosis genes. Most of these genes, including CDD4, ELAVL1, FKBP8, NBR1, RB1CC1, and ZEB1, were correlated

Construction of a Prognostic Risk Model in the TCGA Cohort Based on the Identified Prognostic AF-lncRNAs in HNSC
To investigate the prognostic value of AF-lncRNAs in HNSC patients, we performed univariate Cox analysis to estimate the prognostic relationship between AF-lncRNAs and overall survival (OS) in 493 HNSC samples from TCGA. Our result showed that 17 lncR-NAs, referred to as prognostic AF-lncRNAs (PAF-lncRNAs), were significantly associated with the survival of HNSC patients (p < 0.01, Figure 2b, Table S3). We subsequently performed the multiple stepwise cox regression analysis to investigate the impact of these 17 prognostic-associated AF-lncRNAs on patient survival time and clinical outcomes. PCED1B-AS1, AL512274.1, AL354836.1, MIR9-3HG, MIR4435-2HG and LINC02541 were found to be independent factors in HNSC ( Figure 3a, Table S3). Then, we identified a network between the six key PAF-lncRNAs and autophagy and/or ferroptosis genes. Most of these genes, including CDD4, ELAVL1, FKBP8, NBR1, RB1CC1, and ZEB1, were correlated with patients' overall survival, which again supports that the PAF-lncRNAs were correlated with the process of autophagy and/or ferroptosis (Figure 3b and Figure S1).
Utilizing the aforementioned independent PAF-lncRNAs, we next constructed a prognostic predictive model. The risk-score of each patient was calculated according to the following formula (Table S4): By using univariate and multivariate Cox regression analysis, we compared the prognostic significance of risk-score and different clinical characteristics in HNSC. Our results showed that only age and risk-score calculated by the six PAF-lncRNAs were dramatically correlated with HNSC OS, and the risk-score was shown to be the most significant parameter (p < 0.05, Figure 3c,d). Together, these data suggested the intimated with patients' overall survival, which again supports that the PAF-lncRNAs were correlated with the process of autophagy and/or ferroptosis (Figures 3b and S1). Utilizing the aforementioned independent PAF-lncRNAs, we next constructed a prognostic predictive model. The risk-score of each patient was calculated according to the following formula (Table S4): By using univariate and multivariate Cox regression analysis, we compared the prognostic significance of risk-score and different clinical characteristics in HNSC. Our results

The Risk Model Exhibited Robust Predictive and Discriminative Ability for HNSC Patients
Based on the risk model, HNSC patients from TCGA were divided into high-and lowrisk groups by calculating median risk-score. The PAF-lncRNAs expressions, survival status and risk-scores of these patients were displayed in Figure S2. To assess the accuracy of the stratification, we conducted a set of bioinformatic analyses. First, principal component analysis (PCA) results suggested that risk-score works well in distinguishing high-risk patients with low-risk groups in the TCGA training, testing, and entire sets (Figure 4a). Second, Kaplan-Meier analyses showed the patients' OS rate was dramatically lower in the high-risk group compared to patients in the low-risk group in the TCGA training, testing, and entire sets (Figure 4b). Third, time-dependent receiver operating characteristic (ROC) curves were calculated, and the area under the ROC curves (AUC) of risk-scores at 1, 2, 3 and 4 years for survival prediction were all above 0.6, which suggested the risk-score had moderate prediction accuracy in TCGA training, testing, and entire sets (Figure 4c). In addition, the AUC of the risk-score exhibited moderate performance compared to other measured phenotypes in predicting the prognosis of HNSC patients in TCGA training, testing, and entire sets (Figure 4d). Taken together, these results suggested that the risk model based on the six PAF-lncRNAs has good predictive ability and stratification accuracy for HNSC patients.

Differentially Expressed Genes (DEGs) between High-and Low-Risk Groups Were Dramatically Enriched in Immune-Related Pathways
To study the differences between high-and low-risk groups at the whole genomewide level, we performed differential expression analysis. A total of 437 DEGs (downregulation: 411 and up-regulation: 26) between two groups were identified (Table S5), and the expressions of them were displayed in a heat map (Figure 5a). We also constructed a PAF-lncRNAs ceRNA network through the ENCORI database to explore the interaction between lncRNAs and DEGs (Figure 5b). We found that there are 3 PAF-lncRNAs, AL512274.1, MIR4435-2HG, and MIR9-3HG, which might regulate the 56 DEGs through

Differentially Expressed Genes (DEGs) between High-and Low-Risk Groups Were Dramatically Enriched in Immune-Related Pathways
To study the differences between high-and low-risk groups at the whole genome-wide level, we performed differential expression analysis. A total of 437 DEGs (down-regulation: 411 and up-regulation: 26) between two groups were identified (Table S5), and the expressions of them were displayed in a heat map (Figure 5a). We also constructed a PAF-lncRNAs ceRNA network through the ENCORI database to explore the interaction between lncRNAs and DEGs (Figure 5b). We found that there are 3 PAF-lncRNAs, AL512274.1, MIR4435-2HG, and MIR9-3HG, which might regulate the 56 DEGs through 26 miRNAs. The majority of the DEGs, e.g., SOX11, TBX21, FAM3B, FGF5, HNF1A, MYB, and PLAC8 have been reported to function as promoters or biomarkers in various cancer types [44][45][46][47][48][49][50]. Among these, PLAC8 were found to regulate autophagy-related functions in a variety of cancers [51,52]. In addition, numerous DEGs were reported to be involved in the proliferation and differentiation of immune cells and could regulate immune functions, including CD226, IRF4, LY9, MS4A1, TBX21, TNFRSF17, FCRLA, and SLAMF6 [45,[53][54][55][56][57][58].  To further characterize and understand the biological insights of these DEGs, we performed gene ontology (GO) and KEGG analysis. The DEGs were mainly enriched in immune-related processes, such as T cell differentiation, B cell activation, immune receptor  To further characterize and understand the biological insights of these DEGs, we performed gene ontology (GO) and KEGG analysis. The DEGs were mainly enriched in immune-related processes, such as T cell differentiation, B cell activation, immune receptor activity, cytokine receptor binding, regulation of lymphocyte and CD4-positive, alpha-beta T cell activation (Figure 5c). KEGG also displayed that the DEGs were related to immune pathways (Figure 5d). Together, these results indicated PAF-lncRNAs might mediate the generation of immune microenvironmental differences of HNSC high-and low-risk groups.

Distinct Immune Landscapes between High-and Low-Risk HNSC Patients Were Identified
Following the aforementioned results, we then systematically investigated the immune landscape of the two risk groups. Firstly, we compared the tumor immune microenvironment between high-and low-risk groups (Figure 6a, Table S6). Our data showed that the high-risk group was marked by higher tumor purity and lower immune infiltration levels than the low-risk group (Figure 6b). Secondly, we compared the immune cell composition of the tumor immune microenvironment by using the ImmuCellAI and CIBERSORT algorithm, respectively (Figure 6b). According to the ImmuCellAI calculation results, the largest proportion of infiltrating immune cells in HNSC patients were cytotoxic T cells (Tc), CD4 + T cells, CD8 + T cells, NK cells, Tfh and Th1 ( Figure 6C, Table S7). The relative proportion of these immune cells was significantly higher in low-risk patients than in high-risk patients (p < 0.05; Wilcoxon test, Figure 6c). However, the proportions of dendritic cells, monocytes, macrophages and neutrophils in the high-risk group were significantly higher than the low-risk group (p < 0.05; Wilcoxon test, Figure 6c). Specifically, the CIBERSORT algorithm can distinguish the subtypes of macrophages. By using this algorithm, we found a higher proportion of the M2 macrophages in the high-risk group, while M1 macrophages were higher in the low-risk group (Figure 6d, Table S7). As reported in previous studies, macrophages undergo a switch that leads to differentiation into either inflammatory (M1) or regulatory (M2) subtypes. Among these, M1 is mainly involved in tumor killing, while M2 is mainly involved in supporting tumor growth [59]. Our results suggested that the polarization of macrophages might be a regulatory mechanism for the difference of survival between the high-and the low-risk group. Thirdly, the two groups showed significant differences in immune checkpoint and immune activation gene expression levels (p < 0.05; Figure 6e, Table S8). The high-risk group was associated with relatively lower immune checkpoint and activation signal expression levels, whereas the low-risk group was associated with the higher expression level.
In summary, our results revealed that immune-driven cells are associated with the lowrisk group, while immune-regulatory cells tend to be more common in the high-risk group. The consistency between the immune profile and prognostic profile in the two groups also implied that our classification method is scientific and reasonable. As previous studies have shown that the abundance of T cell subsets, particularly that of tumor infiltrating T cells, could influence clinical curative effects and prognosis [60,61], our risk model has the potential to be applied for predicting the immunotherapy response. In summary, our results revealed that immune-driven cells are associated with the low-risk group, while immune-regulatory cells tend to be more common in the high-risk group. The consistency between the immune profile and prognostic profile in the two groups also implied that our classification method is scientific and reasonable. As previous studies have shown that the abundance of T cell subsets, particularly that of tumor infiltrating T cells, could influence clinical curative effects and prognosis [60,61], our risk model has the potential to be applied for predicting the immunotherapy response.

Discussion
As two closely linked forms of RCD, increasing evidence has shown that autophagy and ferroptosis are intimately associated with tumor progression [62][63][64]. However, the engagement of lncRNAs in HNSC autophagic and ferroptotic processes have not been thoroughly and systematically studied. In the present study, we systematically investigated the expression and prognostic relevance of 363 autophagy-and ferroptosis-related (e) The differential expression of immune checkpoints and immune activation signals between high-and low-risk groups. * p < 0.05, ** p < 0.01, *** p < 0.001.

Discussion
As two closely linked forms of RCD, increasing evidence has shown that autophagy and ferroptosis are intimately associated with tumor progression [62][63][64]. However, the engagement of lncRNAs in HNSC autophagic and ferroptotic processes have not been thoroughly and systematically studied. In the present study, we systematically investigated the expression and prognostic relevance of 363 autophagy-and ferroptosis-related lncRNAs (AF-lncRNAs) in HNSC. Previous reports have shown the connection between lncRNAs and RCD, and their relationship is associated with tumor progression [65][66][67]. In this paper, by performing univariate cox regression analysis, multiple stepwise cox regression analysis, survival analyses and ROC analyses, etc., a total of six key PAF-lncRNAs in HNSC were identified, including MIR4435-2HG, PCED1B-AS1, AL512274.1, LINC02541, AL354836.1, and MIR9-3HG. Previous studies have reported that the expression of MIR4435-2HG, PCED1B-AS1, AL512274.1, AL354836.1, MIR9-3HG were associated with tumorigenesis and regulated cell death in various tumor types, which support the further identification and exploration in HNSC [68][69][70][71][72]. For example, overexpression of MIR4435-2HG will promote tumor cell proliferation while knockdown of MIR4435-2HG will lead to cell death [73][74][75][76]. In addition, researchers have found that inhibition of MIR4435-2HG would downregulate Nrf2, which would alter the resistant status of head and neck cancer cells to a more sensitive status to ferroptosis and eventually promote ferroptotic cell death [13,74].
Functional analyses revealed that these PAF-lncRNAs have a close relationship with several autophagic and/or ferroptotic genes, including FANCD2, CD44, PROM2, ZEB1, IGF1R, AKT1S1, FKBP8, NBR1, RB1CC1, and ELAVL1. Among these, IGF1R, AKT1S1, FKBP8, NBR1, and RB1CC1 have been reported to function as regulatory factors, such as autophagic adaptors or receptors, to regulate the autophagic processes [77][78][79][80][81]. FANCD2, PROM2, and ZEB1 have been reported to regulate ferroptosis through regulating the accumulation of lipid ROS and intracellular iron export, etc. [82][83][84]. As for ELAVL1 and CD44, they have been reported to regulate the interplay of autophagic and ferroptotic processes [85][86][87]. For example, ELAVL1 could promote ferroptosis by regulating autophagy in myocardial ischemia/reperfusion injury [85]. Autophagic flux impairment induced a high expression of CD44 and thus induced mitochondrial dysfunction, oxidative stress and cancer cell death [86]. The interaction between PAF-lncRNAs and these autophagic and/or ferroptotic genes indicated that PAF-lncRNAs might participate in the regulation of autophagy and ferroptosis and thus mediate autophagic and ferroptotic tumor cell death through these genes. However, its detailed molecular mechanism still needs further future genetic and experimental studies.
We then produced a risk model based on these PAF-lncRNAs. The ROC analysis revealed that these PAF-lncRNAs signatures have better diagnostic capability of selecting the high-risk HNSC patients with poor prognosis. Based on the risk model, 493 HNSC patients from TCGA were divided into high-and low-risk groups. The divided high-and lowrisk patients showed distinct gene expression patterns, and the DEGs were dramatically enriched in many immune-related pathways.
Tumor immunity and RCD have been linked together from recent reports [88,89]. Although several findings have supported the importance of immunology in HNSC [90][91][92], the underlying molecular mechanism and potential modulation between tumor immunity and RCD, especially autophagy and ferroptosis, are largely unclear. In this paper, GO and KEGG analyses based on the DEGs enriched many immune-related biological processes and pathways. For instance, T cell differentiation is the process in which precursor cell types acquire characteristics of more mature T cells to achieve immune effects, while B cell activation is defined by the change in morphology and behavior of mature or immature B cells. The function of immune receptor activity is to receive signals and transmit them into cells to initiate the immune response. In addition, we established a connection between 56 DEGs and PAF-lncRNAs from the ENCORI database, which was based on the CLIPseq data. Numerous independent studies have validated the regulatory roles of the six PAF-lncRNA on those DEGs through competitive binding of the miRNA regulators. For example, PCED1B-AS1 could inhibit tumor cell death by cooperating with the miR-194-5p/PCED1B axis in glioma [93], and MIR4435-2HG could regulate cancer proliferation through sponging miR-206, miR-802, miR-128-3 and miR-1224-5p and regulated YAP1, FLOT2, et al., in various cancers [73,[94][95][96]. Among all these DEGs, nearly half (24/56) have been reported to be involved in activation and differentiation of immune cells or regulation of the immune response [45,[53][54][55][56][57][58]. One of them, placenta associated 8 (PLAC8), has been reported to regulate autophagy by suppressing the production of IL-18, which is a pro-inflammatory cytokine that is capable of stimulating interferon gamma production and of regulating T helper cell responses [97]. These pieces of evidence support the idea that the PAF-lncRNAs might participate in regulating HNSC through the tumor immunity process.
HNSC is a disease that was previously characterized by immunosuppression [98]. Recently, considerable progress has been made in immune checkpoint inhibitor (ICI)-based HNSC treatment. However, the response rate of recurrent or metastatic HNSC to PD-1/PD-L1 inhibitors is only 13.3-22%, as per previous clinical trials [98]. Therefore, the selection of patients that can effectively respond to ICIs is crucial. For this reason, we estimated the immune cell infiltration of HNSC tumor samples, and results revealed the discriminated immune microenvironment landscapes of distinct risk groups. The infiltration ratio of effective T cells, NK cells, T helper cells, B cells and M1 macrophages that were related to anti-tumor effects were higher in the low-risk group. Effective T cells and B cells play critical roles in tumor control, and T helper cells can stimulate B cells for an immune response. The natural killer (NK) cells were discovered for their ability to recognize and kill tumor cells [99] and to release a number of cytokines that regulate both innate and adaptive immune responses [100]. As for M2 macrophages, the proportions of dendritic cells (DCs), mast cells and neutrophils that related to pro-tumor effects were significantly higher in the high-risk group than low-risk group by both algorithms. Due to the complex phenotype and cancer heterogeneity, the infiltration of dendritic cells and mast cells have controversial results in predicting clinical outcomes in different tumors [101][102][103]. Because parameters of the immune contexture have been associated with treatment efficacy, it is important to characterize the baseline HNSC immune milieu to clarify the composition and property of tumor-infiltrating immune cells [104]. Together, these results supported the involvement of PAF-lncRNAs in regulating HNSC tumor immunity. Thus, the identification of PAF-lncRNAs not only provides a potential target for anti-cancer immunotherapy but also build a bridge between RCD and immunogenicity of HNSC, which might shed new light on revealing another layer of lncRNA-directed immunogenicity of cancer cells.

Conclusions
In summary, we systematically explored the expression and prognostic value of autophagy-and ferroptosis-related lncRNAs by a series of bioinformatics analyses in HNSC. Our study revealed six prognosis-related AF-lncRNAs and developed a novel prognostic model based on these lncRNAs. This model proved to be an independent prognostic factor, which has a favorable predictive effect on prognosis for HNSC. In addition, we revealed these AF-lncRNAs functioned through multiple critical tumor immune-related processes. Our results would contribute to show the pathogenesis of HNSC and to develop new treatment targets and prognostic molecular markers.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/life11080835/s1, Figure S1: Expression heat map, risk score distribution, and survival status of TCGA-retrieved HNSC patients, Table S1: Clinical characteristics of TCGA-retrieved HNSC patients in the indicated sets, Table S2: The list of autophagic genes, ferroptotic genes, and AF-lncRNAs, Table S3: Univariate /Multivariate Cox proportional hazards analysis, Table S4: List of survival time, status, and risk-score of TCGA-retrieved HNSC patients, Table S5: DEGs between high-and low-risk patients, Table S6: List of stromal score, immune score, estimate score and tumor purity in the indicated groups, Table S7: List of ImmuCellAI score, CIBERSORT score in the indicated groups, Table S8: The expression of immune check point genes in the indicated groups.
Author Contributions: T.S. and Q.G. initiated the project; T.S. and Q.G. designed the analytical process; X.W. organized and supervised the whole project; Q.G., T.S., and X.Z. prepared and performed all bioinformatic analyses; Q.G. drafted the manuscript with input from X.W., T.S. and X.Z., Q.G. and X.Z. contributed equally to this paper. All authors have read and agreed to the published version of the manuscript. Institutional Review Board Statement: All the expression data and clinical information were obtained from publicly available data sets which were free to download and analyze without limitations. Investigators of each study obtained the approval from their local ethics committee and informed patient consent.
Informed Consent Statement: Not applicable.