Clear Cell Renal Cell Carcinoma: A Comprehensive in silico Study in Searching for Therapeutic Targets

Renal cell carcinoma (RCC) is one of the world’s leading causes of death and morbidity. With 431,288 new cases and 179,368 deaths in 2020, this cancer is considered an urgent global problem [1]. So, RCC is an evolving health issue requiring more effective preventive/therapeutic actions. About 70–75% of RCC types are clear cell RCC (ccRCC) that are recognized by their clear appearance in pathology [2]. ccRCC cases are usually resistant to radiotherapy and chemotherapy; therefore, surgery is still the preferred treatment option for patients [3]. However, this approach is only effective in the early stages of cancer (when it is localized). In metastatic forms (systemic spreading of ccRCC), targeted therapies such as tyrosine kinase inhibitors (TKIs), VEGF inhibitors, or mTOR inhibitors are usually prescribed. In addition, immune checkpoint inhibitors were among the recently included therapeutics in the FDA-approved ccRCC treatment list [4]. Although such treatments can improve survival rates in patients, drug responses are not satisfactory, and drug resistance develops in most patients [5]. Therefore, more investigations are required to have a deeper understanding of ccRCC pathogenesis and discover novel, druggable key elements in its pathogenesis.

Although various genetic events like 3p chromosome arm loss, VHL gene mutation, STED2 mutation, and KDM5C mutation have been shown to be linked with the pathogenesis of ccRCC, the exact underlying molecular mechanisms of this cancer are not yet fully understood [3]. In this context, systems biology approaches could be beneficial in dealing with high-throughput data from cancer patients and exploring the complexity behind such convoluted pathological phenomena [6].

Recently, a rising number of researchers have examined high-throughput datasets to understand better the biological roots of various diseases [7, 8]. Analyzing and interpreting transcriptomics datasets with different bioinformatics approaches could provide a holistic view of expressional changes in cancer versus normal states [9].

In this context, the Weighted Gene Co-Expression Network Analysis (WGCNA) algorithm is one of the well-known clustering tools for identifying co-expressed and disease-correlated genes in expression datasets [10‒12]. So far, many studies have used the WGCNA to investigate various biological processes in conditions like genetic disorders, cancers, and other chronic diseases like nephropathy and Alzheimer’s, aiming to discover new biomarkers or therapeutic targets [13, 14].

This study aimed to reveal the cancer-related underlying molecular mechanisms and introduce the principal hub genes in ccRCC pathogenesis by analyzing available ccRCC microarray datasets using network analysis and the WGCNA algorithm. After analyzing the ccRCC-related datasets, the intensity values of the discovered differentially expressed genes (DEGs) were passed to WGCNA in the R environment. Then, by identifying the disease-correlated module, the classified DEGs in the selected module were subjected to functional enrichment analyses. The hub genes were spotted by considering the top genes in the protein-protein interaction (PPI) network, including all the DEGs, and a list ordered by module membership (kME) values of genes in the selected module. The hub genes also were considered for further examinations like expressional validations and drug repositioning. The workflow of the present study is shown in Figure 1.

Gene Expression Omnibus [15] was used to search and download the RCC microarray datasets. Prior to the primary analysis step, principal component analysis was accomplished to distinguish and eliminate probable outliers. Normalization, removal of multiple gene-related probes, and analysis were performed with the help of the Network Analyst online tool [16]. Linear Model for Microarray Analysis (Limma) was selected for the analysis procedure, and a false discovery rate (FDR) cutoff <0.049 was considered to identify the significant DEGs [17, 18].

WGCNA algorithm was applied to identify the co-expression modules among the intensity values of the identified DEGs. A matrix containing the intensity files of the identified DEGs was introduced to the WGCNA. Next to the sample clustering, the best soft threshold power was selected using the “pickSoftThreshold” function in the package. The best soft threshold power (value = 7) was selected based on the mean connectivity and degree of independence values. Afterward, several steps, including adjacency matrix construction, module identification (minModuleSize = 30), topological overlap matrix calculation, construction of modules, and finally, dynamic branch cutting (merging threshold: 0.25), were conducted independently to reach the co-expression modules. Moreover, the most correlated module was identified after assessing the correlation between two states (cancer vs. normal) and the co-expression modules [10].

Functional enrichment analyses, including gene ontology and pathway enrichment, were performed for the extracted DEGs from the module of interest. These analyses (significant enrichment threshold: p < 0.05) were performed using Cytoscape software (version 3.9.1) and the CluGO module (version 2.5.7) [19, 20].

