Immune-related Prognostic Genes of ccRCC

Introduction

Kidney cancer is one of the most commonly diagnosed tumors around the globe.1 According to the statistics from the World Health Organization, annually, there are more than 140,000 RCC-related deaths.2 ccRCC is the most typical subtype of kidney cancer and contributes to the majority of kidney cancer-related deaths.3,4 Until the recent past, the most widely used therapeutic agents were chemotherapy and cytokine therapy using interferon-alpha (INFα) and IL-2. However, these therapies produced effective responses only among a limited proportion of patients and adverse side effects among others.5 Moreover, the overall survival (OS) in advanced ccRCC is less than 30%.6,7 Therefore, the development of new treatment methods with clinical applicability in ccRCC is important.

Tumor microenvironment refers to the highly heterogeneous and dynamic tumor-associated network, consisting of tumor cells and non-tumor components, including blood vessels, immune cells, adipocytes, and cancer-associated fibroblasts.8,9 Numerous previous studies show that the tumor microenvironment plays an important role in the progression and prognosis of multiple cancer types. Generally, the von Hippel Lindau (VHL) gene, a tumor suppressor, is typically inactivated in ccRCC, which results in the overexpression of HIF-2a and VEGF and promotes the progression of ccRCC. Therefore, anti-angiogenic agents (AA), targeting the VEGF pathway, such as the VEGF receptor (VEGFR) and tyrosine kinase inhibitors (TKIs), have been recently used for the therapy of patients with ccRCC.10 Moreover, the immunogenicity of ccRCC has long been acknowledged. Several reports have examined the relationship between the existence of various immune cell types and patient prognosis.11–13 Giraldo et al show that CD8+ T cell infiltration is substantially correlated with poor prognosis in patients with ccRCC, unlike in other cancer types.14 Therefore, some immunotherapies based on immune checkpoint inhibitors (ICI), such as anti-programmed death receptor 1 (anti-PD1), are used for some patients with cancer5 and produce perfect outcomes. However, targeted therapy, such as TKIs and ICI, produce a complete response in a limited number of patients and potentially exerts adverse effects among some patients with cancers.15 In several clinical studies, it has been reported that there remains a considerable proportion of patients with no response or resistance to immune checkpoint inhibitors or anti-angiogenic agents,16 implying that the assessment of several immune checkpoint genes expression alone is insufficient in immunotherapy for most malignancies due to the low accuracy of prediction.17 Roviello et al performed a meta-analysis and found that the efficacy of immune checkpoint inhibitors in PD-L1 negative expression patients is not satisfactory.18 Therefore, it is necessary to detect novel immune-related targets for immunotherapy to elevate the efficacy of treatment.

Most previous studies have examined the role of immune cells in the progression and prognosis of cancers, however, they focus on the association between the immune cells and the prognosis of patients with cancers.9,19 With the advent of gene chips and high-throughput sequencing technologies, it is possible to have an in-depth understanding of the genetic characters of the tumor immune microenvironment through bioinformatic methods. Chen et al have identified prognostic immune genes in endometrioid endometrial adenocarcinoma by analyzing the high-throughput sequencing data using bioinformatic tools.20 Tian et al have analyzed the tumor-infiltrating immune cells (TIICs) and identified prognostic immune cells in adrenocortical carcinoma using bioinformatics.21 In this study, we utilized the high-throughput sequencing data and clinical information of ccRCC specimens extracted from the public databases to conduct immune cell infiltration and survival analyses. The survival-related immune cells and immune-specific genes were selected following the survival analysis. Eventually, four categories of survival-related TIICs and nine immune-specific hub genes were screened; these could be potential immunotherapeutic targets in ccRCC therapy.

Materials and Methods

Data Extraction and Processing

The workflow of this study is shown in Figure 1. The gene expression and corresponding clinical data were downloaded from public databases, Gene Expression Omnibus (GEO, www.ncbi.nlm.nih.gov/geo/) and the Cancer Genome Atlas (TCGA). In total, three ccRCC cohorts (TCGA-KIRC, GSE36895, and GSE53757) were analyzed in this study. TCGA-KIRC RNA sequencing data (count data) were downloaded from UCSC Xena (xenabrowser.net/datapages/), which consisted of 534 tumor- and 73 adjacent normal tissue samples. GSE36895 included 29 ccRCC tumor samples and 23 normal tissue samples.22 GSE53757 included 72 ccRCC tumor samples and 72 normal tissue samples.23 The information on patients from the selected ccRCC cohorts is shown in Table 1. This study design was reviewed by the ethics committee (Lanzhou University Second Hospital Ethics Committee) and deemed exempt for ethical approval.

Table 1 The Information of Clinical Traits

Figure 1 Flowchart of the study flow. CIBERSORT: an algorithm for analyzing the immune cell composition of complex tissues based on gene expression in those tissues.

Calculating the Abundance Ratio Matrix of Immune Cells

To assess immune infiltration in samples, we used the CIBERSORT algorithm, a tool for evaluating the immune cell composition in different organs based on their gene expression patterns.24 CIBERSORT source code and the LM22 gene signature were downloaded from the CIBERSORT website (cibersort.stanford.edu/download.php). The abundance ratios for 22 TIICs were calculated from the gene expression profiles of ccRCC by executing the CIBERSORT source code using the R software. The matrix of abundance ratio was visualized using the “barplot” function of R software and normalized using the “scale” function in R. The correlation between these TIICs was calculated using Pearson correlation through the “cor” function and visualized using the R package “corrplot”.25

Identification of Survival-Related Immune Cells and Grouping

