Comprehensive Characterization of Immune Landscape Based on Tumor Microenvironment for Oral Squamous Cell Carcinoma Prognosis

Objective: This study aims to identify an immune-related signature to predict clinical outcomes of oral squamous cell carcinoma (OSCC) patients. Methods: Gene transcriptome data of both tumor and normal tissues from OSCC and the corresponding clinical information were downloaded from The Cancer Genome Atlas (TCGA). Tumor Immune Estimation Resource algorithm (ESTIMATE) was used to calculate the immune/stromal-related scores. The immune/stromal scores and associated clinical characteristics of OSCC patients were evaluated. Univariate Cox proportional hazards regression analyses, least absolute shrinkage, and selection operator (LASSO) and receiver operating characteristic (ROC) curve analyses were performed to assess the prognostic prediction capacity. Gene Set Enrichment Analysis (GSEA) and Gene Ontology (GO) function annotation were used to analysis the functions of TME-related genes. Results: Eleven predictor genes were identified in the immune-related signature and overall survival (OS) in the high-risk group was significantly shorter than in the low-risk group. An ROC analysis showed the TME-related signature could predict the total OS of OSCC patients. Moreover, GSEA and GO function annotation proved that immunity and immune-related pathways were mainly enriched in the high-risk group. Conclusions: We identified an immune-related signature that was closely correlated with the prognosis and immune response of OSCC patients. This signature may have important implications for improving the clinical survival rate of OSCC patients and provide a potential strategy for cancer immunotherapy.


Background
Oral squamous cell carcinoma (OSCC) is one of the most common malignancies arising in the oral cavity worldwide, and the incidence of OSCC has been increasing by at least 1% annually, especially in China [1,2]. Common treatment options for OSCC include surgery, chemotherapy, radiotherapy, immunotherapy, and these treatment methods have reduced the mortality rate of OSCC partly [3,4]. However, despite the rapid progress in the treatment of OSCC, its mortality and incidence are still rising due to the tumor heterogeneity [5]. Therefore, in order to improve the survival rate of OSCC, the development of new predictive biomarkers to accurately predict the prognosis is of great significance for OSCC patients [6].
Immunotherapy is one of the recent breakthroughs in cancer therapy and becoming a new promising approach in cancer treatment, including OSCC [7,8]. However, due to tumor heterogeneity and multifaceted immunosuppressive signals within the tumor 2 of 16 microenvironment (TME), the efficacy of cancer immunotherapy is limited [9,10]. TME refers to the cellular environment where tumor cells reside, and is composed of tumor cells, immune cells, fibroblasts, stromal cells, endothelial cells, and other components [11][12][13]. TME plays a regulatory role in tumorigenesis and in the interactions between cancer cells and surrounding components (such as cancer-associated fibroblasts, T cells, B cells, macrophages, and lymphocyte), which promotes tumor development and metastasis [8,14,15].
Immune and stromal cells are two most abundant cell types and their degree of infiltration in the TME has a clinical prognostic value in many cancer types [10,16,17]. Therefore, a deeper understanding of the cell components (including immune cells and stromal cells) and the immunosuppressive status of the TME is critical for improving the efficacy of immunotherapy [18]. The Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (i.e., ESTIMATE) algorithm is used to estimate the levels of infiltrating stromal and immune cells by calculating stromal/immune/ESTIMATE scores [19]. At the same time, the therapeutic efficacy of immunotherapy in OSCC is relatively limited. Therefore, it is crucial to identify a more individualized prognostic signature and develop a more accurate immunotherapy for patients with OSCC. Recently, immune-related signatures have been developed and used to predict prognosis in certain cancers [20], and some also have been reported in OSCC (such as those based on epithelialmesenchymal transition, stemness of cancer, or immunosuppression genes) [2,17,21,22]. However, there is no reliable immune-related signature based on TME for predicting the prognosis of OSCC patients currently.
In this study, the stromal and immune scores of 319 patients from the TCGA OSCC database was determined using the ESTIMATE algorithm. A total of 562 upregulated and 31 downregulated immune-related genes were identified. Moreover, a functional enrichment analysis showed that the immune-related genes mainly played a critical role in immune response, activation/proliferation of immune-related cells, and chemokine activity. Finally, a prognostic 11 genes model, i.e., AC103563.1, CCL22 (C-C motif chemokine ligand 22), FLT3 (FMS-like tyrosine 3), GALR2 (galanin receptor 2), IGKV1D.8, IGLV1.36, IGLV4.60, IL10 (interleukin 10), LINC00861, LINC01508, and MS4A2 (membrane-spanning 4 domains, superfamily A, number 2) was constructed and confirmed to be significantly associated with the overall survival (OS) of OSCC patients. In conclusion, we screened the immunerelated signature and explored the relationship between the screened immune-related signature and the prognosis of OSCC patients, and we aimed to use this signature as a potential prognostic biomarker and immune therapeutic target for OSCC patients.