Candidate hub gene identification was performed by considering kME and DEG network centrality (degree) values in a constructed PPI network. In detail, after PPI network construction using all DEGs, the top 100 genes with a high degree of connectivity were listed. On the other hand, the DEGs with kME values of more than 0.8 in the module were extracted and recorded. Lastly, hub genes were selected as the common ones in both lists. Contraction of the PPI network was performed using the STRING database (confidence >0.4), and Cytoscape (v.3.9.1) was utilized for network visualizations and analysis [20, 21].

Gene Expression Profiling Interactive Analysis (GEPIA) [22] was applied for the expression analysis of hub genes [23]. The expression levels of the hub genes in 100 normal kidney tissues and 523 ccRCC tissues were plotted to compare the expression levels between the two states. The log2 transformed transcripts per kilobase million of RNA seq data were used for making the box plots. GEPIA examined the correlation of expressions of hub genes with pathological cancer stages. The log2 transformed transcripts per kilobase million values were used for plotting stage plots. The considered threshold for this analysis was p (>F) < 0.05.

Gene Set Cancer Analysis (GSCA) server was utilized to assess the expression patterns and values of the hub genes in ccRCC and other cancers [24]. The GSCA performs differential expression analysis of paired normal and tumor tissues from TCGA based on normalized RNA seq by expectation-maximization mRNA expression. The data from 11 to 114 paired tumor and normal samples were utilized in the analysis for each cancer type. T test was used to estimate the p value, and FDR was used as the correction method.

To identify potential drugs for the treatment of ccRCC, we performed a drug-gene interaction assessment using Drug-Gene Interaction Database (DGIdb) [25]. DGIdb is an open-source project that freely provides gene-drug interactions from 22 trusted sources [26]. The tool was applied to obtain hub gene-drug interactions. Then, a network of the potential drugs and targeted hub genes was constructed by Cytoscape (version 3.9.1) [20].

Spearman correlation between hub gene GSVA scores and different immune cells was exported as a heatmap using the GSCA server based on the ImmuCellAI method and the Cancer Genome Atlas ccRCC transcriptomics information (FDR ≤0.05) [24, 27].

All the steps of the present work are shown in Figure 1. Three RCC datasets were selected and retrieved from the GEO database, including GSE11151, GSE12606, and GSE36897. Due to the same platform, all three datasets were merged to have a meta dataset, including 58 and 32 cancer and control samples. Principal component analysis, as a valuable tool to check the quality of datasets, was performed before the analysis [17] (Fig. 2a). Likewise, the quantile normalization was performed to ensure the similarity of the expression distributions of each sample across the entire meta dataset. After analysis, 7,411 significant DEGs, including 4,179 downregulated and 3,232 upregulated genes, were detected (FDR cutoff <0.049) and subjected to further analysis. The volcano plot of the analyzed dataset and a heatmap for the top 50 DEGs based on FDR are shown in Figure 2b, c. The results of the DEG analysis are provided in sheet 1 of online supplementary file 1 (for all online suppl. material, see www.karger.com/doi/10.1159/000529861).

Co-expression module identification was carried out using the WGCNA algorithm in the R software environment. Using this package, we performed sample clustering, constructed co-expression networks, and finally recognized the highly correlated module to the RCC state [10]. Sample clustering was performed, and each sample was dedicated to either disease or control state (Fig. 3a). The resulting scale independence and mean connectivity platforms are shown in Figure 3b. According to the resulting platforms, the best soft threshold (soft threshold: 7) was designated to get an approximate scale-free topology.

Further on, after performing hierarchical clustering and module merging, the DEGs were classified into eight different co-expression modules, including blue, black, brown, pink, purple, magenta, red, and salmon (Fig. 4a, b, d). The number and names of the classified DEGs in each module are listed in sheet 2 of online supplementary file 1. The constructed heatmap plot (Fig. 4c) revealed the accuracy of the module division and the topological overlap adjacency among genes in modules.

The correlations between modules and phenotypes (cancer and normal) were checked, and the pink module with 2,441 DEGs was recognized as the most disease-relevant module (r = 0.84; p = 2E-25) (Fig. 4e). Other modules with a positive correlation to the disease state included magenta (r = 0.53; p = 7E-8), black (r = 0.58; p = 2E-9), purple (r = 0.56; p = 8E-9), and red (r = 0.62; p = 9E-11). Other modules including salmon (r = −0.57; p = 5E-9), blue (r = −0.91; p = 3E-35), and brown (r = −0.72; p = 2E-15) had negative correlations with the cancer state. A scatter plot of the gene significance (GS) versus kME values of all DEGs in the pink module is shown in Figure 4f. The high correlation between the values of GS and kME in the pink module (r = 0.9, p < 1e-200), proposing the association of genes with both the module eigengene and the disease state.