To identify the survival-related immune cells from the 21 cell types calculated through CIBERSORT, we performed the survival analysis. The R package, “survival”, was used to select the survival-related tumor-infiltrating cells from 534 samples following the Kaplan–Meier (KM) analysis.26 To filter the immune cells, the LASSO regression analysis was performed to select the significant survival-related TIICs using the R package, “glmnet”.27 Next, the selected survival related TIICs were used to perform multivariate Cox regression analysis and the risk score for each sample was calculated using the relative ratio of TIICs multiplied by the regression coefficient. Samples were divided into high-risk and low-risk groups determined using the median score as a cut-off. The Nomogram based on multiple Cox regression for TIICs and clinical traits was constructed using the R package, “rms”.28 Subsequently, a calibration plot and time-dependent ROC were used to evaluate the prognostic value of the TIICs using the R packages, “rms”28 and “timeROC”,29 respectively.

Examining the Relationship Between Tumor-Infiltrating Cells and Clinical Characteristics

To understand the role of the immune cell types in affecting the clinical characteristics, we compared the immune infiltration score for different clinical parameters. The relationship between the abundance ratio of immune cells and clinical characteristics, including tumor or normal tissue types, tumor grade, and clinical stage were analyzed using the Wilcoxon test. Next, we plotted the immune cells for various clinical characters using the “boxplot” package in R software.

Identification of Differentially Expressed Genes and Enrichment Analysis

To explore the differences between the high-risk group and the low-risk group, we performed the analysis of differentially expressed genes (DEGs) and evaluated their functional enrichment. DEGs between the high-risk and low-risk groups were identified using the R package, “DESeq2”,30 with the cut-off set at P-value < 0.05 and |foldchange| > 1. Next, Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed using the R packages, “clusterProfiler31 and “org.Hs.eg.db”,32 with P-value < 0.05 set as the cut-off. Gene set enrichment analysis (GSEA) was performed using the GSEA function of the R package, “clusterProfiler”.31 The immunological signature gene set (C7 gene set) was selected as the reference for GSEA (www.gsea-msigdb.org/gsea/msigdb/index.jsp).

Construction of Protein–Protein Interaction Network and Identification of Hub Genes

To analyze the correlation among DEGs at the protein level, we constructed a protein–protein interaction (PPI) network for an in-depth understanding of these genes. The PPI network for the DEGs was constructed using the “STRING” plug-in in the Cytoscape software with a mutual score greater than 0.4 set as the threshold.33 Next, the “cytoHubba” plug-in in Cytoscape v3.8.2 was utilized to identify the top 30 genes in the network as the central genes based on four algorithms, namely Edge Percolated Component (EPC), betweenness, closeness, and degree.34 The hub genes were obtained through the intersection of DEGs among the four algorithms.

Multidimensional Validation of Hub Genes

To examine the relationship between hub genes and clinical characteristics, we analyzed their expression in different clinical grades and pathological stages. The outcome was visualized using the R package, “ggpur”.35 The correlation between hub genes and the OS was calculated using the Log rank test using the R package, “survival”.26 Subsequently, a multidimensional validation was performed to reduce the false-positive rates for these hub genes using Gene Expression Profiling Interactive Analysis (GEPIA) (gepia2.cancer-pku.cn/),36 UALCAN (ualcan.path.uab.edu/index.html),37 and the Human Protein Atlas (HPA) (www.proteinatlas.org/).38,39 Additionally, GSE3689522 and GSE5375723 served as the independent external datasets for the transcript-level validation of these hub genes. Finally, univariate and multivariate Cox analyses for these hub genes were used to examine their correlation with the patient prognosis in ccRCC.

Relationship Between Immune Cells and Hub Genes

Herein, we performed a correlation analysis to assess the relationship between the hub genes and immune cells, which can be further understand these genes in immune infiltration. The Spearman correlation between the hub genes and 22 immune cell types was calculated using the R package, “psych”.40 The R package, “pheatmap”, was utilized to visualize the results of the correlation analysis.41 The Tumor Immune Estimation Resource 2.0 (TIMER2.0) (timer.comp-genomics.org/) tool was utilized to validate the association between hub genes and immune cells.42

Statistical Analysis

In this study, groups of boxplots were analyzed using the Wilcoxon test. In the KM survival analysis, the P-value was estimated using Log rank tests, and the hazard ratio (HR) was calculated through the univariate Cox proportional hazards regression analysis. All analyses were performed in the R software version 4.0.3 (The R Foundation for Statistical Computing, 2020). All statistical tests were two-tailed; P-value < 0.05 was considered statistically significant.

Results

Infiltration Landscape of TIICs in ccRCC

The workflow of this research is shown in Figure 1. The information on these datasets is shown in Table 1. The abundance ratio matrix of 22 tumor-infiltrating cell types was calculated for 534 samples using CIBERSORT locally and is shown in Figure 2A and B. The T cells and Macrophages exhibited high infiltration in ccRCC. Moreover, we found that the abundance ratios of several immune cell types had a significant correlation with each other (Figure 2C). Mast cell resting was significantly positively correlated with T cell CD4 memory resting (Figure S1A) and Monocytes (Figure S1B); T cell regulatory (Tregs) were negatively correlated with Mast cells (Figure S1C) and T cell CD4 memory resting (Figure S1D).

Figure 2 (A) The overview of tumor-infiltrating immune cells; (B) Heatmap of tumor-infiltrating immune cell abundance ratio; (C) The coefficient of relationship between the different specific immune cells; (DK) The KM survival analysis for abundance ratios of different tumor-immune infiltrating cells.

Identification of Survival Related TIICs

