Seven Immune-Related Genes’ Prognostic Value and Correlation with Treatment Outcome in Head and Neck Squamous Cell Carcinoma

Abstract

Background. Head and neck squamous cell carcinoma (HNSCC) is a growing concern worldwide, due to its poor prognosis, low responsiveness to treatment, and drug resistance. Since immunotherapy effectively improves HNSCC patients’ survival status, it is important to continuously explore new immune-related predictive factors to accurately predict the immune landscape and clinical outcomes of individuals suffering from HNSCC. Methods. The HNSCC transcriptome profiling of RNA-sequencing data was retrieved from TCGA database, and the microarray of GSE27020 was obtained from the GEO database for validation. The differentially expressed genes (DEGs) between HNSCC and normal samples were identified by multiple test corrections in TCGA database. The univariate and multivariate Cox analyses were performed to identify proper immune-related genes (IRGs) to construct a risk model. The Cox regression coefficient was employed for calculation of the risk score (RS) of IRG signature. The median value of RS was utilized as a basis to classify individuals with HNSCC into high- and low-risk groups. The Kaplan-Meier (K-M) survival analysis and receiver operating characteristic (ROC) curves were employed for the identification of the prognostic significance and precision of the IRG signature. The signature was also evaluated based on clinical variables, predictive nomogram, mutation analysis, infiltrating immune cells, immune-related pathways, and chemotherapeutic efficacy. The protein-protein interaction (PPI) network and functional enrichment pathway investigations were utilized to explore possible potential molecular mechanisms. Finally, the hub gene’s differential mRNA expression levels were evaluated by means of the Gene Expression Profiling Interactive Analysis (GEPIA), and the Human Protein Atlas (HPA) was utilized for the validation of their translational levels. Results. Collectively, 1593 DEGs between HNSCC and normal samples were identified, of which 136 IRGs were differentially expressed. Then, the 136 immune-related DEGs were mostly enriched in the cytokine-related signaling pathways by GO and KEGG analyses. After that, a valuable signature based on seven genes (DKK1, GAST, IGHM, IL12RB2, SLURP1, STC2, and TNFRSF4) was designed. The HNSCC patients into the low-risk group and the high-risk group were divided by using the median RS; the HNSCC patients in the high-risk group had a worse survival than those in the low-risk group. The risk signature was verified to be an independent predictive marker for HNSCC patients. Meanwhile, the RS had the largest contribution to survival of these patients based on the predictive nomogram. In addition, the low-risk HNSCC patients exhibited significantly enriched immune cells, along with an association with high chemosensitivity. Conclusion. The constructed gene signature can independently function as a predictive indicator for the clinical features of HNSCC patients. The low-risk HNSCC subjects might benefit from immunotherapy and chemotherapy.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

As the sixth most prevalent type of malignancy, head and neck squamous cell carcinoma (HNSCC) is the seventh main cause of cancer-related mortalities globally [1]. A study conducted in the United States predicted that by 2022, approximately 66,470 new HNSCC cases would arise and 15,050 HNSCC-related deaths would occur [2]. Despite steady advancements in relevant medical treatments, like surgery, radiotherapy, and chemotherapy, the five-year survival of individuals with HNSCC has not significantly improved [3]. Therefore, finding new and innovative novel prognostic factors for HNSCC patients is an urgent need.

HNSCC is considered an immunodeficiency disease. The main mechanisms underlying the disease include the induction of immune tolerance, local immune escape, and the destruction of T-cell signals [4]. The immune microenvironment of HNSCC has been widely studied [5, 6]. For instance, the human leukocyte antigen (HLA) is responsible for the vital function of transmitting signals between tumor antigen peptides and killer T cells [7]. A previous study demonstrated that more than 50% of HNSCC patients had low HLA expression, with extensive lymph node metastasis and poor prognosis [8]. HNSCC tumor cells could also release chemical factors, to induce many immunosuppressive hematological cells to enter the immune microenvironment, thus suppressing the immune response [9].

Meanwhile, studying how HNSCC survival is linked to infiltrating immune cell proliferation and function could improve the survival of HNSCC patients [10–12]. Immunotherapy’s effect on the clinical outcomes of individuals with HNSCC has also been intensively studied [13–15]. This is why the exploration of immune-related biomarkers to anticipate the clinical features of individuals with HNSCC is imperative. Recent studies demonstrated that immune-related biomarkers could affect the biological behavior of HNSCC as well as the status of patients. For instance, Yao et al.’s model consisted of four immune-related genes (IRGs), including PVR, TNFRSF12A, IL21R, and SOCS1 [16]; Chen et al. constructed predictive model based on three IRGs (SFRP4, CPXM1, and COL5A1) [17]; and Zhang et al. established a model based on six IRGs (PLAU, STC2, TNFRSF4, PDGFA, DKK1, and CHGB) for the prognostic prediction of HNSCC [18]. Although these studies have constructed a proper model to predict the prognosis of patients with HNSCC, the progression of HNSCC is complexity and uncertainty. Therefore, to date, reliable and predictive biomarkers for identifying HNSCC are still limited; it is essential to continuously look for newly representative biomarkers.

During this research, a risk signature of seven immune-related genes was developed for accurately predicting the clinical outcomes of HNSCC subjects, which may provide an effective treatment strategy for these patients.

2. Methods

2.1. Dataset