The results of functional enrichment analyses of the DEGs in the pink module are shown in Table 1 (top 10 terms). Based on the results of Reactome pathway analysis, the DEGs were primarily enriched in “translocation of ZAP-70 to immunological synapse,” “endosomal/vacuolar pathway,” “cell surface interactions at the vascular wall,” and “immune-related pathways.” The enriched terms in the biological process section also confirmed the results of Reactome pathway enrichment, where the primarily enriched terms were related to immunological processes like “regulation of neutrophil activation,” “interleukins, and interferon-mediated signaling pathways,” and “Fc receptor-mediated inhibitory signaling pathway.” In terms of molecular function, the top enriched terms were “transporter associated with antigen processing binding,” “T cell receptor binding,” and “MHC and cytokine receptors activities.” In addition, the top cellular component terms were “complement component C1 complex,” “MHC protein complex,” “NLRP3 inflammasome complex,” “integral component of luminal side of endoplasmic reticulum membrane,” and “phagocytic cup.”

A PPI network with 6,802 nodes and 146,523 edges was constructed, including almost all the classified DEGs in the eight modules. The network is uploaded to the NEDx server and is accessible through the link: https://www.ndexbio.org/#/network/50894774-9c60-11ec-b3be-0ac135e8bacf?accesskey=5fda0e58489fd8a029aad792b07c7200d027a70439915453c3a0d1d6a7d7d431.

The constructed network was analyzed, and the top 100 DEGs were identified based on the degree of centrality. On the other hand, the top DEGs based on kME and GS values (171 DEGs) were identified and listed among the DEGs in the pink module (Table 2). Common DEGs of two lists were 12 genes named hubs in this study. The hubs included “protein tyrosine phosphatase receptor type C (PTPRC),” “integrin subunit alpha M (ITGAM),” “Toll-like receptor 2 (TLR2),” “CD86 molecule (CD86),” “pleckstrin (PLEK),” “TYRO protein tyrosine kinase binding protein (TYROBP),” “integrin subunit beta 2 (ITGB2),” “Rac family small GTPase 2 (RAC2),” “colony-stimulating factor 1 receptor (CSF1R),” “C-C motif chemokine receptor 5 (CCR5),” “C-C motif chemokine ligand 5 (CCL5),” “lymphocyte cytosolic protein 2 (LCP2).”

Based on the expression analysis results of 100 normal kidneys and 523 ccRCC samples by GEPIA, all 12 hub genes showed upregulation patterns in ccRCC samples compared to normal samples (Fig. 5a). Moreover, in the stage analysis step by GEPIA, 4 of 12 hub genes, including CCR5, RAC2, TYROBP, and CCL5, showed significant changes during different ccRCC stages (Fig. 5b). These findings confirmed the present study’s expressional state of the hub DEGs.

GSCA is an integrated database for gene set analysis in cancers. Using this database, we compared the differential expression values of hub genes with other cancers, and the results showed an exclusive alteration pattern in the expressional values of the hub genes in ccRCC samples (Fig. 6a).

In this step, the GSCA database was utilized to obtain the type and patterns of immune cell infiltration in ccRCC. The result of this step is shown in Figure 6b. Based on the findings, there was a positive correlation between the expression of the 12 hub genes and “T helper 1 cells,” “macrophage cells,” “induced Regulatory T cells (iTreg),” “exhausted T cells,” “type 1 regulatory cells,” “CD8-positive T cells,” “central memory T cells,” “T follicular helper cells (Tfh),” natural Treg, dendritic cells, cytotoxic cells, and “CD4 positive T cells.” On the other hand, a negative correlation was observed between the expression of hub genes and “neutrophil cells,” “T helper 17,” “CD8_naive cells,” “T helper 2 cells,” “CD4_naive cells,” “effector memory cells,” and “monocyte immune cells.” In addition, the hub genes positively correlated with the general immune cell infiltration score.

The DGIdb database was searched using the 12 hub DEGs to identify the hubs’ most related drugs. In this step, 60, 15, 18, 10, 12, 4, 6, and 2 drugs were found for CSF1R, CCR5, ITGB2, ITGAM, PTPRC, TLR2, CD86, and RAC2, respectively. A network comprising the hub genes and their related drug molecules is shown in Figure 7.