The results of the KM survival analysis indicated that there were eight immune cell types related to the OS of ccRCC, including Dendritic cells resting, Eosinophils, Monocytes, Mast cells resting, T cell CD4 memory resting, T cell CD4 memory activated, T cell regulatory (Tregs), and Mast cells activated (Figure 2DK). Before multivariate Cox regression analysis, these eight immune cell types were identified through the LASSO regression analysis (Figure 3A and B). Following this, four types of immune cells, including T cells regulatory (Tregs), T cells CD4 memory activated, mast cells resting, and Eosinophils, were found to be possibly related to survival. The risk scores for each sample were calculated using the coefficients of the prior multivariate Cox regression analysis, which was as follows:

Figure 3 Identification of four prognostic immune cells. (A) The LASSO plot determines the number of OS-related genes in survival analysis via LASSO regression. (B) The Lambda plot determines the number of OS-related genes in survival analysis via LASSO regression. (C) The curve of risk score. (D) Survival status of the patients. The higher risk score corresponds to more proportion of dead patients. (E) Heatmap of the four prognostic immune cells relative score in low-risk and high-risk groups. (F) KM survival analysis of the high-risk group and low-risk group divided by the ratio of four types of immune cells. (G) The nomogram integrated ratio of four survival-related immune cells and clinical characteristics. (H) The calibration plot indicates an agreement test between predictions of OS for 1-year, 3-year, and 5-year survival and actual OS outcomes in TCGA dataset. (I) The time-dependent ROC curves for the nomogram based on immune cells and clinical traits in TCGA dataset.

where Frac is the abundance ratio of selected immune cells, and β is the coefficient of the Cox regression.

Subsequently, based on the median risk score as the cut-off, the patients were divided into the high-risk and low-risk groups (Figure 3C). Figure 3D shows the survival status of the patients. Figure 3E shows the heatmap for the fraction of the four immune cell types. As shown, Tregs had high infiltration in the high-risk group, which was in line with the findings of survival analysis. The KM survival analysis revealed that the high-risk group had a worse OS than the low-risk group, which confirmed that this classification was reasonable (Figure 3F). Then the nomogram based on the multivariate Cox regression was constructed for the prediction of overall survival in ccRCC (Figure 3G). The calibration curve plot (Figure 3H) and time-dependent ROC analysis for the nomogram (Figure 3I) showed satisfactory discrimination and calibration of the signature, which further indicated that these four immune cell types were strongly associated with OS in ccRCC.

Relationship Between TIICs and Clinical Characteristics

We examined the relationship between immune cells and clinical traits, including tumor or normal tissue type, clinical stage, and grade. The levels of T cells CD4 memory activated, T cells gamma delta, Macrophages M2, Dendritic cells resting, Mast cells activated, and Eosinophils were not significantly different in their infiltration in TCGA dataset (Figure 4A), and the finding was relatively consistent with those in the GSE53757 dataset (Figure 4B). The relationship between immune cells and grade in TCGA dataset is shown in Figure S2A; the relationship between immune cells and stage in TCGA dataset is shown in Figure S2B. The GSE53757 dataset was utilized to validate the association between immune cells and the clinical stage, and the findings were consistent with the results obtained in TCGA dataset (Figure S2C). The relationship between immune cells and clinical characters showed that the majority of immune cells exhibited differential infiltration in different clinical traits, in particular, several types of T cells showed significantly different infiltration across different stages and grades. These results demonstrated that the tumor immune microenvironment was associated with the progression of ccRCC.

Figure 4 The abundance ratios of immune cells between normal tissue and ccRCC tissues. (A) The comparison of the immune cells between normal and tumor tissues in TCGA dataset. (B) The comparison of the immune cells between normal and tumor tissues in the GSE53757 dataset. ****P < 0.0001, **P < 0.01, *P < 0.05.

Abbreviation: ns, not significant.

Identification and Enrichment of DEGs

The samples were separated into low-risk and high-risk groups to identify the survival-related DEGs in ccRCC, with the median risk score set as a cut-off value. In total, 707 DEGs were identified with P-value < 0.05 and |foldchange| > 1 set as the threshold criteria. The detailed information of these DEGs is shown in Table S1. The volcano plot of DEGs is shown in Figure 5A. The results of GO enrichment and KEGG pathway analyses for these upregulated genes are shown in Figure 5BE (Tables S2 and S3). Some enriched functions related to immune reactions and cytokines included T cell activation, cytokine-cytokine receptor interaction, and humoral immune response; these suggested that immune response and immune microenvironment may contribute substantially to the classification of risk group. For further confirmation of the strong correlation of DEGs with immunological characteristics, gene set enrichment analysis (GSEA) for the immunological signature (C7 gene set) was performed and the immune functions of upregulated DEGs were analyzed; A total of 574 gene sets were significantly enriched (P < 0.05). The top 10 gene sets are shown in Figure 5F (Table S4). These gene sets were mainly correlated with the T cell type infiltration, which implied that these DEGs were associated with immune infiltration.

Figure 5 Functional enrichment analysis of DEGs. (A) Volcano plot of differentially expressed genes between the high-risk group and low-risk group. Red nodes represent significantly upregulated genes with logFC > 1 and P < 0.05. Blue nodes represent significantly downregulated genes with logFC < −1 and P < 0.05. (B) Top 10 GO cellular component enrichment analysis of upregulated genes. (C) Top 10 of GO biological process enrichment analysis of upregulated genes. (D) Top 10 GO molecular function enrichment analysis of upregulated genes. (E) Top 10 KEGG enriched pathways of upregulated genes. (F) Top 10 among GSEA for upregulated genes.

Abbreviation: logFC, log fold change.

Protein–Protein Interaction Network and Hub Genes

