Characterization of Molecular Heterogeneity Associated With Tumor Microenvironment in Clear Cell Renal Cell Carcinoma to Aid Immunotherapy

Clear cell renal cell carcinoma (ccRCC) is the most common type of kidney cancer and has strong immunogenicity. A systematically investigation of the tumor microenvironment (TME) in ccRCC could contribute to help clinicians develop personalized treatment and facilitate clinical decision-making. In this study, we analyzed the immune-related subtype of ccRCC on the basis of immune-related gene expression data in The Cancer Genome Atlas (TCGA, N = 512) and E-MTAB-1980 (N = 101) dataset, respectively. As a result, two subtypes (C1 and C2) were identified by performing non-negative matrix factorization clustering. Subtype C1 was characterized by increased advance ccRCC cases and immune-related pathways. A higher immune score, stromal score, TMB value, Tumor Immune Dysfunction and Exclusion (TIDE) prediction score, and immune checkpoint genes expression level were also observed in C1. In addition, the C1 subtype might benefit from chemotherapy and immunotherapy. The patients in subtype C2 had more metabolism-related pathways, higher tumor purity, and a better prognosis. Moreover, some small molecular compounds for the treatment of ccRCC were identified between the two subtypes by using the Connectivity Map (CMap) database. Finally, we constructed and validated an immune-related (IR) score to evaluate immune modification individually. A high IR score corresponded to a favorable prognosis compared to a low IR score, while more advanced tumor stage and grade cases were enriched in the low IR score group. The two IR score groups also showed a distinct divergence among immune status, TME, and chemotherapy. The external validation dataset (E-MTAB-1980) and another immunotherapy cohort (IMvigor 210) demonstrated that patients in the high IR score group had a significantly prolonged survival time and clinical benefits compared to the low IR score group. Together, characterization of molecular heterogeneity and IR signature may help develop new insights into the TME of ccRCC and provide new strategies for personalized treatment.

Introduction

Clear cell renal cell carcinoma (ccRCC) is a typical subtype of renal cell carcinoma (RCC), which is responsible for 70–80% of RCC cases (Martel and Lara, 2003). It was estimated that, in 2019, there would be 73,820 new cases and 14,770 deaths from ccRCC in America, which seriously affects human health (Siegel et al., 2016; Gao et al., 2020). Currently, the most common treatment strategy is surgical intervention, including laparoscopic partial nephrectomy (LPN) and radical nephrectomy (LN) (Porpiglia et al., 2014; Ahmad and Finelli, 2019). However, about 20∼30% of diagnostic RCC cases are metastatic (Gitlitz and Figlin, 2003; Yang et al., 2012). Moreover, despite early-stage RCC patients achieving a 5-year overall survival (OS) rate of 90%, the 5-year OS rate of metastatic or advance-stage RCC patients was less than 10% (Chaffer and Weinberg, 2011; Hirata et al., 2011). Cumulative studies have shown that targeted therapy is an effective method for the treatment of ccRCC, especially for metastatic cRCC patients (Grüllich, 2018). Targeted agents including axitinib (Escudier and Gore, 2011), pazopanib (Sternberg et al., 2010), sorafenib (Gordon, 2008), and sunitinib (Motzer, 2007) have made some progress in metastatic patients. However, the therapeutics were not effective for all patients and some showed drug resistance (Buchler et al., 2015; Baig et al., 2006). Thus, there is an urgent need to identify some molecular subtypes that may have implications on drug selection.

Tumor immunotherapy is a promising treatment method for tumors and has become a research hotspot in tumor therapy, and its efficacy is tightly associated with the TME (Wang et al., 2020). PD-1/PD-L1 is an important immunotherapy target for the treatment of ccRCC (Motzer et al., 2018). One study discovered that RCC patients with positive PD-L1 expression had a higher overall response rate to PD-1/PD-L1 treatment when compared to negative PD-L1 expression (Topalian et al., 2012). Moreover, dis-regulated genes also had an impact on the TME. Studies have shown that high expression level of HMGB1 in ccRCC can induce IL-10 secretion by T regulatory cells and decrease the anti-tumor activity of CD8+T cells (Li et al., 2017). Tim-3 expression in renal cell carcinoma is associated with the invasion of T cells and was considered as a novel therapeutic target for immunotherapy of ccRCC (Cai et al., 2015).

In the present study, we systematically investigated the immune-related (IR) molecular heterogeneity of ccRCC based on IR gene expression. We identified two subtypes with distinct immune microenvironment characteristics and immunotherapy response. These results lay a foundation for the further study of ccRCC immunotherapy and developing a personalized treatment.

Materials and Methods

Dataset Collection and Processing

The gene expression data and corresponding clinical information of ccRCC were downloaded from The Cancer Genomics Atlas (TCGA¹) and ArrayExpress database (E-MTAB-1980²), respectively. Among them, the TCGA dataset (training dataset) consisted of 539 ccRCC samples, and the E-MTAB-1980 dataset (external validation dataset) included 101 ccRCC samples (Sato et al., 2013). To obtain a reliable result, we excluded the samples with a survival time less than 30 days, and finally a total of 613 ccRCC (512 samples retrieved from TCGA, 101 samples retrieved from the E-MTAB-1980 dataset) samples were included in the downstream analysis. Moreover, the IR gene symbol names were retrieved from the Immunology Database and Analysis Portal (ImmPort³) database (Bhattacharya et al., 2018).

Molecular Subtypes Identification

The IR molecular subtype was identified by performing non-negative matrix factorization (NMF) clustering analysis in the “NMF” R package⁴ (Gillis, 2020). Principal component analysis (PCA) was performed to identify the robustness and reliability of the molecular subtype. Gene Set Variation Analysis (GSVA) was applied to identify the molecular function of each subtype (Hänzelmann et al., 2013). The reproducibility of the subtypes in the E-MTAB-1980 and TCGA datasets was identified by using the subclass mapping algorithm.

Exploration of the Relationship Between Subtype and Immune Cell Infiltration

To assess the immune cell infiltration of ccRCC, we firstly retrieved 28 IR cell gene sets from a previously reported article (Charoentong et al., 2016). We then applied the ssGSEA methods to estimate the relative abundance of each immune cell in ccRCC through using the “GSVA” R package (Hänzelmann et al., 2013). Moreover, the TME including immune score, stromal score, and tumor purity of each patient was calculated using the ESTIMATE algorithm⁵. The divergence between immune cell infiltration level, TME, and molecular subtype was evaluated using Wilcoxon rank sum test analysis.

