# Identification of Candidate Target Genes and Immune Cells in Oral Squamous Cell Carcinoma

Authors: Pengfeng Xie; Shichao Wu; Lijuan Guo; Jun Ren; Kaizhi Cai; Mingyue Zhou; Weiwei Liu; Sen Yang
Journal: Computational and Mathematical Methods in Medicine
 (2021)

---

## Abstract

Background. The advance of new treatment strategies for more effective management of oral cancer requires identification of novel biological targets. Therefore, the purpose of this study is to identify novel biomarkers associated with oral tumorigenesis and prognostic signature by comparing gene expression profile of oral squamous cell carcinomas (OSCCs). Methods. Four datasets including GSE25099, GSE30784, GSE37991, and GSE41613 were collected from Gene Expression Omnibus (GEO) database. Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, Cox model analysis, identification of key genes, and Kaplan-Meier analysis were also performed. The xCell was utilized to analyze the infiltration levels of immune cells. Results. A total of 235 differentially expressed genes (DEGs) were found to be dysregulated in OSCC. These genes were mainly enriched in ECM receptor interaction and focal adhesion. Cox regression analysis identified 10 genes considered as key genes. Kaplan-Meier analysis showed that low expression of SERPINE1 (also known as PAI-1), high expression of CD1C, and C-X3-C motif chemokine receptor 1 (CX3CR1) were associated with well prognostic status in OSCC patients. In addition, we constructed a 3-immune-cell signature (myeloid dendritic cell, T cell CD4+ central memory, and common myeloid progenitor) that may be used to predict the survival status of OSCC patients. Conclusion. Three key genes and 3-immune-cell signature were potential biomarkers for the prognosis of OSCC, and they may serve as potential targets for the treatment of OSCC patients.


---

## Body

## 1. Introduction

Head and neck cancers are mostly derived from oral, pharyngeal, and laryngeal mucosal epithelia, collectively referred to as head and neck squamous cell carcinoma [1]. Oral squamous cell carcinoma (OSCC) is one of the more common 10 malignant tumors in humans, accounting for 40% of head and neck squamous cell cancer (HNSCC) [2]. Approximately 350 000 people are diagnosed with oral squamous cell carcinoma each year [3]. Current treatment of OSCC mainly relied on single or multimodality therapy including surgery, radiation, and/or chemotherapy. Despite rapid advances in these treatment modalities as well as improvements in the early detection of oral cancer, the 5-year survival rate of OSCC remains only 50% in the past decades and 25%–50% of them will suffer from local relapse and distant metastasis after treatment [4]. The incidence of OSCC is high and the prognosis is poor, and an increasing number of investigators have conducted in-depth studies on the pathogenesis and potential therapeutic targets for OSCC [5].Despite there have been several articles reporting the discovery of OSCC biomarkers, only few of them have been validated and successfully applied in routine clinical practice [1–3, 6, 7]. Nonetheless, there were some limitations with regard to most biomarkers in the early detection of OSCC, and sometimes, the prognostic value of them was inaccurate that plagued most people [7]. By using bioinformatics analysis, it is more likely to overcome the defects in the reported biomarkers. Bioinformatics analysis enabled researchers to in-depth study the comprehensive data of a large number of clinical samples from different independent studies, which provided a data infrastructure for unraveling promising biomarkers and updating our understanding of cancer [3, 5, 6].The tumor microenvironment (TME) is a complex and heterogeneous mix of tumor cells and stromal cells, including endothelial cells, cancer-associated fibroblasts (CAFs), and immune cells [8]. The TME assists in tumor initiation and progression [9]. CAFs further assist oral squamous cell carcinoma metastasis to lymph nodes because of the infiltration of immune heterogeneous immune cells [10]. OSCC has lymph node metastasis in more than 40% of patients at the time of diagnosis, and worse, lymph node metastasis has decreased the survival rate of OSCC patients by nearly 50% [11].All malignancies have the special ability to allow continued proliferation, invasion, and metastasis. Pathophysiological processes are driven by some complex molecular signaling pathways [12]. There have been numerous studies showing that aberrant expression of genes in OSCC affects prognosis and is implicated in the alteration of multiple biological functions [13]. However, specific molecular OSCC prophecies have only been partially identified [14]. It is therefore necessary to search for effective diagnostic and prognostic biomarkers for OSCC.Here, we utilized sequencing data of OSCC from public databases and analyzed potential genes for prognosis and identified abnormally regulated biological functions. Additionally, we employed the evaluation of immune cell infiltration to judge the prognosis of OSCC.


