Potential prognostic biomarkers related to immunity in clear cell renal cell carcinoma using bioinformatic strategy.

The clear cell renal cell carcinoma (ccRCC) is the main pathological subtype of renal cell carcinoma. Immune system evasion, one hallmark of cancer, contributes to cancer cells in escaping from the attack of immune cells. In order to identify potential prognostic biomarkers in ccRCC patients and immune cells fraction, we collected and downloaded profiles from The Cancer Genome Atlas (TCGA) database and Gene Expression Omnibus (GEO) database. We obtained 2 modules significantly associated with tumor stage and immune cells; functional enrichment analysis showed that genes in the module 'yellow' were significantly enriched in proteins targeting to membrane and ribosome, as well as the oxidative phosphorylation pathway, while genes in the module 'green' mainly participate in molecular functions associated with immunity like activation of T cells. Four LncRNAs (LINC00472, AL590094.1, AL365203.3, and AC147651.3) and RPL27A and RPL22L1 in the module 'yellow' and two lncRNAs (LINC00426 and AC129507.2) and five protein-coding genes (CSF1, NOD2, ITGAE, CD7, and PDCD1) in the module 'green' represented independent prognostic values in patients with ccRCC. Expression of LINC0042, NOD2, CD7, and PDCD1 were significantly correlated with ratio of immune cells (like T cells CD8 and resting mast cells). LINC00426, with significant correlation with immune cell fraction, shows potential prognostic value in ccRCC patients. Our findings provide a strategy in exploring biomarkers with prognostic significance and significant association with the fraction of immune cells.

Introduction

Renal cell carcinoma (RCC) is the most frequent malignancy in the kidney, accounting for 2–3% of all cancers [1]. Clear cell renal cell carcinoma (ccRCC) is the most common pathological subtype of RCC, originating from the proximal uriniferous tubules and often presenting an aggressive phenotype [2]. ccRCC has high morbidity and mortality, as well as a poor prognosis. Since ccRCC is insensitive to chemotherapy and radiotherapy [3], and the clinicopathological risk factors cannot fully identify ccRCC patients, it is urgent and necessary to mine the potential biomarkers that could help to distinguish patients with ccRCC.

Long non-coding RNAs (lncRNAs) are RNA transcripts longer than 200 nucleotides in length, which cannot be translated into proteins. lncRNAs play an essential role in tumor biology, such as cell proliferation, metabolism, differentiation, and angiogenesis [4]. Previous studies have shown that lncRNA MALAT1 modified the progression of ccRCC by regulating miR-194-5p/ACVR2B signaling [5]. lncRNA MSC-AS1 could activate the Wnt-β-catenin signaling pathway to modulate cell proliferation and migration in ccRCC via miR-3924/WNT5A [6]. A novel lncRNA, LINC00997, contributed to metastasis by regulating S100A11 in ccRCC [7]. Besides, lncRNAs were reported as potential prognostic molecular biomarkers for RRC. A panel of four-lncRNA signature was identified as a potential biomarker for predicting survival in ccRCC [8]. Four lncRNAs were potential therapeutic targets and prognostic biomarkers of ccRCC [9]. A recent publication reported a novel prognostic lncRNA with association of transition of epithelial-mesenchymal in ccRCC [10]. Interactions between cancer cells and immune cells play critical roles in patients in response to cancers. Cancer cells might escape from immune cell attack or influence the dysfunction and apoptosis of immune cells by expressing immune inhibitory molecules [11]. Analysis of lncRNAs associated with immune cells in patients with ccRCC, however, is limited. Therefore, more potential and valuable lncRNA biomarkers, especially immune-related lncRNAs, need to be identified to improve the prognosis.

Weighted gene co-expression network analysis (WGCNA) describes the correlation patterns among genes and clinically relevant traits across microarray samples. Researchers can find modules of highly correlated genes, associate modules with external sample traits, and count module membership measures based on WGCNA [12]. WGCNA has been increasingly used in disease research. Key modules and genes of subtypes of non-small-cell lung cancer were identified by WGCNA [13]. Key pathways and genes were identified in the dynamic progression of HCC based on WGCNA [14]. ILF3-AS1 could regulate PTBP1 by sponging miR-29a in gastric cancer according to WGCNA co-expression network analysis [15]. Hub-methylated differentially expressed genes in patients with gestational diabetes mellitus were identified by multi-omic WGCNA based on epigenome-wide and transcriptome-wide profiling [16].

