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

Abstract

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.

Keywords: immunotherapy; oral squamous cell carcinoma; prognosis; signature; survival.

Figure 1

The relationship between stromal/immune/ESTIMATE scores and survival rate of OSCC patients. (a) KM survival analysis of OSCC patients based on their stromal scores. (b) KM survival analysis of OSCC patients based on their immune scores. (c) KM survival analysis of OSCC patients based on their ESTIMATE scores. Red line represents the survival curve of OSCC patients with higher values of stromal scores, immune scores, and ESTIMATE scores (n = 159). Blue line represents the survival curve of OSCC patients with lower values of stromal scores, immune scores, and ESTIMATE scores (n = 159).

Figure 2

The relationship between stromal/immune/ESTIMATE scores and clinical features of OSCC patients. (a) Relationship between stromal scores and TNM stages of OSCC patients. (b) Relationship between immune scores and TNM stages of OSCC patients. (c) Relationship between ESTIMATE scores and TNM stages of OSCC patients. Sample numbers of each group are as follows: T1, n = 19; T2, n = 54; T3, n = 59; T4, n = 158. (d) Relationship between stromal scores and genders of OSCC patients. (e) Relationship between immune scores and genders of OSCC patients. (f) Relationship between ESTIMATE scores and genders of OSCC patients. Sample numbers of each group are as follows: female, n = 101; male, n = 218. (g) Relationship between stromal scores and tumor grades of OSCC patients. (h) Relationship between immune scores and tumor grades of OSCC patients. (i) Relationship between ESTIMATE scores and tumor grades of OSCC patients. Sample numbers of each group are as follows: G1, n = 51; G2, n = 195; G3, n = 63; G4, n = 2.

Figure 3

Venn diagram and heatmap analysis of the differentially expressed genes (DEGs) based on the immune and stromal scores. (a) Heatmap analysis of the DEGs between the higher stromal scores and lower stromal scores in OSCC patients. (b) Heatmap analysis of the DEGs between the higher immune scores and lower immune scores in OSCC patients. (c) Venn diagrams analysis the number of upregulated genes of higher stromal scores in OSCC patients and upregulated genes of higher immune scores in OSCC patients. (d) Venn diagrams analysis of the number of downregulated genes of higher stromal scores in OSCC patients and downregulated genes of higher immune scores in OSCC patients.

Figure 4

Gene enrichment analysis of the DEGs based on the immune and stromal scores in OSCC patients. (a) Gene ontology (GO) enrichment analysis of DEGs based on the immune and stromal scores in OSCC patients. Top panel is the biological process analysis of GO enrichment. Middle panel is the cellular component analysis of GO enrichment. Bottom panel is the molecular function analysis of GO enrichment. (b) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs based on the immune and stromal scores in OSCC patients.

Figure 5

The protein-protein interaction (PPI) analysis of DEGs based on the immune and stromal scores in OSCC patients. The interaction nodes, scores > 0.95 were set as the cutoff to the PPIs network.

Figure 6

The top 30 key genes analyzed by PPI which associated with the prognosis of OSCC patients, including ATP8B4, CD53, LILRB2, ITGAX, CD3G, IGHV3-11, SUCNR1, P2RY13, GPR183, GALR2, FPR3, CCR8, CCL13, C5AR1, ADORA3, P2RY12, CXCR3, CX3CR1, CCR4, ADRA2A, CCR1, FPR1, CCR2, CCR5, ITGB2, FCER1G, C3AR1, C3 and FPR2.

Figure 7

Univariate Cox and LASSO Cox survival analysis for the DEGs based on PPI analysis in OSCC patients. (a) Univariate Cox survival analysis of the DEGs. (b,c) LASSO Cox analysis identified 23 DEGs with the prognostic value and 23 DEGs were identified as prognostic factors of OSCC patients. Protective factors mean the genes’ HR was <1, while risk factors mean the genes’ HR was >1 in OSCC patients.

Figure 8

Prognostic effect of the immune-related signature. The survival status (a), risk scores (b,c) of OSCC patients, and heatmap of immune-related DEGs expression pattern (d). Kaplan-Meier survival curve (e) for OS of the low-risk and high-risk groups of OSCC patients. (f) Prognostic value evaluation of the 11-gene signature using time-specific ROC curves and dynamic AUC lines analysis. The time-dependent ROC curves are based on 1, 2, and 3 years of follow-up and the dynamic AUC lines of OSCC patients.

Figure 9

Validation of the 11-gene signature in the GEO cohort (GSE65858). The overall survival status (a), risk scores (b), the heatmap of the 11-gene expression pattern (c), Kaplan-Meier curve analysis with risk scores (d), and ROC curve (e) are shown.