## 2. Materials and Methods

### 2.1. Data Collection

GSE25099 included gene expression profile of oral tissue from 57 specimens from patients with OSCC and 22 oral tissues from normal persons. GSE30784 included gene expression profile of oral tissues from 167 OSCC and 45 normal controls. GSE37991 included gene expression profile of biopsies from 40 OSCC patients and 40 nontumor pairwise samples. Survival data for OSCC patients were obtained from GSE41613. These data were all collected from the Gene Expression Omnibus (GEO) database.


### 2.2. Differential Analysis

Differential expression analysis was performed using the limma R package. The genes with an∣log2FoldChange∣>1 and P<0.05 between OSCC and control were assigned as differentially expressed genes (DEGs).


### 2.3. Enrichment Analysis

The enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) for common genes was analyzed using the Enrichr R package. TheP value < 0.05 was considered significantly enriched.


### 2.4. Immune Cell Infiltration

The xCell R package [15] was used to assign and visualize 39 types of immune cells in OSCC and controls.


### 2.5. Identification of Risk Prognostic Factors

Genes with significant prognostic value (P<0.05) in DEGs or immune cells were determined with univariate Cox regression analysis. OSCC samples were divided into high-risk and low-risk groups based on the median level of Cox score. The accuracy of the Cox score was assessed through receiver operating characteristic (ROC) curve with the pROC R package. Then, the prognostic value of high-risk or low-risk groups was assessed using the Kaplan-Meier (KM) curves through the survival R package.


## 2.1. Data Collection

GSE25099 included gene expression profile of oral tissue from 57 specimens from patients with OSCC and 22 oral tissues from normal persons. GSE30784 included gene expression profile of oral tissues from 167 OSCC and 45 normal controls. GSE37991 included gene expression profile of biopsies from 40 OSCC patients and 40 nontumor pairwise samples. Survival data for OSCC patients were obtained from GSE41613. These data were all collected from the Gene Expression Omnibus (GEO) database.


## 2.2. Differential Analysis

Differential expression analysis was performed using the limma R package. The genes with an∣log2FoldChange∣>1 and P<0.05 between OSCC and control were assigned as differentially expressed genes (DEGs).


## 2.3. Enrichment Analysis

The enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) for common genes was analyzed using the Enrichr R package. TheP value < 0.05 was considered significantly enriched.


## 2.4. Immune Cell Infiltration

The xCell R package [15] was used to assign and visualize 39 types of immune cells in OSCC and controls.


## 2.5. Identification of Risk Prognostic Factors

Genes with significant prognostic value (P<0.05) in DEGs or immune cells were determined with univariate Cox regression analysis. OSCC samples were divided into high-risk and low-risk groups based on the median level of Cox score. The accuracy of the Cox score was assessed through receiver operating characteristic (ROC) curve with the pROC R package. Then, the prognostic value of high-risk or low-risk groups was assessed using the Kaplan-Meier (KM) curves through the survival R package.


## 3. Results

### 3.1. Genes Associated with OSCC

To identify OSCC-associated genes, the differential analysis of gene expression between OSCC and controls was performed. We obtained 973 DEGs in GSE25099 (Figure1(a)), 2488 DEGs in GSE30784 (Figure 1(b)), and 2204 DEGs in GSE37991 (Figure 1(c)). We identified 235 common genes among the three sets of DEGs by intersection analysis (Figure 1(d)). The common genes were similarly aberrantly expressed in GSE25099 (Figure 1(e)), GSE30784 (Figure 1(f)), and GSE37991 (Figure 1(g)).Figure 1
Identification of deregulated expressed genes for oral squamous cell carcinoma. Differentially expressed genes between OSCC and controls in (a) GSE25099, (b) GSE30784, and (c) GSE37991. (d) Venny map of differentially expressed genes among the three groups of DEGs. Heatmap of common genes in (e) GSE25099, (f) GSE30784, and (g) GSE37991. OSCC: oral squamous cell carcinoma.
(a)(b)(c)(d)(e)(f)(g)