OSCC Datasets Acquisition and Handing
The RNA-sequence dataset (TCGA, OSCC, n = 319) was obtained from the National Cancer Institute GDC Data Portal (https://portal.gdc.cancer.gov/) (accessed on 3 December 2020). Patients with complete follow-up data and survival status information were included in this analysis. Immune and stromal scores were calculated using the ESTIMATE algorithm [19].

Screening of Differentially Expressed Genes (DEGs)
The package "limma" in R (version 4.0.2) was used to screen for differentially expressed genes (DEGs) between OSCC tissue and normal tissues. FDR < 0.05 and |log2FoldChange| > 1 was considered as threshold values for the identification of DEGs. A heatmap was used to visualize the expression profiles of the DEGs. Venn diagrams were used to show the overlaps among the DEGs [23].

Functional Enrichment Analyses and Functional Annotation of DEGs
Gene Ontology (GO), including biological processes (BPs), molecular functions (MFs), and cellular components (CCs) were determined using the R package "ggplot2", "cluster-Profiler" and "enrichplot", respectively. The Kyoto Encyclopedia of Genes and Genomes Vaccines 2022, 10, 1521 3 of 16 (KEGG) was used for the enrichment analysis of pathways [24]. The Search Tool for Retrieval of Interacting Genes/Proteins (STRING) (https://string-db.org) (accessed on 3 December 2020) was used to predict the functional interaction of proteins. Cytoscape software was used to analyze the interaction scores of the PPI network nodes, and scores > 0.95 were considered the key PPI nodes [25].

Construction of the Immune-Related Prognostic Signatures
OSCC patients (n = 319) were divided into high-risk and low-risk groups based on immune and stromal scores estimated by ESTIMATE. The Cox proportional hazards regression model was used to predict the OS of OSCC patients. p < 0.05 was statistically significant. Univariate Cox proportional hazards regression analyses were used to analyze the survival-related DEGs. The survival-related DEGs were put into the Cox proportional hazards model survival analysis with a least absolute shrinkage and a selection operator (LASSO) penalty. Finally, the TME-related prognostic signature was constructed by weighting the Cox regression coefficients to calculate a risk score for each patient. Using the cut-off value with the median value of the risk score, patients were classified as low-risk and high-risk.

Statistical Analysis
Univariate and multivariate Cox proportional hazards regression analyses were performed using the R package "survival". A LASSO Cox survival analysis was performed using the R package "glmnet", and a 1000-fold cross-validation was used. The R package "time-ROC" and "time AUC" were used to determine the ROC curves and area under the curve (AUC) for the predictive ability of the TME-related signatures. A heatmap was used to visualize the expression profiles of the DEGs. The alluvial diagrams were analyzed using "ggalluvial" in R. p-value < 0.05 was considered to be statistically significant and R software (New Jersey, USA, version 3.6.2) was used to perform the statistical analysis. The R code used in this study was showed in Text S1.

Stromal and Immune Scores Were Correlated to Clinical Features of OSCC Patients
The clinical data of 319 patients with OSCC downloaded from the OSCC-TCGA RNA-seq database were analyzed in the present study. To explore the relationship between stromal and immune scores and the clinical characteristics of OSCC patients, a total of 19,467 genes were extracted from the OSCC-TCGA RNA-seq database. It is well-known that stromal and immune cells are two main types of nontumor components in the TME and have been identified to be of value for the diagnosis and prognostic evaluation of tumor. Therefore, the stromal and immune scores of 319 OSCC patients were determined using the ESTIMATE algorithm, and higher immune and stromal scores in the TME resulted in more stromal and immune components in the TME. Furthermore, the ESTIMATE score is the sum of the immune score and the stromal score and represents the combined ratio of the two components in the TME. Thus, these three scores represent the proportion of stromal cells, the proportion of immune cells, and the purity of tumor in the immune microenvironment, respectively [19,26,27]. In the present study, the stromal scores ranged between −1947.438 and 1969.731, immune scores ranged between −392.987 and 2705.005, and ESTIMATE scores ranged between −2340.425 and 4674.735. The OSCC patients were categorized into high and low scores groups based on the stromal, immune, and ESTIMATE scores firstly, and a Kaplan-Meier (KM) survival analysis was used to calculate the OS (overall survival) of OSCC patients. We found that the OS of lower stromal, immune, and ESTIMATE scores of OSCC patients was relatively shorter (Figure 1a-c). Furthermore, based on the clinical data extracted from the TCGA-OSCC database, we found that the stromal, immune, and ESTIMATE scores in T3 and T4 OSCC patients were relatively lower than those in T1 and T2 OSCC patients (Figure 2a-c). Moreover, we also compared the stromal, immune, and ESTIMATE scores with the gender of OSCC patients, and we found that there was Vaccines 2022, 10, 1521 4 of 16 no significant difference between stromal, immune, and ESTIMATE scores with genders ( Figure 2d-f). We also noticed that there were no significant differences in stromal, immune, and ESTIMATE scores with tumor grades (Figure 2g-i). The above results indicate that stromal, immune, and ESTIMATE scores were relatively associated to the clinical features of OSCC patients.
1a-c). Furthermore, based on the clinical data extracted from the TCGA-OSCC database, we found that the stromal, immune, and ESTIMATE scores in T3 and T4 OSCC patients were relatively lower than those in T1 and T2 OSCC patients (Figure 2a-c). Moreover, we also compared the stromal, immune, and ESTIMATE scores with the gender of OSCC patients, and we found that there was no significant difference between stromal, immune, and ESTIMATE scores with genders (Figure 2d-f). We also noticed that there were no significant differences in stromal, immune, and ESTIMATE scores with tumor grades (Figure 2g-i). The above results indicate that stromal, immune, and ESTIMATE scores were relatively associated to the clinical features of OSCC patients.

Identification of TME-Related Differentially Expressed Genes
In order to determine the exact changes in gene expression profiles related to immune and stromal components in TME, we investigated the TME-related differentially expressed genes (DEGs) with high and low immune and stromal scores in OSCC patients with the R package "limma". A total of 1625 DEGs were obtained from the immune score (high and low score samples), with 1227 genes upregulated and 398 genes downregulated (Figure 3a,c,d). Similarly, 1589 DEGs were obtained from the stromal score, with 1506 upregulated genes and 83 downregulated genes (Figure 3b-d). A Venn diagram analysis showed that there were 562 upregulated genes shared by high scores and 31 downregulated genes shared by low scores in immune and stromal scores, and these DEGs (a total of 593 genes) might be the determinants of TME status (Figure 3c,d). These data indicated that TME-related DEGs had been obtained and identified.  ( Figure 3a,c,d). Similarly, 1589 DEGs were obtained from the stromal score, with 1506 upregulated genes and 83 downregulated genes (Figure 3b-d). A Venn diagram analysis showed that there were 562 upregulated genes shared by high scores and 31 downregulated genes shared by low scores in immune and stromal scores, and these DEGs (a total of 593 genes) might be the determinants of TME status (Figure 3c,d). These data indicated that TME-related DEGs had been obtained and identified.

Functional Annotation of TME-Related DEGs
To explore the biological functions of the above TME-related DEGs, GO and KEGG enrichment analyses were performed using the R package "cluster Profiler". As shown in Figure 4a, the top 10 GO terms of the up-and downregulated DEGs were related to different biological processes (BP), cellular components (CC), and molecular functions (MF). Interestingly, the top five terms related to the biological process included lymphocyte activation, regulation of lymphocyte-mediated immunity, phagocytosis, complement activation, and B cell-mediated immunity (Figure 4a related to molecular functions included antigen binding, immunoglobulin receptor binding, carbohydrate-binding, immune receptor activity, and cytokine receptor activity (Figure 4a, bottom panel). On the other hand, the top 30 KEGG pathways for the DEGs related to the up-and downregulated DEGs included immune/inflammation-related pathways, including B cell receptor signaling pathway, T cell receptor signaling pathway, intestinal immune network for IgA production, natural-killer-cell-mediated cytotoxicity, and Ctype lectin receptor signaling pathway (Figure 4b). The above data indicate that the signaling pathways and biological behaviors enriched by immune-related DEGs are related to the tumor microenvironment or immune functions.

Construction of the Risk Score Prognosis Model
To analyze the survival-related immune-DGEs, univariate Cox proportional hazards regression model survival analyses were performed, and a total of 23 prognostic genes were identified as survival-related DEGs: AC093278.2, AC103563.1, AL365361.1, BTLA,

Construction of the Risk Score Prognosis Model
To analyze the survival-related immune-DGEs, univariate Cox proportional hazards regression model survival analyses were performed, and a total of 23 prognostic genes were identified as survival-related DEGs: AC093278.2, AC103563.1, AL365361.1, BTLA, CCL22, CCR4, CD5, CELF2, F5, FCRL3, FLT3, GALR2, IGKV1D-8, IGLV1-36, IGLV4-60, IL10, LINC00861, LINC01508, MS4A2, P2RY14, P2RY8, PPP1R16B, and RUBCNL (Figure 7). The least absolute shrinkage and selection operator (LASSO) was then used to exclude the overfitting false positives data (Figure 7a,b). Finally, a prognostic model with 11 survivalrelated DEGs was constructed for predicting the OS of OSCC patients. The risk scores were calculated for each patient as follows: risk score = −0.0418 × (expression of AC103563. The median value was taken as the cutoff, the OSCC patients with higher risk scores than the median value were classified in the high-risk group, and the OSCC patients with lower risk scores than the median value were classified in the low-risk group.

Survival Analysis of the 11-Gene Immune-Related Signature
To test the effect of the above model on the evaluation of OSCC patients in the future, we analyzed the outcome survival (OS) of the high-risk group and the low-risk group of OSCC patients. The risk distribution, survival status, and gene expression pattern are shown in Figure 8a-c. The scatter plot (Figure 8a) shows that most of the patients in the high-risk score group died and most of the patients in the low-risk group survived during 15 years of follow-up. Moreover, the gene expression pattern (Figure 8b,c) shows that the

Survival Analysis of the 11-Gene Immune-Related Signature
To test the effect of the above model on the evaluation of OSCC patients in the future, we analyzed the outcome survival (OS) of the high-risk group and the low-risk group of OSCC patients. The risk distribution, survival status, and gene expression pattern are shown in Figure 8a-c. The scatter plot (Figure 8a) shows that most of the patients in the high-risk score group died and most of the patients in the low-risk group survived during 15 years of follow-up. Moreover, the gene expression pattern (Figure 8b,c) shows that the high-and low-risk groups we categorized had corresponding risk scores. The heatmap (Figure 8d) showed that seven immune-related DEGs (AC103563.1, CCL22, FLT3, IGLV4.60, IL10, LINC00861, and MS4A2) were highly expressed in the low-risk group while four immune-related DEGs (GALR2, LINC01508, IGKV1D.8 and IGLV1.36) were highly expressed in the high-risk group. To determine the relationship between immune-related DEGs and OS in OSCC patients, we used the Kaplan-Meier method to plot the survival curves using data obtained from the OSCC-TCGA database. The Kaplan-Meier plots showed that patients in the high-risk score group had a significantly poorer OS than those in the low-risk score group (Figure 8e). Moreover, the AUC for 1-year, 2-year, and 3-year PFS were 0.64, 0.63, and 0.65, respectively (Figure 8f). Thus, these data suggested that the 11-gene immune-related signature performed well for predicting OS in OSCC patients. survival curves using data obtained from the OSCC-TCGA database. The Kaplan-Meier plots showed that patients in the high-risk score group had a significantly poorer OS than those in the low-risk score group (Figure 8e). Moreover, the AUC for 1-year, 2-year, and 3-year PFS were 0.64, 0.63, and 0.65, respectively ( Figure 8f). Thus, these data suggested that the 11-gene immune-related signature performed well for predicting OS in OSCC patients.

Validation of an External GEO Cohort
To verify the clinical value of the 11-gene TME-related signature for predicting OS in OSCC patients, we used an external GEO-OSCC cohort (GSE65858) to validate our study. The 11-gene TME-related signature was constructed in GSE65858, and the risk scores were analyzed by the same methods. High-and low-risk OSCC groups were classified according to the median risk score. A Kaplan-Meier curve analysis showed that the low-risk score was closely associated with a longer overall survival time ( Figure 9D), which was consistent with that the TCGA cohort (Figure 8d). In addition, the ROC curve showed an

Validation of an External GEO Cohort
To verify the clinical value of the 11-gene TME-related signature for predicting OS in OSCC patients, we used an external GEO-OSCC cohort (GSE65858) to validate our study. The 11-gene TME-related signature was constructed in GSE65858, and the risk scores were analyzed by the same methods. High-and low-risk OSCC groups were classified according to the median risk score. A Kaplan-Meier curve analysis showed that the low-risk score was closely associated with a longer overall survival time (Figure 9d), which was consistent with that the TCGA cohort (Figure 8d). In addition, the ROC curve showed an AUC value was 0.55 at 1, 2, and 3 years (Figure 9e), which was consistent with that the TCGA cohort (Figure 8e). Thus, the results showed that the 11-gene TME-related signature could effectively predict the prognosis of OSCC patients partly.

Discussion
In the present study, differentially expressed immune-related genes in oral squamous carcinoma cancer (OSCC) samples were identified and were used to construct a signature with 11 immune-related genes, which performed well in predicting the outcomes of OSCC patients, respectively.
OSCC ranks as the eighth most common form of oral malignancy, with an estimated 500,000 new cases being reported annually, and OSCC is associated with a high grade, rapid progression, poor treatment effects, and bad outcomes [28]. The clinical treatment of OSCC is often a combination of surgical resection with chemotherapy; however, the effect of these treatments is very limited and the 5-year survival rate of OSCC patients remains lower than 50% [29,30]. There are several ways to predict the prognosis of OSCC patients, including but not limited to the following: (1) TNM stages; (2) tumor grade; (3) the number of lymph node metastasis; (4) the expression level of key proteins, such as epidermal growth factor receptor and P53 [31][32][33]. However, a unified and effective standard still lacks to determine the prognosis of OSCC patients.

Discussion
In the present study, differentially expressed immune-related genes in oral squamous carcinoma cancer (OSCC) samples were identified and were used to construct a signature with 11 immune-related genes, which performed well in predicting the outcomes of OSCC patients, respectively.
OSCC ranks as the eighth most common form of oral malignancy, with an estimated 500,000 new cases being reported annually, and OSCC is associated with a high grade, rapid progression, poor treatment effects, and bad outcomes [28]. The clinical treatment of OSCC is often a combination of surgical resection with chemotherapy; however, the effect of these treatments is very limited and the 5-year survival rate of OSCC patients remains lower than 50% [29,30]. There are several ways to predict the prognosis of OSCC patients, including but not limited to the following: (1) TNM stages; (2) tumor grade; (3) the number of lymph node metastasis; (4) the expression level of key proteins, such as epidermal growth factor receptor and P53 [31][32][33]. However, a unified and effective standard still lacks to determine the prognosis of OSCC patients.
The TME is the cellular environment in which tumor cells live, which consists of an extracellular matrix, soluble molecules, and tumor stromal cells [10]. The TME plays a vital role in OSCC progression and consists of tumor cells and other components, such as immune cells, stromal scores, tumor-associated macrophages (TAMs), and carcinomaassociated fibroblasts (CAFs) [34,35]. Moreover, tumor development, progression, and responses to immunotherapies are regulated by cytokines, immune infiltration, TAMs, and CAFs within the TME [36,37]. It is well known that immune cells and stromal cells are two main types of nontumor components in the TME and have been proposed to be of value for the diagnosis and prognostic evaluation of tumors [38]. Recently, ESTIMATE has been used as a tool for inferring immune, stromal, and immune scores in the TME, and the estimation is based on the expression signals of gene sets that characterize stromal and immune cells. Therefore, in this study, we developed an immune-related signature based on stromal and immune scores to predict the prognosis of OSCC patients, which could improve the response rate to immunotherapy. In addition, ESTIMATE is an algorithm that has been used to calculate the immune and stromal scores in several cancers, including pancreatic ductal adenocarcinoma [39], colorectal cancer [40], and OSCC [17,21,22].
In the present study, we conducted a comprehensive analysis between the immune cells and immune-related genes and the clinical indexes and outcomes of OSCC patients. We found the following findings: (1) The clinical indexes (such as TNM stages, grade stages, and genders) of OSCC patients were associated with the stromal, immune, and ESTIMATE scores, but the correlation was not significant (Figure 1). However, high stromal, immune, and ESTIMATE scores correlated with a relatively longer OS of OSCC patients, indicating that the TME components, especially the stromal and immune cells affected the clinical outcomes of OSCC patients ( Figure 2). (2) In total, 593 DEGs were identified from the comparison of high versus low stromal and immune scores groups based on the median value of immune and stromal scores as the standard cutoff to divide patients into high-score and low-score groups. Then, the results of GO terms showed that the DEGs in the high-score group regulated the T cell receptor complex and immunoglobulin complex ( Figure 3a). Moreover, GO and KEGG pathway enrichment data indicated that the DEGs in the high-score group also regulated the T cell receptor signaling pathway and viral myocarditis (Figure 3b). (3) LASSO and random forest (RF) algorithm analyses identified the 11 prognostic immune-related genes (AC103563.1, CCL22, FLT3, GALR2, IGKV1D.8, IGLV1.36, IGLV4.60, IL10, LINC00861, LINC01508, and MS4A2) (summarized in Tables S1 and S2). In particular, high levels of GALR2, LINC01508, IGKV1D.8, and IGLV1.36 were found to be negatively correlated with the OS of OSCC patients, while high levels of AC103563.1, CCL22, FLT3, IGLV4.60, IL10, LINC00861, and MS4A2 were positively correlated with the OS of OSCC patients. (4) To verify the clinical value of the 11-gene immune-related signature for predicting OS in OSCC patients, we used an external GEO-OSCC cohort (GEO65858) for the validation in our study. Altogether, the results showed that the 11-gene TME-related signature could effectively predict the prognosis of OSCC patients.
The clinical implications of these 11 key genes have been reported in different types of cancers. For instance, CCL22 (C-C motif chemokine 22) controls T cell immunity by recruiting regulatory T cells to the tumor tissue and promoting regulatory T cell communication with dendritic cells in TME [41]. Higher infiltration rates of CCL22+ cells are reported to be associated with poor outcomes in cervical cancer patients [42]. Moreover, studies have shown the efficacy of flt3-L and CD40-L combination immunotherapy on prostate tumor growth in TME [43]. Moreover, FLT3 is a potential targeted protein and has been reported as an immune-related prognostic gene in lung adenocarcinoma [44], indicating the clinical value of FLT3. Furthermore, galanin receptor 2 (GALR2) is a G protein-coupled receptor that induces tumor growth and proliferation in SCCHN (squamous cell carcinoma of the head and neck) [45], and GALR2 has been identified as a novel target and biomarker for prostate, colon, and breast cancer screening [46]. In addition, IL10 is an important immune regulatory cytokine in TME [47], and a recent study has shown that IL10 promotes cancer cell metastasis and proliferation via immunosuppression, and also has a predictive value in clear cell renal cell carcinoma [48]. Recently, a study has shown that the MS4A2 (membrane-spanning 4 domains, superfamily A, number 2) gene is located on chromosome 11q13, a region that is linked to asthma-related phenotypes, and it plays an important role in the regulation of human mast cell proliferation and survival [49]. However, very few studies report the role of MAS4A2 in TME. Other immune-related genes are relatively novel and scattered evidence has been reported for their roles in human diseases. For example, IGLV1-36 gene is a potentially therapeutic option for PO-EMS (polyneuropathy, organomegaly, endocrinopathy, monoclonal gammopathy, and skin changes) syndrome [50]. Moreover, the downregulation of LINC01508, a long noncoding RNAs, contributes to cisplatin resistance in ovarian cancer [51]. Additionally, another recent study has shown that the LINC00861/miR-513b-5p axis inhibits cancer cell progression through the PTEN/AKT/mTOR signaling pathway in cervical cancer cells [52]. Furthermore, LINC00861 is also used as a biomarker to predict survival in patients with ovarian cancer [53], the early diagnosis of Parkinson's disease [54], and breast cancer [55]. However, few studies have reported the role of AC103563.1, IGKV1D.8, and IGLV4.60 in TME (especially in OSCC), which will be our future directions. Similar studies have also been shown in OSCC. For instance, Li et al. selected eleven immunosuppression genes (ISGs), including INHBB, BGLAP, CTLA4, CALCA, CXCL8, IL22, FGFR3, HPRT1, ORMDL3, SPHK1, and TLR3, which showed a prognostic potential for OSCC [21]. In this study, a deep-learning-based model showed two subgroups of OSCC samples, while subtype Sub1 displayed a more aggressive phenotype with poor prognosis, with immune-cells-associated pathways and enriched cancer-progression-pertinent pathways. In addition, to test if the stemness and immune-relevant genes were involved in oral cancer, Lin et al. established an eight-gene risk model (ESCO2, CCNA2, COL5A3, RCN3, LMCD1, FMNL3, MMP14, and HEYL), which performed well in predicting overall survival and recurrence-free survival in OSCC patients [17]. Further investigations showed that the eight-gene signature was highly linked to immune suppression. However, we used different methods to independently establish our gene signature, including a proteinprotein interaction analysis, GO and a KEGG enrichment analysis. More importantly, our 11-gene signature involved either LincRNA and a few novel genes, which could potentially contribute to oral cancer progression and treatment resistance in the TME. Altogether, our studies and others have all provided potentially promising ways to obtain a prognosis for OSCC patients.
However, our study still has limitations. For example, the OSCC TME-related signature was not verified by biological experiments both in vivo and in vitro, but in our ongoing projects, these immune-related genes and their potential mechanisms in OSCC are being verified. Moreover, due to the limited samples in the OSCC subgroups, the prognostic value of the immune-related signature in some OSCC subgroups did not show any statistical significance. Therefore, a large cohort of OSCC patient samples is needed for future validations.

Conclusions
In summary, we proposed a signature of 11 immune-related genes based on the TME in OSCC, which could be used as an independent prognostic biomarker for OSCC patients. The signature of 11 immune-related genes had a predominant performance in obtaining the prognosis of OSCC patients, and it could achieve a more personalized and precise immunotherapy effect in OSCC.