The transcriptome profiling of RNA-sequencing (FPKM) was attained from The Cancer Genome Atlas (TCGA) (https://cancergenome.nih.gov/) containing 502 cancerous and 44 healthy tissues samples, along with the relevant clinical information for 528 HNSCC patients. In total, 2483 IRGs were procured from The ImmPort database (https://immport.niaid.nih.gov/) [19, 20]. Besides, the microarray and clinical information of GSE27020 containing 109 HNSCC samples were provided by the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/gds/) for authentication.

2.2. Development of the IRG-Based Signature

The differentially expressed genes (DEGs) between cancerous and healthy samples were identified with the help of TCGA database after multiple test corrections by false discovery rate $\text{FDR}<0.05$ and $\log \text{fold}\text{change}\text{FC}>2$ [21]. Then, screening the intersections of these DEGs with IRGs was carried out. Univariate Cox proportional hazards regression analysis was employed for the determination of the predictive ability of IRGs for the overall survival (OS) of individuals with HNSCC with the aid of the “survival” package in R [22]. The genes with a threshold of $P<0.05$ were subjected to further evaluation using multivariate Cox regression analysis [23]. Then, the expression levels of the hub genes were compared for further exploring their expression features in normal and HNSCC tumor tissues. Subsequently, the calculation of the IRG-based signature-related risk score (RS) was done using the Cox regression coefficient and gene expression formula given below: $$\begin{array}{}\text{(1)}& \text{RS}=\sum _{i=1}^N\text{Expi}\ast \text{Co}‐\text{eff}\end{array}$$

$N$, Expi, and Co-eff indicate signature gene number, gene expression levels, and regression coefficient values, respectively. Using the median value of RS as a criterion, the individuals with HNSCC were classified into the low- and high-risk groups.

2.3. Prognosis Prediction by the IRG-Based Signature

For the verification of the IRG-based signature’s prognostic performance, the signature’s impact on OS in both risk groups was subjected to comparison by Kaplan-Meier (K-M) survival analysis utilizing the “survival” package in R software, followed by the area under the curve (AUC) of the receiver operating characteristic (ROC) curves to assess IRG-based signature’s accuracy through the “survival ROC” package in R software [24]. Subsequently, the GEO database (GSE27020) was used for external validation. It was also subjected to K-M survival analysis and ROC curves to identify the signature’s prognostic value and precision. The RS distribution and survival status of individuals with HNSCC in TCGA were constructed to further understand the prognostic capability of the signature.

2.4. Correlation Analysis

The association of IRG-based signature with clinical variables (age, pathological grade, gender, and tumor and TNM stages) was analyzed. Moreover, the clinical factors and robustness of the signature in predicting the OS were demonstrated by employing univariate and multivariate Cox regression analyses.

2.5. Construction of a Predictive Nomogram

Based on clinical variables, a nomogram was used to establishment a prognostic scoring system for predicting survival in HNSCC patients both in TCGA and GSE27020 databases.

2.6. Somatic Mutation Analysis

We obtained the somatic mutation profiles of all tumor samples from TCGA database and explored the mutation analysis for 528 patients. The R software “maftools” package was utilized to analyzed and visualized for mutation data of the low-risk group and the high-risk group.

2.7. Immune Microenvironment Analysis

Considering the involvement of infiltrating immune cells in the tumor microenvironment, the single-sample gene set enrichment analysis (ssGSEA) algorithm was utilized to evaluate the immune score of each HNSCC sample from TCGA and GSE27020 databases [25]. The different proportions of the infiltrating immune cells between the low- and high-risk groups were assessed by the Wilcoxon test. Moreover, the association between the RS and immune-related biological functions was performed for further exploring the underlining mechanisms. The gene expression profiles corresponding to samples of TCGA and GSE27020 databases were selected to perform the gene set variation analysis (GSVA).

2.8. Prediction of Clinical Application

The calculation of the half inhibitory concentration (IC50) of common chemotherapeutic agents was done, and the differences in the IC50 across the two risk groups were also evaluated for predicting the clinical application of the IRG-based signature both in TCGA and GSE27020 databases.

2.9. Molecular Mechanism Analysis

The STRING biological database (https://string-db.org/) was applied for extraction of the protein-protein interaction (PPI) network [26] as a mathematical representation of the physical contacts among differentially expressed IRGs linked to HNSCC patient survival. Thereafter, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were employed in determining the potential function of the immune-related DEGs [27].

2.10. Investigating the Expression of Hub Genes

The Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/) was utilized for the investigation of the differential mRNA expression profiles of the hub genes in the IRG-based signature. Moreover, the Human Protein Atlas (HPA) (https://www.proteinatlas.org/) was employed for the purpose of validating the translational levels of these hub genes.

2.11. Statistical Analysis

R version 3.6.2 was utilized to conduct statistical analysis procedures. DEGs were compared with multiple test corrections with $\text{FDR}<0.05\text{and}\text{logFC}>2$ were viewed as being dramatically dysregulated. The survival curves were estimated by using the K-M survival analysis and log-rank test between different groups. Clinicopathological features were compared by univariate and multivariate Cox regression analyses. The ssGSEA algorithm and Wilcoxon test were used to compare different proportions of the infiltrating immune cells between different groups. The $t$-test or Wilcoxon test for comparisons of two variables, and a $P<0.05$ (two-side) was taken as a statistically significant value.

3. Results

3.1. Differential Gene Expression Analysis

TCGA database was employed to retrieve the HNSCC RNA-sequencing data comprising 502 tumor samples and 44 healthy samples. Among these patients, 528 HNSCC subjects with gene expression profiles and clinical follow-up data were included. The workflow of this research is illustrated in Figure 1. In the differential gene expression analysis, 1593 DEGs between HNSCC and healthy samples were identified (Figures 2(a) and 2(b)), of which 136 IRGs were differentially expressed (Figures 2(c) and 2(d)). Among these genes, 13 genes were identified to predict survival in the univariate Cox regression analysis. To study their interactions, the STRING biological database was utilized to construct the PPI network, containing 11 nodes and 22 edges. Based on the degree of genes, IL1A, CTLA4, CCR8, IL12RB2, TNFRSF4, CXCL13, and PLAU appeared to be the core genes among these IRGs (Figure 3(a)). As for functional analysis in Figure 3(b), the 136 immune-related DEGs were mostly enriched in immune response/cytokine mediation (BP), immunoglobulin complex/external side of plasma membrane (CC), and cytokine activity/signaling receptor activator activity/receptor ligand activity (MF) by GO analysis. By KEGG pathway analysis, the genes were mostly enriched in cytokine-cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, and chemokine signaling pathway.

[figure(s) omitted; refer to PDF]

3.2. Development and Verification of the IRG-Based Prognostic Signature

The detailed characteristics along with population demographics are given in Table 1. To develop a predictive IRG-based signature, seven IRGs were chosen after univariate and multivariate Cox analyses (Table 2). Meanwhile, the expression levels of the seven IRGs were further investigated. Compared with normal tissues, only SLURP1 was downregulated in tumor tissues, while the expression levels of DKK1, GAST, IGHM, IL12RB2, STC2, and TNFRSF4 were upregulated in tumor tissues (Figure 3(c)). Then, the calculation of the RS of this IRG-based signature was done as follows: $\text{RS}={0.006062}^{\ast }DKK1+{0.010886}^{\ast }GAST+-{0.000928}^{\ast }IGHM+-{0.051088}^{\ast }IL12RB2+-{0.001863}^{\ast }SLURP1+{0.025190}^{\ast }STC2+-{0.089341}^{\ast }TNFRSF4$.

Table 1

Comparison of characteristics between TCGA database and the GEO database (GSE27020).

Characteristic	TCGA		GSE27020		$P^{\ast }$
	Frequency	Deaths (%)	Frequency	Deaths (%)
Total	528	170 (32.20)	109	34 (31.19)
Age (years)					0.063
≤65	345	99 (28.70)	61	19 (31.15)
>65	182	71 (39.01)	48	15 (31.25)
Unknown	1	0
Sex					—
Male	386	114 (29.53)	—	—
Female	142	56 (39.44)	—	—
Grade					0.023
G1-G2	374	120 (32.09)	91	29 (31.87)
G3-G4	132	43 (32.58)	16	5 (31.25)
Unknown	22	7 (31.82)	2	0
Stage					—
Stage I-stage II	120	40 (33.33)	—	—
Stage III-stage IV	394	130 (32.99)	—	—
Unknown	14	0
T					—
T1-T2	189	53 (28.04)	—	—
T3-T4	323	117 (36.22)	—	—
Unknown	16	0
M					—
M0	496	166 (33.47)	—	—
M1	6	3 (50.00)	—	—
Unknown	26	1 (3.85)
N					—
N0	246	88 (35.77)	—	—
N1-N3	260	81 (31.15)	—	—
Unknown	22	1 (4.55)

TCGA: The Cancer Genome Atlas; GEO: Gene Expression Omnibus. ${}^{\ast }$Chi-square test for the comparison of characteristics between TCGA database and the GSE27020 database for each clinical variable.

Table 2

The detailed information of the immune-related gene signature for the survival of HNSCC patients.

Gene name	Coef	HR (95% CI)	$P$ value
DKK1	0.006062	1.006 (0.998-1.014)	0.118
GAST	0.010886	1.011 (1.001-1.021)	0.024
IGHM	-0.000928	0.999 (0.998-1.000)	0.108
IL12RB2	-0.051088	0.950 (0.894-1.010)	0.099
SLURP1	-0.001863	0.998 (0.997-1.000)	0.012
STC2	0.025190	1.026 (1.006-1.046)	0.011
TNFRSF4	-0.089341	0.915 (0.831-1.007)	0.069

HNSCC: head and neck squamous cell carcinoma; Coef: regression coefficient value; HR: hazard ratio; CI: confidence interval.

Moreover, individuals with HNSCC were categorized into the low-risk group ($n=249$) and the high-risk group ($n=249$) as per the median RS. HNSCC patients at high risk showed a worse survival in the K-M analysis, in comparison to the patients at low risk ($P<0.001$) (Figure 4(a)), with the AUC of 0.685 for the 5-year ROC curve (Figure 4(b)), indicating a certain predictive value of the signature in predicting the survival of individuals with HNSCC. Meanwhile, this IRG-based signature was validated in the GEO database (GSE27020) of 109 HNSCC patients who were also grouped into the low-risk group ($n=54$) and the high-risk group ($n=55$). Consistent with TCGA database, the K-M analysis of the GEO data exhibited that high-risk HNSCC individuals presented a worse survival in comparison with the low-risk group ($P<0.05$) (Figure 4(c)), with the AUC of 0.637 for 5-year ROC curve (Figure 4(d)).

[figure(s) omitted; refer to PDF]

Additionally, there were more deaths in HNSCC patients with the elevation in the value of RS (Figures 4(e) and 4(f)), and the seven genes showed differences in mRNA expression across the two groups in the heat map (Figure 4(g)).

3.3. Use of IRG-Based Signature as an Independent Prognostic Marker for HNSCC

Univariate and multivariate Cox regression analyses were conducted for assessing the correlations between the IRG-based signature and clinical variables (age, gender, grade, tumor stage, and TNM stage). The findings indicated that OS of HNSCC patients was significantly associated with age ($\text{HR}=1.022,95\%\text{CI}=1.008–1.037,P=0.003$), M stage ($\text{HR}=3.595,95\%\text{CI}=1.137–11.370,P=0.029$), and the RS calculated from the IRG-based signature ($\text{HR}=1.650,95\%\text{CI}=1.452–1.876,P<0.001$) (Figure 5(a)) in univariate Cox regression analysis and also with age ($\text{HR}=1.020,95\%\text{CI}=1.005–1.036,P=0.010$), M stage ($\text{HR}=4.643,95\%\text{CI}=1.334–16.163,P=0.016$), and the RS ($\text{HR}=1.634,95\%\text{CI}=1.432–1.864,P<0.001$) in multivariate Cox regression analysis (Figure 5(b)).

[figure(s) omitted; refer to PDF]

Furthermore, the associations between clinicopathological parameters and the RS and associations between the seven genes and clinical variables were also evaluated (Table 3, Figure 5(c)). The results revealed that the male gender, high tumor stage, and T stage were linked to a greater value of RS. In addition, mRNA expression levels of IGHM and SLURP1 appeared to be elevated in females in comparison to males. The mRNA expression level of STC2 appeared to be elevated in males when compared with females, and a high pathological grade was correlated with lower mRNA expression of GAST and SLURP1. The results also suggested that higher mRNA expression of IGHM was remarkably linked to a high grade. Greater mRNA expression of GAST was significantly linked to a more advanced tumor stage. The elevated mRNA expression level of DKK1 and GAST was correlated with the advanced T stage. Moreover, lower mRNA expression of IGHM, as well as higher mRNA expression of SLURP1 and STC2, was correlated with the advanced M stage.

Table 3

Correlations between the seven immune-related genes and clinical characteristics.

Gene name	Age ($P$ value)	Sex ($P$ value)	Grade ($P$ value)	Stage ($P$ value)	T ($P$ value)	M ($P$ value)	N ($P$ value)
DKK1	1.314 (0.190)	-0.346 (0.730)	-1.882 (0.062)	-1.725 (0.086)	-2.233 (0.026)	-0.502 (0.649)	-1.181 (0.238)
GAST	-1.550 (0.123)	-0.787 (0.432)	4.759 (<0.001)	-2.389 (0.017)	-3.773 (<0.001)	0.013 (0.990)	-0.107 (0.915)
IGHM	0.207 (0.836)	-2.370 (0.018)	-2.271 (0.025)	0.009 (0.993)	1.094 (0.275)	5.116 (<0.001)	-1.114 (0.266)
IL12RB2	-0.347 (0.729)	1.676 (0.095)	-0.499 (0.618)	1.820 (0.071)	1.411 (0.160)	0.459 (0.677)	1.261 (0.208)
SLURP1	0.071 (0.943)	2.221 (0.028)	2.717 (0.007)	1.622 (0.107)	0.632 (0.528)	5.777 (<0.001)	1.517 (0.130)
STC2	-1.203 (0.230)	-2.498 (0.013)	-1.658 (0.100)	-1.140 (0.256)	-0.647 (0.518)	3.320 (0.015)	0.085 (0.932)
TNFRSF4	1.154 (0.250)	-1.014 (0.311)	-1.484 (0.140)	-0.900 (0.368)	1.407 (0.161)	0.662 (0.535)	-1.381 (0.169)

3.4. Construction of a Predictive Nomogram

The clinical variables and RS were included in the nomogram. As indicated in the nomogram, the RS had the largest contribution to survival of patients with HNSCC both in TCGA and GSE27020 databases (Figures 6(a) and 6(b)).

[figure(s) omitted; refer to PDF]

3.5. Somatic Mutation Analysis

We obtained somatic mutation profiles of 528 patients in TCGA database. Around 241 (97.57%) and 229 (93.47%) samples possessed somatic mutations in the high-risk and low-risk groups, respectively. The top 30 mutated genes for high-risk and low-risk groups are shown in Figures 7(a) and 7(b). The results indicated that the TP53 mutated most frequently approximately accounting for 78% and 62% in the high-risk and low-risk groups, respectively.

[figure(s) omitted; refer to PDF]

3.6. Immune Microenvironment Analysis

The association between 23 immune cells infiltration differences and different risk groups was analyzed in TCGA and GSE27020 databases. Patients in the low-risk group showed higher infiltration levels of these 23 immune cells in TCGA database (Figure 8(a)), and patients in the low-risk group were more correlated with the infiltration of activated CD8 T cell, activated dendritic cell, CD56dim natural killer cell, eosinophil, immature B cell, mast cell, MDSC, monocyte, natural killer cell, natural killer T cell, neutrophil, T follicular helper cell, type 1 T helper cell, and type 17 T helper cell in the GSE27020 database (Figure 8(b)).

[figure(s) omitted; refer to PDF]

In addition, the relationship between immune pathway scores and RS were analyzed in order to better explore the immune-related biological functions. Functions with a correlation greater than 0.2 and $P<0.05$ are shown in Supplementary Figure 1. The results indicated that 14 immune-related pathways were correlated negatively with the RS in TCGA database (Supplementary Figure 1A). In the GSE27020 database, 8 immune-related pathways were correlated negatively with the RS, while 1 was correlated positively (Supplementary Figure 1B). These immune-related pathway scores vary with increasing levels of RS, implying that an imbalance in these pathways is closely related to tumor development.

3.7. Prediction of Clinical Application

The association of risk with the therapeutic efficacy of common chemotherapeutic agents in HNSCC was also studied. The findings exhibited that the low-risk HNSCC patients presented increased sensitivity to Elesclomol, GW843682X, Midostaurin, Pazopanib, QS11, and Salubrinal in TCGA database (Figure 9(a)), and the low-risk group was more likely with higher sensitivity of Bexarotene, BI.2536, MG.132, QS11, Salubrinal, and Thapsigargin in the GSE27020 database (Figure 9(b)). The results indicated that HNSCC patients with low risk represented higher sensitivity to chemotherapy.

[figure(s) omitted; refer to PDF]

3.8. Investigation of the Expression of the Seven IRGs

The expression of the seven IRGs in HNSCC was explored with the help of the GEPIA database. The expression levels of the seven IRGs varied remarkably across cancerous and healthy tissues (Figure 10(a)). However, to validate these findings, more experimental analyses were required. Moreover, the HPA database was employed to investigate the expression of the seven IRGs at the translation level. Among the seven IRGs, expressions of IGHM and SLURP1 were lower in the HNSCC tissues. Moreover, STC2 showed higher expression in HNSCC by immunohistochemistry. No remarkable variations were observed in the expressions of GAST, IL12RB2, and TNFRSF4 across normal and HNSCC tissues, while DKK1 was not detected by immunohistochemistry in the HPA database (Figure 10(b)). However, to further validate the translational relevance of the seven IRGs on HNSCC, more clinical analyses on HNSCC samples are needed.

[figure(s) omitted; refer to PDF]

4. Discussion

During this research, an IRG-based signature was established, which was capable of anticipating the clinical landscapes of HNSCC patients and correlated with clinicopathological characteristics of affected individuals, the numbers of tumor-infiltrating immune cells, and the efficacy of common chemotherapeutics. These findings suggested that this signature may be valuable for predicting HNSCC-related prognosis and provide good clinical application in immunotherapy and chemotherapy.

The IRG-based signature consisted of seven genes (i.e., DKK1, GAST, IGHM, IL12RB2, SLURP1, STC2, and TNFRSF4). Among them, DKK1 is a member of the DKK family and regulates cell proliferation, migration, and apoptosis in various tumor tissues through β-catenin-dependent and β-catenin-independent mechanisms [28]. Moreover, as a tumor suppressor gene, DKK1 causes apoptosis and suppresses cell proliferation [29]. Gao et al. suggested that elevated DKK1 expression levels can predict poor prognosis in HNSCC patients [30]. STC2 regulates tumor cell proliferation, apoptosis, and angiogenesis and is also vital for the invasiveness and metastasis of HNSCC [31]. IL12RB2 is a subunit of the IL-12 receptor, and an increased ratio of IL12RB2-positive tumor-infiltrating lymphocytes is indicative of a good prognosis in laryngeal cancer [32]. SLURP1 belongs to the Ly6/uPAR family that lacks a GPI-anchoring signal sequence and is associated with a poor prognosis of HNSCC [33]. Furthermore, one of the tumor necrosis factor receptors, TNFRSF4, could be a useful target for immunotherapy of HNSCC [34]. Although there are no published reports on GAST and IGHM for HNSCC, these genes may be related to tumorigenesis and development [35, 36]. In general, these previous findings emphasize the importance of these seven genes in HNSCC prognosis prediction. Furthermore, the expression levels of GAST, IL12RB2, and TNFRSF4 in HNSCC samples appeared to be elevated in healthy tissues from the GEPIA database, while no apparent variations were observed between cancerous and healthy tissues from the HPA data. Except for SLURP1 and STC2, IGHM expression in HNSCC tissues was remarkably increased compared to that in healthy samples from the GEPIA database, which was inconsistent with the HPA database. This could be due to abnormal methylation. However, further experimentation is required to confirm this finding.

In the multivariate analysis for the associations between clinicopathological factors and the risk IRG-based signature, a high-immune RS was linked to a high tumor stage and T stage. Also, the signature predicted the possible clinical features of HNSCC subjects, likely by regulating the tumor immune microenvironment. The tumor-infiltrating immune cells are known to be correlated with the progression and OS of HNSCC subjects [37], and a high level of infiltrating immune cells is often a good predictor for the OS of patients [38, 39]. Therefore, the risk IRG-based signature is expected to correlate with infiltrating immune cells. As expected, low-risk HNSCC patients had increased infiltration rates of 23 immune cells, indicating the effectiveness of immunotherapy in the low-risk category compared with that in the high-risk category. Owing to the importance of chemotherapy in HNSCC, the IC₅₀ values of various chemotherapeutic agents were compared in the two groups. A lower RS was linked to a higher IC₅₀ of QS11 and Salubrinal in both TCGA and GSE27020 databases. The QS11 is an inhibitor of ADP-ribosylation factor GTPase-activating protein 1, which modulates Wnt/β-catenin signaling through an effect on protein trafficking [40]; and the Salubrinal is a selective cell complex inhibitor that inhibits endoplasmic reticulum stress-mediated apoptosis. Despite these two drugs are not commonly used as chemotherapy drugs for HNSCC, the findings of this study may be valuable for future research.

This study had some limitations. (i) Owing to limited HNSCC samples in TCGA and GSE27020 databases, an issue of the time period was evident in this study. (ii) The analyses were performed using publicly available data from retrospective studies, and the outcomes must be validated in further research with larger samples and functional experiments. (iii) There is a need for further exploration of other possible predictive factors linked to clinical outcomes in HNSCC. (iv) There is a need for further investigation of the mechanisms underlying the functions of the IRG-based signature in HNSCC. Bioinformatics analysis with a specific reference value was used as a basis to conclude this study. Further corresponding molecular experiments are required to validate these findings.

To conclude, a risk signature based on seven IRGs (DKK1, GAST, IGHM, IL12RB2, SLURP1, STC2, and TNFRSF4) was developed. This signature serves as a potential biological marker and treatment target for immunotherapy and chemotherapy of HNSCC. These findings may facilitate future studies on the molecular mechanisms of HNSCC.

Acknowledgments

This work was supported by the Research Program of Shenzhen Innovation Council (grant number JCYJ20200109140208058), Shenzhen Fund for Guangdong Provincial High-level Clinical Key Specialties (grant number SZGSP008), and Sanming Project of Medicine in Shenzhen (grant number SZSM202111012, Oral and Maxillofacial Surgery Team, Professor Yu Guangyan, Peking University Hospital of Stomatology). The authors thank all data providers, patients, investigators, TCGA database, ImmPort database, GEO database, STRING biological database, GEPIA database, HPA database, and institutions involved in these studies.

References

[1] H. Liu, G. Li, E. M. Sturgis, S. Shete, K. R. Dahlstrom, M. Du, C. I. Amos, D. C. Christiani, P. Lazarus, Q. Wei, "Genetic variants in CYP2B6 and HSD17B12 associated with risk of squamous cell carcinoma of the head and neck," International Journal of Cancer, vol. 151 no. 4, pp. 553-564, DOI: 10.1002/ijc.34023, 2022.

[2] R. L. Siegel, K. D. Miller, H. E. Fuchs, A. Jemal, "Cancer statistics, 2022," CA: a Cancer Journal for Clinicians, vol. 72 no. 1,DOI: 10.3322/caac.21708, 2022.

[3] C. R. Leemans, P. J. F. Snijders, R. H. Brakenhoff, "The molecular landscape of head and neck cancer," Nature Reviews. Cancer, vol. 18 no. 5, pp. 269-282, DOI: 10.1038/nrc.2018.11, 2018.

[4] D. E. Johnson, B. Burtness, C. R. Leemans, V. W. Y. Lui, J. E. Bauman, J. R. Grandis, "Head and neck squamous cell carcinoma," Nature reviews Disease primers, vol. 6 no. 1,DOI: 10.1038/s41572-020-00224-3, 2020.

[5] M. Badr, K. Jöhrens, M. Allgäuer, M. Boxberg, W. Weichert, I. Tinhofer, C. Denkert, P. Schirmacher, A. Stenzinger, J. Budczies, "Morphomolecular analysis of the immune tumor microenvironment in human head and neck cancer," Cancer Immunology, Immunotherapy, vol. 68 no. 9, pp. 1443-1454, DOI: 10.1007/s00262-019-02378-w, 2019.

[6] S. M. Y. Chen, A. L. Krinsky, R. A. Woolaver, X. Wang, Z. Chen, J. H. Wang, "Tumor immune microenvironment in head and neck cancers," Molecular Carcinogenesis, vol. 59 no. 7, pp. 766-774, DOI: 10.1002/mc.23162, 2020.

[7] C. Fu, H. Zhao, Y. Wang, H. Cai, Y. Xiao, Y. Zeng, H. Chen, "Tumor-associated antigens: Tn antigen, sTn antigen, and T antigen," Hla, vol. 88 no. 6, pp. 275-286, DOI: 10.1111/tan.12900, 2016.

[8] R. L. Ferris, J. L. Hunt, S. Ferrone, "Human leukocyte antigen (HLA) class I defects in head and neck cancer: molecular mechanisms and clinical significance," Immunologic Research, vol. 33 no. 2, pp. 113-134, DOI: 10.1385/ir:33:2:113, 2005.

[9] M. Canning, G. Guo, M. Yu, C. Myint, M. W. Groves, J. K. Byrd, Y. Cui, "Heterogeneity of the head and neck squamous cell carcinoma immune landscape and its impact on immunotherapy," Frontiers in Cell and Development Biology, vol. 7,DOI: 10.3389/fcell.2019.00052, 2019.

[10] M. Fiedler, F. Weber, M. G. Hautmann, C. Bohr, T. E. Reichert, T. Ettl, "Infiltrating immune cells are associated with radiosensitivity and favorable survival in head and neck cancer treated with definitive radiotherapy," Oral Surgery, Oral Medicine, Oral Pathology, Oral Radiology, vol. 129 no. 6, pp. 612-620, DOI: 10.1016/j.oooo.2020.02.010, 2020.

[11] R. Duhen, C. Ballesteros-Merino, A. K. Frye, E. Tran, V. Rajamanickam, S. C. Chang, Y. Koguchi, C. B. Bifulco, B. Bernard, R. S. Leidner, B. D. Curti, B. A. Fox, W. J. Urba, R. B. Bell, A. D. Weinberg, "Neoadjuvant anti-OX40 (MEDI6469) therapy in patients with head and neck squamous cell carcinoma activates and expands antigen-specific tumor- infiltrating T cells," Nature Communications, vol. 12 no. 1,DOI: 10.1038/s41467-021-21383-1, 2021.

[12] Y. Peng, L. Xiao, H. Rong, Z. Ou, T. Cai, N. Liu, B. Li, L. Zhang, F. Wu, T. Lan, X. Lin, Q. Li, S. Ren, S. Fan, J. Li, "Single-cell profiling of tumor-infiltrating TCF1/TCF7 + T cells reveals a T lymphocyte subset associated with tertiary lymphoid structures/organs and a superior prognosis in oral cancer," Oral Oncology, vol. 119, article 105348,DOI: 10.1016/j.oraloncology.2021.105348, 2021.

[13] O. Kaidar-Person, Z. Gil, S. Billan, "Precision medicine in head and neck cancer," Drug Resistance Updates, vol. 40, pp. 13-16, DOI: 10.1016/j.drup.2018.09.001, 2018.

[14] J. D. Horton, H. M. Knochelmann, T. A. Day, C. M. Paulos, D. M. Neskey, "Immune evasion by head and neck cancer: foundations for combination therapy," Trends Cancer, vol. 5 no. 4, pp. 208-232, DOI: 10.1016/j.trecan.2019.02.007, 2019.

[15] E. E. W. Cohen, R. B. Bell, C. B. Bifulco, B. Burtness, M. L. Gillison, K. J. Harrington, Q. T. le, N. Y. Lee, R. Leidner, R. L. Lewis, L. Licitra, H. Mehanna, L. K. Mell, A. Raben, A. G. Sikora, R. Uppaluri, F. Whitworth, D. P. Zandberg, R. L. Ferris, "The Society for Immunotherapy of Cancer consensus statement on immunotherapy for the treatment of squamous cell carcinoma of the head and neck (HNSCC)," Journal for Immunotherapy of Cancer, vol. 7 no. 1,DOI: 10.1186/s40425-019-0662-5, 2019.

[16] Y. Yao, Z. Yan, S. Lian, L. Wei, C. Zhou, D. Feng, Y. Zhang, J. Yang, M. Li, Y. Chen, "Prognostic value of novel immune-related genomic biomarkers identified in head and neck squamous cell carcinoma," Journal for Immunotherapy of Cancer, vol. 8 no. 2,DOI: 10.1136/jitc-2019-000444, 2020.

[17] Y. Chen, Z. Y. Li, G. Q. Zhou, Y. Sun, "An immune-related gene prognostic index for head and neck squamous cell carcinoma," Clinical Cancer Research, vol. 27 no. 1, pp. 330-341, DOI: 10.1158/1078-0432.Ccr-20-2166, 2021.

[18] Y. Zhang, P. Chen, Q. Zhou, H. Wang, Q. Hua, J. Wang, H. Zhong, "A novel immune-related prognostic signature in head and neck squamous cell carcinoma," Frontiers in Genetics, vol. 12, article 570336,DOI: 10.3389/fgene.2021.570336, 2021.

[19] T. Liu, H. Wu, J. Qi, C. Qin, Q. Zhu, "Seven immune-related genes prognostic power and correlation with tumor- infiltrating immune cells in hepatocellular carcinoma," Cancer Medicine, vol. 9 no. 20, pp. 7440-7452, DOI: 10.1002/cam4.3406, 2020.

[20] Z. Wang, Q. Song, Z. Yang, J. Chen, J. Shang, W. Ju, "Construction of immune-related risk signature for renal papillary cell carcinoma," Cancer Medicine, vol. 8 no. 1, pp. 289-304, DOI: 10.1002/cam4.1905, 2019.

[21] M. E. Ritchie, B. Phipson, D. Wu, Y. Hu, C. W. Law, W. Shi, G. K. Smyth, "Limma powers differential expression analyses for RNA-sequencing and microarray studies," Nucleic Acids Research, vol. 43 no. 7, article e47,DOI: 10.1093/nar/gkv007, 2015.

[22] C. Y. Li, J.-H. Cai, J. J. P. Tsai, C. C. N. Wang, "Identification of hub genes associated with development of head and neck squamous cell carcinoma by integrated bioinformatics analysis," Frontiers in Oncology, vol. 10,DOI: 10.3389/fonc.2020.00681, 2020.

[23] J. Yang, K. Xie, C. Li, "Immune-related genes have prognostic significance in head and neck squamous cell carcinoma," Life Sciences, vol. 256, article 117906,DOI: 10.1016/j.lfs.2020.117906, 2020.

[24] Z. Hu, J. Zhou, Y. Li, Y. Luan, H. Li, B. Jia, Z. Xie, B. Cheng, T. Wu, "Peripheral immune signature resembles tumor microenvironment and predicts clinical outcomes in head and neck squamous cell carcinoma," Frontiers in Immunology, vol. 13, article 915207,DOI: 10.3389/fimmu.2022.915207, 2022.

[25] J. X. Liang, Q. Chen, W. Gao, D. Chen, X. Y. Qian, J. Q. Bi, X. C. Lin, B. B. Han, J. S. Liu, "A novel glycosylation-related gene signature predicts survival in patients with lung adenocarcinoma," BMC Bioinformatics, vol. 23 no. 1,DOI: 10.1186/s12859-022-05109-8, 2022.

[26] N. T. Doncheva, J. H. Morris, J. Gorodkin, L. J. Jensen, "Cytoscape StringApp: network analysis and visualization of proteomics data," Journal of Proteome Research, vol. 18 no. 2, pp. 623-632, DOI: 10.1021/acs.jproteome.8b00702, 2019.

[27] G. Yu, L. G. Wang, Y. Han, Q. Y. He, "clusterProfiler: an R package for comparing biological themes among gene clusters," OMICS, vol. 16 no. 5, pp. 284-287, DOI: 10.1089/omi.2011.0118, 2012.

[28] G. Zhu, J. Song, W. Chen, D. Yuan, W. Wang, X. Chen, H. Liu, H. Su, J. Zhu, "Expression and role of Dickkopf-1 (Dkk1) in tumors: from the cells to the patients," Cancer Management and Research, vol. Volume 13, pp. 659-675, DOI: 10.2147/cmar.S275172, 2021.

[29] J. Niu, X. M. Li, X. Wang, C. Liang, Y. D. Zhang, H. Y. Li, F. Y. Liu, H. Sun, S. Q. Xie, D. Fang, "DKK1 inhibits breast cancer cell migration and invasion through suppression of β -catenin/MMP7 signaling pathway," Cancer Cell International, vol. 19 no. 1,DOI: 10.1186/s12935-019-0883-1, 2019.

[30] H. Gao, L. Li, M. Xiao, Y. Guo, Y. Shen, L. Cheng, M. Tang, "Elevated DKK1 expression is an independent unfavorable prognostic indicator of survival in head and neck squamous cell carcinoma," Cancer Management and Research, vol. Volume 10, pp. 5083-5089, DOI: 10.2147/cmar.S177043, 2018.

[31] H. F. Ma, G. X. Lv, D. H. Zhang, "miR-381 mediates the development of head and neck squamous cell carcinoma via targeting STC2," Oncotargets and Therapy, vol. Volume 13, pp. 4485-4493, DOI: 10.2147/ott.S246289, 2020.

[32] Y. Tao, N. Gross, Y. Liu, L. Zhang, G. Li, Z. Huang, J. Yang, "A high ratio of IL-12R β 2-positive tumor-infiltrating lymphocytes indicates favorable prognosis in laryngeal cancer," Oral Oncology, vol. 74, pp. 148-156, DOI: 10.1016/j.oraloncology.2017.10.006, 2017.

[33] H. Zheng, H. Liu, Y. Lu, H. Li, "Identification of a novel signature predicting overall survival in head and neck squamous cell carcinoma," Frontiers in Surgery, vol. 8, article 717084,DOI: 10.3389/fsurg.2021.717084, 2021.

[34] S. M. Loick, A. Fröhlich, J. Gabrielpillai, A. Franzen, T. J. Vogt, J. Dietrich, C. Wiek, K. Scheckenbach, S. Strieth, J. Landsberg, D. Dietrich, "DNA methylation and mRNA expression of OX40 (TNFRSF4) and GITR (TNFRSF18, AITR) in head and neck squamous cell carcinoma correlates with HPV status, mutational load, an interferon- γ signature, signatures of immune infiltrates, and survival," Journal of Immunotherapy, vol. 45 no. 4, pp. 194-206, DOI: 10.1097/cji.0000000000000407, 2022.

[35] E. Ness-Jensen, E. A. Bringeland, F. Mattsson, P. Mjønes, J. Lagergren, J. E. Grønbech, H. L. Waldum, R. Fossmark, "Hypergastrinemia is associated with an increased risk of gastric adenocarcinoma with proximal location: a prospective population-based nested case-control study," International Journal of Cancer, vol. 148 no. 8, pp. 1879-1886, DOI: 10.1002/ijc.33354, 2020.

[36] Z. Lu, Y. Gao, "Screening differentially expressed genes between endometriosis and ovarian cancer to find new biomarkers for endometriosis," Annals of Medicine, vol. 53 no. 1, pp. 1377-1389, DOI: 10.1080/07853890.2021.1966087, 2021.

[37] A. Almangush, I. Leivo, A. A. Mäkitie, "Overall assessment of tumor-infiltrating lymphocytes in head and neck squamous cell carcinoma: time to take notice," Acta Oto-Laryngologica, vol. 140 no. 3, pp. 246-248, DOI: 10.1080/00016489.2020.1720284, 2020.

[38] M. E. Spector, E. Bellile, L. Amlani, K. Zarins, J. Smith, J. C. Brenner, L. Rozek, A. Nguyen, D. Thomas, J. B. McHugh, J. M. G. Taylor, G. T. Wolf, for the University of Michigan Head and Neck SPORE Program, "Prognostic value of tumor-infiltrating lymphocytes in head and neck squamous cell carcinoma," JAMA Otolaryngology. Head & Neck Surgery, vol. 145 no. 11, pp. 1012-1019, DOI: 10.1001/jamaoto.2019.2427, 2019.

[39] Q. Xu, C. Wang, X. Yuan, Z. Feng, Z. Han, "Prognostic value of tumor-infiltrating lymphocytes for patients with head and neck squamous cell carcinoma," Translational Oncology, vol. 10 no. 1, pp. 10-16, DOI: 10.1016/j.tranon.2016.10.005, 2017.

[40] M. K. Singh, H. Gao, W. Sun, Z. Song, R. Schmalzigaug, R. T. Premont, Q. Zhang, "Structure-activity relationship studies of QS11, a small molecule Wnt synergistic agonist," Bioorganic & Medicinal Chemistry Letters, vol. 25 no. 21, pp. 4838-4842, DOI: 10.1016/j.bmcl.2015.06.062, 2015.