### 3.2. Molecular Mechanisms of OSCC

To identify the molecular mechanisms associated with OSCC, we performed an enrichment analysis of common genes. GO results (Figure2(a)) showed that these genes were mainly enriched in extracellular matrix organization, cytokine-mediated signaling pathway, and response to type I interferon of biological processes (BP). Endoplasmic reticulum lumen and platelet alpha granule were significantly enriched in cellular components (CC). And metalloendopeptidase activity and metallopeptidase activity of molecular functions (MF) were enriched by common genes. In addition, ECM receptor interaction, cytokine receptor interaction, and focal adhesion were significantly enriched by common genes in KEGG enrichment results (Figure 2(b)).Figure 2
Enrichment analysis of common genes. (a) The top significant 10 terms of biological processes, cellular components, and molecular functions of common gene enrichment. (b) KEGG pathways significantly enriched in common genes.
(a)(b)


### 3.3. Candidate Genes for OSCC

To identify candidate target genes for OSCC, we performed Cox regression analysis on common genes. The top 10 most significant genes were selected as key genes (Table1). The OSCC samples were divided into high-risk and low-risk groups by the median value of the risk score (Figure 3(a)). In the high-risk group, there were fewer OSCC patients ultimately surviving, compared with the low-risk group. The expression of key genes differed in the high-risk and low-risk groups. In addition, the risk score had a good ability to predict the prognosis of OSCC patients (AUC>0.8) (Figure 3(b)). The low-risk group had a significantly higher survival probability than the high-risk group (Figure 3(c)). The nomogram results showed that low expression of SERPINE1, high expression of CD1C, and CX3CR1 may benefit the survival probability of OSCC patients (Figure 4).Table 1
Top 10 genes with smallestP values in Cox proportional hazards regression analysis.
SymbolHR (exp(coef))Coef95% CI lower95% CI upperZP valueSERPINE11.4335680.3601660.1707670.5495653.727121.94E−04SCG51.2971250.260150.119180.401123.6169732.98E−04MUC150.770907-0.26019-0.40435-0.11602-3.537314.04E−04CD1C0.528052-0.63856-0.99697-0.28015-3.491994.79E−04NEK61.9764340.6812940.2925061.0700833.4345465.94E−04INHBA1.437770.3630930.1527530.5734333.3833347.16E−04LAMC21.4164730.348170.1462060.5501343.3788217.28E−04CX3CR10.573106-0.55668-0.91193-0.20144-3.071360.002131SPINK70.882734-0.12473-0.20696-0.0425-2.972930.00295TGM30.885544-0.12155-0.20225-0.04086-2.952260.003155Figure 3
Risk score analysis of gene signature in OSCC. (a) Risk score and survival status distribution according to the median of the risk score for each signature. (b) ROC curves of risk scores to predict patients’ 1-, 3-, and 5-year survival. (c) Kaplan-Meier survival curves of high- and low-risk groups for OSCC patients. OSCC: oral squamous cell carcinoma.
(a)(b)(c)Figure 4
Nomogram of risk gene signature.


### 3.4. Immune Infiltration of OSCC

Among the signaling pathways enriched in common genes, we found a large number of immune inflammation-related pathways. Therefore, we calculated the infiltration of immune cells in OSCC patients (Figure5(a)). Cox regression analysis identified three immune cells that may be significant for the prognosis of OSCC patients (P<0.05). Nomogram results showed that high infiltration levels of myeloid dendritic cell, T cell CD4+ central memory, and common myeloid progenitor were beneficial to patient prognosis (Figure 5(b)). OSCC patients were divided into high-risk and low-risk groups by the median values of the risk scores for the three immune cells (Figure 6(a)). In the high-risk group, there were fewer OSCC patients ultimately surviving, compared with the low-risk group. The infiltration of immune cells differed in the high-risk and low-risk groups. The risk score had a potential ability to predict the prognosis of OSCC patients (AUC>0.6) (Figure 6(b)). In addition, the low-risk group had a significantly higher survival probability than the high-risk group (Figure 6(c)).Figure 5
Immune cell infiltration in OSCC. (a) Heatmap of immune cell infiltration in OSCC and control of GSE25099, GSE30784, and GSE37991. (b) Nomogram of risk immune cell signature. OSCC: oral squamous cell carcinoma.
(a)(b)Figure 6
Risk score analysis of immune cell signature in OSCC. (a) Risk score and survival status distribution according to the median of the risk score for immune cell signatures. (b) ROC curves of risk scores to predict patients’ 1-, 3-, and 5-year survival. (c) Kaplan-Meier survival curves of high- and low-risk groups for OSCC patients. OSCC: oral squamous cell carcinoma.
(a)(b)(c)


