Single-cell RNA sequencing identifies cancer-associated fibroblast marker genes for determining cancer subtypes in oral squamous cell carcinoma and predicting patient prognosis

Introduction

The incidence and mortality rates of oral cancer, one of the most common malignancies worldwide, have increased in recent years. According to surveys, 377,713 novel cases of oral cancer and 177,757 oral cancer-related fatalities were recorded globally in 2020 (1). Oral squamous cell carcinoma (OSCC) is a typical subtype of oral cancer that accounts for approximately 90% of all oral malignancies. Although primary tumor resection is currently the conventional treatment for OSCC (2), its highly aggressive nature leads to postoperative recurrence and metastasis in most post-surgically treated patients (3), resulting in a suboptimal 5-year overall survival (OS) in patients with advanced OSCC (4). Although targeted therapy and immunotherapy have evolved as crucial modalities for cancer treatment in recent years, their application in OSCC has been impeded by individual patient variance and acquired resistance (5,6). Hence, it is crucial to investigate new targets and therapies that can effectively impede tumor progression and increase the survival time of OSCC patients.

OSCC is a solid tumor primarily composed of tumor cells and intricate tumor microenvironment (TME) components. Among these, cancer-associated fibroblasts (CAFs) represent the most dominant TME component and exert a crucial influence in promoting tumor advancement (7). Studies have demonstrated that CAFs remodel stromal components through influencing the production and deposition of the extracellular matrix (ECM). CAFs also facilitate the formation of an immunosuppressive and invasive microenvironment for OSCC cells (8-10). Moreover, CAFs mediate resistance to chemotherapy for tumor cells. For instance, CAFs can transfer miR-196a into head and neck squamous cell carcinoma (HNSCC) cells via exosomes, thereby downregulating CDKN1B and ING5 expression and conferring cisplatin resistance in tumor cells (11). Given the pivotal role of CAFs in the tumor immune microenvironment (TIME) and chemotherapy resistance, developing therapies targeting CAFs appears promising (12). CAFs are a heterogeneous stromal cell population (13,14). Therefore, novel therapeutic strategies targeting CAFs require comprehensive understanding of cell biology features of CAF and identification of specific CAF marker genes.

With the emergence of single-cell RNA sequencing (scRNA-seq) technology, the intricate and heterogeneous nature of tumor tissue composition can be probed at the cellular level. This state-of-the-art technology detects gene expression profiles of cellular subpopulations and enables the examination of intercellular communication at single-cell resolution (15,16). In addition, it permits the identification of stromal cell phenotype-specific marker genes (17,18), providing a new direction for promising cancer therapies.

Through scRNA-seq analysis, we have identified CAF marker genes in this study. We integrated data on OSCC from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases and subsequently classified 338 samples into two separate CAF subtypes based on prognosis-associated CAF marker genes. We also analyzed the complex TME of different subtypes and explored potential regulatory mechanisms. Furthermore, we constructed a predictive model for patients’ OS by utilizing the differentially expressed genes (DEGs) between the two CAF subtypes. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-2129/rc).

Methods

Data acquisition