Evaluating the Benefit of Each Subtype From Immunotherapy and Four Chemo Drugs

The chemotherapeutic response of each ccRCC sample was evaluated through the online database Genomics of Drug Sensitivity in Cancer (GDSC⁶) (Yang et al., 2016). Four common ccRCC treatment chemo drugs including Sorafenib, Sunitinib, Pazopanib, and Axitinib were selected (Liu et al., 2019). The prediction process was completed by using the “pRRophetic” R package and the half-maximal inhibitory concentration (IC50) of each ccRCC patient was calculated by using the ridge regression algorithm (Geeleher et al., 2014). The parameters were set as default. In addition, we also applied the Tumor Immune Dysfunction and Exclusion (TIDE)⁷ algorithm to evaluate the clinical response of ccRCC (Jiang et al., 2018).

Small Molecular Compounds Prediction for the Treatment of Clear Cell Renal Cell Carcinoma

To predict the small molecule compounds of ccRCC, we first identified the differentially expressed genes (DEGs) between subtypes C1 and C2 using the “limma” R package with | logFC| > 1 and adjusted P-value < 0.05 (Ritchie et al., 2015). Then the DEGs were uploaded to the Connective Map (CMap) database, respectively (Lamb et al., 2006). The small molecule compounds were further predicted based on the enrichment value and P-value.

Gene Ontology Terms and Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analysis

To investigate the molecular function of DEGs, we performed gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis to explore the molecular functions and pathways through the “ClusterProfiler” R package. Significant GO terms or pathways were screened with the cutoff adjusted P-value < 0.05 (Yu et al., 2012).

Construction of Immune-Related Score

Firstly, we used univariate Cox regression analysis to select the prognostic DEGs with a P-value < 0.05, and a total of 274 prognostic genes were selected for further analysis. Then, we performed PCA analysis to construct an IR relevant score. PC1 and PC2 were selected as the signature scores. The advantage of this method is that the score is concentrated on the set with the largest block of related (or anti-related) genes, while the weight is not combined with the contribution of the genes tracked by other set members (Zhang et al., 2020):

IR = Σ (PC1_i + PC2_i), where i represent the expression of each prognostic signature gene. To validate the prognostic value and clinical benefit, RNAseq datasets of 348 patients with bladder cancer were retrieved from the IMvigor 210 dataset (Mariathasan et al., 2018).

Statistical Analysis

All the analyses was carried out in the R/3.6.1 environment. The continuous data were evaluated using the Wilcoxon rank sum test. A log-rank test Kaplan–Meier (KM) curve was applied to evaluate the survival divergence of different subtypes. In addition, the Fisher test was employed to calculated the difference for the categorical data. For all tests, a P-value < 0.05 was regarded as statistically significant.

Results

Identification of Molecular Subtypes of Clear Cell Renal Cell Carcinoma

To investigate the IR subtype of ccRCC, we firstly retrieved the reads count data from the TCGA database and transformed it into TPM values. We then used the NMF algorithm to cluster the patients based on the IR gene expression data. To ensure a robust clustering result, we first removed the low expression level genes and subsequently included the selected genes in univariate Cox regression analysis. The genes with a P-value less than 0.05 were subjected to NMF clustering analysis in the TCGA dataset. The cophenetic correlation coefficients were applied to determine the optimal k value, and k = 2 was then selected as the optimal cluster number after a comprehensive consideration (Figure 1A). As shown in Figure 1B, when k = 2, we observed that the two subtypes (named C1 and C2) had clear boundaries, suggesting a stable and reliable clustering for the ccRCC sample. PCA was further applied to validate the assignments of the two subtypes. As shown in Figure 1C, the two-dimensional PCA distribution and the two subtypes had similar consistency. The survival curve revealed that C2 had a significantly better survival outcome when compared to C1 in overall survival (P-value < 0.001) (Figure 1D). In addition, we also performed NMF clustering analysis in the E-MTAB-1980 dataset, and the clustering result showed two distinct molecular subtypes, which was consistent with the training dataset result (Supplementary Figures 1A–C). The survival difference in the external validation dataset also showed a significant survival difference between the two subtypes (Supplementary Figure 1D) (P-value < 0.05). To validate the reproducibility of the molecular subtypes in the TCGA and E-MTAB-1980 datasets, we used the subclass mapping algorithm to compare the expression profiles of the subtypes from the two datasets. As Supplementary Figure 2 shows, we discovered that subtype C1 (N = 162) from TCGA was significantly associated with subtype C2 (N = 35) from the E-MTAB-1980 dataset, while subtype C2 (N = 350) from TCGA was significantly related to subtype C1 (N = 66) from the E-MTAB-1980 dataset (Bonferroni-corrected P-value = 0.003996).

FIGURE 1

Figure 1. Non-negative matrix factorization clustering analysis to identify potential molecular subtypes of ccRCC based on the immune-related gene expression in TCGA dataset. (A) The cophenetic correlation coefficient for the cluster number from 2 to 7. (B) Consensus matrix heatmap when k = 2. (C) Principal component analysis (PCA) for the 512 ccRCC patients, each dot represents a single sample. (D) KM survival curve analysis for the overall survival of the two subtypes in ccRCC.

Function Enrichment of the Molecular Subtype

To investigate the function of the two subtypes, we performed Gene Set Variation Analysis (GSVA) for the all gene expression datasets based on the subtypes. We then calculated the pathway activity for each sample. We discovered that metabolism-related pathways including fatty acid metabolism, beta alanine metabolism, tryptophan metabolism, and the propanoate metabolism pathway were upregulated in C2, while the p53 signaling pathway, natural killer cell-mediated cytotoxicity, cytokine cytokine receptor interaction, intestinal immune network for IGA production, and primary immunodeficiency pathway were enriched in C1 (Supplementary Figure 3).

The Relationship Between Molecular Subtype and Tumor Microenvironment

The tumor microenvironment (TME) is essential for immune function and has diverse clinical implications for immunotherapy. The divergence of immune cell infiltration and the TME were then explored in the two subtypes. The 28 immune cells’ infiltration level of each ccRCC patient was evaluated by applying the ssGSEA algorithm. As shown in Figure 2A, we observed that the level of activated B cells, CD4 T cells, CD8 T cells, central memory T cells, immature B cells, and T helper cells were significantly high in C1, while CD56dim natural killer cells and neutrophil cells presented at a lower level in C1. Moreover, we also calculated the immune score, stromal score, and tumor purity of each ccRCC sample using the ESTIMATE algorithm. We found that the immune and stromal score were significantly higher in subtype C1 when compared to subtype C2 (Figures 2B,C). However, the tumor purity showed a reverse trend that was lower in C1 (Figure 2D).

