A Novel Immune-Related Gene Signature to Identify the Tumor Microenvironment and Prognose Disease Among Patients With Oral Squamous Cell Carcinoma Patients Using ssGSEA: A Bioinformatics and Biological Validation Study

Abstract

Oral squamous cell carcinoma (OSCC) is the most invasive oral malignancy in adults and is associated with a poor prognosis. Accurate prognostic models are urgently needed, however, knowledge of the probable mechanisms behind OSCC tumorigenesis and prognosis remain limited. The clinical importance of the interplay between the immune system and tumor microenvironment has become increasingly evident. This study explored immune-related alterations at the multi-omics level to extract accurate prognostic markers linked to the immune response and presents a more accurate landscape of the immune genomic map during OSCC. The Cancer Genome Atlas (TCGA) OSCC cohort (n = 329) was used to detect the immune infiltration pattern of OSCC and categorize patients into two immunity groups using single-sample gene set enrichment analysis (ssGSEA) and hierarchical clustering analysis. Multiple strategies, including lasso regression (LASSO), Cox proportional hazards regression, and principal component analysis (PCA) were used to screen clinically significant signatures and identify an incorporated prognosis model with robust discriminative power on the survival status of both the training and testing set. We identified two OSCC subtypes based on immunological characteristics: Immunity-high and immunity low, and verified that the categorization was accurate and repeatable. Immunity_ high cluster with a higher immunological and stromal score. 1047 differential genes (DEGs) integrate with immune genes to obtain 319 immue-related DEGs. A robust model with five signatures for OSCC patient prognosis was established. The GEO cohort (n = 97) were used to validate the risk model's predictive value. The low-risk group had a better overall survival (OS) than the high-risk group. Significant prognostic potential for OSCC patients was found using ROC analysis and immune checkpoint gene expression was lower in the low-risk group. We also investigated at the therapeutic sensitivity of a number of frequently used chemotherapeutic drugs in patients with various risk factors. The underlying biological behavior of the OSCC cell line was preliminarily validated. This study characterizes a reliable marker of OSCC disease progression and provides a new potential target for immunotherapy against this disease.

Keywords: immune infiltration; immune-related gene; oral squamous cell carcinoma; prognostic biomarker; single-sample gene set enrichment analysis.

Figure 1

Construction and verification of oral squamous cell carcinoma clustering. (A) All 22 invading immune cells are represented by a correlation matrix. Immune cells were shown to be favorably associated and are represented in red, while others were found to be negatively related and are represented in blue. The threshold was set at P < 0.05. (B) Using ssGSEA analysis, gene expression data from OSCC patients were divided into two clusters. (C) The PCA plot of the distribution status of the two OSCC clusters. (D) The heatmap showed that the 29 immune-related cell types had high expression in the high-immune cell infiltration group (Immunity-high), and low expression in the low immune cell infiltration group (Immunity-low). The tumor purity and ESTIMATE, Immune, and Stromal Scores of each patient are shown with clustering information using the ESTIMATE algorithm. (E) The violin plot shows the difference in the ESTIMATE, Immune, and Stromal Scores between the two clusters. (F) The box plot shows a statistically significant difference in HLA family expression. (G) The box plot shows a statistical difference in immune cell infiltration between the two clusters. *P < 0.05, **P < 0.01, ***P < 0.001.

Figure 2

Gene functional enrichment analysis of the immunity-high and immunity-low clusters.

Figure 3

Analysis of differentially expressed immune-related genes in oral squamous cell carcinoma patients. (A) The volcano plot shows 761up-regulated genes and 286 down-regulated genes in the immunity-high and immunity-low clusters. Up-regulated and down-regulated genes are represented by red and green bars, respectively. The gene names of genes with |log2FC| > 6 are displayed. (B) The heatmap shows the degree of DEG expression in the immune-high and immune-low clusters. (C) The heatmap shows immune-related gene expression in the immunity-high and immunity-low clusters. (D) The Venn diagram shows 319 genes from both gene sets. DEGs, differentially expressed genes.

Figure 4