Biological systems are typically shown as networks, which are collections of binary interactions between various components [28]. The nodes of such networks might be the biological entities that are differently expressed in a situation like cancer (for example, genes, proteins, metabolites, or noncoding RNAs). In this situation, creating, analyzing, and interpreting such networks may improve our comprehension of cancer biology and point us in the direction of the primary pathogens involved in the disease [12, 29]. Numerous methods have been created and used to analyze biological networks [30, 31]. One of the well-known methods for co-expressional network analysis of high-throughput expressional data, for instance, is WGCNA [10]. The program may group the co-expressed genes into different modules, find the disease most associated module, and then identify the key genes in that module. So far, this method has been applied to several different illness situations, introducing multiple possible treatment targets [11, 14, 32]. Another popular method for locating the important nodes (potential treatment targets) in diseased networks is the topological analysis of PPI networks [33]. The centrality characteristics of nodes, such as degree, betweenness, proximity, eigenvector, etc., are often taken into account in this method to identify the important nodes in PPI networks [34]. In this strategy, due to their higher number of interactions with other nodes in the network, the top centrality-based nodes are considered potential therapeutic targets in this respect.

In this integrative in silico study, the selection of hubs was based on a novel strategy considering the extracted information by analyzing both the PPI and the co-expression networks of DEGs in ccRCC. The 12 hub genes were identified as essential genes in the pathogenesis of ccRCC. Notably, some of the identified hub genes in the present study were previously recognized as key players in the pathogenesis of ccRCC [35‒39]. A brief description of the identified hub genes is provided in online supplementary file 2. The identified hub genes had high interaction degrees in the PPI network and high module membership values in the disease most correlated co-expressed module. All 12 hub genes were upregulated in ccRCC versus normal conditions, and the expressional validation results using the TCGA database also were in line with our analysis.

Moreover, the expressional pattern of the hub genes was different in other cancers, which might indicate the exclusive role of these 12 genes in ccRCC pathogenesis (Fig. 6a). In addition, based on the stage analysis results, 4 of the hub genes showed increasing patterns in their expression in different stages of ccRCC. Accordingly, these 4 genes, including TYROBP, CCL5, CCR5, and RAC2, might have a part in the progression of ccRCC, and as a biomarker, the expression levels of these genes can indicate the stages of the disease.

Based on gene ontology results, immunological pathways were the most enriched terms for the DEGs in the disease most correlated module. Also, all the identified hub genes were involved in immune response processes. Such findings were in line with previous histopathologic and transcriptomic studies in which ccRCC was recognized as one of the most immune cell infiltrated neoplasms [40, 41]. Based on our results, the expression of the identified hub genes is associated with increased Treg and decreased Th17 and CD8+ T immune cell infiltration in the tumor microenvironment. According to some experiments, infiltration of these cells predicts poor prognosis for ccRCC patients [42]. Accordingly, it is suggested that the identified hub genes are involved in the ccRCC immune microenvironment dynamics. Specifically, T cell immune checkpoint inhibitors have created a new promising era in the treatment of ccRCC, and current checkpoint inhibitors (alone or in combination with TKIs) are becoming the standard treatment for metastatic ccRCC [41]. In addition, our gene ontology analysis showed that the innate immune responses (like IFN pathways, regulation of neutrophil activation and degranulation, antigen processing and presentation, Toll-like receptor cascades, and Fc receptor signaling) have a principal value in the ccRCC pathogenesis. Such findings are consistent with other studies indicating the importance of innate immune cells in ccRCC immunology, such as myeloid-derived suppressor cells [40, 43]. Therefore, the outlook of ccRCC immunotherapy might need to pass through from only focusing on the adaptive immune system.

In the next step of this work, we focused on drug-gene interaction identification to propose candidates for more evaluations in drug repurposing studies. Based on the results, 127 drugs were identified to interact with the 8 hub genes. The results of this part indicated the importance and therapeutic potential of the identified hub genes in ccRCC. Interestingly, some of the identified drugs (sunitinib, pazopanib, sorafenib) that target CSF1R are currently applied for ccRCC treatment as TKIs [4]. Notably, CSF1R inhibition has recently gained more attention for targeting immune suppression in solid tumors [44]. Furthermore, ipilimumab (anti-CTLA4), which is responsible for suppressing the interaction of CD28 and CD86, is another approved drug for ccRCC immunotherapy [41].

