Identification of Survival Risk and Immune-Related Characteristics of Kidney Renal Clear Cell Carcinoma

Wu, Xiaobin; Liang, Yonghui; Chen, Xian; Long, Xiangyang; Xu, Wujun; Liu, Li; Wang, Binhui; Zou, Xiong

doi:https://doi.org/10.1155/2022/6149369

Journal of Immunology Research

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusion Data Availability Ethical Approval Conflicts of Interest Authors’ Contributions Acknowledgments Supplementary Materials References Copyright Related Articles

Research Article | Open Access

Volume 2022 | Article ID 6149369 | https://doi.org/10.1155/2022/6149369

Identification of Survival Risk and Immune-Related Characteristics of Kidney Renal Clear Cell Carcinoma

Xiaobin Wu,¹Yonghui Liang,²Xian Chen,¹Xiangyang Long,¹Wujun Xu,¹Li Liu,¹Binhui Wang,¹and Xiong Zou^3,4,5,6

Academic Editor: Arif Siddiqui

Received30 Apr 2022

Accepted14 Jun 2022

Published04 Jul 2022

Abstract

Background. Immunity exerts momentous functions in the progression and treatment of kidney renal clear cell carcinoma (KIRC). A better understanding of the relationship between KIRC and immunity may make a great contribution to evaluating the prognosis and immune-related therapeutic response of KIRC. Methods. A series of information such as RNA sequence, clinical data, and tumor mutation burden (TMB) of KIRC patients were downloaded through The Cancer Genome Atlas (TCGA). Next, combining the survival information and gene expression data of TCGA and Gene Expression Omnibus (GEO), we established an immune gene-related prognosis model (IGRPM) and analyzed it. Then we constructed a nomogram which was convenient for clinicians to judge the prognosis of KIRC. Last but not the least, the expressions of some genes used to construct IGRPM in early KIRC, and adjacent normal tissues were verified through real-time fluorescence quantitative polymerase chain reaction (RT-qPCR). Perl (strawberry-perl-5.30.0.1-64bit), R software (4.0.3), and GraphPad Prism 7 were used to process the relevant data. Results. The single-sample gene set enrichment analysis (ssGSEA) showed that there were significant differences in StromalScore, ImmuneScore, ESTIMATEScore, TumorPurity, 22 kinds of human immune cells infiltration, and HLA genes expression between high immunity group (Immunity_H) and low immunity group (Immunity_L). The Immunity_H expressed more immune-related genes and enriched more immune-related functions than the Immunity_L. In addition, compared with the low-risk group, the high-risk group had worse survival outcome and higher TMB. Combining IGRPM-based risk characteristic and TMB, we found that low-TMB + low-risk was the most beneficial to the survival outcome of KIRC patients. The risk characteristic based on IGRPM could be used as an independent prognostic factor for KIRC, and the nomogram constructed for evaluating the prognosis of KIRC showed excellent predictive potential. The RT-qPCR results suggested that not all the genes used to construct IGRPM showed differential expression in early KIRC compared with adjacent normal tissues, but all these genes had significant influence on the prognosis of KIRC. Conclusion. These comprehensive immune assessments and survival predictions, integrating multiple aspects of data and clinical information, can provide additional value to the current Tumor Node Metastasis staging system for risk stratification of KIRC and may facilitate the development of KIRC immunotherapy.

1. Introduction

Renal cell carcinoma (RCC), which originates from nephron [1], is one of the 10 most common malignancies in humans [2]. However, kidney renal clear cell carcinoma (KIRC) is one of the most invasive pathological types of RCC and the main pathological type [3], accounting for about 75% of all RCC [4]. Moreover, KIRC is a crucial cause of cancer-related death worldwide [5]. It is important to analyze the clinicopathological features of KIRC patients according to the WHO classification and Tumor Node Metastasis (TNM) staging for selecting appropriate treatment methods [6]. However, there may be significant differences in prognosis even in patients with the same TNM staging and pathological grade after similar treatment [7, 8], indicating that KIRC is a heterogeneous tumor, and the existing clinicopathological features and classification are not sufficient to evaluate the prognosis and risk stratification of KIRC patients. Therefore, finding a new method with high predictive value is of great significance to improve the prognosis of KIRC.

The development and occurrence of tumor and the effectiveness of immunotherapy are closely related to tumor immune microenvironment (TIM) [9, 10]. More and more studies have emphasized that there is heterogeneity in the proportion of immune-related genes and immune cell populations in tumors [11, 12], which may be one of the key reasons for different prognosis of patients. Oncologists are committed to exploring the immune responses and its regulatory mechanisms and screening immune-related biomarkers in order to coordinate the immune microenvironment of tumors and improve survival prognosis [13]. Immune cells are related to tumor metastasis and invasion and are potentially important forces to slow down or inhibit tumor growth [14]. However, although the research on the relationship between KIRC and immunity has increased in recent years, there is still a lack of research that can accurately predict the survival of KIRC patients in terms of immunity.

With the sharing of massive data, bioinformatics analysis has been widely used in the study of tumor prognosis [15]. In this study, we constructed an immune gene-related prognosis model (IGRPM) for patients with KIRC using data from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO). In addition, we systematically understood and distinguished the immune-related characteristics of KIRC, and looked into the correlation of tumor immune-related genes with TIM in depth. Additionally, we constructed a nomogram which was convenient for clinical application and verified its similarity with the ideal model. Finally, we verified the expressions of some of these genes used to construct IGRPM in early KIRC and adjacent normal tissues through real-time fluorescence quantitative polymerase chain reaction (RT-qPCR). Moreover, we confirmed that all the genes used to construct IGRPM had a significant effect on the prognosis of KIRC through UALCAN database.

2. Materials and Methods

2.1. Data Source