Development of a prognostic signature for oral squamous cell carcinoma based on immune-related genes. (A) Using univariable Cox regression analysis, the HR and p-value for the chosen genes in the immune terms. (B) The interaction between the 14 immune-related genes and transcription factors is depicted as an alluvial diagram. (C) The interaction network of the immune-related prognostic genes. (D) The 14 immune-related gene LASSO coefficient profiles. (E) 1,000-round cross-validation was used to find the best values for the penalty parameter. (F) The heatmap shows the correlation between immune-related prognostic genes and immune cell infiltration. (G) PCA results for prognostic genes in two clusters of immunity level in training set. (H) PCA results for prognostic genes in high- and low-risk groups in training set.

Figure 5

Construction and validation of the immune-related gene prognostic signature in the training and testing sets. The survival status of patients in the high-risk and low-risk groups in the training (A) and testing sets (B). Kaplan-Meier survival curves for OSCC patients in the training (C) and testing sets (D). The prognostic signature’s time-independent ROC curve at 1-, 3-, and 5-years in the training (E) and testing sets (F). Each OSCC sample’s risk curve is reordered by the riskscore in the training (G) and testing sets (H). A scatter plot depicts the survival of OSCC samples in the training (I) and testing sets (J). Interaction analysis of the immune-related prognostic genes in the training (K) and testing sets (L).

Figure 6

Clinical response to anti-tumor therapy as well as immune checkpoint-related gene expression in the high-risk and low-risk groups. (A) The chemotherapy and molecular drugs prediction of risk groups. (B) The difference of immune checkpoints expression in risk groups of traing set. (C) The difference of immune checkpoints expression in risk groups of testing set. (D) Correlations between riskscore and immune infiltration cells. (E) Correlations between riskscore and the enrichment scores of immunotherapy-predicted pathways. *P < 0.05, **P < 0.01, ***P < 0.001.

Figure 7

Construction of a nomogram and verification that the immune-related gene prognostic signature is an independent prognostic factor. Univariate (A) and multivariate (B) Cox regression analysis of the immune-related gene prognostic signature in OSCC patients to determine independent risk variables. (C) The development of a nomogram based on the immune-related gene prognostic signature in the TCGA training cohort. (D) The calibration curve of the nomogram. (E) The combined ROC for riskscore, nomogram, gender, stage, age, and TMN at 1-, 3-, and 5-years. (F) Time-independent ROC curves of overall survival for immune-related gene prognostic model, Lv geneSig, Ribeiro geneSig, Zhao geneSig, and Zhang geneSig at 1-, 3-, and 5-years.

Figure 8

The correlation between the immune-related gene prognostic signature and TMB. (A) The box plot for TMB levels among patients in the high- and low-risk groups. (B) Kaplan–Meier curves for the high- and low-TMB of OSCC patients. (C) Kaplan–Meier curves for OSCC patients by TMB status in the high-risk and low-risk groups. (D) The oncoPrint was constructed based on CNV profile in the high-risk scores of OSCC patients. (E) The oncoPrint was constructed based on CNV profile in the low-risk scores of OSCC patients. Individual patients are represented in each column. TMB, tumor mutational burden; CNV, copy number variation.

Figure 9

CTSG and TNFRSF4 overexpression using pEXP-RB-Mam-EGFP transfection inhibits OSCC cell line viability and clonogenicity. CTSG (A) and TNFRSF4 (B) protein expression in OSCC cell lines and normal human oral cavity epithelial cells. CTSG (C) and TNFRSF4 (D) mRNA expression in OSCC cell lines and normal human oral cavity epithelial cells. CTSG (E) and TNFRSF4 (F) protein and mRNA expression in OSCC cell lines transfected with pEXP-RB-Mam-EGFP. Colony formation assay of OSCC cell lines treated with specific pEXP-RB-Mam-EGFP and the negative control of CTSG and TNFRSF4 (G). Transwell migration (H) and invasion (I) assay of OSCC cell lines treated with specific pEXP-RB-Mam-EGFP and the negative control of CTSG and TNFRSF4. Wound healing assays of OSCC cell lines treated with specific pEXP-RB-Mam-EGFP and the negative control of CTSG and TNFRSF4 (J). *P < 0.05