To further understand these DEGs, we constructed a PPI network using the“STRING” app in Cytoscape software. From the four algorithms in the cytoHubba plug-in, the top 30 DEGs in each algorithm were used to create the protein interaction network (Figure 6AD). These genes in the protein interaction network played essential roles in the occurrence and development of the disease. The intersection of genes obtained from the four algorithms consisted of nine hub genes (Figure 6E), including APOE, CASR, CTLA4, CXCL8, EGF, F2, IL6, KNG1, and MMP9. The information of these nine hub genes is summarized in Table S5. The detailed results of the four algorithms are shown in Table S6. The hub genes consisted of several immune checkpoint inhibitor genes and cytokine genes, such as CTLA4, CXCL8, and IL6, which suggested that the method of obtaining the immune-related hub genes was reasonable and these hub genes warranted further evaluation.

Figure 6 Establishment of PPI network. (AD) The top 30 genes in the PPI network were identified as central genes using four algorithms: betweenness, EPC, closeness, and degree. (E) A Venn plot of a shared gene was obtained using four algorithms.

Abbreviations: EPC, edge percolated component; PPI, protein–protein interaction.

Validation of Hub Genes in Public Datasets

To further verify the correlation between the hub genes and clinical characteristics, their expressions according to different clinical traits were analyzed. CASR, CTLA4, EGF, F2, KNG1, and MMP9 were significantly differentially expressed between normal and cancer tissue types (Figure 7A), which was validated in both the GEO datasets, GSE36895 and GSE53757 (Figure 7B and Figure S3C). The expressions of hub genes in GSE36895 were consistent with the results obtained in TCGA cohort. The expression levels of APOE, CTLA4, F2, IL6, and MMP9 increased with the increase in clinical grade and pathological stage (Figure S3A and B), and this was validated in the GSE55375 dataset (Figure S3D). The Kaplan–Meier survival analysis showed that five hub genes, APOE, CTLA4, F2, IL6, and MMP9, were significantly associated with the OS (Figure 7C). GEPIA indicated that the expressions of these hub genes at the genomic level were consistent with our results (Figure S4). The Human Protein Atlas showed immunohistochemically verified expressions of APOE, CXCL8, F2, IL6, and KNG1 in renal carcinoma (Figure 7D). UALCAN analysis for these hub genes at the protein level is shown in Figure S5, and the findings were consistent at the mRNA level. These results indicated that hub genes were associated with clinical characters and prognosis both at the mRNA and protein levels, which suggested that these genes contributed to the development or progression of ccRCC.

Figure 7 The expression of nine hub genes between normal tissue and tumor tissues in TCGA dataset (A) and GSE36895 dataset (B). (C) KM survival analysis of the nine hub genes. The expression of CTLA4, CXCL8, F2, IL6, and MMP9 is significantly correlated with overall survival. (D) Immunohistochemistry (IHC) of five genes in renal cell carcinoma and normal kidney tissue. ****P < 0.0001; *P < 0.05.

Abbreviation: ns, not significant.

Validation of Hub Genes in TCGA Dataset Using the Prognostic Signature

To verify the relationship between these candidate genes and OS in patients with ccRCC, both univariate and multivariate Cox analyses were performed. The univariate Cox regression analysis showed that all gene expressions were significantly correlated with the OS of patients with ccRCC, except for two genes, EGF and KNG1 (Figure 8A). The multivariate Cox regression analysis for the seven genes showed that these genes were correlated with OS, and CASR (P < 0.05, HR = 0.9322) and CTLA4 (P < 0.01, HR = 1.1692) were independent prognostic factors for ccRCC patients (Figure 8B). The nomogram for the candidate genes and clinical traits, based on the multivariate Cox regression, was developed for prognostic prediction of ccRCC (Figure 8C), and the calibration curve showed strong discrimination and calibration power of the model (Figure 8D). Time-dependent ROC analysis demonstrated that the AUC values of the nomogram were 0.853, 0.806, and 0.777 for 1-year, 3-year, and 5-year, respectively (Figure 8EG); this indicated that the nomogram based on genes and clinical traits had a higher accuracy than the stage- and grade-based models. Decision curve analysis showed that the nomogram based on genes and clinical traits had greater benefits than the stage- and grade-based models when the threshold probabilities were set greater than 0.2 (Figure 8HJ). The evaluation of the Nomogram based on hub genes showed that they were well capable of predicting the prognosis in ccRCC, which indicated that these hub genes had a prognostic role in the ccRCC and deserve further analysis.

Figure 8 The results of prognostic signature based on hub genes and clinical traits. (A) The forest plot of univariate Cox regression analysis for hub genes. (B) The forest plot of multivariate Cox regression analysis for hub genes. (C) Nomogram integrated nine immune-related genes and clinical traits. (D) The calibration plot indicates an agreement test between predictions of OS for 1, 3, and 5 years and actual OS outcomes in the TCGA dataset. (EG) The time-dependent ROC curves for the nomogram and clinical trait model in TCGA dataset. (HJ) Decision curve analysis for the prediction of OS in 1-year, 3-year, and 5-year using the nomogram. The blue line indicates that there are no high-risk patients. The green line indicates that all patients are considered to be high-risk. The red line indicates that the nomogram can provide higher net benefits for accurately forecasting patient survival. **P < 0.01, *P < 0.05.

Relationship Between TIICs and Hub Genes

The relationship between hub genes and the TIICs is shown in Figure 9A, which demonstrated that the proportion of TIICs was linked to several hub genes. For instance, the expression of CTLA4 was positively related to the infiltration level of T cell CD8 and T cell follicular helper, while negatively related to the abundance ratio of mast cell resting. The expression of MMP9 was positively related to the proportion of Macrophages M0. These findings suggested that the hub genes may play a significant role in the function of the TIICs and are potential immunotherapy targets. Next, TIMER 2.0 was used to validate the correlation between hub genes and immune cells based on the XCELL algorithm (Figure 9B). The correlation estimated by XCELL was consistent with the result obtained using the CIBERSORT algorithm. The relationship between immune cells and hub genes in the GSE53757 dataset was evaluated and it was consistent with the findings in TCGA dataset (Figure S6).