## 3.1. Genes Associated with OSCC

To identify OSCC-associated genes, the differential analysis of gene expression between OSCC and controls was performed. We obtained 973 DEGs in GSE25099 (Figure1(a)), 2488 DEGs in GSE30784 (Figure 1(b)), and 2204 DEGs in GSE37991 (Figure 1(c)). We identified 235 common genes among the three sets of DEGs by intersection analysis (Figure 1(d)). The common genes were similarly aberrantly expressed in GSE25099 (Figure 1(e)), GSE30784 (Figure 1(f)), and GSE37991 (Figure 1(g)).Figure 1
Identification of deregulated expressed genes for oral squamous cell carcinoma. Differentially expressed genes between OSCC and controls in (a) GSE25099, (b) GSE30784, and (c) GSE37991. (d) Venny map of differentially expressed genes among the three groups of DEGs. Heatmap of common genes in (e) GSE25099, (f) GSE30784, and (g) GSE37991. OSCC: oral squamous cell carcinoma.
(a)(b)(c)(d)(e)(f)(g)


## 3.2. Molecular Mechanisms of OSCC

To identify the molecular mechanisms associated with OSCC, we performed an enrichment analysis of common genes. GO results (Figure2(a)) showed that these genes were mainly enriched in extracellular matrix organization, cytokine-mediated signaling pathway, and response to type I interferon of biological processes (BP). Endoplasmic reticulum lumen and platelet alpha granule were significantly enriched in cellular components (CC). And metalloendopeptidase activity and metallopeptidase activity of molecular functions (MF) were enriched by common genes. In addition, ECM receptor interaction, cytokine receptor interaction, and focal adhesion were significantly enriched by common genes in KEGG enrichment results (Figure 2(b)).Figure 2
Enrichment analysis of common genes. (a) The top significant 10 terms of biological processes, cellular components, and molecular functions of common gene enrichment. (b) KEGG pathways significantly enriched in common genes.
(a)(b)


## 3.3. Candidate Genes for OSCC

To identify candidate target genes for OSCC, we performed Cox regression analysis on common genes. The top 10 most significant genes were selected as key genes (Table1). The OSCC samples were divided into high-risk and low-risk groups by the median value of the risk score (Figure 3(a)). In the high-risk group, there were fewer OSCC patients ultimately surviving, compared with the low-risk group. The expression of key genes differed in the high-risk and low-risk groups. In addition, the risk score had a good ability to predict the prognosis of OSCC patients (AUC>0.8) (Figure 3(b)). The low-risk group had a significantly higher survival probability than the high-risk group (Figure 3(c)). The nomogram results showed that low expression of SERPINE1, high expression of CD1C, and CX3CR1 may benefit the survival probability of OSCC patients (Figure 4).Table 1
Top 10 genes with smallestP values in Cox proportional hazards regression analysis.
SymbolHR (exp(coef))Coef95% CI lower95% CI upperZP valueSERPINE11.4335680.3601660.1707670.5495653.727121.94E−04SCG51.2971250.260150.119180.401123.6169732.98E−04MUC150.770907-0.26019-0.40435-0.11602-3.537314.04E−04CD1C0.528052-0.63856-0.99697-0.28015-3.491994.79E−04NEK61.9764340.6812940.2925061.0700833.4345465.94E−04INHBA1.437770.3630930.1527530.5734333.3833347.16E−04LAMC21.4164730.348170.1462060.5501343.3788217.28E−04CX3CR10.573106-0.55668-0.91193-0.20144-3.071360.002131SPINK70.882734-0.12473-0.20696-0.0425-2.972930.00295TGM30.885544-0.12155-0.20225-0.04086-2.952260.003155Figure 3
Risk score analysis of gene signature in OSCC. (a) Risk score and survival status distribution according to the median of the risk score for each signature. (b) ROC curves of risk scores to predict patients’ 1-, 3-, and 5-year survival. (c) Kaplan-Meier survival curves of high- and low-risk groups for OSCC patients. OSCC: oral squamous cell carcinoma.
(a)(b)(c)Figure 4
Nomogram of risk gene signature.