FIGURE 2

Figure 2. The landscape of immune cell infiltration in the two subtypes. (A) Twenty-eight immune cells’ infiltration levels that were calculated by the ssGSEA algorithm in the two subtypes. The comparisons of stromal score (B), immune score (C) and tumor purity (D) between the two subtypes. ns: not significant, *P < 0.05, **P < 0.01, ***P < 0.001.

Evaluation of Immunotherapeutic Response in the Clear Cell Renal Cell Carcinoma Subtypes

Recently, the immune checkpoint blockade has become a benchmark in the treatment of many tumors. Therefore, we subsequently explored the expression level of some known immune checkpoints including CD274 (PD-L1), PDCD1LG2 (PD-L2), PDCD1 (PD-1), LAG3, TIGIT, and CTLA4. As shown in Figures 3A–F, we found that the PD-L2, PD-1, LAG3, TIGIT, CTLA4 expression levels in subtype C1 were higher than in C2, indicating that patients in subtype C1 with a high expression of immune checkpoint inhibitors were more likely to form an immunosuppressive microenvironment, further leading to tumor immune escape (Dunn et al., 2002). The tumor mutational burden (TMB) has been demonstrated to be an effective predictor of immunotherapy response. In our result, we identified that the TMB level in subtype C2 was significantly lower than C1, which is consistent with a previous report where ccRCC patients’ overall survival corresponded to a low TMB level (Figure 3G; Huang et al., 2021). We also used the TIDE algorithm to estimate immunotherapy response and discovered that the TIDE value in subtype C1 was significantly higher than in C2, indicating that C2 was more likely to respond to immunotherapy (Figure 3H).

FIGURE 3

Figure 3. Evaluation of the immunotherapeutic response of the two subtypes in ccRCC. (A–F) The expression level of six immune checkpoint inhibitors including CD274 (PD-L1), PDCD1LG2 (PD-L2), PDCD1 (PD-1), LAG3, TIGIT, and CTLA-4 in the two molecular subtypes. (G) The comparisons of TMB value between the two subtypes. (H) The Tumor Immune Dysfunction and Exclusion (TIDE) prediction score of the two subtypes.

Chemo Drugs for the Treatment of Different Subtypes

Despite the fact that the curative effect of chemotherapy is limited in ccRCC, especially in advanced ccRCC, and considering that chemotherapy is the conventional therapy method for ccRCC, we thus attempted to evaluate the response of the two subtypes to four chemo drugs (Axitinib, Pazopanib, Sorafenib, Sunitinib). Therefore, we applied the ridge regression algorithm to train the GDSC cell line dataset, and obtained satisfactory prediction accuracy through 10-fold cross-validation. The IC50 value was estimated for each ccRCC sample according to the predictive model. As shown in Figures 4A–D, apart from Axitinib (P-value = 0.058), most of the chemo drugs including Pazopanib (P-value = 1.09e-14), Sorafenib (P-value = 2.14e-07), and Sunitinib (P-value = 4.82e-14) showed a lower IC₅₀ value in subtype C1, which indicated that patients in C1 were more likely to respond well to chemo drugs.

FIGURE 4

Figure 4. The IC₅₀ values of four chemo drugs in the two molecular subtypes of ccRCC including Axitinib (A), Pazopanib (B), Sorafenib (C), and Sunitinib (D).

Immunotherapy Response Prediction and Small Molecular Compounds Identification Between the Two Subtypes

Currently, immune checkpoint inhibitors are regarded as the routine drugs for the treatment of advanced ccRCC cases (Parikh and Bajwa, 2020). Here, we investigated the likelihood of response to immunotherapy by using the TIDE algorithm, and it demonstrated that patients in subtype C2 (108/350 = 0.309) were more likely to respond to immunotherapy than those in C1 (25/162 = 0.154) (Fisher exact test P-value = 0.0002). In addition to the TIDE prediction, we also compared the expression level of the two subtypes and identified another published dataset that included 47 patients with melanoma that responded to immune checkpoint inhibitors therapy by performing submap mapping analysis (Roh et al., 2017). Considering the worse prognosis of C1, we were delighted to discover that patients in subtype C1 were more sensitive to anti-PD-L1 therapy (Figure 5A). To identify the small molecular compounds of the ccRCC subtypes, we first performed differentially expressed analysis between C1 and C2. A total of 349 DEGs were identified, including 99 upregulated genes and 250 downregulated genes (Supplementary Table 1). These genes were then uploaded to the CMap database and the compounds with an absolute enrichment value more than 0.5 were retained. As a result, CMap mode-of-action (MoA) analysis of the 19 compounds indicated 16 mechanisms of action shared by the above drugs. As shown in Figure 5B, we can observed that scoulerine and suloctidil shared the MoA of adrenergic receptor antagonist, and alclometasone and piretanide shared the MoA of alucocorticoid receptor agonist. These results might provide potential therapeutic targets for the subtypes.

FIGURE 5

Figure 5. Evaluation of potential compounds and immunotherapy response to the two subtypes. (A) Submap analysis demonstrated that cluster 1 had more sensitivity to the programmed cell death protein 1 inhibitor (Bonferroni-corrected, P = 0.024). (B) CMap database analysis identified candidate drugs targeting the two molecular subtypes based on the DEGs.

Generation of Immune-Related Gene Signatures and Functional Annotation