Figure 9 The correlation between hub genes and immune cells. (A) The Spearman correlation coefficients between hub genes and immune cells in TCGA database are presented on a heatmap. (B) The relationship between hub genes and immune cells via XCELL algorithm.

Discussion

ccRCC is the most common type of kidney renal cell carcinoma and accounts for the majority of kidney cancer-related deaths. In recent years, several studies show that immunotherapy is beneficial for patients with ccRCC,43–45 however, the OS of patients with ccRCC remains unsatisfactory. Therefore, pursuing novel immune-related genes as therapeutic targets is critical for the development of effective therapy for ccRCC. The purpose of this study was to identify the immune cell types and immune-related genes correlated with the prognosis of ccRCC. These findings provided potential novel immunotherapeutic targets for ccRCC.

In this study, four hub immune cells were identified, including T cells regulatory (Tregs), T cells CD4 memory activated, Mast cells resting, and Eosinophils. Tregs are responsible for regulating the basic balance of the adaptive immune response, adjusting function, and maintaining homeostasis in response to diverse tissue contexts and immunological circumstances.20,46 Giraldo et al show that ICOS+ Tregs can be used as a prognostic marker in localized ccRCC.47 Some previous studies also show that CD4+ T cells are upregulated and are significantly correlated with OS in RCC.48–50 These previous results were consistent with our findings as a higher abundance ratio of Treg cells was identified in patients with worse OS. Nakanishi et al showed that mast cells played an important role in tumor growth, progression, and malignant aggressiveness of RCC tissues,51 which was consistent with our findings that the larger quantity of mast cells resting was associated with a better prognosis. These studies show, to some extent, that these four immune cell subgroups may be predictors of ccRCC prognosis. However, the abundance ratio of eosinophils was low in the majority of samples, and therefore the conclusion for eosinophils should be analyzed with caution.

Functional enrichment was performed for the DEGs between the high-risk and low-risk groups. The biological processes associated with upregulated DEGs were enriched in T cell activation and humoral immune response. These results implied that immune-related upregulated genes may be associated with the function of several T cell subtypes, which was also consistent with the findings in risk grouping based on immune cells. KEGG enrichment analysis demonstrated that most upregulated DEGs were enriched in cytokine-cytokine receptor interaction, chemokine signaling pathway, IL-17 signaling pathway, T cell receptor signaling pathway, and NF-κB signaling pathway. According to the findings, these signaling pathways played a crucial role in ccRCC. Previous studies have shown that cytokine-cytokine receptor interaction, chemokine signaling pathway, and IL-17 signaling pathway play significant roles in the proliferation, epithelial-mesenchymal transition, immune evasion, and metastases of several cancer types.52,53 Therefore, the regulation of these DEGs and tumor microenvironment could impact the tumor progression and prognosis in ccRCC.

Nine hub genes were selected through four algorithms in the cytoHubba plug-in in Cytoscape, including IL6, EGF, APOE, KNG1, MMP9, F2, CTLA4, CASR, and CXCL8. In the early stages of the immune responses, CTLA4 is expressed on the surface of CD4+ and CD8+ cells.54 CTLA4 binding reduces IL2 production and T cell proliferation.55 Some studies show that tumor cells can use the pathways regulated by CTLA4, including checkpoint protein signaling, to create a microenvironment that permits tumor cell growth.56,57 Anti-CTLA4 antibodies (ipilimumab, tremelimumab) are now being tested in different clinical trials as putative immune checkpoint inhibitors.55 In 2007, the first clinical study on anti-CTLA4 agents in ccRCC was undertaken, which demonstrated the efficacy of CTLA-4 inhibition in some patients with metastatic RCC.58 Thus, to a certain extent, the analysis of identification of hub genes is justified and other hub genes are also valuable to an in-depth investigation of the tumor immune microenvironment. IL6 is a cytokine produced by the immune and some tumor cells and is related to tumor proliferation and metastasis.59,60 These results were consistent with our findings as IL6 was highly expressed in grade 4 and stage 4. IL6 is a key factor in the differentiation of CD4+ T cell subsets and is required for the formation of T follicular helper cells.61 The IL6 receptor antibody, tocilizumab, is used for the treatment of rheumatoid arthritis.55 These studies indicate that IL6 correlates with immune cells and processes, which were also validated by IHC of renal cancer samples. The Epidermal growth factor, EGF, belongs to the epidermal growth factor superfamily. This protein is a robust mitogenic factor that plays a significant role in the growth, proliferation, and differentiation of numerous cell types.62 Polimeno et al show that EGF is a marker of host immunity in patients with renal cell carcinoma,50 which was also consistent with our findings. APOE is essential for the normal catabolism of triglyceride-rich lipoprotein constituents. Ostendorf et al demonstrate that APOE genotype is a biomarker for the outcome and therapeutic response in patients with melanoma.63 Tavazoie et al show that the APOE-associated axis regulates innate immune response suppression and plays an important role in enhancing the efficacy of cancer immunotherapy.64 Taken together, APOE influences the prognosis of cancer patients and correlates with immune therapy. However, there have been only a few studies on the association of APOE as a novel therapeutic target and renal cell carcinoma. A previous study also shows that MMP9 and CASR are hub genes significantly correlated with OS in ccRCC.6 In summary, among these hub genes, several genes have already been investigated in clinical trials, for example CTLA4, indicating these genes are associated with cancer therapy. Therefore, some of hub genes were novel genes that deserve in-depth investigation in the tumor immune environment.