## 3.4. Immune Infiltration of OSCC

Among the signaling pathways enriched in common genes, we found a large number of immune inflammation-related pathways. Therefore, we calculated the infiltration of immune cells in OSCC patients (Figure5(a)). Cox regression analysis identified three immune cells that may be significant for the prognosis of OSCC patients (P<0.05). Nomogram results showed that high infiltration levels of myeloid dendritic cell, T cell CD4+ central memory, and common myeloid progenitor were beneficial to patient prognosis (Figure 5(b)). OSCC patients were divided into high-risk and low-risk groups by the median values of the risk scores for the three immune cells (Figure 6(a)). In the high-risk group, there were fewer OSCC patients ultimately surviving, compared with the low-risk group. The infiltration of immune cells differed in the high-risk and low-risk groups. The risk score had a potential ability to predict the prognosis of OSCC patients (AUC>0.6) (Figure 6(b)). In addition, the low-risk group had a significantly higher survival probability than the high-risk group (Figure 6(c)).Figure 5
Immune cell infiltration in OSCC. (a) Heatmap of immune cell infiltration in OSCC and control of GSE25099, GSE30784, and GSE37991. (b) Nomogram of risk immune cell signature. OSCC: oral squamous cell carcinoma.
(a)(b)Figure 6
Risk score analysis of immune cell signature in OSCC. (a) Risk score and survival status distribution according to the median of the risk score for immune cell signatures. (b) ROC curves of risk scores to predict patients’ 1-, 3-, and 5-year survival. (c) Kaplan-Meier survival curves of high- and low-risk groups for OSCC patients. OSCC: oral squamous cell carcinoma.
(a)(b)(c)


## 4. Discussion

