Basic information
Finally, a total of 522 patients with HNC were included for analysis. The basic character of the included patients was summarized in Table 1. The mean age was 60.87 years (range 18–90), approximately two-thirds of the patients were male, and about 66.67% and 8.93% of the patients had histories of drinking and smoking, respectively. Approximately one-third of the patients were HPV-positive. Pathological tumor node metastasis staging results were as follows: 33 (6.32%) cases in stage I, 98 (18.77%) cases in stage II, 105 (20.11%) cases in stage III, and 286 (54.79%) cases in stage IV. According to the ESTIMATE algorithm, stromal and immune scores were − 436.34 (interquartile range (IRQ), -904.04-136.41), and 361.16 (IQR, -160.20-1038.52), respectively. Average immune scores were higher than stromal scores in all pathology stages (p < 0.05) (Fig. 1). A previous study has shown that TP53 mutation was associated with decreased overall survival in patients with HNC 19. We then calculated immune scores and stromal scores in HNC patients with or without p53 mutation. Compared with p53 wildtype patients, the immune score in p53 mutation patients was significantly higher (p < 0.001), while no significant difference was noticed in stromal scores between these two groups of patients (Fig. 2A and 2B). Moreover, to determine the relationship between immune scores/stromal scores and all-cause mortality, patients with HNC were classified into high immune/stromal scores and low immune/stromal score groups based on the median value of immune scores or stromal scores. There was no significant difference in age, sex, HPV infection, or tumor pathology stage between the two groups (Table 1). Kaplan-Meier survival analysis revealed that the median overall survival time of patients in the immune low-score group was shorter than in the high-score group (Fig. 2C, 625 vs. 680 days, log-rank test, p = 0.1716). While compared to patients with higher stromal scores, patients with lower stromal scores also had shorter median overall survival times (Fig. 2D, 623 vs. 701 days, log-rank test, p = 0.8528).
Table 1
The basic information of the included head and neck cancer
|
Total patients
(n = 522)
|
Low immune score group (n = 261)
|
High immune score group (n = 261)
|
P value
|
Age
|
60.87(19–90)
|
60.3(19–88)
|
61.2(24–90)
|
0.776
|
Sex (male, %)
|
385(73.7%)
|
202(77.39%)
|
183(70.1%)
|
0.059
|
Alcohol
|
348(66.67%)
|
173(66.28%)
|
174(66.67%)
|
0.926
|
Smoking
|
20/224(8.93%)
|
11/116(9.48)
|
9/108(7.63%)
|
0.763
|
HPV
|
39/112(34.82%)
|
12/46(26.09%)
|
27/66(40.90%)
|
0.105
|
Stage
|
|
|
|
0.294
|
Stage I
|
35(6.71%)
|
21(4.02%)
|
14(2.68%)
|
|
Stage II
|
98(18.77%)
|
48(9.20%)
|
50(9.58%)
|
|
Stage III
|
105(20.11%)
|
58(11.11%)
|
47(9.00%)
|
|
Stage IV
|
286(54.79%)
|
134(25.67%)
|
152(29.12%)
|
|
Association between gene expression and ESTIMATE scores
To further understand the relationship between the gene expression and immune scores and/or stromal scores, we generated distinct gene expression profiles of HNC based on high vs. low immune/stromal score groups. The heat maps in Fig. 3A and 3B show the expression patterns of genes differentially changed in the immune and stromal score groups, respectively. Differentially expressed genes (DEGs) were then identified from the comparison of high vs. low immune score and stromal score groups. According to the immune scores, a total of 925 genes were upregulated, while a total of 72 genes were downregulated (FC > 2, p < 0.05). For the stromal scores, there were about 1414 genes that were upregulated, and 26 genes were downregulated (FC > 2, p < 0.05). The Venn diagrams were then generated and indicated that 425 genes were commonly upregulated, while only 2 genes were commonly downregulated in the high-scoring groups (Fig. 3C and 3D). A volcano plot of the acquired data highlights that more genes are upregulated (right side), and fewer genes are downregulated (left side) (Fig. 4A and 4B). To describe the potential function of the identified DEGs, the commonly upregulated genes (The genes are listed in Supplementary Table S1) in the high-immune scores group were selected for Gene Ontology (GO) enrichment analysis. The results indicated that the top GO terms associated with these genes included T cell costimulation, regulation of immune response, external side of the plasma membrane, integral component of the plasma membrane, transmembrane signaling receptor activity, receptor activity for biological processes (Fig. 3E), cellular component (Fig. 3F), and molecular function (Fig. 3G).
Association between DEGs and overall mortality
Kaplan-Meier survival curves were then obtained to further analyze the association between the individual DEGs and all-cause mortality. During an average of 1.77 years (IQR 1.04–3.23) follow-up, a total of 222 patients died. Kaplan-Meier survival analysis showed that from the 480 upregulated DEGs in the high-immune scores group, about 176 DEGs (Supplementary Table S2) were significantly associated with all-cause mortality (p < 0.05; representative genes are shown in Fig. 5).
Protein-protein interactions
To further search the hub genes of DEGs, we used the STRING tool to construct protein-protein interaction (PPI) networks and the MCODE app to identify gene clusters, respectively. The constructed PPI networks and the top five clusters are shown in Figs. 6 and 7, respectively. Then, we named these modules PTPRC, CD247, CD4, FCGR3A, and FCGR1A (Table 2) for ease of description. In the PTPRC module, the network forms 461 edges, involving 36 nodes; CTLA4, ITGAL, CD40LG, IL10RA, CD28, and SLAMF1 were the prominent nodes, which have the most connections with other members of the module. In the CD247 module, 26 nodes and 124 edges were observed. Among the clusters, CD247, CD3D, IKZF1, CD3G, PLEK, and HLA-DRB1 had higher degree values. For the CD4 module, CD4, IFNG, CXCR3, CD48, CD5, and CXCR5 had higher degree values. In addition, we performed a functional enrichment analysis for the five hub genes, and the results showed that 19 GO terms were significantly enriched in the biological process (BP), 7 GO terms were significantly enriched in the cellular compartment (CC), and 4 GO terms were significantly enriched in molecular function (MF) (FDR < 0.05, -log FDR > 1.30). Consistent with previous results, most of the genes are mainly involved in immune function. The top GO terms included T cell costimulation, regulation of immune response (Fig. 8A), external side of the plasma membrane, integral component of the plasma membrane (Fig. 8B), transmembrane signaling receptor activity, and receptor activity (Fig. 8C). Besides, the KEGG pathways analysis was also performed (Fig. 8D). In addition, A circos graph was constructed to visualize the prognosis-related genes involving in GO analysis according to the logFC value (Fig. 9).
Table 2
Functional roles of five hub genes with degree ≥ 8
Cluster
|
Nodes
|
Edges
|
Gene
|
1
|
36
|
461
|
SLAMF1, CD28, CD40LG, IL10RA, ITGAL, CTLA4, IL2RG, TLR7, TLR8, VCAM1, PDCD1, SPN, PRF1, GZMB, CD3E, CCL5, CCR5, KLRK1, CXCL10, CD8A, GZMA, CXCL9, LAG3, SELL, TBX21, CXCL13, KLRD1, CD52, CCL19, IDO1, CD2, KLRB1, CCR2, FASLG, PTPRC, CD27
|
2
|
26
|
124
|
TRAT1, PLEK, MPEG1, CYBB, HLA-DRB1, C1QA, C1QB, HLA-DMB, HLA-DMA, SIGLEC1, CD3D, HLA-DRB5, HLA-DQB1, HLA-DQA1, CD247, HLA-DPB1, HLA-DPA1, CD3G, IKZF1, GRAP2, GZMH, NKG7, HLA-DOA, HLA-DQA2,
|
3
|
24
|
82
|
SH2D1A, GBP5, ICOS, HCST, CD5, IFNG, GBP4, CIITA, CXCR3, FCGR1B, TNFSF13B, CD8B, GZMK, CXCL11, BTLA, CXCR5, CXCR6, IL21R, CD226, TIGIT, CD4, CD7, ITK, CD48
|
4
|
15
|
40
|
FCGR3A, EVI2B, MS4A6A, SASH3, DOCK2, NCKAP1L, CD84, FGL2, WDFY4, MNDA, GPR174, P2RY10, CD53, AOAH, CD74
|
5
|
9
|
16
|
CD163, CLEC10A, CD37, FCGR1A, IGSF6, IRF8, HLA-DRA, MS4A1, C1QC
|
The relationship between the hub genes and immune infiltration and survival
To further analyze the association of the five hub genes with immune infiltration. We found that all the 5 hub genes are significantly associated with tumor purity, B cell, CD8+ T cell, CD4+ T cell, macrophage, neutrophil, and dendritic cell (Fig. 10). Kaplan-Meier survival analysis showed that low PTPRC, CD247, and CD4 groups had poor overall survival except for FCGR3A and FCGR1A (Fig. 11). Also, multivariate Cox regression analysis further confirmed that PTPRC, CD247, and CD4 are associated with all-cause mortality (Table 3).
Table 3
Univariate and multivariate cox regression analysis of five hub genes
|
Overall survival
|
Univariate analysis
|
Multivariate analysis
|
Gene
|
HR (95% CI)
|
P Value
|
HR (95% CI)
|
P Value
|
PTPRC
|
0.884(0.809–0.996)
|
0.006
|
0.854(0.761–0.958)
|
0.007
|
CD247
|
0.809(0.72–0.909)
|
༜0.001
|
0.791(0.681–0.919)
|
0.002
|
CD4
|
0.903(0.818–0.996)
|
0.041
|
0.861(0.760–0.976)
|
0.019
|
FCGR3A
|
1.017(0.932–1.11)
|
0.703
|
0.979(0.882–1.087)
|
0.694
|
FCGR1A
|
0.997(0.879–1.13)
|
0.959
|
0.947(0.815-1.100)
|
0.474
|