However, the present study has certain limitations. First, the data of high-throughput sequencing of samples were downloaded from UCSC Xena, which is a depository for TCGA. Immune infiltration investigations were not considered at the time that the samples were acquired. Consequently, patients with immune system illnesses were enrolled, which may have resulted in selection bias at the beginning of the trial. Second, these survival-related immune cells and hub genes were identified from TCGA database, where the majority of patients are Americans or white. Considering the genetic heterogeneity, these immune cells and hub genes need to be validated in more public databases. Finally, this study lacks relevant in vivo and in vitro experiments to validate the functional roles of the identified hub genes.

Conclusion

In conclusion, four survival-related immune cell types and nine immune-related hub genes were identified using the CIBERSORT algorithm and regression analyses. Furthermore, these candidate hub genes were validated in clinical datasets and may also have implications as novel potential targets for ccRCC immunotherapy.

Abbreviations

ccRCC, clear cell renal cell carcinoma; KIRC, kidney renal clear cell carcinoma;KM, Kaplan–Meier; TIICs, tumor-infiltrating immune cells; IC, immune checkpoint; LASSO, Least Absolute Shrinkage and Selection Operator; HR, hazard ratio; DEGs, differentially expressed genes; GEPIA, Gene Expression Profiling Interactive Analysis; HPA, Human Protein Atlas; OS, overall survival; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Data Sharing Statement

All high-throughput sequencing data and clinical data can be obtained from UCSC Xena (xenabrowser.net/datapages/). The datasets, GSE36895 and GSE53757, can be obtained from the Gene Expression Omnibus (www.ncbi.nlm.nih.gov/geo/) database. The Immunohistochemistry (IHC) images of hub genes can be obtained from Human Protein Atlas (HPA) (www.proteinatlas.org/). The protein expression of genes can be obtained from UALCAN (ualcan.path.uab.edu/index.html).

Ethics Approval

Not applicable. All datasets used in this work are obtained from public databases and are freely available. This work did not include any experiments on humans or animals. The patients involved in the public databases have been enrolled after ethical approval. Additionally, this study was reviewed by the ethics committee (Lanzhou University Second Hospital Ethics Committee) and deemed exempt for ethical approval.

Acknowledgments

The information for this study was obtained from UCSC Xena, Gene Expression Omnibus (GEO). We are grateful to them for providing the source of data used in our study. We thank Bullet Edits Limited for the language editing assistance in this manuscript. We also thank Guozi Learn Bioinformatics (WeChat Official Accounts).

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

This study was supported by grants from the (1) Science and Technology Project of Chengguan District, Lanzhou City, Gansu Province Science and Technology Bureau (Project number: 2017KJGG0052); (2) the Fundamental Research Funds for the Central Universities (Project number: 561219007); (3) Cuiying Graduate Supervisor Applicant Training Program of Lanzhou University Second Hospital (Project number: 201704); (4) Cuiying Scientific and Technological Innovation Program of Lanzhou University Second Hospital (Project number: CY2017-BJ16); (5) Lanzhou City Talent Innovation and Entrepreneurship Project (Project number: 2019-RC-37), (6) Gansu Provincial Department of Education: Outstanding Postgraduate “Innovation Star” Project (Project number: 2021CXZX-154), and (7) Gansu Health Industry Research Project(Project number: GSWSKY2017-10).

Disclosure

The authors report no potential conflicts of interest.

References

1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA a Cancer J Clin. 2020;70(1):7–30. doi:10.3322/caac.21590

2. Ferlay J, Soerjomataram I, Dikshit R, et al. Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012: globocan 2012. Int J Cancer. 2015;136(5):E359–E386. doi:10.1002/ijc.29210

3. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68(6):394–424. doi:10.3322/caac.21492

4. Zhang Z, Lin E, Zhuang H, et al. Construction of a novel gene-based model for prognosis prediction of clear cell renal cell carcinoma. Cancer Cell Int. 2020;20(1):27. doi:10.1186/s12935-020-1113-6

5. Hammers H. Immunotherapy in kidney cancer: the past, present, and future. Curr Opin Urol. 2016;26(6):543–547. doi:10.1097/MOU.0000000000000338

6. Du B, Zhou Y, Yi X, et al. Identification of immune-related cells and genes in tumor microenvironment of clear cell renal cell carcinoma. Front Oncol. 2020;10:1770. doi:10.3389/fonc.2020.01770

7. Motzer RJ, Escudier B, McDermott DF, et al. Nivolumab versus everolimus in advanced renal-cell carcinoma. N Engl J Med. 2015;373(19):1803–1813. doi:10.1056/NEJMoa1510665

8. Clark DJ, Dhanasekaran SM, Petralia F, et al. Integrated proteogenomic characterization of clear cell renal cell carcinoma. Cell. 2019;179(4):964–983.e31. doi:10.1016/j.cell.2019.10.007

9. Şenbabaoğlu Y, Gejman RS, Winer AG, et al. Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol. 2016;17(1):231. doi:10.1186/s13059-016-1092-z

10. Heidegger I, Pircher A, Pichler R. Targeting the tumor microenvironment in renal cell cancer biology and therapy. Front Oncol. 2019;9:490. doi:10.3389/fonc.2019.00490

11. Giraldo NA, Sanchez-Salas R, Peske JD, et al. The clinical role of the TME in solid cancer. Br J Cancer. 2019;120(1):45–53. doi:10.1038/s41416-018-0327-z

