The DKK1 mRNA Expression in Different Cancers
The mRNA expression level of DKK1 in cancers and normal tissues of various tumors were analyzed on Oncomine database. As shown in Figure 1A, higher mRNA expression level of DKK1 was presented in brain and CNS, head and neck, liver, esophageal, lung and pancreatic, while DKK1 had a lower expression in bladder and prostate cancers. The mRNA-seq data of TCGA analyzed by the TIMER database was used to validate the expression levels of DKK1 in cancers. DKK1 with a higher expression level existed in cholangiocarcinoma(CHOL),
colon adenocarcinoma(COAD), esophageal carcinoma(ESCA), head and neck squamous cell carcinoma(HNSC), liver hepatocellular carcinoma(LIHC), lung squamous cell carcinoma
(LUSC), thyroid carcinoma(THCA) and stomach adenocarcinoma(STAD). In contrast, the level of DKK1 was significantly declined in bladder urothelial carcinoma(BLCA), breast invasive carcinoma(BRCA), kidney chromophobe(KICH), kidney renal papillary cell carcinoma(KIRP) and prostate adenocarcinoma(PRAD). These data indicated that DKK1 had different expression levels in different cancers, suggesting that DKK1 exerted diverse functions in various cancers (Figure 1B). Besides, we used the UALCAN database to verify the findings in Oncomine and TIMER website and reported higher expression of DKK1 in HNSC tissues than in normal tissues(Figure 1C). Notably, DKK1 expression was associated with individual cancer stages, tumor grade, HPV status, nodal metastasis status in HNSC(Figures 1D–H).
Prognostic Significance of DKK1 mRNA Expression in Human Cancers
We explored the Kaplan-Meier plotter databases for the prognostic significance of DKK1 expression in human cancers. High level of DKK1 exhibited a poor survival period for patients with HNSC, KIRP, BRCA, LUAD, PAAD, STAD and LIHC(Figures 2A-G).While as for SARC and OV, highly expressed DKK1 has little effect on overall survival(Figures 2H, I). As Kaplan-Meier plotter analyzes only OS value, We evaluated the multiple clinical prognostic indicators of DKK1 in a variety of cancers by R project using “survival” and “ggplot2” packages. Forest plot showed DKK1 as a risk factor of HNSC, STAD, ACC, LUAD and PADD in OS, DSS and PFI.(Figure 3A-C).These findings demonstrated that DKK1 is a hazard for predicting worse prognostic in HNSC.
Correlation of DKK1 Expression With Clinical subtypes of Head and Neck Squamous Cell Cancer
Then, we evaluated the association of DKK1 expression with different clinicopathological factors of HNSC using the Kaplan-Meier Plotter database (Table 1).Upregulated levels of DKK1 was related to worse prognostic outcomes in females (OS: HR=2.49, P=0.0033), males (OS: HR=2.3,P=2.1e−06), stage 2 (OS: HR=2.94, P=0.0097), stage 3 (OS: HR=5.52, p=0.00022), stage 4 (OS: HR=1.73,P=0.0026) , grade1(OS: HR=3.81, P=0.0049), grade2(OS: HR=2.02, P=6.8e−05; RFS:HR=3.36, P=0.019) and grade3(OS: HR=2.69, P=0.00027) in HNSC. Similarly, High DKK1 level was correlated with poorer OS and RFS in high mutation burden (OS: HR=2.48, P=1.3e−06; RFS: HR=0.25, p=0.027) and low mutation burden(OS: HR=1.91, P=0.0026). From above mentioned clinicopathological factors, high expression of
DKK1 could be considered to be the worse prognostic indictor for HNSC.
DKK1 Expression is Correlated With Immune Infiltration in Head and Neck Squamous Cell Carcinoma
Although DKK1 was associated with prognosis and immune infiltration of hepatocellular carcinoma, esophageal cancer and pancreatic cancer, the relationship between DKK1 immune infiltration and head and neck squamous cell carcinoma is still unclear. The TIMER database was utilized to analyze the correlation between DKK1 level with immune infiltration levels in different cancer types. The results exhibited that DKK1 expression is significantly negative correlated with B cells (r=-0.198, p=1.37e-05) and CD8+T cells (r=-0.181, p=7.19e-05) in HNSC (Figure 4A). In addition, we investigated prognostic value of DKK1 level and tumor infiltrating immune cells in HNSC, using cox proportional hazard model by TIMER. The result reveal that DKK1 expression (p<0.001) was significantly correlated with clinical prognosis in HNSC (Table 2).Then, we utilized TISIDB database to furtherly survey the relationship between DKK1 level and 28 tumor immune infiltrating cell subtypes. These results indicated that DKK1 is associated with eleven immune cell subtypes in HNSC (Figure 4B, Table3). Especially, activated CD8 T cell (r=-0.143, p=0.00105), central memory CD8 T cell (r=0.179, p=4.1e-05), effector memory CD8 T cell (r=-0.137, p=0.00169), type 17 T helper cell (r=-0.115, p=0.00877), activated B cell (r=-0.187, p=1.78e-05), immature B cell(r=-0.156, p=0.000338), and DKK1 are moderately correlated (Figure 4C). These results strongly imply that DKK1 could serve as a major tumor immune infiltration regulator in HNSC.
DKK1 Expression Was Correlated With Markers of Immune Cells
We evaluated the correlation between DKK1 expression and tumor-infiltrating immune cell gene marker in HNSC tissues by checking the TIMER database.Our results showed that the DKK1 level in HNSC tissues negatively correlated with the marker genes of B cells and CD8+T cells. Notably, we found that the DKK1 was correlated with various subtypes of B cells markers, CD79A, CD19, T cells marker levels, including CD8+T markers, CD8A, CD8B, cell T (general) markers, CD3D, CD3E, CD2, T cell exhaustion marker, PDCD1, Dendritic cell markers, NRP1, CD1C, M1 Macrophage markers, PTGS2, IRF5, M2 Macrophage markers, VSIG4, CD163, Monocyte marker, CD86, Natural killer cell marker, KIR3DL1, Neutrophils markers, CCR7, ITGAM, TAM markers, CCL2, CD68, Th1 markers, TBX21, STAT4, IFNG , Th2 markers, GATA3, STAT6, STAT5A, and Treg markers, STAT5B, FOXP3, TGF-b, in HNSC (Table4). These findings indicate that DKK1 is involved in regulating tumor immune infiltration in HNSC.
Prognostic Value of DKK1 Expressions in Head and Neck Carcinoma Based on Immune Cells
This study illustrated that the DKK1 level was associated with the immune infiltration of HNSC. Also, upregulated DKK1 has a worse prognosis in HNSC patients.Therefore, we state a hypothesis that DKK1 may affect the prognosis of HNSC patients partly through immune infiltration.We use Kaplan-Meier plotter to analyze DKK1 expression in HNSC following B cells, CD4+ memory T-cells, CD8+ T-cells, regulatory T-cells, type 1 T-helper cells, type 2 T-helper cells, natural killer T-cells, mesenchymal stem cells and macrophage cells.We found that high DKK1 levels in HNSC in enriched B cells (p=2.7e−06), CD4+memory T cells (p=2e−08), CD8+T cells (1.6e−06), regulatory T cells (p=1.5e−05), type 1 T helper cells (p=3.3e-05), type 2 T-helper cells (p=7.6e-08 ), natural killer T cells (p= 0.011), mesenchymal stem cells(p=0.001), macrophages (p=1.7e−05) had a worse prognosis (Figures 5A-I). The above analysis suggested that immune infiltration may affect the prognosis of HNSC patients with high DKK1 expression to a certain extent.
Gene expression, Copy Number Variation, and Mutation feature analysis of DKK1 in HNSC
DKK1 expression was significantly rised in HNSC. We assessed the source of elevated DKK1 levels. DNA methylation, gene mutation, copy number variation were evidently related to genetic and epigenetic regulation, and were critically associated with the occurrence and development of tumors.We analyzed the genetic alterations of DKK1 in HNSC patients by using the cBioPortal online tool. We observed the genetic alteration status of DKK1 in different tumor samples in the TCGA cohort, which mainly include “mutation”, “amplification”, and “copy number deletion”. The “amplification” type of CNA was the primary type in the head and neck cancer cases, which show an alteration frequency of ~2.2% (Figure 6A-B). The types of the DKK1 genetic alteration were further presented in Figure 6B. Besides, we confirmed the copy number variation, somatic mutation, DNA methylation levels of the DKK1 in HNSC via the UCSC Xena database. The heatmap shows that the expression of DKK1 mRNA was involved in copy number variation and DNA methylation, but not with a somatic mutation in HNSC(Figure 6C). The cBioPortal results show that DKK1 mRNA is increased in the samples in which DKK1 is amplified(Figure 6D). So, we suggested that CNV and DNA methylation might contribute to the elevated level of DKK1 in HNSC.
Table 1
Correlation of DKK1 mRNA expression and prognosis in HNSC with different clinical subtypes by Kaplan-Meier plotter
Clinical subtypes
|
HNSC
|
OS
|
RFS
|
HR P
|
HR P
|
Gender
|
|
|
|
|
Female (n=4305)
|
2.49 (1.33-4.67)
|
0.0033
|
0.34 (0.11-1.02)
|
0.043
|
Male (n=3184)
|
2.3 (1.61-3.27)
|
2.1e−06
|
2.76 (0.96 − 7.99)
|
0.051
|
Stage
|
|
|
|
|
1 (n=1790)
|
6.71 (0.69-65.63)
|
0.06
|
0.35 (0.08-1.59)
|
0.16
|
2 (n=1707)
|
2.94 (1.25-6.93)
|
0.0097
|
0 (0-Inf)
|
0.11
|
3 (n=1269)
|
5.52(2.03-14.96)
|
0.00022
|
2.34 (0.6-9.08)
|
0.21
|
4 (n=676)
|
1.73 (1.2-2.47)
|
0.0026
|
-
|
-
|
Grade
|
|
|
|
|
1 (n=304)
|
3.81(1.41-10.29)
|
0.0049
|
2.36(0.47-11.94)
|
0.29
|
2 (n=1297)
|
2.02 (1.42-2.87)
|
6.8e−05
|
3.36 (1.15-9.87)
|
0.019
|
3 (n=1510)
|
2.69 (1.54-4.67)
|
0.00027
|
0.3 (0.06-1.51)
|
0.12
|
4 (n=93)
|
-
|
-
|
-
|
-
|
Mutation burden
|
|
|
|
|
High(n=3510)
|
2.48 (1.7-3.63)
|
1.3e−06
|
0.25 (0.07-0.94)
|
0.027
|
Low(n=3390)
|
1.91 (1.25-2.94)
|
0.0026
|
2.3 (0.91-5.85)
|
0.071
|
Table 2
The cox proportional hazard model of DKK1 and six tumor-infiltrating immune cells in HNSC (TIMER)
|
HNSC
|
coef
|
HR
|
95%CI_l
|
95%CI_u
|
p.value
|
B cell
|
-1.769
|
0.170
|
0.013
|
2.155
|
0.172
|
CD8 T cell
|
-0.892
|
0.410
|
0.064
|
2.638
|
0.348
|
CD4 T cell
|
-2.990
|
0.050
|
0.002
|
1.136
|
0.060
|
Macrophage
|
2.269
|
9.668
|
0.800
|
116.786
|
0.074
|
Neutrophil
|
-0.976
|
0.377
|
0.021
|
6.831
|
0.509
|
Dendritic
|
0.863
|
2.371
|
0.534
|
10.520
|
0.256
|
DKK1
|
0.133
|
1.142
|
1.060
|
1.230
|
0.000
|
Table 3
The correlation between DKK1 expression and tumor lymphocyte infiltration in HNSC (TISIDB)
|
HNSC
|
r
|
p
|
Activated CD8 T cell (Act _CD8)
|
-0.143
|
0.00105
|
Central memory CD8 T cell (Tcm _CD8)
|
0.179
|
4.1e-05
|
Effector memory CD8 T cell (Tem _CD8)
|
-0.137
|
0.00169
|
Activated CD4 T cell (Act _CD4)
|
-0.011
|
0.795
|
Central memory CD4T cell (Tcm _CD4)
|
0.1
|
0.0227
|
Effector memory CD4 T cell (Tem _CD4)
|
-0.078
|
0.0751
|
T follicular helper cell(Tfh)
|
-0.073
|
0.0956
|
Gamma delta T cell (Tgd)
|
0.056
|
0.204
|
Type 1 T helper cell (Th1)
|
-0.074
|
0.0908
|
Type 17 T helper cell (Th17)
|
-0.115
|
0.00877
|
Type 2 T helper cell (Th2)
|
0.067
|
0.126
|
Regulatory T cell (Treg)
|
0.019
|
0.671
|
Activated B cell (Act _B)
|
-0.187
|
1.78e-05
|
Immature B cell (Imm _B)
|
-0.156
|
0.000338
|
Memory B cell (Mem _B)
|
-0.026
|
0.561
|
natural killer cell (NK)
|
-0.066
|
0.129
|
CD56bright natural killer cell (CD56bright)
|
0.104
|
0.017
|
CD56dim natural killer cell (CD56dim)
|
-0.068
|
0.119
|
Myeloid derived suppressor cell (MDSC)
|
-0.1
|
0.0225
|
Natural killer T cell (NKT)
|
0.038
|
0.386
|
Activated dendtritic cell (Act _DC)
|
-0.023
|
0.594
|
Plasmacytoid dendtritic cell (pDC)
|
0.045
|
0.306
|
Immature dendtritic cell (iDC)
|
-0.044
|
0.314
|
Macrophage (Macrophage)
|
-0.011
|
0.807
|
Eosinophi (Eosinophil)
|
-0.144
|
0.001
|
Mast (Mast)
|
-0.099
|
0.0235
|
Monocyte(Monocyte)
|
0.005
|
0.914
|
Neutrophil (Neutrophil)
|
0.052
|
0.239
|
Table 4
Correlation analysis between DKK1 and related genes and markers of immune cells in TIMER.
Description
|
Gene markers
|
HNSC
|
None
|
Purity
|
cor
|
p
|
cor
|
p
|
B cell
|
CD19
|
-0.162
|
1.95e-04
|
-0.171
|
1.42e-04
|
|
CD79A
|
-0.15
|
5.88e-04
|
-0.155
|
5.77e-04
|
CD8+T cell
|
CD8A
|
-0.146
|
8.15e-04
|
-0.138
|
2.09e-03
|
|
CD8B
|
-0.181
|
3.06e−05
|
-0.175
|
9.13e−05
|
T cell (general)
|
CD3D
|
-0.179
|
4.05e−05
|
-0.172
|
1.26e−04
|
|
CD3E
|
-0.148
|
6.94e−04
|
-0.139
|
2.08e−03
|
|
CD2
|
-0.153
|
4.53e−04
|
-0.145
|
1.29e−-03
|
T cell exhaustion
|
CTLA4
|
0.065
|
1.35e−01
|
-0.06
|
1.85e−01
|
|
LAG3
|
-0.08
|
6.69e−02
|
-0.071
|
1.15e−01
|
|
HAVCR2
|
0.001
|
8.7e−01
|
0.031
|
4.98e−01
|
|
GZMB
|
-0.082
|
6.01e−02
|
-0.074
|
1.03e−01
|
|
PDCD1
|
-0.117
|
7.49e−03
|
-0.109
|
1.52e−02
|
Dendritic cell
|
ITGAX
|
-0.01
|
8.19e−01
|
0.004
|
9.30e−01
|
|
NRP1
|
0.185
|
2.15e−05
|
0.208
|
3.38e−06
|
|
CD1C
|
-0.126
|
4.04e−03
|
-0.108
|
1.64e−02
|
|
HLA-DPA1
|
-0.07
|
1.09e−01
|
-0.051
|
2.58e−01
|
|
HLA-DRA
|
-0.07
|
1.13e−01
|
-0.054
|
2.30e−01
|
|
HLA-DQB1
|
-0.044
|
3.21e−01
|
-0.031
|
4.95e−01
|
|
HLA-DPB1
|
-0.107
|
1.46e−02
|
-0.095
|
3.48e−02
|
M1 Macrophage
|
PTGS2
|
0.124
|
4.66e−03
|
0.123
|
6.27e−03
|
|
IRF5
|
-0.128
|
3.37e−03
|
-0.124
|
6.07e−03
|
|
NOS2
|
-0.005
|
9.17e−01
|
0.003
|
9.44e−01
|
M2 Macrophage
|
MS4A4A
|
0.07
|
1.08e−01
|
0.095
|
3.49e−02
|
|
VSIG4
|
0.098
|
2.56e−02
|
0.123
|
6.27e−03
|
|
CD163
|
0.076
|
8.43e−02
|
0.098
|
3.01e−02
|
Monocyte
|
CSF1R
|
-0.015
|
7.35e−01
|
0.006
|
8.89e−01
|
|
CD86
|
0.115
|
8.52e−03
|
0.145
|
1.28e−03
|
Natural killer cell
|
KIR2DS4
|
-0.078
|
7.67e−02
|
-0.075
|
9.77e−02
|
|
KIR3DL3
|
-0.051
|
2.47e−01
|
-0.034
|
4.56e−01
|
|
KIR3DL2
|
-0.054
|
2.18e−01
|
-0.033
|
4.67e−01
|
|
KIR3DL1
|
-0.096
|
2.84e−02
|
-0.105
|
1.95e−02
|
|
KIR2DL4
|
-0.032
|
4.59e−01
|
-0.023
|
6.17e−01
|
|
KIR2DL3
|
-0.066
|
1.31e−01
|
-0.052
|
2.52e−01
|
|
KIR2DL1
|
0.028
|
5.22e−01
|
0.042
|
3.49e−01
|
Neutrophils
|
CCR7
|
-0.138
|
1.59e−03
|
-0.137
|
2.33e−03
|
|
ITGAM
|
0.1
|
2.22e−02
|
-0.08
|
7.48e−02
|
|
CEACAM8
|
−0.016
|
7.14e−01
|
-0.022
|
6.30e−01
|
TAM
|
CCL2
|
0.075
|
8.8e−02
|
0.094
|
3.76e−02
|
|
IL10
|
-0.005
|
9.06e−01
|
-0.012
|
7.88e−01
|
|
CD68
|
0.11
|
1.16e−02
|
0.135
|
2.77e−03
|
Tfh
|
BCL6
|
-0.037
|
4e−01
|
-0.035
|
4.35e−01
|
|
IL21
|
-0.064
|
1.44e−01
|
-0.042
|
3.51e−01
|
Th1
|
TBX21
|
-0.12
|
6.22e−03
|
-0.11
|
1.42e−02
|
|
STAT4
|
0.098
|
2.54e−02
|
0.126
|
5.23e−03
|
|
STAT1
|
0.047
|
2.82e−01
|
0.065
|
1.51e−01
|
|
IFNG
|
-0.11
|
1.19e−02
|
-0.101
|
2.44e−02
|
|
IL13
|
-0.067
|
1.25e−01
|
-0.049
|
2.82e−01
|
Th2
|
GATA3
|
-0.073
|
9.56e−02
|
-0.076
|
9.17e−02
|
|
STAT6
|
0.065
|
1.39e−01
|
0.093
|
4.01e−02
|
|
STAT5A
|
-0.103
|
1.86e−02
|
0.095
|
2.49e−02
|
Th17
|
STAT3
|
-0.063
|
1.49e−01
|
-0.042
|
3.52e−01
|
|
IL17A
|
-0.095
|
3.03e−02
|
-0.081
|
7.30e−02
|
Treg
|
FOXP3
|
-0.098
|
2.54e−02
|
-0.088
|
5.16e−02
|
|
CCR8
|
-0.028
|
5.28e−01
|
-0.008
|
8.52e−01
|
|
STAT5B
|
0.108
|
1.37e−02
|
0.121
|
7.07e−03
|
|
TGFB1
|
0.199
|
4.63e−06
|
0.209
|
1.84e−06
|