To further uncover the function role of each molecular subtype, we firstly performed univariate Cox regression analysis on the 349 DEGs, and a total of 274 prognostic DEGs were identified. We then using the “clusterProfiler” R package to conduct GO and KEGG enrichment analysis on the 274 DEGs. As shown in Figure 6A, we discovered that the genes were mainly enriched in complement activation, B cell-mediated immunity, and humoral immune response in the biological processes (Supplementary Table 2). Moreover, the complement and coagulation cascades, NF-kappa B signaling pathway, and cytokine-cytokine receptor interaction were observed in the KEGG pathway enrichment result (Figure 6B and Supplementary Table 3). Interestingly, these GO terms and pathway results were similar to the GSEA results. We further conducted NMF clustering analysis based on the 274 DEGs to categorize patients into different immune genomic subtypes. Interestingly, two immune genomic subtypes (named gene clusters A and B) were identified (Figures 7A,B). The KM curve analysis result suggested that patients in gene cluster B had a better prognosis than A (P-value < 0.001) (Figure 7C). In addition, the PCA analysis could completely distinguish ccRCC samples between A and B (Figure 7D). Interestingly, we observed that most of the cases from subtype C2 were in gene cluster A (291/350 = 83.14%), while 96.30% (156/162) of cases from subtype C1 were classified into gene cluster B (Supplementary Figure 4). These results indicated that these DEGs were immune phenotype-related genes and might explain why gene cluster A was relevant to worse survival outcome. We further assessed the relationship between the gene clusters and clinical trait. As shown in Figure 8, we observed that most low-grade cases were enriched in gene cluster B, while advanced stage and grade cases were mainly focused on gene cluster A. Surprisingly, the patients in gene cluster B had large consistency with subtype C1. We also explored the difference of immune checkpoint expression between gene clusters A and B. As expected, except PD-L1 (Figure 9A), we found all immune checkpoints including PD-L2, PD-1, LAG3, TIGIT, and CTLA4 had a lower expression level in gene cluster B (Figures 9B–F). However, the infiltration level of activated B cells, activated CD4 T cells, activated CD8 T cells, central memory T cells, immature B cells, and T helper cells were higher in gene cluster A (Figure 9G).

FIGURE 6

Figure 6. GO (A) and KEGG (B) enrichment analysis was performed to identify potential function and pathway of DEGs.

FIGURE 7

Figure 7. The non-negative matrix factorization (NMF) clustering analysis to identify the genomic subtype of ccRCC based on the 274 DEGs. (A) Consensus matrix heatmap when k = 2. (B) The cophenetic correlation coefficient for the cluster number from 2 to 7. (C) KM survival curve analysis for the overall survival of the two subtypes in ccRCC. (D) Principal component analysis (PCA) for the 512 ccRCC patients, each dot represents a single sample.

FIGURE 8

Figure 8. The NMF clustering of DEGs in TCGA cohorts to categorize patients into two genomic subtypes (A and B). The gene clusters, immune-related subtype, tumor stage, age, survival status, gender, grade, and age were termed as patient annotations.

FIGURE 9

Figure 9. The evaluation of the immune checkpoint and immune cell infiltration in the two genomic clusters. The expression level of six immune checkpoint including CD274 (A), PDCD1LG2 (B), PDCD1 (C), LAG3 (D), TIGIT (E), CTLA4 (F) were showed in the two genomic clusters. (G) Twenty eight immune cells infiltration level were compared in the two genomic clusters. ns: not significant, *P < 0.05, **P < 0.01, ***P < 0.001.

Immune Score Construction

Considering the individual heterogeneity and complexity of IR subtypes, we developed a scoring system to quantify the immune pattern of individual patients with ccRCC based on 274 DEGs, named as the IR score. The KM curve analysis revealed that those with a high IR score had a better survival outcome (Figure 10A). To determine whether IR score can serve as a prognostic predictor, we performed ROC analysis and found that the IR score had the potential capability to ascertain the prognosis of ccRCC (Figure 10B). We further investigated the relationship between IR and clinical trait (grade and stage), and demonstrated that IRS was negatively correlated with grade and stage (Figures 10C,D). We also evaluated whether IR can serve as an independent prognostic factor via Cox regression analysis. As shown in Supplementary Figure 5A, we found that IR was a prognostic factor among the clinical factors in the univariate Cox regression analysis (P-value < 0.001), while there was no significant divergence in the multivariate Cox regression analysis (Supplementary Figure 5B). We further evaluated the immune cell infiltration and immune microenvironment, and observed that a high infiltration level of activated B cells, CD4 T cells, CD8 T cells, central memory T cells, immature B cells, and T helper cells was enriched in the low score group, while the level of CD56dim natural killer cells and neutrophils was higher in the high score group (Supplementary Figure 6). In addition, the immune and stromal scores were higher in subtype C1 when compared to subtype C2, whereas the level of tumor purity was low in the IR score group (Supplementary Figures 7A–C). Considering that the curative effect of chemotherapy is limited in ccRCC, we subsequently examined the sensitivity of the four chemo drugs to the IR score group. As shown in Figure 11, we found that the patients in the low IR score group showed more sensitivity to chemo drugs, indicating that the low IR score group may benefit more from chemo drugs. Finally, we identified several small molecular compounds between the low and high IR score groups. Among these compounds, ciprofibrate and clofibrate shared the MoA of PPAR receptor agonist, tiabendazole shared the MoA of angiogenesis inhibitor, and digitoxigenin shared the MoA of ATPase inhibitor (Supplementary Figure 8).

FIGURE 10

Figure 10. Construction of immune-related (IR) score via principal component analysis. (A) Kaplan–Meier curves for patients with high and low IR score subgroups in TCGA cohort. (B) The receiver operator curve analysis for the IR score. The relationship between IR score and tumor grade (C) and stage (D).

FIGURE 11

Figure 11. The IC₅₀ values of four chemo drugs in the IR score group of ccRCC including Axitinib (A), Pazopanib (B), Sorafenib (C), and Sunitinib (D).

Validation of the Immune-Related Score and the Role in Predicting Immunotherapeutic Benefits

To validate the prognostic value of the IR score, we retrieved the disease-free survival (DFS), disease-specific survival (DSS), and progression-free survival (PFS) datasets of ccRCC. As shown in Figures 12A–C, we discovered that a significant survival difference between the high IR score and low IR score groups existed in these datasets. We also validated it in the external validation dataset, and found that patients in the high IR score group had a better survival outcome (Figure 12D). The immunotherapy represented by PD-L1 and PD-1 blockades has achieved unprecedented success in cancer treatment. We then explored whether the IR score could predict patient response to PD-L1 and CTLA-4 blockade therapy based on IMvigor210. As expected, we were delighted to see that patients with a high IR score showed significant clinical benefits and prolonged survival (Figure 12E). The immune response and positive therapeutic benefits of patients in the distinct IR score group to immune checkpoint blockade treatment was further validated (Figure 12F).

FIGURE 12

Figure 12. Validation of the prognostic value of IR in disease-free survival (A), disease-special survival (B), and progression-free survival (C) of ccRCC, E-MTAB-1980 dataset (D), and IMvigor210 cohort (E), respectively. (F) The proportion of patients with response to PD-L1 blockade immunotherapy in low or high IR score groups. Responder/non-responder: 10%/90% in the low IR score groups and 25%/75% in the high IR score groups.

Discussion