12. Sadeghi Rad H, Monkman J, Warkiani ME, et al. Understanding the tumor microenvironment for effective immunotherapy. Med Res Rev. 2021;41(3):1474–1498. doi:10.1002/med.21765

13. Binnewies M, Roberts EW, Kersten K, et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med. 2018;24(5):541–550. doi:10.1038/s41591-018-0014-x

14. Giraldo NA, Becht E, Pagès F, et al. Orchestration and prognostic significance of immune checkpoints in the microenvironment of primary and metastatic renal cell cancer. Clin Cancer Res. 2015;21(13):3031–3040. doi:10.1158/1078-0432.CCR-14-2926

15. Lai Y, Tang F, Huang Y, et al. The tumour microenvironment and metabolism in renal cell carcinoma targeted or immune therapy. J Cell Physiol. 2021;236(3):1616–1627. doi:10.1002/jcp.29969

16. Zhong J, Liu Z, Cai C, Duan X, Deng T, Zeng G. A modification patterns and tumor immune landscape in clear cell renal carcinoma. J Immunother Cancer. 2021;9(2):e001646. doi:10.1136/jitc-2020-001646

17. Zhao Q, Xie R, Lin S, You X, Weng X. Anti-PD-1/PD-L1 antibody therapy for pretreated advanced or metastatic nonsmall cell lung carcinomas and the Correlation between PD-L1 expression and treatment effectiveness: an update meta-analysis of randomized clinical trials. Biomed Res Int. 2018;2018:3820956. doi:10.1155/2018/3820956

18. Roviello G, Corona SP, Nesi G, Mini E. Results from a meta-analysis of immune checkpoint inhibitors in first-line renal cancer patients: does PD-L1 matter? Ther Adv Med Oncol. 2019;11:1758835919861905. doi:10.1177/1758835919861905

19. George DJ, Martini JF, Staehler M, et al. Immune biomarkers predictive for disease-free survival with adjuvant sunitinib in high-risk locoregional renal cell carcinoma: from randomized phase III S-TRAC study. Clin Cancer Res. 2018;24(7):1554–1561. doi:10.1158/1078-0432.CCR-17-2822

20. Chen B, Wang D, Li J, Hou Y, Qiao C. Screening and Identification of prognostic tumor-infiltrating immune cells and genes of endometrioid endometrial adenocarcinoma: based on the cancer genome atlas database and bioinformatics. Front Oncol. 2020;10:554214. doi:10.3389/fonc.2020.554214

21. Tian X, Xu W, Wang Y, et al. Identification of tumor-infiltrating immune cells and prognostic validation of tumor-infiltrating mast cells in adrenocortical carcinoma: results from bioinformatics and real-world data. Oncoimmunology. 2020;9(1):1784529. doi:10.1080/2162402X.2020.1784529

22. Peña-Llopis S, Brugarolas J. Simultaneous isolation of high-quality DNA, RNA, miRNA and proteins from tissues for genomic applications. Nat Protoc. 2013;8(11):2240–2255. doi:10.1038/nprot.2013.141

23. von Roemeling CA, Radisky DC, Marlow LA, et al. Neuronal pentraxin 2 supports clear cell renal cell carcinoma by activating the AMPA-selective glutamate receptor-4. Cancer Res. 2014;74(17):4796–4810. doi:10.1158/0008-5472.CAN-14-0210

24. Newman AM, Steen CB, Liu CL, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. 2019;37(7):773–782. doi:10.1038/s41587-019-0114-2

25. Wei T, Viliam S. R package “corrplot”: visualization of a correlation matrix (version 0.84). 2017.

26. Therneau TM. A package for survival analysis in R. R package version 32–7. 2020.

27. Friedman JH, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33(1):1–22. doi:10.18637/jss.v033.i01

28. Frank EHJ. rms: regression modeling strategies. R package version 51–4. 2019.

29. Jacqmin-Gadda PB, Dartigues J-F, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med. 2013;32:5381–5397. doi:10.1002/sim.5958

30. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. doi:10.1186/s13059-014-0550-8

31. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–287. doi:10.1089/omi.2011.0118

32. Carlson M. org.Hs.eg.db: genome wide annotation for human. R package version 3120. 2020.

33. Doncheva NT, Morris JH, Gorodkin J, Jensen LJ. Cytoscape stringapp: network analysis and visualization of proteomics data. J Proteome Res. 2019;18(2):623–632. doi:10.1021/acs.jproteome.8b00702

34. Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8(Suppl 4):S11. doi:10.1186/1752-0509-8-S4-S11

35. Kassambara A. ggpubr: “ggplot2” based publication ready plots. R package version 040. 2020.

36. Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 2017;45(W1):W98–W102. doi:10.1093/nar/gkx247

37. Chandrashekar DS, Bashel B, Balasubramanya SAH, et al. UALCAN: a portal for facilitating tumor subgroup gene expression and survival analyses. Neoplasia. 2017;19(8):649–658. doi:10.1016/j.neo.2017.05.002

38. Uhlen M, Zhang C, Lee S, et al. A pathology atlas of the human cancer transcriptome. Science. 2017;357(6352):eaan2507. doi:10.1126/science.aan2507

39. Thul PJ, Åkesson L, Wiking M, et al. A subcellular map of the human proteome. Science. 2017;356(6340):eaal3321. doi:10.1126/science.aal3321

40. Revelle W. Psych: procedures for psychological, psychometric, and personality research. R package version 2012. 2020.

41. Kolde R. Pheatmap: pretty Heatmaps. R package version 1012. 2019.

42. Ju A, Tang J, Chen S, Fu Y, Luo Y. Pyroptosis-related gene signatures can robustly diagnose skin cutaneous melanoma and predict the prognosis. Front Oncol. 2021;11:709077. doi:10.3389/fonc.2021.709077