Among the predicted drugs, there were also some immunosuppressive ones like hydrocortisone, prednisone, dexamethasone, methylprednisolone (glucocorticoids), adalimumab, infliximab, etanercept (anti-TNF alpha antibodies), and indomethacin (NSAIDs) that target CD86, ITGB2, ITGAM, PTPRC, and CD86. Interestingly, there are some clinical case reports of spontaneous remission following glucocorticoid therapy in ccRCC patients [45‒48]. Despite such reports, we believe glucocorticoids could be potential drug candidates for the treatment of ccRCC [49]. Likewise, different in vitro and in vivo studies have shown the inhibitory effects of these drugs on ccRCC [50, 51]. Although immunosuppressive, like several approved TKIs for the treatment of ccRCC, glucocorticoids have antiangiogenesis functions [52, 53]. In addition, even though some TKIs have immunomodulatory properties [54], different studies have shown their ability to reprogram the tumor immune microenvironment and promote tumor immunity [43]. Also, some evidence indicates the antiproliferative/immunostimulatory effect of anti-TNF alpha antibodies and indomethacin in ccRCC cases [55‒58].

Other identified drugs were chemotaxis-targeting molecules, including rovelizumab, which is known as a safe and tolerable anti-CD18/CD11b (ITGB2/ITGAM), and maraviroc that target the CCR5/CCL5 axis [59, 60]. It was shown that the CD18/CD11b complex on macrophages and neutrophils has a role in the recruitment of myeloid-derived suppressor cells to the tumor microenvironment. Currently, the efficiency of several CD18/CD11b inhibitors for solid tumor therapy is under evaluation [59]. Likewise, the role of the CCR5/CCL5 axis in immunosuppressive cell recruitment in tumor progression was shown by previous studies. As an inhibitor of this axis, maravirocs approved for treating HIV-1 infection [60‒62]. Resveratrol hexanoic acid, a resveratrol derivate with improved pharmacokinetic functions, was another identified drug for the hub genes in ccRCC condition (targeting TLR2) [63]. Several studies have indicated this FDA-approved antioxidant agent’s immunostimulatory and tumor-suppressing effect in ccRCC. Therefore, considering the excellent safety profile of resveratrol, it can be a good target for further investigations in the treatment of the ccRCC [64‒67].

By utilizing an integrative bioinformatics approach, this in silico study tried to shed light on the underlying mechanisms in the pathogenesis of ccRCC. As far as we know, there is no such study considering the information in the PPI network and the co-expression network of DEGs in ccRCC to identify key molecules in this cancer. By applying this strategy, 12 genes were introduced as therapeutic target candidates to hinder the progression of ccRCC. Most of the introduced hubs were previously identified as key players in ccRCC, which could confirm the applied integrative in silico strategy in this study. Moreover, the identified hubs could be the target of future investigations to explore their exact roles in the progression of this cancer. Recently, single-cell RNA sequencing has been recognized as an excellent solution for unraveling the complex aspect of the tumor microenvironment. Therefore, our next step is to identify the gene expression profiles of these 12 genes in different ccRCC tumor microenvironment cell populations, especially immune cells.

We thank Regenerative Medicine Research Center members for their help in some parts of the bioinformatics analysis steps.

An ethics statement was not required for this study type, and no human or animal subjects or materials were used. However, this experiment’s protocol was reviewed by the Ethics Committee of the Isfahan University of Medical Sciences.

The authors declare that they have no competing interests.

This work was supported by the Isfahan University of Medical Sciences [Grant No. 2400109].

Amir Roointan participated in the study design, dataset analysis using the WGCNA algorithm, data interpretation, and manuscript drafting. Mohammadjavad Naghdibadi participated in the design, dataset analysis using the Limma algorithm, data interpretation, and manuscript drafting. Maryam Momeni participated in further exploring the hub genes, including expressional validation, and also participated in manuscript drafting. Parvin Yavari participated in dataset analysis, hub gene selection strategy, immune cell infiltration prediction, drug-gene interaction analysis steps, and manuscript drafting. Alieh Gholaminejad contributed to dataset selection, interpretation of the obtained results, and drafting and editing of the manuscript. All authors participated in preparing and reviewing the draft of the manuscript.

The analyzed datasets by the present study are available in the GEO repository. Further inquiries can be directed to the corresponding author Amir Roointan ([email protected]).

[https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE11151]

[https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE12606]

[https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36897]