Clear cell renal cell carcinoma is a heterogeneous disease that is associated with a poor prognosis in the advanced stage. How to prolong patient survival time and reduce death in advanced ccRCC remains a difficult clinical problem (Gao et al., 2020). Moreover, ccRCC has strong immunogenicity. ccRCC is characterized by a substantial number of inflammatory cells including NK cells, macrophages, T cells, and NK cells. Immune suppression is usually generated in the TME in advanced ccRCC through multiple mechanisms to avoid immune attack. For example, mechanisms weaken the presentation of effective antigens, reduce the impact of T cells, and promote immune tolerance, which are vital to tumor occurrence and development (Noessner et al., 2012). Thus, the TME plays a crucial role in the development of ccRCC. In recent years, immunotherapy has emerged as a new pillar for the treatment of tumors (Rini et al., 2019). The immune checkpoint blockade has attracted a lot of attention and has achieved remarkable results in immunotherapy, such as the use of PD-1/PD-L1 inhibitors and CTLA-4 inhibitors (Chae et al., 2018). However, the curative effects for some patients are still not satisfactory (Gulati and Vaishampayan, 2020). Several studies have demonstrated that individualized variation may contribute to the phenomenon (Takagi et al., 2020). Thus, exploring and understanding the subtypes of ccRCC will contribute to a better understanding of the heterogeneity of this disease (Jiang et al., 2020).

The ImmPort database collected a large number of known IR genes, which facilitated a comprehensive analysis of the role of IR genes in ccRCC (Bhattacharya et al., 2018). Wan et al. (2019) constructed a risk model to predict the prognosis of ccRCC based on seven IR genes including LAU, ISG15, IRF9, ARG2, RNASE2, SEMA3G, and UCN. Liao et al. (2021) used four prognostic IR genes (CRABP2, LTB4R, PTGER1, and TEK) to establish a gene signature to predict the prognosis of ccRCC. Zou and Hu (2020) developed a 14 IR gene signature to predict clinical outcome of kidney renal clear cell carcinoma. However, these studies mainly focused on the prognostic gene value and did not consider the heterogeneity and complexity of ccRCC. Thus, making a comprehensive analysis of IR subtypes will not only provide new insights into the molecular mechanism of ccRCC but also help us to better develop a personalized treatment in ccRCC. In the present study, we identified two IR molecular subtypes of ccRCC according to the gene expression profile. We found that patients in subtype C2 had a better survival outcome than those in C1. We then compared several known immune checkpoint’s (PD-L1, PD-L2, PD-1, LAG3, TIGIT, IDO1, CTLA-4) expression levels between the two subtypes. We found that most of the immune checkpoint inhibitors had a higher expression level in C1 when compared to C2. The previously studies reported that patients with a high expression level of immune checkpoint genes were more likely to form an immunosuppressive microenvironment and promote tumor immune escape (Dunn et al., 2002). These results suggest that patients in C1 may benefit from immune checkpoint inhibitor therapies. The TME (immune score, stromal score, and immune purity) and immune cell infiltration have been reported to be tightly associated with the prognosis and immunotherapy of cancers (Luo et al., 2020). Thus, we also investigated the relationship between subtype and TME and immune cell infiltration. Interestingly, we found that patients in subtype C1 tended to have a higher expression level of CD8⁺ and CD4⁺ T cell infiltration than those in C2. This indicates that the anti-tumor effect of high T cell infiltration is offset by the strong immunosuppressive pathway activated by over-expressed immune checkpoint proteins (Matsushita et al., 2016; Hua et al., 2020). However, further research is needed to demonstrate the potential molecular interactions between the molecular subtype and the immune cells of tumor immune status in ccRCC. In addition, we also discovered that the immune score and stromal score were higher in C1 when compared to C2. Previous studies have reported that high immune score and stromal score corresponded to poor survival (Luo et al., 2020). These findings were consistent with our results. Similarly, TMB value also showed a similar conclusion with the previous studies that presented a high level in C1, but a lower level in the C2 subtype (Huang et al., 2021).

Cumulative studies have demonstrated that chemo drugs including axitinib, pazopanib, sorafenib, and sunitinib have achieved some progress in advanced ccRCC cases (Koudijs et al., 2018). Here, we investigated the curative effect of four common chemo drugs for the subtypes C1 and C2. Interestingly, most of the IC₅₀ values of chemo drugs in the C1 subtype were lower than subtype C2, suggesting that patients in subtype C1 were more sensitive to chemo drugs and may benefit from these drugs. In addition, we also identified 18 small molecular compounds and 16 mechanisms of action for the treatment of the molecular subtypes. These drugs included adrenergic receptor antagonist (scoulerine and suloctidil), glucocorticoid receptor agonist (alclometasone and piretanide), adrenergic receptor agonist (oxymetazoline), protein synthesis inhibitor (puromycin), and HMGCR inhibitor (lovastatin). We also discovered other candidate drugs which might pave the way for implementation of subtype treatments for ccRCC patients.

To further reveal the molecular function and pathway that are involved in the subtypes, we conducted GSVA enrichment analysis for the two subtypes. We found that C2 corresponded to metabolism-related pathway, while C1 was associated with p53 signaling and cytokine-cytokine receptor interaction pathways. Studies have revealed that cytokine-cytokine receptor interaction is an important IR pathway that regulates the interaction between cytokines in tumors and is involved in the development and occurrence of tumors (Lee and Rhee, 2017). Considering that subtype C1 has a poor prognosis and more advanced ccRCC patients, it is reasonable that tumor-related pathways were enriched in subtype C1. In addition, we also performed limma analysis and obtained 274 DEGs between C1 and C2, and pathway enrichment results showed that DEGs mainly focused on the complement and coagulation cascades, NF-kappa B signaling pathway, and cytokine-cytokine receptor interaction, which is similar with our previously result, indicating that these DEGs were regarded as immune phenotype-related gene signatures. Similar to the IR gene clustering results, we also identified two genomic clusters based on these DEGs. The two genomic clusters were significantly related to distinct survival outcomes and TME landscapes. Due to the poor prognosis of ccRCC, more potential and valuable biomarkers are urgently needed. Zhang et al. (2021) developed an IR lncRNA-based model (AC012236.1, AC078778.1, AC078950.1, AC087318.1, and AC092535.4) for survival prediction in ccRCC. Liao et al. (2021) developed and validated the prognostic value of IR genes (CRABP2, LTB4R, PTGER1, and TEK) in ccRCC. Xiang et al. (2021) developed a four IR lncRNA signature for the prognosis of ccRCC through WGCNA analysis and Cox regression analysis. However, these risk models may be unstable and easily influenced by the gene expression level, especially considering the strong immunogenicity and heterozygosity of ccRCC. Thus, to guide therapeutic strategies for individual patients more precisely, we further constructed a IR score to quantify the IR patterns of individual tumors. The IR scores were closely associated with clinical information and can also serve as a prognostic biomarker for ccRCC survival. In addition, the prognostic value and therapeutic benefits of IR score were further validated in an external validation dataset and the IMvigor 210 cohort (Mariathasan et al., 2018).