43. Chevrier S, Levine JH, Zanotelli VRT, et al. An immune atlas of clear cell renal cell carcinoma. Cell. 2017;169(4):736–749.e18. doi:10.1016/j.cell.2017.04.016

44. The Cancer Genome Atlas Research Network. Comprehensive molecular characterization of clear cell renal cell carcinoma. Nature. 2013;499(7456):43–49. doi:10.1038/nature12222

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

46. Liston A, Gray DHD. Homeostatic control of regulatory T cell diversity. Nat Rev Immunol. 2014;14(3):154–165. doi:10.1038/nri3605

47. Giraldo NA, Becht E, Vano Y, et al. Tumor-infiltrating and peripheral blood T-cell immunophenotypes predict early relapse in localized clear cell renal cell carcinoma. Clin Cancer Res. 2017;23(15):4416–4428. doi:10.1158/1078-0432.CCR-16-2848

48. Nishida K, Kawashima A, Kanazawa T, et al. Clinical importance of the expression of CD4+CD8+ T cells in renal cell carcinoma. Int Immunol. 2020;32(5):347–357. doi:10.1093/intimm/dxaa004

49. Siddiqui SA, Frigola X, Bonne-Annee S, et al. Tumor-infiltrating Foxp3 CD4 + CD25 + T cells predict poor survival in renal cell carcinoma. Clin Cancer Res. 2007;13(7):2075–2081. doi:10.1158/1078-0432.CCR-06-2139

50. Polimeno M, Napolitano M, Costantini S, et al. Regulatory T cells, interleukin (IL)-6, IL-8, vascular endothelial growth factor (VEGF), CXCL10, CXCL11, epidermal growth factor (EGF) and hepatocyte growth factor (HGF) as surrogate markers of host immunity in patients with renal cell carcinoma: surrogate markers of host immunity in RCC. BJU Int. 2013;112(5):686–696. doi:10.1111/bju.12068

51. Nakanishi H, Miyata Y, Mochizuki Y, et al. Pathological significance and prognostic roles of densities of CD57+ cells, CD68+ cells, and mast cells, and their ratios in clear cell renal cell carcinoma. Hum Pathol. 2018;79:102–108. doi:10.1016/j.humpath.2018.05.007

52. Agrawal R, Garg A, Benny Malgulwar P, Sharma V, Sarkar C, Kulshreshtha R. p53 and miR-210 regulated NeuroD2, a neuronal basic helix-loop-helix transcription factor, is downregulated in glioblastoma patients and functions as a tumor suppressor under hypoxic microenvironment: regulatory and functional dissection of NeuroD2. Int J Cancer. 2018;142(9):1817–1828. doi:10.1002/ijc.31209

53. Zeng Q, Sun S, Li Y, Li X, Li Z, Liang H. Identification of therapeutic targets and prognostic biomarkers among CXC chemokines in the renal cell carcinoma microenvironment. Front Oncol. 2020;9:1555. doi:10.3389/fonc.2019.01555

54. Li X, Shao C, Shi Y, Han W. Lessons learned from the blockade of immune checkpoints in cancer immunotherapy. J Hematol Oncol. 2018;11(1):31. doi:10.1186/s13045-018-0578-4

55. Kennedy LB, Salama AKS. A review of cancer immunotherapy toxicity. CA A Cancer J Clin. 2020;70(2):86–104. doi:10.3322/caac.21596

56. Bhatia A, Kumar Y. Cellular and molecular mechanisms in cancer immune escape: a comprehensive review. Expert Rev Clin Immunol. 2014;10(1):41–62. doi:10.1586/1744666X.2014.865519

57. Vinay DS, Ryan EP, Pawelec G, et al. Immune evasion in cancer: mechanistic basis and therapeutic strategies. Semin Cancer Biol. 2015;35:S185–S198. doi:10.1016/j.semcancer.2015.03.004

58. Yang JC, Hughes M, Kammula U, et al. Ipilimumab (anti-CTLA4 antibody) causes regression of metastatic renal cell cancer associated with enteritis and hypophysitis. J Immunother. 2007;30(8):825–830. doi:10.1097/CJI.0b013e318156e47e

59. Mihara M, Hashizume M, Yoshida H, Suzuki M, Shiina M. IL-6/IL-6 receptor system and its role in physiological and pathological conditions. Clin Sci. 2012;122(4):143–159. doi:10.1042/CS20110340

60. Gudbrandsdottir G, Aarstad HH, Bostad L, et al. Serum levels of the IL-6 family of cytokines predict prognosis in renal cell carcinoma (RCC). Cancer Immunol Immunother. 2021;70(1):19–30. doi:10.1007/s00262-020-02655-z

61. Kang S, Tanaka T, Narazaki M, Kishimoto T. Targeting Interleukin-6 signaling in clinic. Immunity. 2019;50(4):1007–1023. doi:10.1016/j.immuni.2019.03.026

62. Akbari B, Farajnia S, Ahdi Khosroshahi S, et al. Immunotoxins in cancer therapy: review and update. Int Rev Immunol. 2017;36(4):207–219. doi:10.1080/08830185.2017.1284211

63. Ostendorf BN, Bilanovic J, Adaku N, et al. Common germline variants of the human APOE gene modulate melanoma progression and survival. Nat Med. 2020;26(7):1048–1053. doi:10.1038/s41591-020-0879-3

64. Tavazoie MF, Pollack I, Tanqueco R, et al. LXR/ApoE activation restricts innate immune suppression in cancer. Cell. 2018;172(4):825–840.e18. doi:10.1016/j.cell.2017.12.026

Read more here: Source link