We acquired the scRNA-seq dataset (GSE103322) for OSCC from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database. Similarly, we downloaded the transcriptomic data, somatic mutation data, and clinical data of 243 OSCC samples from the TCGA (https://portal.gdc.cancer.gov/) database and converted the normalized expression data from TCGA to transcripts per million (TPM). Additionally, we acquired the transcriptomic data and clinical data for 97 OSCC samples from GSE41613 in the GEO database. The microarray data (GSE41613) were normalized using the limma R package with the RMA method for background correction and log₂ transformation. The TCGA RNA sequencing (RNA-seq) data [fragments per kilobase of transcript per million mapped reads (FPKM)] were converted to log₂(TPM +1) values. To integrate the datasets, we retained only genes common to both platforms. Batch effects, primarily due to the different experimental platforms (RNA-seq vs. microarray), were removed using the ComBat algorithm from the “sva” R package (19), with the data source (TCGA or GEO) specified as a known batch variable. Finally, we obtained 338 OSCC samples with complete transcriptome and survival data for further analysis (table available at https://cdn.amegroups.cn/static/public/tcr-2025-2129-1.docx).

Analysis of scRNA-seq data

The analysis focused on primary tumor samples from all 18 patients in the GSE103322 dataset. Cells annotated as originating from lymph node metastases were excluded. Data processing was performed using R and the Seurat software package (20). Initial quality control filtered out cells expressing fewer than 200 genes and genes detected in fewer than three cells. Cells with a mitochondrial gene percentage exceeding 5% were also excluded. Gene expression counts were log-normalized using the LogNormalize method (scale factor =10,000). The top 1,500 highly variable genes were identified for downstream analysis. Dimensionality reduction was performed using principal component analysis (PCA), and significant principal components were selected using the JackStraw tool (19,21). Cell clusters were identified based on a shared nearest neighbor (SNN) modularity optimization algorithm. Marker genes for each cluster were identified using the FindAllMarkers function. A parent fibroblast cluster was identified based on canonical markers (e.g., FAP, PDGFRA), and its marker genes were extracted for subsequent analyses. Cluster annotation was performed using the SingleR package with reference to cell-type-specific marker genes (22). Interactions among distinct cell clusters were analyzed using the CellChat package (23).

Enrichment analysis based on CAF markers

Functional enrichment analysis of CAF gene markers was evaluated using the clusterprofiler R package, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, and the gene set variation analysis (GSVA) R program was utilized to identify pathways most closely associated with CAF markers (24), with P values <0.05 being considered statistically significant.

Identification of CAF subtypes

By conducting univariate Cox regression analysis, we determined 24 CAF marker genes relevant for prognosis. We utilized the expression profiles of these 24 genes to classify 338 patients with OSCC with the assistance of the “ConsensusClusterPlus” software package (25). t-distributed stochastic neighbor embedding (t-SNE) and PCA analyses were performed to validate the classification results. Our study compared survival differences among subtypes using the Kaplan-Meier analysis. Additionally, we examined TME heterogeneity and the relative enrichment scores of immune cell signatures between subtypes using the ESTIMATE method and single-sample gene set enrichment analysis (ssGSEA), respectively (26,27). To assess the level of CAF infiltration in each subtype, we used the Tumor Immune Dysfunction and Exclusion (TIDE) computational method (http://tide.dfci.harvard.edu/) (28). Furthermore, we collected immune checkpoint scores for CTLA4 and PD1 in OSCC samples that were retrieved from The Cancer Immunome Atlas (TCIA) ‘Head and Neck Squamous Cell Carcinoma’ dataset to assess the variance in immunotherapy efficacy across diverse subtypes (29). Finally, we explored the potential regulatory mechanisms of each subtype using the “GSVA” package.

Construction and evaluation of prognostic models

The 338 patients with OSCC were randomly divided into training and test groups, with equal sample sizes in each group (1:1). Prognostic DEGs associated with subtypes were analyzed using the least absolute shrinkage and selection operator (LASSO) and multivariate Cox regression methods. Eight crucial genes were selected to create the prognostic models for the training group. The OSCC patients were classified into high- and low-risk subgroups and scored according to the risk scoring formula. Riskscore=∑i=1n(Coefi×Expi) or risk score = coef gene₁ × exprgene₁ + coef gene₂ ×exprgene₂ + ... + coef gene_n × exprgene_n (30). Kaplan-Meier analysis was used to assess the survival differences between the two subgroups. The model’s accuracy for predicting patients’ OS was evaluated using the area under the receiver operating characteristic (ROC) curve (AUROC), and the risk scores among various CAF subtypes were compared. Survival analysis was performed to further validate the model’s applicability by combining risk scores with clinical features and evaluating the relationship between clinical characteristics and risk scores. Both univariate and multivariate Cox regression analyses were conducted to determine whether the risk score had an independent prognostic value for OSCC. Finally, a nomogram was created by integrating the risk score to forecast patient survival rates for up to 5 years, and the reliability of the results obtained from the nomogram was assessed by using calibration curves.

Differences in immune microenvironment between risk subgroups

We used ssGSEA to analyze the variance in the enrichment scores of immune cell signatures within the TIME between the two risk subcategories of patients. Furthermore, we employed the Spearman method to assess the correlation of risk scores with immune cells. To evaluate the anti-tumor immune function score of each patient, we used the Tracking Tumor Immunophenotype (TIP) meta-server (31) and compared variations in anti-tumor immune function between the two patient groups.

Somatic mutation analysis and sensitivity analysis of immunotherapy

We employed the “Maftools” package (32) to create somatic mutation landscapes for both risk subgroups. The tmb() function was used to calculate the tumor mutation burden (TMB) scores for every patient with OSCC, and the differences in TMB were compared between the two subgroups. The investigation of the correlation between the TMB and risk scores utilized the Spearman approach. Additionally, we used the TIDE score to predict patients’ responses to immunological checkpoint blockade therapy. The TIDE score was calculated using the publicly available TIDE web tool, which integrates signatures of T-cell dysfunction and exclusion to predict immune checkpoint blockade response. A higher TIDE score indicates a greater likelihood of immune evasion and thus a poorer predicted response to immunotherapy. We conducted a Chi-squared test to explore the variance in response rates between patients in the two risk subgroups (33).

Clinical samples

We collected 30 pairs of cancer and paracancerous tissues from OSCC patients who underwent surgery at the Oral and Maxillofacial Surgery Department of Peking University Shenzhen Hospital, Shenzhen, China. Patients’ clinicopathologic features were shown in Table S1. None of the patients had received any chemotherapy prior to surgery. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethical Committee of Peking University Shenzhen Hospital (No. 2022-117) and all patients provided written informed consent.

RNA extraction and quantitative real-time polymerase chain reaction (qRT-PCR)

Total RNAs were purified by applying 1 mL TRIzol reagent (Takara, Dalian, China) for 50 mg tissues, which were reverse transcribed into complementary DNA (cDNA) using a reverse transcription kit (Accurate Biology, Changsha, China). The expression levels of eight key genes were detected by qRT-PCR with SYBR Green PCR Master Mixes (Accurate Biology). The 2^−ΔΔCt method was applied to calculate the relative quantification, and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was selected as the internal control.

Statistical analysis

We utilized R version 4.2.1 and its affiliated packages for all analyses. For comparison between two groups, if the normal distribution is met, a t-test is utilized. If it is a non-normal distribution, the nonparametric test is used. We conducted correlation analysis using Spearman’s rank correlation method. The log-rank test, Cox regression, and Kaplan-Meier curves were applied to evaluate prognostic value. The Wilcoxon rank sum test was applied to evaluate the relationship between immune checkpoints, TMB, and immune cell infiltration between these two groups. All statistical analysis was two-sided, and P<0.05 was considered statistically significant.

Results

CAFs involved in regulating the immune microenvironment in OSCC

After strict quality control, we obtained 4,540 high-quality single-cell samples from GSE103322 and performed PCA. We then used the t-SNE algorithm to cluster and visualize the first 20 principal components and successfully classified these cells into 19 subclasses (Figure 1A). These cell subclasses were subsequently annotated using the “SingleR” R package, The result showed that these 19 cell subclasses were annotated into eight major cell types, including T cells, tissue stem cells, epithelial cells, fibroblasts, endothelial cells, B cells, common myeloid progenitor (CMP), and monocytes which confirmed the high heterogeneity of OSCC tissues (Figure 1B). Moreover, we validated the reliability of cluster 7 as fibroblasts based on the expression levels of CAF marker genes (22,26). FAP, PDGFRA, PDGFRB, SPARC, FGF7, MME, and PDPN (Figure 1C), and a bubble plot showing the marker genes for all annotated cell types is provided in Figure S1. We then conducted enrichment analyses using GO and KEGG methods on 320 markers from cell cluster 7 (Table S2) to investigate their biological functions and pathways. Those marker genes were demonstrated playing significant roles in production and remodeling of components of ECM (Figure 1D). Furthermore, they were highly focused on the PI3K-Akt signaling pathway, focal adhesion, ECM-receptor interactions, and regulation of the actin cytoskeleton (Figure 1E). Importantly, these biological processes are closely associated with fibroblast properties, further confirming the appropriate annotation of this cluster. Finally, we analyzed intercellular communication using the CellChat package. The overall communication network among eight major cell clusters, based on ligand-receptor interaction strength, is shown in Figure 1F. Specifically, interactions where fibroblasts act as ligand cells with monocytes and T cells exhibited higher frequency and intensity compared to other interactions.

Figure 1 scRNA-seq analysis of 4,540 cells from OSCC tissues. (A) A total of 4,540 cells were classified into 19 subclasses using the t-SNE method. (B) Nineteen cell subclasses were annotated into eight major cell types. (C) Expression levels of recognized CAF marker genes in each cell type. Functional enrichment analysis of 320 marker genes in fibroblast clusters. (D) GO enrichment analysis. (E) KEGG pathway enrichment analysis. (F) Strength of interactions between cell types. BP, biological process; CAF, cancer-associated fibroblast; CC, cellular component; CMP, common myeloid progenitor; ECM, extracellular matrix; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; OSCC, oral squamous cell carcinoma; scRNA-seq, single-cell RNA sequencing; t-SNE, t-distributed stochastic neighbor embedding.

Consensus clustering identified CAF-high and CAF-low OSCC subtypes

We performed univariate Cox regression analysis on the 320 CAF marker genes and recognized 24 CAF marker genes associated with the prognosis of OSCC (Table S3). Based on the expression profiles of these 24 genes, we conducted a consensus clustering analysis of 338 OSCC samples to gain comprehensive insights into the significance of CAF marker genes in driving OSCC progression. Our analysis identified two distinct OSCC subtypes: cluster A (n=203) and cluster B (n=135) (Figure 2A). We further validated the clustering results using the PCA and t-SNE algorithms and confirmed that the two clusters were well separated and distinct (Figure 2B,2C).

Figure 2 Identification of two CAF subtypes in OSCC. (A) Consensus clustering analysis of two clusters (k=2). (B) PCA and (C) t-SNE indicated distinct transcriptional profiles between the two subtypes. (D) Heat map of 320 CAF marker genes expression between the two clusters; yellow represents high expression, and blue represents low expression. (E) Comparison of the abundance of CAF infiltration between the two subtypes. (F) Kaplan-Meier curve for OS differences between the CAF-high subtype and the CAF-low subtype. CAF, cancer-associated fibroblast; OS, overall survival; OSCC, oral squamous cell carcinoma; PC, principal component; PCA, principal component analysis; t-SNE, t-distributed stochastic neighbor embedding; TCGA, The Cancer Genome Atlas.

Furthermore, we conducted a differential expression analysis of CAF marker genes for these two subtypes (Figure 2D) and found that cluster B exhibited significantly higher level of CAF marker gene expression than cluster A. Next, we evaluated and compared the CAF infiltration scores of the two subtypes using the TIDE method and found that cluster B had a high CAF infiltration score (Figure 2E). Consequently, cluster B was classified as the CAF-high subtype and cluster A as the CAF-low subtype. We compared the survival outcomes of the two subtypes using the Kaplan-Meier method and discovered that patients with the CAF-high subtype demonstrated a notably worse prognosis (Figure 2F), and increased CAF infiltration was strongly linked to worse OS among the 338 patients with OSCC (Figure S2). These results revealed that CAF-related OSCC can be distinguished into the CAF-high and CAF-low subtypes, and CAF-high OSCC patients had worse prognosis outcomes than the CAF-low group.

CAF-high subtype presented immunosuppressive microenvironment and resulted in poor prognosis in OSCC patients

We utilized the ESTIMATE to analyze the TME heterogeneity among the subtypes. We found that the CAF-high subtype has significantly elevated immune and stromal scores than the CAF-low subtype (Figure 3A,3B). While the CAF-high subtype has significantly reduced tumor purity than the CAF-low subtype (Figure 3C). Additionally, we utilized ssGSEA to analyze the variance in the enrichment scores of 23 immune cell signatures between the two subtypes. The results revealed that the CAF-high subtype had significantly higher enrichment scores for most immune cell signatures than the CAF-low subtype. Furthermore, the enrichment scores for signatures of immunosuppressive cells, such as myeloid-derived suppressor cells (MDSCs) and regulatory T cells (Tregs) were also notably increased in the CAF-high subtype, whereas the enrichment score of neutrophils significantly decreased (Figure 3D). We identified many immune checkpoint genes that were notably upregulated in the CAF-high subtype, as visualized by the boxplots (where points denote outliers). Conversely, the opposite pattern was observed in the CAF-low subtype (Figure 3E). Moreover, GSVA revealed that pathways linked to immunosuppression and cancer promotion were notably enriched in the CAF-high subtype, including the inflammatory response, IL6-JAK-STAT3 and TGF-β signaling pathways, epithelial-mesenchymal transition (EMT), and angiogenesis. Pathways associated with metabolisms, such as the reactive oxygen pathway, fatty acid metabolism, oxidative phosphorylation, and the P53 pathway, were primarily enriched in the CAF-low subtype (Figure 3F).

Figure 3 Heterogeneity of TME between two subtypes. Comparison of immune scores (A), stromal scores (B), and tumor purity (C) between CAF-high and CAF-low subtypes. (D) Differences in the infiltration of 23 immune cells between the two subtypes. (E) Differential expression of immune checkpoint genes between CAF-high and CAF-low subtypes (points represent outliers). (F) GSVA of differentially enriched pathways between CAF-high and CAF-low subtypes. CAF, cancer-associated fibroblast; GSVA, gene set variation analysis; MDSC, myeloid-derived suppressor cell; TCGA, The Cancer Genome Atlas; TME, tumor microenvironment.

We used TCIA and TIDE approaches to assess the sensitivity of immune therapy in the two subtypes for better understanding the impact of the immunosuppressive microenvironment on immune therapy in patients. We compared the effects of the two subtypes on immune checkpoint blockade therapy by predicting immunophenotype score (IPS) under four scenarios reflecting different expression states of CTLA4 and PD-1. Based on our results (Figure 4A-4D). It can be inferred that the CAF-high subtype patients had a worse response to the CTLA4 therapy approach than the CAF-low subtype patients. Additionally, the TIDE method results (Figure 4E,4F) indicated that the CAF-high subtype exhibited a higher TIDE score than the CAF-low subtype, and a less proportion of CAF-high patients positively responded to immunotherapy than the CAF-low subtype patients.

Figure 4 Sensitivity of immune therapy in the two subtypes. (A-D) Differential sensitivity (IPS) of the two subtypes to immune checkpoint blockade therapies under four simulated conditions: targeting neither CTLA4 nor PD-1 (CTLA4_neg/PD-1_neg), targeting PD-1 only (CTLA4_neg/PD-1_pos), targeting CTLA4 only (CTLA4_pos/PD-1_neg), and targeting both (CTLA4_pos/PD-1_pos). (E,F) Differences in TIDE scores and response rates between the two subtypes. IPS, immunophenotype score; neg, negative; pos, positive; TIDE, Tumor Immune Dysfunction and Exclusion.

Prognostic model constructed and verified for OSCC patients based on CAF subtypes

Based on the two subtypes, we identified 1,782 DEGs. Subsequently, 240 prognostic DEGs were found by univariate Cox regression analysis. Using LASSO and multivariate Cox regression analyses, we selected eight crucial genes to construct OSCC prognostic model (Table 1). While we validated the expression of these eight genes in 30 pairs of OSCC tissues and paracancerous tissues collected from Peking University Shenzhen Hospital. The results showed that ADAMTS15, THBS1, CCR7, and IL13RA2 were significantly up-regulated in tumor tissues. While SERPINB7, SLFN11, and TSPAN11 showed an upward trend and AQP1 showed a downward trend, consistent with the bioinformatics analysis direction, these differences did not reach statistical significance in this validation cohort (Figure 5A-5H).

Figure 5 Eight OSCC prognostic genes expression validation. (A-H) The relative expression of eight OSCC prognostic key genes, ADAMTS15, THBS1, AQP1, CCR7, IL13RA2, SERPINB7, SLFN11, and TSPAN11, in 30 pairs of cancer and paracancerous tissues of OSCC patients. GAPDH, glyceraldehyde-3-phosphate dehydrogenase; mRNA, messenger RNA; OSCC, oral squamous cell carcinoma.

Using half of OSCC samples (169/338) as the training group, a risk signature was created utilizing these eight crucial genes. The patients were categorized into high-risk or low-risk groups based on the median risk score of the training group. Kaplan-Meier analysis revealed that the high-risk group had notably reduced OS than the low-risk group (Figure 6A). And the number of deaths was increasing as the risk score increased (Figure 6B). The model’s accuracy in predicting OS was assessed using the ROC curve. The results indicated that the AUC at 1, 3, and 5 years were 0.783, 0.799, and 0.793 correspondingly (Figure 6C), demonstrating the model’s good predictive power. This analysis yielded consistent results in the test group (another half of OSCC samples, n=169) (Figure 6D-6F). To assess the model’s applicability, we conducted a survival analysis by merging the risk score with the clinical characteristics. The findings showed that all clinical subgroups of patients with high risk had worse OS than the low-risk groups, demonstrating the model’s good applicability (Figure S3). Furthermore, using the combined cohort, we compared the risk scores across the clinical feature subgroups and found that patients with a higher tumor stage and more lymph node metastases had significantly higher risk scores (Figure 6G-6I). Subsequently, by comparing the risk scores for both CAF subtypes, we found that the CAF-high subtype patients had significantly higher risk scores than the CAF-low subtype (Figure 6J). And a notably positive relationship existed between the risk score and the level of CAF infiltration in OSCC (Figure 6K,6L). Finally, we combined CAF subtypes with risk scores for survival analysis and found that patients with characteristics of both the high-risk score and the CAF-high subtype had the poorest prognosis (Figure S4).

Figure 6 Construction and evaluation of eight gene signature prognostic models in OSCC. (A-C) Kaplan-Meier survival analysis, risk-score distribution map, and time-dependent ROC analysis in the training cohort. (D-F) Kaplan-Meier survival analysis, risk-score distribution map, and time-dependent ROC analysis in the test cohort. (G-I) Risk scores of patients with OSCC in different clinical feature groups (combined cohort). (J) Comparison of risk scores between the two CAF subtypes. (K) Differences in CAF infiltration scores between the two risk groups. (L) Positive correlation map between the risk score and CAF infiltration level. AUC, area under the curve; CAF, cancer-associated fibroblast; CI, confidence interval; HR, hazard ratio; N, node; OSCC, oral squamous cell carcinoma; ROC, receiver operating characteristic; T, tumor.

Based on our analysis of the complete sample group, both univariate and multivariate regression analyses revealed that the tumor stage and risk score were independent prognostic indicators for OSCC (Figure 7A,7B). To enhance the model’s clinical application, we devised a nomogram that integrates tumor stage and risk score to forecast the OS of patients with OSCC over periods of 1-, 3-, and 5-year (Figure 7C). Furthermore, the nomogram-predicted outcomes and actual survival results demonstrated a significant level of consistency, as indicated by the calibration curve (Figure 7D).

Figure 7 Performance evaluation of prognostic models. (A,B) Forest plot results of univariate and multivariate regression analyses of risk score. (C) Risk score combined with tumor staging to create a nomogram. (D) The calibration chart demonstrates the nomogram’s predictive capacity. (E-I) ROC curves comparing the performance of our prognostic model with four CAF-derived gene feature models previously reported in literature studies at 1, 3, and 5 years. (J) Concordance index comparing the accuracy of our prognostic model with the four CAF-derived gene feature models in predicting patient OS. ***, P<0.001. AUC, area under the curve; C-index, concordance index; CAF, cancer-associated fibroblast; CI, confidence interval; OS, overall survival; ROC, receiver operating characteristic.

Furthermore, we compared our prognostic model with four models based on CAF-related gene signatures previously reported in the literature for HNSCC and OSCC (34-37). Based on the complete sample group, the ROC curves (Figure 7E-7I) indicated that our model demonstrated a greater AUC value than the other four models at 1, 3, and 5 years, implying that our model possesses a stronger capacity to forecast patient OS than other methods. The concordance index demonstrated that our signature was the most accurate in predicting the patient’s OS (Figure 7J). The above results, based on internal comparative analysis within the cohort, suggest that our prognostic model displays promising performance. However, validation in an independent external cohort is warranted to further confirm its generalizability and comparative advantage. Nonetheless, these findings indicate its potential to function as an independent prognostic indicator for patients with OSCC.

The high- and low-risk OSCC patients exhibited notable variations in TIME and anti-tumor immune function

The enrichment scores of 23 different immune cell signatures and 29 immune-related pathways were assessed in patients using the ssGSEA algorithm (Figure 8A,8B). According to our analysis, the low-risk subgroup exhibited higher enrichment scores for signatures of neutrophils, mast cells, activated B cells, and CD8⁺ T cells, whereas the high-risk subgroup exhibited a notably higher enrichment score for the natural killer T cell signature. Furthermore, the high-risk subgroup displayed significantly lower scores for immune-related pathways such as T cell co-stimulation and type II interferon (IFN) responses. Correlation analysis revealed that the risk score was negatively connected with the enrichment scores of several immune cell signatures (Figure 8C). Since the TIME was highly associated with anti-tumor immune status, we quantified the anti-tumor immune function scores of each patient using the TIP meta-server and compared the activities of seven continuous anti-tumor immune processes between the two different risk groups (Figure 8D). According to our findings, the high-risk subgroup had considerably lower activity levels in the third steps (priming and activation) and the fifth steps (infiltration of immune cells into tumors) than the low-risk subgroup.

Figure 8 Analysis of TIME and anti-tumor immune function in high- and low-risk groups. (A,B) Differences in immune cell types and immune-related pathways between the high- and low-risk groups. (C) Heatmap of the correlation between risk score and immune cell enrichment score. (D) Comparison of anti-tumor immune process activities between high- and low-risk groups. aDCs, activated dendritic cells; APC, antigen-presenting cell; CCR, C-C chemokine receptor; DCs, dendritic cells; HLA, human leukocyte antigen; iDCs, immature dendritic cells; IFN, interferon; MDSC, myeloid-derived suppressor cell; MHC, major histocompatibility complex; NK, natural killer; pDCs, plasmacytoid dendritic cells; Tfh, T follicular helper; Th, T helper; TIL, tumor-infiltrating lymphocyte; TIME, tumor immune microenvironment; Treg, regulatory T cell.

In many types of cancer, the TMB can serve as an effective predictive biomarker for immunotherapy. Therefore, we scrutinized the somatic mutation scenario of the TCGA dataset among the two risk groups. Our findings showed that 97.6% of the high-risk patients had somatic mutations, whereas only 91.82% of the low-risk patients had somatic mutations (Figure 9A,9B). Although the high-risk subgroup had a slightly higher TMB than the low-risk subgroup, the variance was insignificant (Figure 9C). According to the survival analysis, patients with high TMB experienced notably reduced OS compared with those with low TMB (Figure 9D). TMB and risk score showed a positive link in correlation analysis (Figure 9E). When TMB was combined with the risk score for Kaplan-Meier analysis, patients with characteristics of both high-risk score and high TMB had the lowest OS (Figure 9F).

Figure 9 TMB analysis in OSCC between high- and low-risk groups. (A) The somatic mutation landscape of the high-risk group in the TCGA dataset. (B) The somatic mutation landscape of the low-risk group in the TCGA dataset. (C) Comparison of TMB between high- and low-risk groups. (D) Survival analysis of patients with high TMB vs. low TMB. (E) Correlation analysis between the risk score and TMB. (F) Kaplan-Meier analysis of the OS of OSCC patients stratified by risk score and TMB. (G) ROC curve analysis of TMB for predicting the response to immunotherapy in patients with OSCC using the TIDE method. (H-J) Disparities in immunotherapy sensitivity between patients in high- and low-risk groups. (K) ROC curve analysis of risk score for predicting the response to immunotherapy in patients with OSCC using the TIDE method. AUC, area under the curve; CAF, cancer-associated fibroblast; CI, confidence interval; H, high-risk; L, low-risk; OS, overall survival; OSCC, oral squamous cell carcinoma; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas; TIDE, Tumor Immune Dysfunction and Exclusion; TMB, tumor mutation burden.

The efficacy of immune checkpoint blockade can be predicted using TMB. Therefore, we used the TIDE method to evaluate whether TMB could predict the response to immunotherapy in patients with OSCC. The results showed that the AUC value was 0.543 [95% confidence interval (CI): 0.465–0.623] (Figure 9G), indicating that TMB may have a certain predictive ability for the immune therapy response in patients. In addition, we evaluated the difference in immunotherapy sensitivity between the two risk groups and assessed the capability of the risk score to predict immunotherapy response in patients with OSCC. The results showed that the risk score of the immunotherapy-non-responsive subgroup (n=230) was notably greater than that of the responsive subgroup (n=105) and that patients categorized as low-risk exhibited a greater response to immunotherapy than those categorized as high-risk (Figure 9H-9J). Moreover, as depicted in Figure 9K, the AUROC curve was 0.611, indicating that the predictive ability of the risk score for the immunotherapy response in patients was superior to that of TMB, and we also validated the immunotherapy response in the TCGA-OSCC dataset (Figure S5). In all, the predictive ability of the risk score OSCC prognostic model to immunotherapy response in patients was superior to that in TMB.

Discussion

As the main component of tumor stroma in cancerous tissue, CAFs have been demonstrated to enhance the growth, invasion, and metastasis of OSCC (38-40) and modulate the TIME to attenuate host immune response (41,42). Furthermore, studies have shown a negative association between the OS of patients with OSCC and the proportion of CAFs (43). Consistent with this, our evaluation of CAF infiltration scores in 338 OSCC patients revealed a significant association between more higher CAF infiltration scores and more poorer OS. Therefore, targeting CAFs may be an important strategy for treating OSCC and improving the TIME.

Increasing evidence suggests that CAF can impede the anti-tumor immune response by elevating the proportion of immunosuppressive cells and the expression of immune checkpoint molecules within the TME (44-47). This CAF-driven immunosuppressive microenvironment has deleterious effects on the effectiveness of tumor immunotherapy (48). Consistent with this, we found that patients with the CAF-high subtypes exhibited reduced responsiveness to immunotherapy. Moreover, the results of GSVA revealed significant enrichment of inflammatory response, IL6-JAK-STAT3 signaling, EMT, angiogenesis, and TGF-β signaling pathway in the CAF-high subtype. Previous research has shown that CAF abnormally activates these signaling pathways to promote OSCC progression and immune escape (49-52). These results imply that targeting CAFs may help restore impaired immune responses in patients with OSCC and that concomitant co-treatment with checkpoint inhibitors may yield more desirable therapeutic outcomes (53).

This study also successfully created a risk model of eight crucial genes based on prognostic DEGs between the two subtypes. The model was used to forecast the prognosis of OSCC and was found to have good predictive ability. The risk of this model was compared with that of other models based on CAF-related gene signatures and was found to outperform the other models. Furthermore, compared with TMB, the risk score exhibited a higher predictive ability for the immunotherapy response in patients with OSCC.

Notably, patients in the high-risk subgroup demonstrated diminished anti-tumor immune function levels, particularly in the priming and activation of immune cells and the infiltration of immune cells into tumors. Moreover, the ratio of activated B-cells to CD8⁺ T-cells was significantly lower in the high-risk subgroup. A previous study found that CD8⁺ T cells are the most crucial component of anti-tumor immunity and are significantly associated with a good prognosis in many cancers (54). Although the abundance of CD8⁺ T cells is not an independent predictor of OSCC survival, it is closely associated with patient prognosis and immunotherapy efficacy (55,56). We found that the decrease in T cell co-stimulation and type II IFN response activity may be important reasons for weakening the activation of immune cells and killing tumor cells in high-risk patients (57,58). These findings suggest that enhancing T cell co-stimulation and the type II IFN response to improve the priming and activation of immune cells, as well as promoting immune cell infiltration into tumor tissues, may be a potential method to enhance the anti-tumor immune function of high-risk patients.

In previous studies, some prognostic models that can predict OSCC have been verified (59,60). In our study, we further expanded the sample size and enrolled 338 OSCC samples from TCGA and GEO databases, eliminating inter-individual differences and greatly improving the accuracy of prognosis prediction. Using different Cox analysis and bioinformatics analysis methods, we established a scoring CAF-related prognostic signature that could well reflect risk and positive factors for clinical prognosis in patients with OSCC. Based on the characteristics of TIME, it might be applied to guide individualized immune-therapy for OSCC patients. Therefore, it has a good survival prediction ability and provides an important reference for clinical treatment.

However, there are several limitations in this work. The clinical information of OSCC patients was incomplete, which will affect the predictability of the risk score as an independent prognostic factor, and the clinical characteristics of patients need to be considered comprehensively in the future. Our prognostic model is limited by the public databases. More in vitro and in vivo experiments and clinical research were needed to validate the significance of CAF-related prognostic signature in clinical practice. Additionally, in our preliminary clinical validation, the expression differences for four of the eight signature genes (SERPINB7, SLFN11, TSPAN11, and AQP1) did not reach statistical significance, which may be attributed to the limited sample size and high heterogeneity of OSCC. Future studies with larger, well-annotated cohorts are warranted to further validate their expression patterns and clinical relevance. Although the direct correlation between the continuous risk score and CAF infiltration level was modest, the consistent associations observed through group comparisons and survival analysis support the biological relevance of CAF enrichment in the high-risk phenotype. Future studies are warranted to elucidate the precise molecular interactions between these signature genes and CAF activation. At last, further experiments are required to explore the potential mechanisms of CAF genes in OSCC patients.

Conclusions

We revealed that CAF genes were associated with OSCC survival and established a novel risk score prognostic model based on a signature of eight differentially expressed CAF-related genes, which may also act as potential immune-related therapeutic targets for OSCC patients.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-2129/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-2129/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-2129/prf

Funding: This study was supported by the Basic Research Program of Shenzhen Innovation Council (No. JCYJ20220818102816036), the Shenzhen Fund for Guangdong Provincial High-level Clinical Key Specialties (No. SZGSP008), the Sanming Project of Medicine in Shenzhen (Oral and Maxillofacial Surgery Team, Professor Yu Guangyan, Peking University Hospital of Stomatology) (No. SZSM 202111012), and the Shenzhen Clinical Research Center for Oral Diseases (No. 20210617170745001).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-2129/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethical Committee of Peking University Shenzhen Hospital (No. 2022-117) and all patients provided written informed consent.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
2. Tan Y, Wang Z, Xu M, et al. Oral squamous cell carcinomas: state of the field and emerging directions. Int J Oral Sci 2023;15:44. [Crossref] [PubMed]
3. Yanamoto S, Yamada S, Takahashi H, et al. Clinicopathological risk factors for local recurrence in oral squamous cell carcinoma. Int J Oral Maxillofac Surg 2012;41:1195-200. [Crossref] [PubMed]
4. Ferreira AK, Carvalho SH, Granville-Garcia AF, et al. Survival and prognostic factors in patients with oral squamous cell carcinoma. Med Oral Patol Oral Cir Bucal 2021;26:e387-92. [Crossref] [PubMed]
5. Li H, Zhang Y, Xu M, et al. Current trends of targeted therapy for oral squamous cell carcinoma. J Cancer Res Clin Oncol 2022;148:2169-86. [Crossref] [PubMed]
6. von Witzleben A, Wang C, Laban S, et al. HNSCC: Tumour Antigens and Their Targeting by Immunotherapy. Cells 2020;9:2103. [Crossref] [PubMed]
7. Wu F, Yang J, Liu J, et al. Signaling pathways in cancer-associated fibroblasts and targeted therapy for cancer. Signal Transduct Target Ther 2021;6:218. [Crossref] [PubMed]
8. Takahashi H, Sakakura K, Kawabata-Iwakawa R, et al. Immunosuppressive activity of cancer-associated fibroblasts in head and neck squamous cell carcinoma. Cancer Immunol Immunother 2015;64:1407-17. [Crossref] [PubMed]
9. Li YY, Tao YW, Gao S, et al. Cancer-associated fibroblasts contribute to oral cancer cells proliferation and metastasis via exosome-mediated paracrine miR-34a-5p. EBioMedicine 2018;36:209-20. [Crossref] [PubMed]
10. Dourado MR, Korvala J, Åström P, et al. Extracellular vesicles derived from cancer-associated fibroblasts induce the migration and invasion of oral squamous cell carcinoma. J Extracell Vesicles 2019;8:1578525. [Crossref] [PubMed]
11. Qin X, Guo H, Wang X, et al. Exosomal miR-196a derived from cancer-associated fibroblasts confers cisplatin resistance in head and neck cancer through targeting CDKN1B and ING5. Genome Biol 2019;20:12. [Crossref] [PubMed]
12. Custódio M, Biddle A, Tavassoli M. Portrait of a CAF: The story of cancer-associated fibroblasts in head and neck cancer. Oral Oncol 2020;110:104972. [Crossref] [PubMed]
13. Costea DE, Hills A, Osman AH, et al. Identification of two distinct carcinoma-associated fibroblast subtypes with differential tumor-promoting abilities in oral squamous cell carcinoma. Cancer Res 2013;73:3888-901. [Crossref] [PubMed]
14. Puram SV, Tirosh I, Parikh AS, et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell 2017;171:1611-1624.e24. [Crossref] [PubMed]
15. Katzenelenbogen Y, Sheban F, Yalin A, et al. Coupled scRNA-Seq and Intracellular Protein Activity Reveal an Immunosuppressive Role of TREM2 in Cancer. Cell 2020;182:872-885.e19. [Crossref] [PubMed]
16. Lähnemann D, Köster J, Szczurek E, et al. Eleven grand challenges in single-cell data science. Genome Biol 2020;21:31. [Crossref] [PubMed]
17. Hughes TK, Wadsworth MH 2nd, Gierahn TM, et al. Second-Strand Synthesis-Based Massively Parallel scRNA-Seq Reveals Cellular States and Molecular Features of Human Inflammatory Skin Pathologies. Immunity 2020;53:878-894.e7. [Crossref] [PubMed]
18. Sebastian A, Hum NR, Martin KA, et al. Single-Cell Transcriptomic Analysis of Tumor-Derived Fibroblasts and Normal Tissue-Resident Fibroblasts Reveals Fibroblast Heterogeneity in Breast Cancer. Cancers (Basel) 2020;12:1307. [Crossref] [PubMed]
19. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
20. Mangiola S, Doyle MA, Papenfuss AT. Interfacing Seurat with the R tidy universe. Bioinformatics 2021;37:4100-7. [Crossref] [PubMed]
21. Zheng H, Liu H, Li H, et al. Weighted Gene Co-expression Network Analysis Identifies a Cancer-Associated Fibroblast Signature for Predicting Prognosis and Therapeutic Responses in Gastric Cancer. Front Mol Biosci 2021;8:744677. [Crossref] [PubMed]
22. Chen LJ, Li JY, Nguyen P, et al. Single-cell RNA sequencing unveils unique transcriptomic signatures of endothelial cells and role of ENO1 in response to disturbed flow. Proc Natl Acad Sci U S A 2024;121:e2318904121. [Crossref] [PubMed]
23. Jin S, Plikus MV, Nie Q. CellChat for systematic analysis of cell-cell communication from single-cell transcriptomics. Nat Protoc 2025;20:180-219. [Crossref] [PubMed]
24. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
25. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 2010;26:1572-3. [Crossref] [PubMed]
26. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
27. Qin Y, Liu Y, Xiang X, et al. Cuproptosis correlates with immunosuppressive tumor microenvironment based on pan-cancer multiomics and single-cell sequencing analysis. Mol Cancer 2023;22:59. [Crossref] [PubMed]
28. Fu J, Li K, Zhang W, et al. Large-scale public data reuse to model immunotherapy response and resistance. Genome Med 2020;12:21. [Crossref] [PubMed]
29. Konishi D, Umeda Y, Yoshida K, et al. Regulatory T cells induce a suppressive immune milieu and promote lymph node metastasis in intrahepatic cholangiocarcinoma. Br J Cancer 2022;127:757-65. [Crossref] [PubMed]
30. Yue T, Chen S, Zhu J, et al. The aging-related risk signature in colorectal cancer. Aging (Albany NY) 2021;13:7330-49. [Crossref] [PubMed]
31. Xu L, Deng C, Pang B, et al. TIP: A Web Server for Resolving Tumor Immunophenotype Profiling. Cancer Res 2018;78:6575-80. [Crossref] [PubMed]
32. Mayakonda A, Lin DC, Assenov Y, et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018;28:1747-56. [Crossref] [PubMed]
33. Zheng H, Liu H, Ge Y, et al. Integrated single-cell and bulk RNA sequencing analysis identifies a cancer associated fibroblast-related signature for predicting prognosis and therapeutic responses in colorectal cancer. Cancer Cell Int 2021;21:552. [Crossref] [PubMed]
34. Ko YC, Lai TY, Hsu SC, et al. Index of Cancer-Associated Fibroblasts Is Superior to the Epithelial-Mesenchymal Transition Score in Prognosis Prediction. Cancers (Basel) 2020;12:1718. [Crossref] [PubMed]
35. Yang Y, Ma B, Han L, et al. Integrated single-cell and bulk RNA sequencing analyses reveal a prognostic signature of cancer-associated fibroblasts in head and neck squamous cell carcinoma. Front Genet 2022;13:1028469. [Crossref] [PubMed]
36. Dong L, Sun Q, Song F, et al. Identification and verification of eight cancer-associated fibroblasts related genes as a prognostic signature for head and neck squamous cell carcinoma. Heliyon 2023;9:e14003. [Crossref] [PubMed]
37. Wu S, Huang C, Su L, et al. Cancer associated fibroblast derived gene signature determines cancer subtypes and prognostic model construction in head and neck squamous cell carcinomas. Cancer Med 2023;12:6388-400. [Crossref] [PubMed]
38. Elmusrati AA, Pilborough AE, Khurram SA, et al. Cancer-associated fibroblasts promote bone invasion in oral squamous cell carcinoma. Br J Cancer 2017;117:867-75. [Crossref] [PubMed]
39. Kato K, Miyazawa H, Kawashiri S, et al. Tumour: Fibroblast Interactions Promote Invadopodia-Mediated Migration and Invasion in Oral Squamous Cell Carcinoma. J Oncol 2022;2022:5277440. [Crossref] [PubMed]
40. Wang D, Tian M, Fu Y, et al. Halofuginone inhibits tumor migration and invasion by affecting cancer-associated fibroblasts in oral squamous cell carcinoma. Front Pharmacol 2022;13:1056337. [Crossref] [PubMed]
41. Barrett RL, Puré E. Cancer-associated fibroblasts and their influence on tumor immunity and immunotherapy. Elife 2020;9:e57243. [Crossref] [PubMed]
42. Mao X, Xu J, Wang W, et al. Crosstalk between cancer-associated fibroblasts and immune cells in the tumor microenvironment: new findings and future perspectives. Mol Cancer 2021;20:131. [Crossref] [PubMed]
43. Dourado MR, Guerra ENS, Salo T, et al. Prognostic value of the immunohistochemical detection of cancer-associated fibroblasts in oral cancer: A systematic review and meta-analysis. J Oral Pathol Med 2018;47:443-53. [Crossref] [PubMed]
44. Kinoshita T, Ishii G, Hiraoka N, et al. Forkhead box P3 regulatory T cells coexisting with cancer associated fibroblasts are correlated with a poor outcome in lung adenocarcinoma. Cancer Sci 2013;104:409-15. [Crossref] [PubMed]
45. Harper J, Sainson RC. Regulation of the anti-tumour immune response by cancer-associated fibroblasts. Semin Cancer Biol 2014;25:69-77. [Crossref] [PubMed]
46. Sun Q, Zhang B, Hu Q, et al. The impact of cancer-associated fibroblasts on major hallmarks of pancreatic cancer. Theranostics 2018;8:5072-87. [Crossref] [PubMed]
47. Xiang H, Ramil CP, Hai J, et al. Cancer-Associated Fibroblasts Promote Immunosuppression by Inducing ROS-Generating Monocytic MDSCs in Lung Squamous Cell Carcinoma. Cancer Immunol Res 2020;8:436-50. [Crossref] [PubMed]
48. Monteran L, Erez N. The Dark Side of Fibroblasts: Cancer-Associated Fibroblasts as Mediators of Immunosuppression in the Tumor Microenvironment. Front Immunol 2019;10:1835. [Crossref] [PubMed]
49. Mirkeshavarz M, Ganjibakhsh M, Aminishakib P, et al. Interleukin-6 secreted by oral cancer- associated fibroblast accelerated VEGF expression in tumor and stroma cells. Cell Mol Biol (Noisy-le-grand) 2017;63:131-6. [Crossref] [PubMed]
50. Li YY, Zhou CX, Gao Y. Interaction between oral squamous cell carcinoma cells and fibroblasts through TGF-β1 mediated by podoplanin. Exp Cell Res 2018;369:43-53. [Crossref] [PubMed]
51. Wang Y, Jing Y, Ding L, et al. Epiregulin reprograms cancer-associated fibroblasts and facilitates oral squamous cell carcinoma invasion via JAK2-STAT3 pathway. J Exp Clin Cancer Res 2019;38:274. [Crossref] [PubMed]
52. Kennel KB, Bozlar M, De Valk AF, et al. Cancer-Associated Fibroblasts in Inflammation and Antitumor Immunity. Clin Cancer Res 2023;29:1009-16. [Crossref] [PubMed]
53. Ford K, Hanley CJ, Mellone M, et al. NOX4 Inhibition Potentiates Immunotherapy by Overcoming Cancer-Associated Fibroblast-Mediated CD8 T-cell Exclusion from Tumors. Cancer Res 2020;80:1846-60. [Crossref] [PubMed]
54. Borsetto D, Tomasoni M, Payne K, et al. Prognostic Significance of CD4+ and CD8+ Tumor-Infiltrating Lymphocytes in Head and Neck Squamous Cell Carcinoma: A Meta-Analysis. Cancers (Basel) 2021;13:781. [Crossref] [PubMed]
55. Gao A, Pan X, Yang X, et al. Predictive factors in the treatment of oral squamous cell carcinoma using PD-1/PD-L1 inhibitors. Invest New Drugs 2021;39:1132-8. [Crossref] [PubMed]
56. Lequerica-Fernández P, Suárez-Canto J, Rodriguez-Santamarta T, et al. Prognostic Relevance of CD4(+), CD8(+) and FOXP3(+) TILs in Oral Squamous Cell Carcinoma and Correlations with PD-L1 and Cancer Stem Cell Markers. Biomedicines 2021;9:653. [Crossref] [PubMed]
57. Castro F, Cardoso AP, Gonçalves RM, et al. Interferon-Gamma at the Crossroads of Tumor Immune Surveillance or Evasion. Front Immunol 2018;9:847. [Crossref] [PubMed]
58. Sundqvist KG. T Cell Co-Stimulation: Inhibition of Immunosuppression? Front Immunol 2018;9:974. [Crossref] [PubMed]
59. Liu Y, Wang T, Li R. A prognostic Risk Score model for oral squamous cell carcinoma constructed by 6 glycolysis-immune-related genes. BMC Oral Health 2022;22:324. [Crossref] [PubMed]
60. Yue K, Yao X. Prognostic model based on telomere-related genes predicts the risk of oral squamous cell carcinoma. BMC Oral Health 2023;23:484. [Crossref] [PubMed]