Although we identified two molecular subtypes and constructed an IR score based on expression profile data of ccRCC, several limitations should be noted. Firstly, the samples with clinical information in the external validation dataset were relatively few and need to be expanded. Besides, due to the absence of an appropriate ICI-based ccRCC dataset, we used the different immunotherapy dataset across different malignancies (urothelial cancer) to validate the effects of IR score. Finally, all the discoveries were based on the dataset analysis and in vitro and in vivo experiments need to be performed to further verify these results.

In conclusion, we comprehensively investigated IR subtypes among 613 ccRCC samples based on the IR gene expression data, and systematically linked these subtypes with TME characteristics. We also evaluated the difference of chemo drug sensitivity and immunotherapy response between the subtypes. Finally, we developed an IR score to quantify the immune pattern of individual patients with ccRCC. These results may contribute to promoting our understanding of the characteristics of TME infiltration and provide new strategies for personalized treatment.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

JwH, YaL, and JyH designed the study. CH, YiL, and YY collected the clinical information and gene expression data. WZ and YY analyzed the data. WZ, YiL, and HZ wrote and revised the manuscript. All the authors contributed to the article and approved the submitted version.

Funding

This study was supported by the Health Science Research Personnel Training Program of Fujian Province (2017-CXB-22), Fujian Science and Technology Plan Guiding Projects (2019D026), and Xiamen Medical Advantage Subspecialty Vascular Access Construction Fund ([2018] 296). Moreover, we also thank the support of the research program fund of the fifth hospital of Xiamen.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We appreciate TCGA and the ArrayExpress database for providing the original study data. Moreover, our deepest gratitude goes to the reviewers and editors for their careful work and thoughtful suggestions that have helped improve this article substantially.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.736540/full#supplementary-material

Supplementary Figure 1 | The non-negative matrix factorization (NMF) clustering analysis to identify potential molecular subtypes of ccRCC based on the immune-related gene expression in the E-MTAB-1980 dataset. (A) The cophenetic correlation coefficient for the cluster number from 2 to 7. (B) Consensus matrix heatmap when k = 2. (C) Principal component analysis (PCA) for the 101 ccRCC patients, each dot represents a single sample. (D) KM survival curve analysis for the two subtypes.

Supplementary Figure 2 | Subclass mapping analysis revealed that subtype cluster 1 from TCGA was significantly associated with subtype cluster 2 from the E-MTAB-1980 dataset, while subtype cluster 2 from TCGA was significantly related to subtype cluster 1 from the E-MTAB-1980 dataset (Bonferroni-corrected P-value = 0.003996).

Supplementary Figure 3 | Gene Set Variant Analysis (GSVA) for the two subtypes on the basis of the 512 samples.

Supplementary Figure 4 | Sankey plot showing the changes of immune-related molecular subtypes, gene cluster, IR score, and patient survival status.

Supplementary Figure 5 | Identification of the independence of the IR score and clinical information through univariate Cox regression analysis (A) and multivariate Cox regression analysis (B).

Supplementary Figure 6 | The landscape of immune cells in the two IR score groups.

Supplementary Figure 7 | The comparisons of the stromal score (A), immune score (B), and tumor purity (C) in the IR score groups.

Supplementary Figure 8 | Identification of small molecular compounds by using the CMap drug database between the high IR score and low IR score groups.

Footnotes

1. ^ https://portal.gdc.cancer.gov/
2. ^ https://www.ebi.ac.uk/arrayexpress/
3. ^ https://www.immport.org/
4. ^ https://cran.r-project.org/web/packages/NMF/index.html
5. ^ http://r-forge.r-project.org
6. ^ https://www.cancerrxgene.org/
7. ^ http://tide.dfci.harvard.edu/login/

References

Ahmad, A., and Finelli, A. (2019). Renal function outcomes following radical or partial nephrectomy for localized renal cell carcinoma: should urologists rely on preoperative variables to predict renal function in the long term? Eur. Urol. 75, 773–774. doi: 10.1016/j.eururo.2018.12.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Baig, M. A., Lin, Y. S., Rasheed, J., and Mittman, N. (2006). Renal medullary carcinoma. J. Natl. Med. Assoc. 98, 1171–1174.

Google Scholar

Bhattacharya, S., Dunn, P., Thomas, C., Smith, B., Schaefer, H., Chen, J., et al. (2018). ImmPort, toward repurposing of open access immunological assay data for translational and clinical research. Sci. Data 5:180015. doi: 10.1038/sdata.2018.15

PubMed Abstract | CrossRef Full Text | Google Scholar

Buchler, T., Bortlicek, Z., Poprach, A., Pavlík, T., Veskrnova, V., Honzirkova, M., et al. (2015). Outcomes for patients with metastatic renal cell carcinoma achieving a complete response on targeted therapy: a registry-based analysis. Eur. Urol. 70, 469–475. doi: 10.1016/j.eururo.2015.12.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, C., Xu, Y.-F., Wu, Z.-J., Dong, Q., Li M-y, Olson, J., et al. (2015). Tim-3 expression represents dysfunctional tumor infiltrating T cells in renal cell carcinoma. World J. Urol. 34, 561–567. doi: 10.1007/s00345-015-1656-1657

CrossRef Full Text | Google Scholar

Chae, Y., Arya, A., Iams, W., Cruz, M., Chandra, S., Choi, J., et al. (2018). Current landscape and future of dual anti-CTLA4 and PD-1/PD-L1 blockade immunotherapy in cancer; lessons learned from clinical trials with melanoma and non-small cell lung cancer (NSCLC). J. ImmunoTherapy Cancer 6:39. doi: 10.1186/s40425-018-0349-343

CrossRef Full Text | Google Scholar

Chaffer, C., and Weinberg, R. (2011). A perspective on cancer cell metastasis. Science (New York, NY) 331, 1559–1564. doi: 10.1126/science.1203543