We downloaded transcriptome and clinical data of KIRC from TCGA (https://portal.gdc.cancer.gov/; December, 2021) and GEO (https://www.ncbi.nlm.nih.gov/geo/; December, 2021; GSE29609). The transcriptome data (HTseq-FPKM) of TCGA included 539 tumor samples and 72 paracancerous normal samples. 537 KIRC samples were downloaded from TCGA including clinical data. Masked somatic mutation data for 339 KIRC samples were downloaded from TCGA, of which tumor mutation burden (TMB) was successfully calculated for 336 samples. The GSE29609 cohort including 39 KIRC samples from GEO database and 530 KIRC samples from TCGA were used to construct prognostic models. During the data processing, we deleted some samples with incomplete information. The genetic data of related samples were converted into matrix files by Perl software (strawberry-perl-5.30.0.1-64bit). We analyzed the data through R software (4.0.3).

2.2. Single-Sample Gene Set Enrichment Analysis (ssGSEA) and Immune Classification Analysis

Using the ssGSEA algorithm of R software package (“GSEABase,” “limma,” and “GSVA”), the immune-related characteristics of each sample in TCGA-KIRC were evaluated comprehensively, based on 29 immune gene sets, including different immune cells, immune-related functions, and immune-related genes [16]. According to the immunological characteristics of the samples, 539 tumor samples from TCGA were analyzed by Ward’s linkage and Euclidean distance, which were divided into high immunity group (Immunity_H) and low immunity group (Immunity_L) [17]. The accuracy of immunotyping of KIRC patients was verified by R package (“Rtsne”).

2.3. Estimation of Tumor Purity and Immune Characteristics

We estimated the infiltration of stromal cells and immune cells in 539 KIRC tumor samples from TCGA by R software package (“estimate”), thus determining the StromalScore, ImmuneScore, ESTIMATEScore, and TumorPurity of the samples [18]. We also calculated the abundance of 22 types of human immune cells in each sample by CIBERSORT method [19] and compared the cell composition of the two subtypes in 539 KIRC samples of TCGA.

2.4. Differentially Expressed Immune Genes in Immunity_H and Immunity_L

We identified the differentially expressed genes (DEGs) between Immunity_H and Immunity_L by “limma” package [20]. Then, the DEGs were intersected with the immunologically relevant list of genes downloaded from ImmPort to obtain the differentially expressed immune genes (DEIGs) [21]. In addition, we also compared the differences in the expressions of HLA-related genes in Immunity_H and Immunity_L.

2.5. Functional Enrichment Analysis

We performed gene set enrichment analysis of 539 KIRC patients from TCGA by GSEA (Version4.1.0), in order to reveal the key signal pathways and functions involved in the two subtypes [22]. The screening thresholds of KEGG pathway and GO analysis in the Immunity_H and Immunity_L were . The enriched pathways were visualized with R packages (“reshape” and “ggplot2”).

2.6. Analyses of Prognosis Associated Immune Genes and Regulatory Networks

We obtained the DEIGs shared by TCGA and GEO through R packages (“limma” and “sva”). Next, we merged the shared DEIGs with survival data from TCGA and GEO into a single file and carried out univariate Cox proportional hazard regression analysis with R package (“survival”). It was considered that the survival risk gene of was the prognosis associated immune genes (PIGs) of KIRC patients. Next, we used the CISTROME (http://www.cistrome.org/) to find transcription factors (TFs) correlated to cancer occurrence and progression [23, 24], retrieved the differentially expressed transcription factors (DETFs) from DEGs, and constructed the regulatory network of DETFs and PIGs by Pearson’s correlation coefficient analysis. Correlation coefficient (|cor|) >0.4 and false discovery rate (FDR) <0.001 were the cut-off point of significant correlation. For coexpressed PIGs and DETFs, protein-protein interaction (PPI) analysis was also performed using STRING (https://cn.string-db.org/).

2.7. Construction of Immune Gene-Related Prognostic Model Based on TCGA and GEO Data

Using the LASSO Cox regression model based on R packages (“Matrix,” “glmnet,” and “survival”) and taking PIGs as the research object, an optimal immune gene-related prognosis model (IGRPM) was constructed for KIRC. The formula for calculating the risk score is as follows: the . In this formula, Expression_(k) denotes the expression level of gene in the TCGA and GEO samples, and coef_(k) denotes the regression coefficient of characteristic gene . All patients were separated into two groups based on IGRPM’s median risk score: low-risk and high-risk; and R packages (“survival” and “survminer”) were used to conduct the survival analysis. In order to reflect the predictive ability of risk characteristic based on IGRPM, we created a receiver operating characteristic curve (ROC) depending time through R packages (“survminer,” “survival,” “rms,” and “timeROC”) and calculated the area under the curve (AUC) of different survival times. Spearman correlation analysis of R software was used to explore the correlations of IGRPM-based risk characteristic with immune cells infiltration and TMB. Finally, we performed multivariate and univariate Cox regression analyses of IGRPM-based risk characteristic as well as other clinicopathological parameters to judge if IGRPM-based risk characteristic could be applied as an individual prognostic indicator for KIRC patients. Last but not the least, we constructed a nomogram through the risk characteristic based on IGRPM and the other clinicopathological features; what is more, the proposed nomogram was tested using ROC and calibration curve analyses to see if it was suitable for clinical application.

2.8. Real-Time Fluorescence Quantitative Polymerase Chain Reaction and UALCAN

Three early KIRC samples and adjacent normal tissues were collected from the Second Affiliated Hospital of University of South China. The clinicopathological characteristics of the patients were shown in Supplementary table 1. Total RNA was extracted by the FastPure Cell/Tissue Total RNA Isolation Kit V2 (RC112-01, Vazyme). PrimeScript RT reagent kit (RR037A, Takara) and Bio-Bad T100 Thermal Cycler according to the manufacturer’s instructions were used to perform reverse transcription to synthesize cDNAs. Reverse transcription was performed with a final quality of 1000 ng of RNA. For quantification of gene expression levels, we performed RT-qPCR using FastStart Universal SYBR Green Master (04913914001, Roche) and Light Cycler96 Instrument (Roche). The comparative cycle threshold (Ct) method was used for relative quantification of each target [25–27]. . The GAPDH was used as an internal reference. The 2^-ΔΔCt method was applied to calculate the expression of target genes in tumor tissues compared with adjacent normal tissues. All groups set up three biological replicates, and experiments were repeated three times for each sample. Meanwhile, we explored the effect of all genes used to construct IGRPM on the overall survival (OS) of KIRC through the UALCAN database. The following primers were used for RT-qPCR.

GAPDH:

5-ACAACTTTGGTATCGTGGAAGG-3(Forward),

5-GCCATCACGCCACAGTTTC-3(Reverse);

OASL (2-5-Oligoadenylate Synthetase Like):

5-CTCCCACACTCACATCTATCTG-3(Forward),

5-ACTGTCTTGGATGCCATAGATC-3(Reverse);

THRB (Thyroid Hormone Receptor Beta):

5-CTCATCAAAACTGTCACCGAAG-3(Forward),

5-CAACTTTTTGGCAAAATCCACC-3(Reverse);

CLDN4 (Claudin 4):

5-CCACAACATCATCCAAGACTTC-3(Forward),

5-CAGAATACTTGGCGGAGTAAGG-3(Reverse);

SAA1 (Serum Amyloid A1):

5-CAGACAAATACTTCCATGCTCG-3(Forward),

5-TCTCTGGATATTCTCTCTGGCA-3(Reverse);

NR3C2 (Nuclear Receptor Subfamily 3 Group C Member 2):

5-GCTGGAAGAAATGATTGCATCA-3(Forward),

5-AACTTCTTTGACTTTCGTGCTC-3(Reverse);

LGR4 (Leucine-Rich Repeat-Containing G Protein-Coupled Receptor 4):

5-CACACTTGGGCCAATAACTAAC-3(Forward);

5-ACAAAGTCTTTTGCTGCTAAGG-3(Reverse).

2.9. Statistical Analysis

The results of RT-qPCR were statistically analyzed by -test using GraphPad Prism 7. indicated statistically significant difference.

3. Results

3.1. Identification of Two Subtypes of KIRC By Immune Gene-Related Characteristics of Samples

To systematically characterize the immunological features of KIRC, the ssGSEA algorithm was used to evaluate 539 tumor samples from the TCGA-KIRC cohort through 29 immune cells or immune-related functions. According to ssGSEA analysis, the immune score and hierarchical clustering algorithm were performed for each sample, and the 539 tumor samples of TCGA were divided into two clusters: Immunity_L and Immunity_H (Figure 1(a)). Moreover, we further confirmed the cluster analysis of KIRC patients by using tSNE algorithm and obtained similar results (Figure 1(b)). In addition, we visually showed the scores of 29 immune cells or immune-related functions as well as immune microenvironment scores in Immunity_H and Immunity_L through heat map and violin plot (Figures 1(c) and 1(d)). According to the two pictures, it can be concluded that the StromalScore, ImmuneScore, and ESTIMATEScore in Immunity_H were higher than those in Immunity_L, while the Immunity_H’s TumorPurity was lower than the Immunity_L’s TumorPurity (wilcox.test).

(a)

(b)

(c)

(d)

3.2. Comparison of Immune-Related Characteristics in the Two Subtypes

CIBERSORT algorithm was used to compare the infiltration degree of 22 kinds of human immune cells in Immunity_H and Immunity_L. We discovered that the fractions of memory B cells, gamma delta T cells, follicular helper T cells, CD8 + T cells, regulatory T cells (Tregs), activated memory CD4 + T cells, resting memory CD4 + T cells, resting dendritic cells, and M1 macrophages in Immunity_H were significantly higher than that in Immunity_L; however, the fractions of resting Mast cells, M2 Macrophages, and M0 Macrophages in Immunity_H were significantly lower than that in Immunity_L (wilcox.test) (Figure 2(a)). Next, we also explored the expressions of HLA-related genes between Immunity_H and Immunity_L. To our surprise, the expressions of all HLA-related genes in Immunity_H were significantly higher than that in Immunity_L (Figure 2(b)). All these results indicated that it was very reasonable and of great significance for us to classify KIRC patients into the Immunity_H and Immunity_L.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2

The levels of immune cells infiltration, HLA genes expression, DEGs, and DEIGs analyses in Immunity_H and Immunity_L. (a) Comparison of the levels of immune cells infiltration between Immunity_H and Immunity_L. (b) Comparison of HLA genes expression levels between Immunity_H and Immunity_L. (c) Volcano plot of all DEGs, Immunity_H vs. Immunity_L, |LogFC| >0.585, and FDR <0.05. (d) Venn diagram of DEGs and immune genes. (e) Landscape of DEGs between Immunity_H vs. Immunity_L. (f) Landscape of DEIGs between Immunity_H and Immunity_L. , , and .

Furthermore, we analyzed the DEGs between subtypes and obtained the DEIGs based on ImmPort database. Compared with Immunity_L, Immunity_H had more DEGs upregulation (Figure 2(c)). 1811 genes in the two subtypes were authenticated as DEGs, of which 402 were DEIGs (Figures 2(d)–2(f)).

3.3. The KEGG and GO Analyses of Two Subtypes

We carried out KEGG and GO analyses according to DEGs in order to further understand the functional differences between Immunity_H and Immunity_L (Immunity_H vs. Immunity_L). We found that the KEGG pathways enriched in Immunity_H were mainly antigen processing and presentation, cytokine-cytokine receptor interaction, T cell receptor signaling pathway, natural killer cell-mediated cytotoxicity, JAK-STAT signaling pathway, cell adhesion molecules cams, toll-like receptor signaling pathway, chemokine signaling pathway, B cell receptor signaling pathway, NOD-like receptor signaling pathway, complement and coagulation cascades, and leukocyte transendothelial migration; all these pathways are closely related to immunity (Figure 3(a)). However, there was no significant correlation of the enriched pathways with immunity in Immunity_L (Figure 3(a)). Biological process (BP), cellular component (CC), and molecular function (MF) of GO analysis also showed similar results as KEGG (Figures 3(b)–3(d)). These results demonstrated that the immunotyping of KIRC was accurate and that the immune activity was significantly improved in the Immunity_H.

(a)

(b)

(c)

(d)

3.4. Screening of PIGs and Construction of Regulatory Networks

Next, combining the survival and gene expression data of TCGA and GEO, we identified 47 PIGs and calculated the hazard ratio (HR) by univariate Cox proportional hazard regression analysis (Figure 4(a)). There were 35 high-risk PIGs (HR >1) and 12 low-risk PIGs (HR <1). To better understand the multidimensional and complex regulatory network of PIGs involved in the onset and progression of KIRC, we first investigated the coexpression analysis of PIGs and DETFs. Surprisingly, all PIGs and DETFs that met our filter conditions showed a positive correlation (Supplementary table 2). The regulatory network between PIGs and DETFs was shown in Figure 4(b). In order to further confirm that there were significant correlations of DETFs with PIGs, we used STRING to analyze PPI and observed that there were complicated connections within them (Figure 4(c)).

(a)

(b)

(c)

3.5. Combining TCGA and GEO Samples to Construct an IGRPM of KIRC through Lasso Cox Regression Analysis

To better utilize of the PIGs of KIRC, the dimension of PIGs was reduced by lasso Cox regression analysis. Combining with the KIRC cohorts of TCGA and GEO, we constructed an IGRPM containing 17 key genes to predict OS (Figure 5(a)). The genes and coef used to determine the risk score for each subject were shown in Supplementary table 3. The corrected gene expression level, survival time, and survival status of TCGA and GEO samples were shown in Supplementary file1. We evaluated the risk score for KIRC patients of GSE29609 and TCGA and separated them into high-risk and low-risk groups based on the median risk score (Figure 5(b)). It is not difficult to find that the high-risk group had considerably more deaths than the low-risk group (Figure 5(c)). In addition, our further analysis showed that the risk score was negatively correlated with OS (Figure 5(d)), suggesting that the OS of KIRC patients gradually decreased with the elevation of the risk score. We also drew the Kaplan-Meier curve of KIRC patients with R packages (“survival” and “survminer”). The results demonstrated that patients in the high-risk group had a significantly shorter survival time than those in the low-risk group (, Figure 5(e)). What is more, we conducted a time-dependent ROC analysis to better examine the prediction accuracy of IGRPM-based risk characteristic, and the predicted AUC at 1, 3, and 5 years was 0.751, 0.717, and 0.728, respectively (Figure 5(f)). Additionally, the calibration diagram revealed the anticipated value of IGRPM-based risk characteristic was basically consistent with the real value (Figure 5(g)).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 5

Construction and analysis of IGRPM. (a) Cross-validation for tuning parameter selection in the lasso model. (b) Risk score distribution of KIRC patients in the combined TCGA and GEO cohorts. (c) Distribution of the survival time and survival status in the combined TCGA and GEO cohorts. (d) Correlation analysis of OS and risk score in the combined TCGA and GEO cohorts. (e) Kaplan-Meier survival analysis for different risk groups in the combined TCGA and GEO cohorts. (f) The 1-, 3-, and 5-year ROC and AUC predicted based on the risk characteristic of IGRPM in the combined TCGA and GEO cohorts. (g) The calibration diagram of IGRPM-based risk characteristic for the combined TCGA and GEO cohorts. OS: overall survival.

3.6. The Correlations between IGRPM-Based Risk Characteristic and Immune Cells Infiltration and TMB

Due to the absence of some data from GEO, we performed further analysis based on TCGA data. We analyzed the correlations between 17 genes used to construct IGRPM and 22 kinds of human immune cells by “spearman” method. Our results revealed that 17 key genes were closely associated with a variety of immune cells (Figure 6(a)), which once again demonstrated the importance of the immune microenvironment in the prognosis of KIRC patients. We further explored whether TMB differed in high-risk and low-risk groups. The results demonstrated that the TMB of high-risk group was substantially elevated compared with that of low-risk group (Figure 6(b)). Moreover, the survival time of patients with high TMB (H-TMB) was significantly shorter than that of patients with low TMB (L-TMB) (Figure 6(c)). Considering that both TMB and IGRPM had significant effects on the survival of patients with KIRC, we further investigated if the combination of IGRPM and TMB could better predict the survival outcome of KIRC patients. According to the risk score and TMB, we separated all samples from TCGA into 4 groups, including L-TMB + high-risk group, L-TMB + high-risk group, H-TMB + high-risk group, and H-TMB + low-risk group. The result revealed that the survival outcome of KIRC patients with L-TMB + low-risk group was the best in four groups (Figure 6(d)).

(a)

(b)

(c)

(d)

Figure 6

The correlations between IGRPM-based risk characteristic and immune cells infiltration and TMB. (a) Heat map showing correlation of 17 genes used to construct IGRPM with 22 kinds of human immune cells. (b) Comparison of TMB between the low-risk and high-risk groups in TCGA-KIRC cohort. (c) Kaplan-Meier survival analysis based on the TMB in the TCGA-KIRC cohort. (d) Combining with the risk characteristic based on IGRPM and the TMB in the TCGA-KIRC cohort, the survival analysis of the four groups stratified was carried out.

3.7. Independent Prognostic Analysis and Establishment and Verification of Nomogram

Since the risk score was substantially associated with the progression of KIRC, univariate and multivariate Cox regression analyses were performed to verify whether IGRPM-based risk characteristic could be used as an independent prognostic factor for KIRC patients. The results showed that the value of the risk score was less than 0.001 in both univariate and multivariate regression analyses (Figures 7(a) and 7(b)), which confirmed that the risk score based on IGRPM was an independent variable affecting the survival of KIRC patients. Combining with the patient’s age, sex, TNM staging, grade, and risk score, we designed a nomogram to expand clinical applicability of the risk characteristic based on IGRPM (Figure 7(c)). By scoring the comprehensive characteristics of each patient, in a certain range, the higher the total scored, the worse the prognosis of the patient would be. In addition, the ROC proved that nomogram could accurately judge the survival status of KIRC patients (Figure 7(d)). Moreover, the calibration diagram showed a good similarity between the performance of the nomogram and the ideal model (Figure 7(e)).

(a)

(b)

(c)

(d)

(e)

3.8. Validation of the Expression of Genes Used to Construct IGRPM in Early KIRC and Adjacent Normal Tissues

In order to further deepen the understanding of the genes used to construct IGRPM, we verified the expressions of some of these genes in early KIRC and adjacent normal tissues. The results of RT-qPCR indicated that the expressions of CLDN4, NR3C2, OASL, SAA1, and THRB in early KIRC were considerably suppressed than those in adjacent normal tissues (Figures 8(a)–8(e)); however, the expression of LGR4 in early KIRC was no statistical significance compared with that in adjacent normal tissues (Figure 8(f)). These results demonstrated that not all genes used to construct IGRPM were significantly differentially expressed in early KIRC compared with adjacent normal tissues. Last but not least, we explored the impact of all genes used to construct IGRPM on the prognosis of KIRC through the UALCAN database. To our surprising, the expressions of all these genes significantly affected the OS of KIRC (Figures 9(a)–9(q)), which showed that all these genes used to construct IGRPM had a significant impact on the prognosis of KIRC, although not all these genes showed differential expression in early KIRC compared with adjacent normal tissues. Therefore, it was more accurate to judge the prognosis of KIRC patients by integrating the expressions of these genes, which indicated that the IGRPM constructed by us was also more reliable and meaningful than the expression of a single or few genes to judge the prognosis of KIRC patients.

(a)

(b)

(c)

(d)

(e)

(f)

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

4. Discussion

As the prognosis of KIRC patients varies greatly, it is particularly important to establish a stable classifier for risk stratification of different KIRC patients, which can maximize the individualized treatment and timely follow-up of KIRC patients. Therefore, scientists have invested a lot of time and energy to explore the complex mechanism of KIRC, but the current understanding of KIRC, especially the understanding of TIM and immune-related prognostic factors of KIRC, is far from satisfactory. In this study, the multifaceted mining of transcriptome data, tumor immune-related characteristics, and TMB was designed to comprehensively understand the immune characteristics of KIRC and to construct a simple and practical method to help clinicians determine the prognosis of patients with KIRC, so as to develop personalized treatment options for KIRC patients.

First of all, we divided KIRC into two immune-related subtypes, Immunity_L and Immunity_H, and analyzed the molecular characteristics of subtype specificity, such as immune-related gene expression, function, and pathway enrichment. Immunity_H had lower tumor purity and higher immune score than Immunity_L, indicating that immune cells exerted critical role in inhibiting tumor progression. These findings supported previous researches that the activation of immune-related responses can effectively control the progression of KIRC [28, 29]. Additionally, Immunity_H had higher levels of HLA genes expression and immune cells infiltration than Immunity_L, indicating that Immunity_H had stronger immunogenicity. It has been reported that the expressions of HLA genes play an important role in antitumor activity and in the recognition and elimination of cancer cells [30, 31]. Then we visualized the DEGs and DEIGs in Immunity_H and Immunity_L, and it could be seen that the expressions of immune-related genes were found to be significantly different in the two groups. GO and KEGG analyses showed that Immunity_H had more abundant immune-related signals than Immunity_L. These results demonstrated the accuracy and necessity of dividing KIRC into two subtypes associated with immunity. Then we combined the gene expression data and survival data of TCGA and GEO to identify the genes that may have a significant impact on the prognosis of KIRC patients. We screened out 47 genes that had a significant impact on the prognosis of KIRC patients. Most of these genes have been reported in previous studies whose results were similar to ours. For example, it has been reported that NR3C2 expression could be employed as an independent prognostic factor of KIRC [32], and overexpression of NR3C2 have potential to restrict the proliferation and invasion of KIRC [33]. For HNF4A (Hepatocyte Nuclear Factor 4, Alpha), a tumor suppressor, its expression could affect the prognosis of KIRC [34]. Overexpression of IKBKE (Inhibitor of Kappa Light Polypeptide Gene Enhancer in B-cells, Kinase Epsilon) is associated with an increased risk of death in KIRC [35]. These results fully showed that the PIGs screened out by us had great reliability and that our method of screening prognostic genes had certain superiority. Then coexpression analysis of KIRC-related transcription factors and PIGs was performed, and their interactions were again demonstrated by PPI. It has been reported that PPARG (peroxisome proliferation-activated receptor gamma) promotes apoptosis and inhibits the migration and proliferation of KIRC by inhibiting SIX2 (SIX homeobox 2) [36]. BATF (Basic Leucine Zipper ATF-Like Transcription Factor) affects the prognosis of many kinds of cancer by combining with HHLA2 (HERV-H LTR-associating 2) DNA, especially the prognosis of KIRC [37]. CIITA (Class II Major Histocompatibility Complex Transactivator) methylation is associated with immune heterogeneity in hematological tumors [38]. These transcription factors are associated with tumors through different mechanisms. All these results provided a basis for us to further understand the occurrence and development mechanism of KIRC. Based on the above information, we could have a more comprehensive understanding of the relationship between KIRC and immunity, and dividing KIRC patients into different subtypes through immune-related genes would lay a foundation for us to estimate the prognosis of KIRC patients from a new perspective.

Based on the expressions of 17 immune genes, we constructed an IGRPM to assess the prognosis of patients with KIRC and evaluated the reliability of the model through ROC and calibration diagram. The results demonstrated that the risk score was negatively correlated with the OS of KIRC patients, and there was a substantial difference in survival time between the high-risk score group and low-risk score group. Moreover, most of the 17 genes used to construct IGRPM are related to the onset and progression of cancer and have the potential to become therapeutic targets for cancer, such as OASL [39], NR3C2 [32, 33], SAA1 [40], CXCL5 (C-X-C Motif Chemokine Ligand 5) [41], TLR3(toll-like receptor 3) [42], CLDN4 [43, 44], THRB [45], and LGR4 [46]. However, the specific molecular mechanisms of these genes affecting the prognosis of tumor patients are not clear, so the molecular mechanisms of the effect of these genes on KIRC patients need to be further studied in the future. Nevertheless, we speculated that the 17 genes affecting the prognosis of KIRC patients were closely related to the TIM, because these 17 genes were significantly associated with 22 types of human immune cells.

For the sake of further confirming the scientificity of classifying KIRC patients into the high- and low-risk groups, we compared the TMB of the two groups and found that the TMB of the high-risk group was substantially elevated compared with that of the low-risk group, and the survival time of KIRC patients with L-TMB was significantly longer than that of KIRC patients with H-TMB. This was consistent with the conclusion that TMB could be used as a biomarker to judge the prognosis of tumor patients reported in other researches [47, 48]. Similar results have been obtained from the research of lung squamous cell carcinoma [49]. Moreover, combining the risk score and TMB, we found that L-TMB + low-risk was most beneficial to the survival outcome of patients with KIRC, which provided a more accurate method for us to evaluate the prognosis of KIRC patients.

Finally, through univariate and multivariate Cox regression analyses, we proved that the risk characteristic based on IGRPM could be employed as an independent prognostic factor for patients with KIRC. Meanwhile, in order to make it convenient for clinicians to judge the prognosis of KIRC patients simply and effectively, we constructed a nomogram combining the clinicopathological information and the risk score of patients and verified its reliability by ROC and calibration diagram. Last but not the least, we verified the expressions of some of these genes used to construct IGRPM in early KIRC and adjacent normal tissue by RT-qPCR. The RT-qPCR results demonstrated that not all these genes used to construct IGRPM were significantly differentially expressed in early KIRC compared with adjacent normal tissue, but all these genes had significant influence on the prognosis of KIRC. Because of this, it further indicated that it was more accurate and effective to judge the prognosis of patients with KIRC by risk stratification based on the expression of these genes than by the expression of a single or few genes.

Although our study could make an important contribution to differentiating the immune status and predicting the prognosis of patients with KIRC, admittedly, there are some limitations to our study. On the one hand, the immune-related features of KIRC are constructed through public datasets, but these features have not yet been verified by clinical cohorts. On the other hand, due to the different information contained in TCGA and GEO databases, some data can only be obtained from TCGA, which makes some of the data less representative. Nevertheless, our results still have important guiding significance for survival prediction and immune assessment of patients with KIRC and may promote the development of immunotherapy for KIRC.

5. Conclusion

In a word, we elucidated the information of the interaction between tumor and immunity in KIRC from different aspects and used a variety of algorithms to establish IGRPM-based risk characteristic to evaluate the prognosis of KIRC patients. In addition, we constructed a nomogram which was convenient for clinicians to combine the clinicopathological information of patients to judge the prognosis of patients with KIRC. These comprehensive immune assessments and survival predictions, integrating multiple aspects of data and clinical information, can provide additional value to the current Tumor Node Metastasis staging system for risk stratification of KIRC and may facilitate the development of KIRC immunotherapy.

Data Availability

The data used to support the findings of this study are available in TCGA (https://gdc.cancer.gov/) and GEO databases. The raw data of RT-qPCR will be made available by the authors, without undue reservation, to any qualified researcher.

Ethical Approval

This research was approved by the Ethics Committee of the Second Affiliated Hospital, Hengyang Medical School, University of South China. The patients provided their written informed consent to participate in this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Xiaobin Wu and Xiong Zou designed the study. Xiong Zou, Xiaobin Wu, and Yonghui Liang collected the data. Xiaobin Wu and Xiong Zou performed RT-qPCR. Xiong Zou, Xiaobin Wu, Yonghui Liang, Xian Chen, Xiangyang Long, Wujun Xu, Li Liu, and Binhui Wang completed the data analysis and interpretation. Xiaobin Wu and Xiong Zou drafted and revised the article for important intellectual content. All authors have read and approved the manuscript.

Acknowledgments

This work was supported by the grant from the Key Guidance Project of Hunan Provincial Health Commission (20201908).

Supplementary Materials

table1: the characteristics of patients for RT-qPCR. Supplementary table 2: the coexpression relationship between PIGs and DETFs. Supplementary table 3: the gene and coef used to calculate the risk score for each sample. Supplementary file1: the corrected gene expression level, survival time, and survival status of TCGA and GEO samples. (Supplementary Materials)

References

Y. Luo, D. Shen, L. Chen et al., “Identification of 9 key genes and small molecule drugs in clear cell renal cell carcinoma,” Aging (Albany NY), vol. 11, no. 16, pp. 6029–6052, 2019.
View at: Publisher Site | Google Scholar
M. Akhtar, I. A. Al-Bozom, and H. T. Al, “Molecular and metabolic basis of clear cell carcinoma of the kidney,” Advances in Anatomic Pathology, vol. 25, no. 3, pp. 189–196, 2018.
View at: Publisher Site | Google Scholar
L. Vuong, R. R. Kotecha, M. H. Voss, and A. A. Hakimi, “Tumor microenvironment dynamics in clear-cell renal cell carcinoma,” Cancer Discovery, vol. 9, no. 10, pp. 1349–1357, 2019.
View at: Publisher Site | Google Scholar
R. E. Gray and G. T. Harris, “Renal cell carcinoma: diagnosis and management,” American Family Physician, vol. 99, no. 3, pp. 179–184, 2019.
View at: Google Scholar
R. L. Siegel, K. D. Miller, and A. Jemal, “Cancer statistics, 2020,” CA: a Cancer Journal for Clinicians, vol. 70, no. 1, pp. 7–30, 2020.
View at: Publisher Site | Google Scholar
R. J. Motzer, E. Jonasch, N. Agarwal et al., “Kidney cancer, version 2.2017, NCCN clinical practice guidelines in oncology,” Journal of the National Comprehensive Cancer Network, vol. 15, no. 6, pp. 804–834, 2017.
View at: Publisher Site | Google Scholar
J. H. Wei, Z. H. Feng, Y. Cao et al., “Predictive value of single-nucleotide polymorphism signature for recurrence in localised renal cell carcinoma: a retrospective analysis and multicentre validation study,” The Lancet Oncology, vol. 20, no. 4, pp. 591–600, 2019.
View at: Publisher Site | Google Scholar
X. Zou and Z. Mo, “CYP2J2 is a diagnostic and prognostic biomarker associated with immune infiltration in kidney renal clear cell carcinoma,” BioMed Research International, vol. 2021, Article ID 3771866, 15 pages, 2021.
View at: Publisher Site | Google Scholar
W. H. Fridman, L. Zitvogel, C. Sautes-Fridman, and G. Kroemer, “The immune contexture in cancer prognosis and treatment,” Nature Reviews. Clinical Oncology, vol. 14, no. 12, pp. 717–734, 2017.
View at: Publisher Site | Google Scholar
D. F. Quail and J. A. Joyce, “Microenvironmental regulation of tumor progression and metastasis,” Nature Medicine, vol. 19, no. 11, pp. 1423–1437, 2013.
View at: Publisher Site | Google Scholar
D. Zhang, W. He, C. Wu et al., “Scoring system for tumor-infiltrating lymphocytes and its prognostic value for gastric cancer,” Frontiers in Immunology, vol. 10, p. 71, 2019.
View at: Publisher Site | Google Scholar
D. Zeng, M. Li, R. Zhou et al., “Tumor microenvironment characterization in gastric cancer identifies prognostic and immunotherapeutically relevant gene signatures,” Cancer Immunology Research, vol. 7, no. 5, pp. 737–750, 2019.
View at: Publisher Site | Google Scholar
Y. Eso and H. Seno, “Current status of treatment with immune checkpoint inhibitors for gastrointestinal, hepatobiliary, and pancreatic cancers,” Therapeutic Advances in Gastroenterology, vol. 13, article 1756284820948773, p. 175628482094877, 2020.
View at: Publisher Site | Google Scholar
A. Hajiran, N. Chakiryan, A. M. Aydin et al., “Reconnaissance of tumor immune microenvironment spatial heterogeneity in metastatic renal cell carcinoma and correlation with immunotherapy response,” Clinical and Experimental Immunology, vol. 204, no. 1, pp. 96–106, 2021.
View at: Publisher Site | Google Scholar
Y. Jiang, Q. Zhang, Y. Hu et al., “ImmunoScore signature: a prognostic and predictive tool in gastric cancer,” Annals of Surgery, vol. 267, no. 3, pp. 504–513, 2018.
View at: Publisher Site | Google Scholar
S. Hanzelmann, R. Castelo, and J. Guinney, “GSVA: gene set variation analysis for microarray and RNA-seq data,” BMC Bioinformatics, vol. 14, no. 1, p. 7, 2013.
View at: Publisher Site | Google Scholar
X. Guan, Z. Y. Xu, R. Chen, J. J. Qin, and X. D. Cheng, “Identification of an immune gene-associated prognostic signature and its association with a poor prognosis in gastric cancer patients,” Frontiers in Oncology, vol. 10, article 629909, 2020.
View at: Publisher Site | Google Scholar
K. Yoshihara, M. Shahmoradgoli, E. Martínez et al., “Inferring tumour purity and stromal and immune cell admixture from expression data,” Nature Communications, vol. 4, no. 1, p. 2612, 2013.
View at: Publisher Site | Google Scholar
A. M. Newman, C. L. Liu, M. R. Green et al., “Robust enumeration of cell subsets from tissue expression profiles,” Nature Methods, vol. 12, no. 5, pp. 453–457, 2015.
View at: Publisher Site | Google Scholar
M. E. Ritchie, B. Phipson, D. Wu et al., “limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Research, vol. 43, no. 7, article e47, 2015.
View at: Publisher Site | Google Scholar
S. Bhattacharya, P. Dunn, C. G. Thomas et al., “ImmPort, toward repurposing of open access immunological assay data for translational and clinical research,” Sci Data., vol. 5, no. 1, article 180015, 2018.
View at: Publisher Site | Google Scholar
A. Subramanian, P. Tamayo, V. K. Mootha et al., “Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles,” Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 43, pp. 15545–15550, 2005.
View at: Publisher Site | Google Scholar
S. Mei, Q. Qin, Q. Wu et al., “Cistrome Data Browser: a data portal for ChIP-Seq and chromatin accessibility data in human and mouse,” Nucleic Acids Research, vol. 45, no. D1, pp. D658–D662, 2017.
View at: Publisher Site | Google Scholar
R. Zheng, C. Wan, S. Mei et al., “Cistrome Data Browser: expanded datasets and new tools for gene regulatory analysis,” Nucleic Acids Research, vol. 47, no. D1, pp. D729–D735, 2019.
View at: Publisher Site | Google Scholar
A. J. Siddiqui, J. Bhardwaj, M. Goyal et al., “Immune responses in liver and spleen against Plasmodium yoelii pre-erythrocytic stages in Swiss mice model,” Journal of Advanced Research, vol. 24, pp. 29–41, 2020.
View at: Publisher Site | Google Scholar
A. J. Siddiqui, J. Bhardwaj, M. Goyal et al., “Assessment of real-time method to detect liver parasite burden under different experimental conditions in mice infected with Plasmodium yoelii sporozoites,” Microbial Pathogenesis, vol. 89, pp. 35–42, 2015.
View at: Publisher Site | Google Scholar
A. J. Siddiqui, J. Bhardwaj, and S. K. Puri, “mRNA expression of cytokines and its impact on outcomes after infection with lethal and nonlethal Plasmodium vinckei parasites,” Parasitology Research, vol. 110, no. 4, pp. 1517–1524, 2012.
View at: Publisher Site | Google Scholar
N. Klümper, D. J. Ralser, R. Zarbl et al., “CTLA4promoter hypomethylation is a negative prognostic biomarker at initial diagnosis but predicts response and favorable outcome to anti-PD-1 based immunotherapy in clear cell renal cell carcinoma,” Journal for Immunotherapy of Cancer, vol. 9, no. 8, p. e002949, 2021.
View at: Publisher Site | Google Scholar
W. Zhong, H. Zhong, F. Zhang, C. Huang, Y. Lin, and J. Huang, “Characterization of hypoxia-related molecular subtypes in clear cell renal cell carcinoma to aid immunotherapy and targeted therapy via multi-omics analysis,” Frontiers in Molecular Biosciences, vol. 8, article 684050, 2021.
View at: Publisher Site | Google Scholar
A. Sadagopan, T. Michelakos, G. Boyiadzis, C. Ferrone, and S. Ferrone, “Human leukocyte antigen class I antigen-processing machinery upregulation by anticancer therapies in the era of checkpoint inhibitors: a review,” JAMA Oncology, vol. 8, no. 3, pp. 462–473, 2022.
View at: Publisher Site | Google Scholar
M. Kawazu, T. Ueno, K. Saeki et al., “HLA class I analysis provides insight into the genetic and epigenetic background of immune evasion in colorectal cancer with high microsatellite instability,” Gastroenterology, vol. 162, no. 3, pp. 799–812, 2022.
View at: Publisher Site | Google Scholar
Y. Liu, Z. Huang, G. Cheng et al., “Development of a four-gene prognostic model for clear cell renal cell carcinoma based on transcriptome analysis,” Genomics, vol. 113, no. 4, pp. 1816–1827, 2021.
View at: Publisher Site | Google Scholar
Z. Zhao, M. Zhang, X. Duan, T. Deng, H. Qiu, and G. Zeng, “Low NR3C2 levels correlate with aggressive features and poor prognosis in non-distant metastatic clear-cell renal cell carcinoma,” Journal of Cellular Physiology, vol. 233, no. 10, pp. 6825–6838, 2018.
View at: Publisher Site | Google Scholar
J. Lu, Z. Chen, H. Zhao et al., “ABAT and ALDH6A1, regulated by transcription factor HNF4A, suppress tumorigenic capability in clear cell renal cell carcinoma,” Journal of Translational Medicine, vol. 18, no. 1, p. 101, 2020.
View at: Publisher Site | Google Scholar
M. A. Hildebrandt, W. Tan, P. Tamboli et al., “Kinome expression profiling identifies IKBKE as a predictor of overall survival in clear cell renal cell carcinoma patients,” Carcinogenesis, vol. 33, no. 4, pp. 799–803, 2012.
View at: Publisher Site | Google Scholar
Y. Wu, T. Song, M. Liu et al., “PPARG negatively modulates Six2 in tumor formation of clear cell renal cell carcinoma,” DNA and Cell Biology, vol. 38, no. 7, pp. 700–707, 2019.
View at: Publisher Site | Google Scholar
B. Wang, Z. Ran, M. Liu, and Y. Ou, “Prognostic significance of potential immune checkpoint member HHLA2 in human tumors: a comprehensive analysis,” Frontiers in Immunology, vol. 10, p. 1573, 2019.
View at: Publisher Site | Google Scholar
O. Dufva, P. Pölönen, O. Brück et al., “Immunogenomic landscape of hematological malignancies,” Cancer Cell, vol. 38, no. 3, pp. 380–399.e13, 2020, e13.
View at: Publisher Site | Google Scholar
W. Zhong, C. Huang, J. Lin et al., “Development and validation of nine-RNA binding protein signature predicting overall survival for kidney renal clear cell carcinoma,” Frontiers in Genetics, vol. 11, article 568192, 2020.
View at: Publisher Site | Google Scholar
M. Djurec, O. Graña, A. Lee et al., “Saa3 is a key mediator of the protumorigenic properties of cancer-associated fibroblasts in pancreatic tumors,” Proceedings of the National Academy of Sciences of the United States of America, vol. 115, no. 6, pp. E1147–E1156, 2018.
View at: Publisher Site | Google Scholar
A. N. Hosein, G. Dangol, T. Okumura et al., “Loss of Rnf43 accelerates Kras-mediated neoplasia and remodels the tumor immune microenvironment in pancreatic adenocarcinoma,” Gastroenterology, vol. 162, no. 4, pp. 1303–1318.e18, 2022.
View at: Publisher Site | Google Scholar
X. Zou, B. Guo, Q. Ling, and Z. Mo, “Toll-like receptors serve as biomarkers for early diagnosis and prognosis assessment of kidney renal clear cell carcinoma by influencing the immune microenvironment: comprehensive bioinformatics analysis combined with experimental validation,” Frontiers in Molecular Biosciences, vol. 9, article 832238, 2022.
View at: Publisher Site | Google Scholar
P. Wang, T. Liu, Z. Zhao, Z. Wang, S. Liu, and X. Yang, “SPTBN2 regulated by miR-424-5p promotes endometrial cancer progression via CLDN4/PI3K/AKT axis,” Cell Death Discov., vol. 7, no. 1, p. 382, 2021.
View at: Publisher Site | Google Scholar
H. Wang, H. Xu, F. Ma et al., “Zinc finger protein 703 induces EMT and sorafenib resistance in hepatocellular carcinoma by transactivating CLDN4 expression,” Cell Death & Disease, vol. 11, no. 4, p. 225, 2020.
View at: Publisher Site | Google Scholar
J. Schnoell, U. Kotowski, B. J. Jank et al., “Prognostic relevance of thyroid-hormone-associated proteins in adenoid cystic carcinoma of the head and neck,” J Pers Med., vol. 11, no. 12, p. 1352, 2021.
View at: Publisher Site | Google Scholar
Z. Yue, Z. Yuan, L. Zeng et al., “LGR4 modulates breast cancer initiation, metastasis, and cancer stem cells,” The FASEB Journal, vol. 32, no. 5, pp. 2422–2437, 2018.
View at: Publisher Site | Google Scholar
T. A. Chan, M. Yarchoan, E. Jaffee et al., “Development of tumor mutation burden as an immunotherapy biomarker: utility for the oncology clinic,” Annals of Oncology, vol. 30, no. 1, pp. 44–56, 2019.
View at: Publisher Site | Google Scholar
M. Palmeri, J. Mehnert, A. W. Silk et al., “Real-world application of tumor mutational burden-high (TMB-high) and microsatellite instability (MSI) confirms their utility as immunotherapy biomarkers,” ESMO Open., vol. 7, no. 1, p. 100336, 2022.
View at: Publisher Site | Google Scholar
T. Jiang, J. Chen, X. Xu et al., “On-treatment blood TMB as predictors for camrelizumab plus chemotherapy in advanced lung squamous cell carcinoma: biomarker analysis of a phase III trial,” Molecular Cancer, vol. 21, no. 1, p. 4, 2022.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2022 Xiaobin Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

314

Downloads

447

Citations