Identification of DEGs in GEO dataset
With |log2FC|>1 and P-value <0.05 as the critical values, 115 DEGs were identified in GSE105261, 424 DEGs in GSE22541, and 265 DEGs in GSE85258 using R limma packet(Fig. 1a-c).Detailed information for the three datasets were listed in Table1.In addition,20 overlapping DEGs were found in GSE105261, GSE85258 and GSE22541 datasets (Fig. 1d and Additional file 2: Table S2).
Table1.Data set information for GSE105261, GSE22541 and GSE85258
Data Set
|
Platform
|
|
Sample size
|
Total differentially expressed genes
|
Up-regulated
|
Down-regulated
|
GSE85258
|
GPL570
|
|
Total: 29
Mets: 15
Primary: 14
|
265
|
183
|
82
|
GSE22541
|
GPL570
|
|
Total: 48
Mets: 24
Primary: 24
|
424
|
263
|
161
|
GSE105261
|
GPL10558
|
|
Total: 35
Mets: 26
Primary: 9
|
115
|
79
|
36
|
KEGG and GO Enrichment Analyses of DEGs
To understand the biological functions and pathways of overlapping DEGs, GO and KEGG analyses were performed. Gene Ontology includes biological process (BP), cellular component (CC) and molecular function (MF). The results of biological process showed that DEGs were enriched inextracellular matrix organization. The results of cell component showed that DEGs were enriched incollagen-containing extracellular matrix. The results of molecular function showed that DEGs were enriched inextracellular matrix structural constituent(Fig.2a and Additional file 3: Table S3). Furthermore, theresult of KEGG pathway showed that DEGs were enriched in protein digestion and absorption, aldosterone-regulated sodium reabsorption, proteoglycans in cancer, ECM-receptor interaction, platelet activation, renin-angiotensinsystem, focal adhesion, diabetic cardiomyopathy (Fig. 2b and Additional file 3: Table S3). These results suggest that overlapping DEGs may be associated with tumor invasion and metastasis in ccRCC.
Screening and validation ofthekey gene
In order to further identify the key gene involved in ccRCC metastasis, survival analysis of overlapping genes was performed.The result showed that the high expression of KDELR3(Logrank p = 5.7e-07) and PDGFRL(Logrank p = 8.3e-09)in ccRCC patients suggested poor prognosis. Low expression of EMX2 (Logrank p= 5.5e-09) and HSD11B2 (Logrank p =4.2e-05) were associated with poorprognosis(P < 0.001)(Fig. 3a-3t).Meanwhile,we also used GEPIA2 database to detect mRNA expression levels of prognostic genesand found that the expression level of HSD11B2 in the primary group was higher than that in the metastatic group. (Fig.4a-d).Therefore, these results suggest that HSD11B2 may be the key gene in the process of ccRCC metastasis.
Low expression of HSD11B2 in ccRCC and qRT-PCR validation
Firstly, we found that the expression level of HSD11B2 mRNA in ccRCC tissues was significantly lower than that in normal tissues in TCGA-KIRC dataset(p = 3.68e-38, Fig. 5a) Then, 15 pairs of ccRCC tissues and adjacent normal kidney tissues were used to verify that HSD11B2 was significantly down-regulated in ccRCC by qRT-PCR(p<0.0001)(Fig. 5b and Additional file 4: Table S4).In addition,the expression of HSD11B2 protein in normal kidney tissue and ccRCC tissue was verified using Human Protein Atlas (HPA) database.The results showed that the expression of HSD11B2 protein in tumor tissue was significantly lower than that in normal tissue(Fig. 5c).
Correlation between HSD11B2 expression level and clinical features
Through the analysis of the correlation between HSD11B2 expression level and clinical features, it was further proved that HSD11B2 is the key gene in the metastasis process of ccRCC.The association between HSD11B2 expression and clinical features in ccRCC patients is summarized in Table 2.In addition, TCGA-KIRC dataset analysis revealed significant differences in HSD11B2 among different clinical subgroups.The low expression of HSD11B2 was significantly correlated with tumor stage(p=6.08e-06), T stage(p=9.43e-05), N stage(p=0.0001) and M stage(p=2.24e-05)(Figure 6a,b, c, d). These results suggest that HSD11B2 expression can be used as an indicator to assess the risk of metastatic renal clear cell carcinoma.
HSD11B2 expression may be an independent prognostic factor
Univariate Cox analysis showed that HSD11B2 expression was significantly correlated with OS(p=1.05e-05, HR=0.801)(Fig.7a andTable 3).Multivariate Cox analysis showed that HSD11B2 expression could be an independent prognostic factor for OS of ccRCC patients(p=0.009, HR=0.873)(Fig.7b and Table 3).
Table2.Univariate and multivariate analyses of HSD11B2 and clinical factors associated with OS
Parameter
|
Univariate analysis
|
Multivariate analysis
|
HR
|
95% CI
|
p
|
HR
|
95% CI
|
P
|
Overall survival
|
Age
|
1.03
|
1.02–1.04
|
8.50E-06
|
1.03
|
1.019–1.048
|
4.95E-06
|
Gender
|
1.06
|
0.77–1.45
|
0.706
|
1.1
|
0.799–1.522
|
0.548
|
T
|
1.92
|
1.63–2.27
|
8.61E-15
|
0.945
|
0.627–1.421
|
0.785
|
M
|
4.39
|
3.22–5.99
|
1.04E-20
|
1.54
|
0.796–2.968
|
0.200
|
Tumor Stage
|
1.89
|
1.66–2.16
|
3.98E-21
|
1.67
|
1.068–2.620
|
0.0245
|
HSD11B2
|
0.801
|
0.72–0.88
|
1.05E-05
|
0.873
|
0.788–0.967
|
0.00955
|
Establishment and validation of a prognostic Nomogram
Based on univariate and multivariate COX regression analysis, HSD11B2 expression, age and tumor stage can be used as independent prognostic factors to construct a Nomogram to predict 1-year, 3-year and 5-year survival(Fig.8a). The total scoreforall patients provided the estimated 1-year, 3-year and 5-year survival time. The results show that the actual OS of 1-year, 3-year and 5-yearin accord with the nomogram-predicted OS (Fig.8b).
HSD11B2 enriches multiple cancer pathways in ccRCC
GSEA analysis showed that the low expression of HSD11B2 was enriched in cancer-related pathways such as nod like receptor signaling pathway, P53 signaling pathway and cytoplasmic DNA sensing pathway. (Fig.9a-k and Additional file 5: Table S5). These results provide further evidence that low expression of HSD11B2 may contribute to the development of cancer.
Correlation between HSD11B2 expression and immune cells
In order to investigate the relationship between HSD11B2 expression and tumor immune invasion, we used the TISIDB database to analyze the correlation between HSD11B2 expression and tumor infiltrating lymphocytes (TILs) and three kinds of immunomodulators. For TILs, HSD11B2 expression was negatively correlated with most of the tumor infiltrating lymphocytes of CCRCC, such as ACT_ DC_ abundance (Spearman: r = -0.485, P < 2.2e-16), MDSC_abundance (Spearman: r = -0.469, P < 2.2e-16), Act_CD4_abundance (Spearman: r = -0.459, P < 2.2e-16), and Macrophage_abundance(Spearman: r= -0.438, P < 2.2e-16) (Figure10A). Immunomodulators include immunoinhibitors, immunostimulators and major histocompatibility complex (MHC) molecules.For immunoinhibitors, HSD11B2 expression was negatively correlated with CSF1R (Spearman: r = -0.411, P < 2.2e-16), CTLA4 (Spearman: r = -0.366, P < 2.2e-16), LGALS9(Spearman: r = -0.346, P = 2.11e-16), and BTLA (Spearman: r = -0.341, P = 6.76e-16) in ccRCC(Figure 7B).Forimmunostimulators, HSD11B2 expression was negatively correlated with TNFSF13B (Spearman: r = -0.472, P < 2.2e-16), CD86(Spearman: r = -0.454, P < 2.2e-16), TNFSF14 (Spearman: r = -0.403, P < 2.2e-16), and CD80 (Spearman: r = -0.395, P < 2.2e-16)(Figure 8A). For MHC molecules, HSD11B2 expression was negative correlated with HLA-DOB (Spearman: r = -0.305, P = 8.18e-13), HLA-DMB (Spearman: r = -0.292, P = 7.83e-12), HLA-DRA (Spearman: r = -0.289, P = 1.16e-11), and HLA-DBP1 (Spearman: r = -0.263, P = 7.64e-10) (Figure 10B).These results suggest that the expression of HSD11B2 may regulate tumor immune invasion by interacting with these immune cells.
Table 3.Association of HSD11B2 expression and the clinicopathological characters in ccRCC
Clinicopathological
|
Taylor cohort
|
N
|
mean±SD
|
P value
|
Age(years)
|
|
|
|
<60
|
259
|
4.82±1.67
|
0.5644
|
≥60
|
262
|
4.83±1.78
|
Gende
|
|
|
|
Female
|
182
|
5.08±1.69
|
0.01444
|
Male
|
339
|
4.69±1.73
|
Tumor stage
|
|
|
Stage I-II
|
318
|
5.09±1.78
|
5.755e-06
|
Stage III-IV
|
200
|
4.42±1.55
|
T
|
|
|
|
T1-T2
|
335
|
5.06±1.76
|
3.013e-05
|
T3-T4
|
186
|
4.42±1.59
|
N
|
|
|
|
N0
|
233
|
4.89±1.78
|
0.00130
|
N1
|
16
|
3.05±1.85
|
M
|
|
|
|
M0
|
423
|
4.99±1.69
|
5.238e-05
|
M1
|
79
|
4.11±1.70
|