PubMed Abstract | CrossRef Full Text | Google Scholar

Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2016). Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. 18, 248–262.

Google Scholar

Dunn, G., Bruce, A., Ikeda, H., and Old, L. (2002). Cancer immunoediting: from immunosurveillance to tumor escape. Nat. Immunol. 3, 991–998. doi: 10.1038/ni1102-991

PubMed Abstract | CrossRef Full Text | Google Scholar

Escudier, B., and Gore, M. (2011). Axitinib for the management of metastatic renal cell carcinoma. Drugs RD 11, 113–126. doi: 10.2165/11591240-000000000-00000

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, X., Yang, J., and Chen, Y. (2020). Identification of a four immune-related genes signature based on an immunogenomic landscape analysis of clear cell renal cell carcinoma. J. Cell. Physiol. 235, 9834–9850. doi: 10.1002/jcp.29796

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N., and Huang, R. (2014). pRRophetic: an R Package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 9:e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Gillis, N. (2020). Nonnegative Matrix Factorization. Philadelphia, PA: SIAM.

Google Scholar

Gitlitz, B., and Figlin, R. (2003). Cytokine-based therapy for metastatic renal cell cancer. Urol. Clin. North Am. 30, 589–600. doi: 10.1016/S0094-0143(03)00027-22

CrossRef Full Text | Google Scholar

Gordon, M. S. (2008). Sorafenib in advanced clear-cell renal-cell carcinoma. Yearbook Oncol. 2008, 95–96. doi: 10.1016/S1040-1741(08)79114-79110

CrossRef Full Text | Google Scholar

Grüllich, C. (2018). Prospective single center trial of next-generation sequencing analysis in metastatic renal cell cancer: the MORE-TRIAL. Future Sci. OA 4:FSO299. doi: 10.4155/fsoa-2017-2150

CrossRef Full Text | Google Scholar

Gulati, S., and Vaishampayan, U. (2020). Current state of systemic therapies for advanced renal cell carcinoma. Curr. Oncol. Rep. 22:26. doi: 10.1007/s11912-020-0892-891

CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-Seq data. BMC Bioinformatics 14:7. doi: 10.1186/1471-2105-14-17

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirata, H., Hinoda, Y., Ueno, K., Nakajima, K., Ishii, N., and Dahiya, R. (2011). MicroRNA-1826 directly targets beta-catenin (CTNNB1) and MEK1 (MAP2K1) in VHL-inactivated renal cancer. Carcinogenesis 33, 501–508. doi: 10.1093/carcin/bgr302

PubMed Abstract | CrossRef Full Text | Google Scholar

Hua, X., Juan, C., Su, Y., and Liang, C. (2020). Identification of an immune-related risk signature for predicting prognosis in clear cell renal cell carcinoma. Aging 12, 2302–2332. doi: 10.18632/aging.102746

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, J., Li, Z., Fu, L., Lin, D., Wang, C., Wang, X., et al. (2021). Comprehensive characterization of tumor mutation burden in clear cell renal cell carcinoma based on the three independent cohorts. J. Cancer Res. Clin. Oncol. 147:1745. doi: 10.1007/s00432-020-03299-3293

CrossRef Full Text | Google Scholar

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24, 1550–1558. doi: 10.1038/s41591-018-0136-131

CrossRef Full Text | Google Scholar

Jiang, Y., Gou, X., Wei, Z., Tan, J., Yu, H., Zhou, X., et al. (2020). Bioinformatics profiling integrating a three immune-related long non-coding RNA signature as a prognostic model for clear cell renal cell carcinoma. Cancer Cell Int. 20:116. doi: 10.1186/s12935-020-01242-1247

CrossRef Full Text | Google Scholar

Koudijs, K., Scheltinga, A., Böhringer, S., Schimmel, K., and Guchelaar, H.-J. (2018). Personalised drug repositioning for clear cell renal cell carcinoma using gene expression. Sci. Rep. 8:5250. doi: 10.1038/s41598-018-23195-23198

CrossRef Full Text | Google Scholar

Lamb, J., Crawford, E., Peck, D., Modell, J., Blat, I., Wrobel, M., et al. (2006). The connectivity map: using gene-expression signatures to connect small molecules, genes, and disease. Science (New York, NY) 313, 1929–1935. doi: 10.1126/science.1132939

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, M., and Rhee, I. (2017). Cytokine signaling in tumor progression. Immune Network 17:214. doi: 10.4110/in.2017.17.4.214

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Sun, J., Rong, R., Li, L., Shang, W., Song, D., et al. (2017). HMGB1 promotes myeloid-derived suppressor cells and renal cell carcinoma immune escape. Oncotarget 8, 63290–63298. doi: 10.18632/oncotarget.18796

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, Z., Yao, H., Wei, J., Feng, Z., Chen, W., Luo, J., et al. (2021). Development and validation of the prognostic value of the immune-related genes in clear cell renal cell carcinoma. Trans. Androl. Urol. 10, 1607–1619. doi: 10.21037/tau-20-1348

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Zhou, L., Chen, Y., Liao, B., Ye, D., Wang, K., et al. (2019). Hypertension as a prognostic factor in metastatic renal cell carcinoma treated with tyrosine kinase inhibitors: a systematic review and meta-analysis. BMC Urol. 19:49. doi: 10.1186/s12894-019-0481-485

CrossRef Full Text | Google Scholar

Luo, J., Xie, Y., Zheng, Y., Wang, C., Qi, F., Hu, J., et al. (2020). Comprehensive insights on pivotal prognostic signature involved in clear cell renal cell carcinoma microenvironment using the ESTIMATE algorithm. Cancer Med. 9, 4310–4323. doi: 10.1002/cam4.2983

PubMed Abstract | CrossRef Full Text | Google Scholar

Mariathasan, S., Turley, S., Nickles, D., Castiglioni, A., Yuen, K., Wang, Y., et al. (2018). TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 554, 544–554. doi: 10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

Martel, C., and Lara, P. (2003). Renal cell carcinoma: current status and future directions. Crit. Rev. Oncology/Hematol. 45, 177–190. doi: 10.1016/S1040-8428(02)00076-78

CrossRef Full Text | Google Scholar

Matsushita, H., Sato, Y., Karasaki, T., Nakagawa, T., Kume, H., Ogawa, S., et al. (2016). Neoantigen load, antigen presentation machinery, and immune signatures determine prognosis in clear cell renal cell carcinoma. Cancer Immunol. Res. 4, 463–471. doi: 10.1158/2326-6066.CIR-15-0225

