Screening out of Two Prognosis-Related Modules
After weeding out 21 outlier samples, 509 ccRCCs were included for WGCNA (Fig. 2). The beta (β) = 5 (scale free R2 = 0.84) was further set as soft-thresholding for further adjacencies calculation (Fig. 3). Then we split IRGs into gene modules. High-related modules were further merged for avoiding over-fitting. Eventually, as shown in Fig. 4A, totally 18 modules were identified. In this method, no-significant genes were stored in grey module, which were excluded from subsequent analysis. These genes were regarded as biomarkers with no-significance. Furthermore, brown module was associated with OS time (P = 0.001, R2 = -0.14, Fig. 4B) and OS (P = 0.35, R2 = 0.35, Fig. 4B). Turquoise module was associated with not only OS time (P = 2e-05, R2 = -0.19, Fig. 4B) but also OS (P = 5e-05, R2 = 0.18, Fig. 4B). The MS (of OS time) of turquoise module was higher than MSs of any other modules meanwhile the MS (of OS) of brown module was the highest compared to others (Fig. 4D). Moreover, the associations between MM and GS for OS time (cor = 0.4, P = 2.3e-27)/ OS (cor = 0.66, P = 1.3e-120) in brown module were also significant (Fig. 4C). The same trends existed in turquoise module (OS time: cor = 0.65, P = 5.1e-142; OS: cor = 0.42, P = 1.4e-78). Both the two modules (brown, turquoise) reached the standards previous mentioned, which were considered as prognosis-related modules in the present study. We also created a network heatmap in the present study (Fig. 5A). As shown in Fig. 5B, the classical MDS plot concluded that the 18 modules were independent of each other.
Function Exploring of Genes in Prognosis-Related Modules
GO and KEGG pathway enrichment analysis were conducted to explore the potential of genes in prognosis-related modules. As shown in Fig. 6A, GO analysis indicated that genes in brown module were involved in mitotic nuclear division, nuclear division, organelle fission, chromosome segregation, mitotic sister chromatid segregation, sister chromatid segregation, microtubule cytoskeleton organization involved in mitosis, regulation of mitotic nuclear division, regulation of nuclear division and nulcear chromosome segregation (Fig. 6A). Moreover, genes in turquoise module were mainly enriched in epithelial cilium movement, cilium assembly, cilium organization and cilium movement (Fig. 6C). The KEGG enrichment analysis indicated that genes in brown module significantly enriched in cell cycle, progesterone-mediated oocyte maturation, p53 signaling pathway and oocyte meiosis (Fig. 6B) meanwhile genes in turquoise module were mainly enriched in herpes simplex virus1 infection (Fig. 6D).
Two Genes Were Regarded as Novel Immune-Related Prognostic Biomarkers
10 genes from brown module (Fig. 7A) and 80 genes from turquoise module (Fig. 7B) were selected by using the cut-off criterion of |cor.geneModuleMembership| > 0.8 and |cor.geneTraitSignificance| > 0.2. Finally, two genes including GNRH1 (gonadotropin releasing hormone 1) and LTB4R (leukotriene B4 receptor) overlapping in hub genes in WGCNA and IRGs were screened out (Fig. 7C).
Internal Validation of the Two Immune-Related Prognostic Biomarkers
After screening out four hub genes though comprehensive bioinformatics analysis, we validated these potential prognostic biomarkers. Firstly, based on the GEPIA webtool, Higher expression of GNRH1 was significantly related to worse OS (Hazard ratio [HR] = 1.9, P = 3.4e-05, Fig. 8A). Moreover, we concluded that patients in LTB4R high-expression group obviously occupied worse OS (HR = 1.9, P = 3.2e-05, Fig. 8D), accurately. Then the mRNA expression of prognostic biomarkers between tumor tissues and normal tissues was compared, the results indicated that both the genes including GNRH1 (P < 0.05, Fig. 8B), and LTB4R (P < 0.05, Fig. 8E) were significantly higher expressed in ccRCC samples compared to normal samples. Moreover, the results also suggested that high expression of GNRH1 (F = 2.63, P = 0.0497; Fig. 8C) and LTB4R (F = 3.16, P = 0.0243; Fig. 8f) were significantly related to higher tumor stage.
External Validation of the Two Immune-Related Prognostic Biomarkers
Expression levels of prognostic biomarkers between normal samples and ccRCC samples were compared based on Oncomine database. We indicated the same conclusion that GNRH1 (P = 0.049, Fig. 9A) and LTB4R (P = 0.003, Fig. 9D) had higher expression in ccRCC tissues than in normal tissues. By using HPA database, we investigated the translational level expression of immune-related prognostic biomarkers. GNRH1 obtained medium staining in ccRCC samples compared to normal samples (not detected) (Fig. 10A). Unfortunately, there was no significant difference between ccRCC and normal tissues on the level of translation (Fig. 10B). The prognostic values of immune-related prognostic biomarkers were further validated after verification of expressions. As the results suggested, GNRH1 expression was significantly associated with not only OS time (r = -0.19, P < 0.001, Fig. 11A) but also OS status (r = 0.21, P < 0.001, Fig. 11B). Moreover, LTB4R expression was significantly related to OS time (r = -0.17, P < 0.001, Fig. 11C) and OS status (r = 0.26, P < 0.001, Fig. 11D). Also, the survival analysis demonstrated that there was a trend that higher expression of GNRH1 was significantly correlated to worse OS (GSE29609), CSS (GSE29609), PFS (GSE29609), PFS (E-MTAB-3267); suggested by Fig. 12A (P = 0.310), Fig. 12B (P = 0.450), Fig. 12C (P = 0.660), Fig. 12D (P = 0.063), respectively. Meanwhile ccRCC patients with higher LTB4R expression was related to short OS time, accurately (P = 0.034, GSE29609, Fig. 12E). The trend that higher expression of LTB4R was associated with worse CSS was verified by using GSE29609 (P = 0.660, Fig. 12F). Furthermore, ccRCC patients with higher expression of LTB4R was significantly correlated to worse PFS, as Fig. 12G (P = 0.035, GSE29609) and Fig. 12H (P = 0.0028, E-MATB-3267) suggested.
In addition, we also conducted ROC analysis. The result demonstrated that GNRH1 could significantly distinguish ccRCC samples from normal samples, suggested by Fig. 13A (AUC = 0.770). Interestingly, LTB4R showed strong potential for ccRCC diagnosis (AUC = 0.829, Fig. 13B). All the results above indicated that the two immune-related prognostic biomarkers we screened were credible.
Experimental validation
We further validated the expression of GNRH1 and LTB4R in ccRCC cell line by qRT-PCR and western blotting. qRT-PCR analysis showed that GNRH1 and LTB4R mRNA levels were significantly higher in ccRCC tissues than normal kidney tissues (Fig. 14A-B). In other hand, compared with the normal renal epithelial cell line HK2, expression of GNRH1 and LTB4R is significantly increased in the ccRCC cell line Caki1 (Fig. 14C-D).
Identification of Immune-Related Prognostic Biomarker Related Pathways
The results of GSEA indicated that GNRH1 was significantly enriched in one KEGG signaling pathway including olfactory transduction (nominal P = 0.014, n = 368, FDR = 5.164%, Fig. 15C). Meanwhile, we found that LTB4R was significantly related to adipogenesis (nominal P = 0.035, n = 191, FDR = 22.685%).
Association between Immune-Related Prognostic Biomarker Expression and Immune Infiltration Level in ccRCC
Immune infiltration was reported to be associated with survival and progression of cancers. Thus, the association between prognostic biomarker and immune infiltration level was obtained by applying ssGSEA. GNRH1 expression was significantly associated with several immune cell types, such as central memory CD8 T cell, CD56dim natural killer cell, and immature dendritic cell (Fig. 15A). Moreover, LTB4R was significantly related to activated CD4 T cell, CD56dim natural killer cell, and MDSC. These results indicated that expression of the two prognostic biomarkers might affect the immune infiltration levels of some tumor-infiltrating immune cell types.