Our study utilized different public data to analyze the genes differentially expressed in OSCC and investigate the related signaling pathways to identify possible biological roles and target genes in the process of OSCC. Immune cell infiltration levels in OSCC were also analyzed utilizing public data. Cox regression analysis was further utilized to identify candidate genes and immune cells that may influence the prognosis of OSCC. In addition, we also found abnormal immunoinflammatory signaling pathways in OSCC.Utilizing three datasets, we identified common genes that were differentially expressed between OSCC and controls. The enrichment results of these genes revealed that OSCC-related genes were mainly involved in extracellular matrix organization, response to type I interferon, and cell proliferation. We identified important pathways in KEGG results. This is consistent with previous studies. The extracellular matrix has been reported to be potentially involved in the progression and spread of OSCC [16, 17]. Type I interferon, as an autocrine or paracrine pathway underlying cancer immunosurveillance, is involved in the immune response process of OSCC [18, 19]. ECM remodeling can promote tumor growth and invasiveness, with a key role in oral squamous cell carcinogenesis and progression [20, 21]. The formation of ECM receptor interaction and focal adhesion signaling pathways may be associated with integrin gene expression and contribute to OSCC progression [22, 23]. Focal adhesion is a cell substratum adhesion structure mediated by integrins, whose functions include immobilizing the ends of actin filaments, facilitating strong attachment to substrates, and functioning as an integrin signaling platform [24].The results of Cox regression analysis and nomogram showed that the abnormal expression of SERPINE1, CD1C, and CX3CR1 might affect the prognosis of OSCC patients. SERPINE1 is an important gene associated with metastasis and promotes OSCC invasion and metastasis [25]. SERPINE1, also known as plasminogen activator inhibitor-1, is a regulator of the urokinase and tissue-type plasminogen activators. These serine proteases in turn activate the proenzyme plasminogen to plasmin that promotes invasion by degradation of the extracellular matrix, as well as activation of matrix metalloproteinases. SERPINE1 is a potential prognostic tool in OSCC through cell adhesion, migration, invasion, and cell proliferation [26]. CD1C is a relevant gene for dendritic cells and is involved in the immune microenvironment of OSCC [27]. Dendritic cell subtype which highly expresses CD1C can elicit a wide range of responses, including activation of Th1 and Th2 [28]. The gene coding for C-X3-C motif chemokine receptor 1 (CX3CR1) is located on chromosome 3p22.2. CX3CR1 is a 7-transmembrane receptor coupled to heterotrimeric G proteins (GPCRs). The CX3CL1-CX3CR1 axis mediates the adhesion of leukocytes and is also involved in cell survival and recruitment of immune cell subpopulations. Downregulation of CX3CR1 can significantly inhibit lung cancer and pancreatic ductal adenocarcinoma cell proliferation and increase apoptosis [29, 30]. Therefore, we conclude that SERPINE1, CD1C, and CX3CR1 are therapeutic targets for OSCC and have future clinical implications. As previously described, a set of prognostic signatures including PLAU, CLDN8, and CDKN2A was identified from the database to predict overall survival in patients with OSCC. PSMA7, ITGA6, ITGB4, and APP were overexpressed in HNSCC, and they exhibited a significant correlation with poor overall survival in patients. In addition, it was found that LGALS1 was linked to oral cancer progression and metastasis and potentially regarded as a prognostic biomarker for oral cancer therapy [30].Recent studies have shown that tumor-infiltrating immune cells can regulate tumor development and progression [31]. Previous studies have shown that AUNIP, a novel prognostic biomarker, has been shown to be related to stromal and immune scores in OSCC. It was demonstrated that AUNIP was associated with TME, human papillomavirus infection, and cell cycle in OSCC. The inhibited AUNIP led to the suppression of OSCC cell proliferation and brought about G0/G1 phase arrest in OSCC cells. The survival analysis showed that the overexpression of AUNIP predicted poor prognosis of OSCC patients [27–31]. Here, utilizing Cox regression analysis, we identified 3 immune cells with predictive significance in the prognosis of OSCC. The higher the risk score, the lower the level of immune cells in OSCC patients, and the worse the overall viability, suggesting that 3 immune cell function defects are closely related to the prognosis of oral squamous cell carcinoma patients. We observed that the infiltration level of myeloid dendritic cell, T cell CD4+ central memory, and common myeloid progenitor was strongly associated with better survival. Myeloid dendritic cells (MDCs), one of the major professional antigen-presenting cells, play a crucial determinant role in the immune responses to antigens [32]. The role of MDCs in OSCC has not been extensively studied. Recently, it is reported that memory T cells in tumors are an independent prognostic factor for several human malignancies [33, 34]. Memory T cells in the tumors may influence the development of lymphatic metastasis of OSCC [35]. Induction of memory T cells may be beneficial leading to cancer therapy. The common myeloid progenitor is a critical part in the immune regulation of tumors, where they are able to suppress both innate and adaptive immunity, including immunity to cancer [36]. While in the results of our analysis, high infiltration level of myeloid progenitor was associated with good prognosis in OSCC patients. Therefore, the clinical importance of these immune cells in OSCC may be ubiquitous.Our study also has some limitations. First, our study mainly focused on the key genes and immune cells derived from Cox regression analysis, and other genes were not studied in depth. Second, we lacked clinical samples to validate the prognostic value of key genes and immune cells. More studies are needed to explore the novel therapeutic targets for OSCC for key genes and immune cells we identified. Finally, the Cox model we constructed will be validated in a large clinical cohort in the further.


## 5. Conclusion

Some of the genes found in this study already have prognostic significance and may be repositioned as targets for the treatment of OSCC. We used a bioinformatics approach to screen core genes and explore the impact of immune cell infiltration on the prognosis of OSCC patients. This study broadens the existing findings and reveals genes and immune cells associated with OSCC, opening the way for the development of novel targeted therapies.



---

*Source: 5802110-2021-12-30.xml*