PubMed Abstract | CrossRef Full Text | Google Scholar

Motzer, R. J., Hutson, T. E., Tomczak, P., Michaelson, M. D., Bukowski, R. M., Rixe, O., et al. (2007). Sunitinib versus interferon alfa in metastatic renal-cell carcinoma. N. Engl. J. Med. 356, 115–124. doi: 10.1056/NEJMoa065044

PubMed Abstract | CrossRef Full Text | Google Scholar

Motzer, R., Tannir, N., McDermott, D., Aren, O., Melichar, B., Choueiri, T., et al. (2018). Nivolumab plus ipilimumab versus sunitinib in advanced renal-cell carcinoma. N. Engl. J. Med. 378, 1277–1290. doi: 10.1056/nejmoa1712126

PubMed Abstract | CrossRef Full Text | Google Scholar

Noessner, E., Brech, D., Mendler, A., Masouris, I., Schlenker, R., and Prinz, P. (2012). Intratumoral alterations of dendritic-cell differentiation and CD8+ T-cell anergy are immune escape mechanisms of clear cell renal cell carcinoma. Oncoimmunology 1, 1451–1453. doi: 10.4161/onci.21356

PubMed Abstract | CrossRef Full Text | Google Scholar

Parikh, M., and Bajwa, P. (2020). Immune checkpoint inhibitors in the treatment of renal cell carcinoma. Sem. Nephrol. 40, 76–85. doi: 10.1016/j.semnephrol.2019.12.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Porpiglia, F., Bertolo, R., Amparore, D., Cattaneo, G., and Fiori, C. (2014). Mini-Retroperitoneoscopic clampless partial nephrectomy for “low-complexity” renal tumours (PADUA Score <= 8). Eur. Urol. 66, 778–783. doi: 10.1016/j.eururo.2014.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Rini, B., Battle, D., Figlin, R., George, D., Hammers, H., Hutson, T., et al. (2019). The society for immunotherapy of cancer consensus statement on immunotherapy for the treatment of advanced renal cell carcinoma (RCC). J. ImmunoTherapy Cancer 7:81. doi: 10.1186/s40425-019-0813-818

CrossRef Full Text | Google Scholar

Ritchie, M., Phipson, B., Wu, D., Hu, Y., Law, C., Shi, W., et al. (2015). LIMMA powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Roh, W., Chen, P.-L., Reuben, A., Spencer, C. N., Prieto, P., Miller, J., et al. (2017). Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD-1 blockade reveals markers of response and resistance. Sci. Trans. Med. 9:eaah3560. doi: 10.1126/scitranslmed.aah3560

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, Y., Yoshizato, T., Shiraishi, Y., Maekawa, S., Okuno, Y., Kamura, T., et al. (2013). Integrated molecular analysis of clear-cell renal cell carcinoma. Nat. Genet. 45, 860–867. doi: 10.1038/ng.2699

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R., Miller, K., and Jemal, A. (2016). Cancer statistics, 2017. CA Cancer J. Clin. 67, 7–30. doi: 10.3322/caac.21387

PubMed Abstract | CrossRef Full Text | Google Scholar

Sternberg, C., Davis, I., Mardiak, J., Lee, E., Wagstaff, J., Barrios, C., et al. (2010). Pazopanib in locally advanced or metastatic renal cell carcinoma: results of a randomized Phase III trial. J. Clin. Oncol. 28, 1061–1068. doi: 10.1200/JCO.2009.23.9764

PubMed Abstract | CrossRef Full Text | Google Scholar

Takagi, T., Yoshida, K., Wada, A., Kondo, T., Fukuda, H., Ishihara, H., et al. (2020). Predictive factors for recurrence after partial nephrectomy for clinical T1 renal cell carcinoma: a retrospective study of 1227 cases from a single institution. Int. J. Clin. Oncol. 25, 892–898. doi: 10.1007/s10147-020-01632-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Topalian, S., Hodi, F., Brahmer, J., Gettinger, S., Smith, D., McDermott, D., et al. (2012). Safety, activity, and immune correlates of Anti–PD-1 antibody in cancer. N. Engl. J. Med. 366, 2443–2454. doi: 10.1056/NEJMoa1200690

PubMed Abstract | CrossRef Full Text | Google Scholar

Wan, B., Liu, B., Huang, Y., Yu, G., and Lv, C. (2019). Prognostic value of immune-related genes in clear cell renal cell carcinoma. Aging 11, 11474–11489. doi: 10.18632/aging.102548

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Xu, H., Zhu, L., He, T., Lv, W., and Wu, Z. (2020). Establishment and evaluation of a 6-Gene survival risk assessment model related to lung adenocarcinoma microenvironment. BioMed Res. Int. 2020, 1–8. doi: 10.1155/2020/6472153

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiang, Z., Shen, E., Li, M., Hu, D., Zhang, Z., and Yu, S. (2021). Potential prognostic biomarkers related to immunity in clear cell renal cell carcinoma using bioinformatic strategy. Bioengineered 12, 1773–1790. doi: 10.1080/21655979.2021.1924546

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Shi, L., Fu, Q., Xiong, H., Zhang, M., and Yu, S. (2012). Efficacy and safety of sorafenib in advanced renal cell carcinoma patients: results from a long-term study. Oncol. Lett. 3, 935–939. doi: 10.3892/ol.2012.585

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Lightfoot, H., Bignell, G., Behan, F., Cokelear, T., Haber, D., et al. (2016). Genomics of Drug Sensitivity in Cancer (GDSC): a resource for biomarker discovery in cancer cells. Eur. J. Cancer 69:S82. doi: 10.1016/S0959-8049(16)32839-32838

CrossRef Full Text | Google Scholar

Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for comparing biological themes among gene clusters. Omics 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y. (2020). M6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol. Cancer 19:53. doi: 10.1186/s12943-020-01170-1170

CrossRef Full Text | Google Scholar

Zhang, Z., Tang, Y., Liu, Y., Zhuang, H., Lin, E., Xie, L., et al. (2021). A novel immune-related lncRNA-Based model for survival prediction in clear cell renal cell carcinoma. J. Immunol. Res. 2021, 1–37. doi: 10.1155/2021/9921466

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, Y., and Hu, C. (2020). A 14 immune-related gene signature predicts clinical outcomes of kidney renal clear cell carcinoma. PeerJ 8:e10183. doi: 10.7717/peerj.10183

PubMed Abstract | CrossRef Full Text | Google Scholar