In this study, our study aims to identify potential prognostic biomarkers related to immunity in distinguishing ccRCC patients at early and advanced stages. We identified modules related to ccRCC tumor-stage and immune infiltration using WGCNA; functional enrichment and protein–protein interactions analyses were conducted to explore relevant signaling pathways and gene clusters; univariate and multivariate Cox regression analyses were conducted to identify potential prognostic biomarkers for ccRCC; correlation among immune cells and candidate genes was conducted to analyze genes with significant association with immune infiltration.

Materials and methods

Data collection, processing, and quality control

TCGA-ccRCC (containing 530 RCC samples and 72 control tissue samples) was downloaded from the TCGA database (https://portal.gdc.cancer.gov/) for analysis. To investigate genes associated with tumor stage, gender, and tumor immune microenvironment (TIME) [17], 530 tumor samples were taken for analysis. Due to the lack of cancer stage information, 3 data were excluded, leaving 527 data for further analysis.

Immune cell distribution analysis

To explore the immune cell distributions, an analytical tool CIBERSORTx (https://cibersortx.stanford.edu/) was used to examine the fractions of 22 immune cells in tumor samples. The samples were separated into two groups: early (I/II) and advanced (III/IV) ccRCC according to the clinical information. The abundance of immune cells was depicted using a heatmap.

Differential gene analysis

Patients were divided into two groups: early (I/II) and advanced (III/IV) ccRCC as per tumor stages. All data from TCGA-ccRCC were normalized using the limma package. The threshold to identify differentially expressed genes (DEGs) was set to |fold change| >>1.2 P < 0.05 to ensure sufficient genes for downstream analysis. Volcano plots and heatmaps were used to visualize the DEGs using the ggplot2 package and pheatmap software, respectively.

WGCNA network construction

The gene expression matrix and clinical information (tumor stage (TS), gender (G), and immune cells fraction (ICF)) matrix network were constructed by setting the ‘minModuleSize’ parameter to 80. The ‘mergeCloseModules’ function was used to merge modules with the correlation of 0.75 and other default parameters were added. The correlation between clinical traits and sample expression was calculated with the ‘cor’ function; the P-value was calculated with the ‘corPvalueStudent’ function. Genes with significant relationships (P < 0.05) with clinical traits were identified for further analysis.

Functional enrichment analysis and protein–protein interaction (PPI)

DEGs associated with immunity and tumor stages in identified modules using WGCNA were used for functional enrichment analysis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using clusterProfiler [18]. GO terms can be classified into 3 categories: biological processes (BP), cellular components (CC), and molecular functions (MF). GO terms (top ten of each category) and KEGG pathways were visualized with bar plots and bubble plots, respectively. PPI was analyzed using an online tool STRING (https://string-db.org/), which provides functional protein association networks. PPI was visualized using Cytoscape software to predict gene clusters.

Cox regression analysis and Kaplan-Meier (KM) survival analysis

Univariate Cox regression analysis was performed to explore potential prognosis-associated genes. Genes with clinical prognostic value were visualized with forest trees for further analysis. Then, multivariate Cox regression analysis was conducted using genes identified to evaluate the independent prognostic value of genes as biomarkers. The Cox regression analysis was visualized by forest plot.

To compare the expression level of genes between patients with ccRCC at early and advanced stages, we calculated the difference of two groups with the T-test method; the expression levels were depicted with box plots. The patients with ccRCC were divided into two groups according to the expression of genes, KM survival plots are depicted to explore the association between expression levels and survival probability.

Correlation analysis of genes and immune cells

To explore the impact of genes on immunity, we analyzed the association between gene expression and fraction of immune cells. In detail, we conducted a matrix containing the expression of genes with significant prognostic value and proportion of immune cells and calculated the correlation scores among them with the Pearson method; significant value (P-value) was calculated and set to 0.01 for significance.

Validation of candidate genes

To validate the expression of genes obtained in this analysis, we collected independent profiles related to different stages of ccRCC from Gene Expression Omnibus (GEO) database. GSE150404, with 60 tumor samples (30 in early stage, 30 in advanced stage), was downloaded and normalized. We annotated the genes according to the annotation profile, GPL17692[HuGene-2_1-st] Affymetrix Human Gene 2.1 ST Array [transcript (gene) version]. Boxplots were depicted with method described above.

Results

Immune cell distribution

The fraction of immune cells between two tumor stages was analyzed to explore the variety of immune infiltration; the analysis revealed differences in multiple immune cells (T cells CD4 memory resting, Macrophages M2 and T cells CD8) among different samples (Figure 1). Representative immune cells (B cells naïve, T cells CD8, T cells CD4 memory resting, NK cells resting, Monocytes, Macrophages M0, Macrophages M1, Macrophages M2, and Mast cells resting) were selected for immune-related gene analysis.

Figure 1.

Immune cell distribution analysis between two groups

DEGs at different stages of ccRCC

DEGs between two tumor stages of ccRCC were analyzed; a total of 2902 DEGs were detected, including 2556 protein-coding genes (PC) and 218 lncRNAs according to the threshold. The DEGs were visualized using volcano plots (Figure 2a) and heatmaps (Figure 2b).

Figure 2.

Visualization of differentially expressed genes (DEGs) at different stages of ccRCC

Weighted co-expression network construction

Topology trees were constructed based on expression levels and the outlier samples were removed (Figure 3a). Then, we matched the samples with clinical information provided (Figure 3b).

Figure 3.

Sample clustering dendrogram of TCGA-ccRCC

After removing the outlier samples, 448 samples (270 early stages, 178 advanced stages) with 2092 genes were finally selected to build the network. The soft-threshold power β = 7 was selected to construct a network based on the scale-free topology (Figure 4(a, b)), and 7 modules were obtained (Figure 4(d, e), Table 1). Module ‘green’ with 554 genes indicated a significant correlation with multiple immune cells including T cells CD8, Monocytes, and Mast cells resting (Figure 4e). Given that the module ‘yellow’ correlated with survival, cancer stage, and multiple immune cell distribution (Figure 4e), this module (containing 343 genes) was also selected for subsequent analysis.

Figure 4.

Weighted gene co-expression network analysis (WGCNA) of TCGA-ccRCC

Table 1.

Modules identified using weighted gene co-expression network analysis

Module colors	Number of genes
black	188
blue	454
brown	353
green	554
grey	204
pink	150
turquoise	656
yellow	343

Functional enrichment analysis and PPI network construction

Interactions of proteins, forming a machine, play a crucial role in regulating ccRCC. Genes in the most significant modules related to traits (module ‘yellow’) were selected for functional enrichment analysis. These genes were enriched in 59 BP, 6 CC, and 2 MF. These genes were significantly enriched in BP terms related to membranes, such as the establishment of protein localization to the membrane, co-translational protein targeting to membrane, and protein targeting to membrane; CC terms related to the ribosome (Figure 5a). Genes in the module ‘green’ were enriched in 935 terms including 826 terms in BP, 45 in CC, and 64 in MF. Multiple molecular functions associated with immunity such as T cell activation, immunological synapse, immune receptor activity, and regulation of cytokine were enriched (Figure 5c). KEGG enrichment analysis showed that genes in module ‘yellow’ were enriched in pathways like prion disease, oxidative phosphorylation, and retrograde endocannabinoid signaling (Figure 5b). KEGG analysis of genes in the module ‘green’ indicated the enrichment of pathways associated with immunity such as the interaction of cytokine-cytokine receptor, signaling pathway of T cell, and primary immunodeficiency (Figure 5d).

Figure 5.

Functional enrichment analysis of genes in yellow modules

PPI network was constructed using the STRING database (Figure 6a for module ‘yellow’, 6c for ‘green’). Genes in the gene cluster with the highest score predicted by Cytoscape based on interactions of genes in the module ‘yellow’ were selected for further analysis (Figure 6b). Given the complex interactions among genes in the module ‘green’, we selected the top 2 gene clusters for further analysis (Figure 6d).

Figure 6.

Protein-protein interaction (PPI) and gene cluster prediction

Univariate and multivariate Cox regression analysis

For the ‘yellow’ module, univariate regression analysis identified 271 genes of prognostic significance out of 343 genes, including 21 lncRNAs (as in Figure 7a) and 17 genes in gene clusters (Figure 7b). All three tests (likelihood ratio test, Wald test, and Logrank test) suggested reliable regression models (P < 0.05). Multivariate Cox regression analysis identified four LncRNAs (LINC00472, AL590094.1, AL365203.3, and AC147651.3) (Figure 7c) and RPL27A and RPL22L1 (Figure 7d) in the gene cluster with independent prognostic value in distinguishing patients with ccRCC in early stages from advanced stages. We obtained 19 lncRNAs and 26 protein-coding genes in top 2 gene clusters with prognostic significance (Figure 7(e, f)); Multivariate Cox regression analysis revealed that two lncRNAs (LINC00426 and AC129507.2) and five protein-coding genes (CSF1, NOD2, ITGAE, CD7, and PDCD1) had potential independent prognostic value (Figure 7(g, h)).

Figure 7.

Prognostic analysis visualized with forest plots

Expression level and KM analysis

The expression of four LncRNAs (LINC00472, AL590094.1, AL365203.3, and AC147651.3) and RPL27A and RPL22L1 was visualized with box plots (Figure 8a). LINC00472 was significantly decreased in the advanced stages of ccRCC compared with that in early stages; AL590094.1, AL365203.3, AC147651.3, RPL27A, and RPL22L1 were significantly increased in advanced stages compared with that in early stages (Figure 8a). The KM plots showed that high expression of LINC00472 had a greater survival probability, while the other five genes showed the opposite trend (Figure 8b).

Figure 8.

Expression level and Kaplan-Meier (KM) survival of genes in module ‘yellow’

For module ‘green’, gene expression was significantly stimulated in patients with ccRCC in advanced stages, compared with those in early stages, except lncRNA, AC129507.2, showed a decrease in patients with ccRCC in advanced stages (Figure 9a). LnRNA, LINC00426, did not present a significant prognostic value in the KM method; AC129507.2 expression indicated a positive correlation with survival probability, while the other five protein-coding genes represented a negative association (Figure 9b).

Figure 9.

Expression level and Kaplan-Meier (KM) survival of genes in module ‘green’

Immune association analysis

Six genes (4 lncRNAs and 2 protein-coding genes) showed weak or no relation to the fraction of immune cells; LINC00472 indicated a negative correlation with the expression of the other threelncRNAs (Figure 10a). We observed an obvious association of genes in the module ‘green’, the immune cells proportion; expression of LINC00426, NOD2, CD7, and PDCD1 indicated a positive correlation with a fraction of T cells CD8; the increase of LINC00426 decreased the accumulation of monocytes and mast cells (resting); negative correlation of NOD2 and mast cells (resting) fraction was observed; PDCD1 expression limited the accumulation of NK cells (resting), monocytes, and mast cells (resting) (Figure 10b).

Figure 10.

Correlation of genes and immune cells

Validation of crucial DEGs

The expression of seven protein coding genes (RPL27A, RPL22L1, CSF1, NOD2, ITGAE, CD7, and PDCD1) between the early and advanced stage of ccRCC were validated with GSE150404. We failed to examine the expression of RPL22L1 due to the absence of expression in the matrix. Five genes (CSF1, NOD2, ITGAE, CD7, and PDCD1) showed significantly increased expression level in patients in advanced stage compared with those in early stage; expression of RRL27A was increased in advanced stage though not significant.

Discussion

WGCNA was used to find the most significant module related to ccRCC tumor stage and immunity, and we found that genes in the ‘yellow’ and ‘green’ modules were significantly related to tumor stage and immunity (Figure 4). Previous publications reported the histological tumor types in TCGA-ccRCC project [19]; similarly, we observed multiple outliers in TCGA-ccRCC (Figure 3a). GO enrichment analysis showed that genes in the ‘yellow’ module were significantly enriched in CC terms related to the ribosome (Figure 5a). The previous research showed that knockdown of ribosomal protein S15A inhibited human kidney cancer cell growth in vitro and in vivo [20]. Ribosomal s6 protein kinase 4 was a prognostic factor for RCC [21]. Genes in module ‘green’ were enriched in immune-related terms like T cell activation and cytokine activity, which have been widely reported in the process of RCC [22-25]. KEGG enrichment analysis showed that identified genes in module ‘yellow’ were enriched in pathways like oxidative phosphorylation (Figure 5b). Papillary renal cell carcinomas rewired glutathione metabolism and were deficient in both anabolic glucose synthesis and oxidative phosphorylation [26]. System biology analysis showed that eight proteins were mainly distributed in oxidative phosphorylation in RCC with bone metastasis via mitochondrial dysfunction [27]. Some co-expressed genes in RCC were also reported to be enriched in oxidative phosphorylation [28,29], which is consistent with our research. Cytokine-cytokine receptor interaction and chemokine signaling pathway were significantly enriched by genes in the ‘green’ module; cytokines play a critical role in regulating the RCC microenvironment and mediating immune cell infiltration [30]. C-X-C Motif Chemokine Ligand 10 (CXCL10), a protein-coding gene engaged in chemokine activity, has been revealed to be associated with disease-free survival (DFS) in patients remedied with adjuvant VEGFR TKIs in RCC [31].

The Cox regression analysis showed that LINC00472 had a independent diagnostic value in distinguishing patients with ccRCC in early stages from advanced stages (Figure 7(a, c)). The KM plot showed that ccRCC patients with high expression of LINC00472 had a greater survival probability (Figure 8b). LINC00472 acted as a tumor suppressor in NSCLC through the KLLN-Mediated p53-signaling pathway via microRNA-149-3p and microRNA-4270 [32]. Downregulation of LINC00472 promotes osteosarcoma tumorigenesis by reducing FOXO1 expressions via miR-300 [33]. LINC00472 regulated cell stiffness and inhibited the migration and invasion of lung adenocarcinoma by binding to YBX1 [34]. LINC00472 suppressed proliferation and promoted apoptosis through elevating PDCD4 expression by sponging miR-196a in colorectal cancer [35]. It was reported that LINC00472 may serve as a potential diagnostic marker for diabetic kidney disease [36]. LINC00472 was identified as an RNA transcriptional marker associated with ccRCC prognosis by a competing endogenous RNA network analysis [37]. However, that study differed from our study in terms of the study object. Our research aimed to identify biomarkers to distinguish patients with ccRCC at early and advanced stages, while the previous research aimed to distinguish patients with ccRCC from healthy individuals. We found that RPL22L1 and RPL27A had the potential in distinguishing ccRCC patients at early and advanced stages (Figure 7). The KM plots showed that patients with lower expression of RPL22L1 and RPL27A had higher survival probability (Figure 8b). The previous research showed that ribosomal subunit protein RPL27A is a target of miR-595 and may contribute to the myelodysplastic phenotype through ribosomal dysgenesis [38]. Heavy-ion radiation-induced DNA damage mediated apoptosis via RPL27A-RPL5-MDM2-p53/E2F1 signaling pathway in mouse spermatogonia [39].

Recent works revealed that LINC00426 expression represented significantly prognostic value in patients with lung cancer and hepatocellular carcinoma [40,41]. Zhang, et al. found the potential of immune-related lncRNAs (including LINC00426) in predicting immune checkpoint blockade [41]. In this study, we observed a significant correlation between LINC00426 expression and fraction of T cells CD8 and mast cells (resting) (Figure 10b). LINC00426 showed independent prognostic significance in patients with ccRCC (Figure 7(e, g)). However, we did not obtain a prognostic value for LINC00426 consistent with the KM approach (Figure 8d). Kaplan–Meier curves are powerful as categorical predictors; however, they might represent instability in analyzing quantitative predictors like gene expression [42]. More expression data are required to assess the prognostic significance of LINC00426. Expression of NOD2, CD7, and PDCD1 was significantly increased in advanced stage and indicated a significant positive correlation with T cells CD8 fraction (Figures 9a, 10b, 11). NOD2 was reported in regulating innate and adaptive immunity in patients with Crohn’s disease [43]. Researchers found that the NOD2 expression showed potential as a biomarker for kidney cancer patients [44]. CD7 plays a crucial role in mediating the apoptosis of galectin-3-induced type II T-cells and was considered as a critical target for acute lymphoblastic leukemia [45,46]. Programmed cell death 1 (PDCD1, also called PD1), a critical inhibitory molecule mainly expressed in pro-B-cells, participates in T cell receptor signaling and innate immune system; one ligand of PDCD1, programmed death-ligand-1 (PD-L1), was utilized as a critical target in treating patients with cancers [47,48]. In our analysis, we observed the differential expression of NOD2, CD7, and PDCD1 between two tumor stages of ccRCC and the potential of NOD2, CD7, and PDCD1 as independent prognostic biomarkers in patients with ccRCC.

Validation of candidate genes

In conclusion, based on systematic analysis of the correlation of genes expression, tumor stage, and immune cells fraction, we observed 6 lncRNAs and 7 protein coding genes with independent prognostic significance in patients with ccRCC; one lncRNA (LINC00426) and three genes (NOD2, CD7, and PDCD1) indicated significant correlation with immune cells fraction, especially T cells CD8.

Conclusions

Six lncRNAs (LINC00472, AL590094.1, AL365203.3, AC147651.3, AC129507.2, and LINC00426) and 7 genes (RPL27A, RPL22L1, CSF1, NOD2, ITGAE, CD7, and PDCD1) show expression differences in patients with ccRCC at early stages compared to those at advanced stages, representing independent prognostic significance; LINC00426 is one potential biomarker with prognostic value and indicates a significant correlation with immune cell fraction.

Click here for additional data file.