Prediction and analysis of novel key genes ITGAX, LAPTM5, SERPINE1 in clear cell renal cell carcinoma through bioinformatics analysis

Data collection

Four gene expression profiles ; (GSE12606, GSE53000, GSE68417, and GSE66272) were downloaded from the GEO (Barrett et al., 2013) database (http://www.ncbi.nlm. nih.gov/geo/). The TCGA KIRC data was downloaded from TCGA (https://www.cancer.gov/about-nci /organization /ccg/research/structural-genomics/tcga). The GSE12606 (Stickel et al., 2009) data set was obtained by GPL570 platform, including three adjacent normal kidney specimens, and three CCRCC samples. The GSE53000 (Gerlinger et al., 2014) and GSE68417 (Thibodeau et al., 2016) data sets were based on GPL6244 Platforms. GSE53000 data set was comprised of 60 samples including 6 adjacent normal kidney specimens and 54 CCRCC samples and GSE68417 data set was comprised of 43 samples including 14 adjacent normal kidney specimens and 29 CCRCC samples. GSE66272 (Wotschofsky et al., 2016; Liep et al., 2016) data set was obtained by GPL5029 platform and contained 27 CCRCC tumor samples and 26 normal kidney samples. The TCGA KIRC data contained 539 CCRCC tumor samples and 72 normal kidney samples.

Study design and data processing

In order to clarify our study, we designed a flow chart to demonstrate data collection, processing, analysis and verification (Fig. 1). The online tool GEO2R (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga) is an online analysis of GEO series based on the R programming language, which is used to screen DEGs between normal kidney and CCRCC samples from the GSE data sets. The data was standardized and filtered to select significant DEGs. P-value < 0.05, log Fold Change (—log FC—) ≥1 were chosen as the cutoff criteria. log Fold Change (— log FC —) ≥ 1 means that the multiple of change is greater than or equal to 2. It is generally considered that there is a difference between 2 times and more. Then, DEGs were further screened according to cutoff criteria: P-value < 0.05, log FC ≥1 as up-regulated genes, P-value < 0.05, log FC ≤1 as down-regulated genes. Finally, importing all the up-regulated genes or down-regulated genes in the 4 datasets into Funrich 3.1.1 software, and taking the intersection of the up-regulated genes or down-regulated genes in the 4 datasets respectively. TCGA RNA-seq simple converter was used to standardize and log2 conversion of TCGA KIRC data. TCGA KIRC data were used for expression verification and ROC curve analysis of Key genes.

Enrichment analysis of DEGs

DAVID 6.8 (Database for Annotation, Visualization and Integrated Discovery, https://david.ncifcrf.gov/)  (Huang da, Sherman & Lempicki, 2009b; Huang da, Sherman & Lempicki, 2009a) database was used to analyze Gene Ontology (GO) (The Gene Ontology Consortium, 2017) such as the biological process (BP) (Berchtold, Csaba & Zimmer, 2017), cellular component (CC) (Borg & Baudino, 2011) andmolecular function (MF), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa & Goto, 2000) pathway enrichment analysis of important DEGs, which promoted the visualization of gene and protein function (Dennis Jr et al., 2003). Among them, the count value represents the number of genes enriched in the pathway. The cutoff value was P < 0.05.

Protein-Protein Interaction (PPI) Network construction and hub gene screening

The online tool STRING (Von Mering et al., 2003) (http://string-db.org) was utilized for the analysis of protein-protein interactions . Importing all DEGs (including up-regulated genes and down-regulated genes) into the STRING database for analysis, and the confidence level ≥0.4 is considered to be significant for PPI. STRING and Cytoscape software (version 3.6.1) were used to construct the PPI network, which is regarded as a network software platform for visualizing protein molecular interactions. Using the degree algorithm in the cytohubba plug-in of Cytoscape, the degree value could be calculated and thus the hub genes could also be screened from PPI network. The importance of genes is directly related to the degree of protein. According to the degree value (degree ≥37), the first 35 genes were selected as hub genes. Three key genes ITGAX, LAPTM5 and SERPINE1with high degree were screened through survival and prognosis analysis as well as the novelty of genes.

Key genes analysis and verification

Three key genes were analyzed and verified comprehensively by using TCGA, GEPIA, Oncomine, Kaplan–Meier plotter and UALCAN databases. Similarly, TCGA database also applied to ROC curve analysis. ROC curve was plotted by R software, which could be used to analyze the diagnostic value of key genes in CCRCC. GEPIA (http://gepia.cancer-pku.cn/) analyzed the expression of key genes in CCRCC. Oncomine (Rhodes et al., 2007) (http://www.on comine.org) was a database consisting of microarray data of various tumors, which verified the expression of key genes in CCRCC. Kaplan-Meier Plotter (https://kmplot.com/analysis/) was used to analyze the survival of key genes in CCRCC. UALCAN (Chandrashekar et al., 2017) (http://ualcan.path. uab.edu/index.html) online database was used to analyze the relationship between key genes and Clinical pathology in CCRCC, and the promoter methylation level of key genes in CCRCC was analyzed.

Statistical analysis

Statistically significant differences between the normal tissues group and tumor tissues group were determined using Student’s t tests. Also, Kaplan–Meier analysis was used to assess OS. ROC curve is a graphical plot (Cao & López-de Ullibarri, 2019) that reflects the sensitivity and specificity of continuous variables. P < 0.05 was considered to be statistically significant for all the tests. All the statistical analyses applied GraphPad prism 6.0 or R software.

Volcano plots of the differentially expressed genes in four datasets

The data sets of CCRCC, which were GSE12606, GSE53000, GSE68417 and GSE66272, were downloaded from the GEO database (Table 1) and analyzed by GEO2R separately. A total of 9,560 DEGs were screened from the GSE12606 data sets, among which 1568 genes were up-regulated and 5436 genes were down-regulated instead. A total of 1,286 DEGs were screened from the GSE53000 data sets, among which 583 were up-regulated and 708 were down-regulated. A total of 1,890 DEGs were screened from the GSE68417 data sets among them 726 up-regulated and 1,164 down-regulated genes were selected. There were 3,627 up-regulated genes and 5,170 down-regulated genes among the 8,797 DEGs screened from the GSE66272 data sets. The screening criteria were P-value 0.05, log Fold Change (— log FC —) ≥ 1. All of the DEGs from the four data sets were presented in Volcano plots (Figs. 2A–2D). Among them, red represents high-expressed genes, green represents low-expressed genes, and black represents genes whose expression levels are not significant in each data set. We gained the intersection of four independent data sets through Funrich 3.1.1 software and made a visualization analysis of Venn Diagram. After being overlapped, the common 738 genes (Fig. S1) were identified, including 289 up-regulated and 449 down-regulated genes (Figs. 2E–2F).

GO and KEGG enrichment analysis of DEGs

To further explore the function of DEGs, enrichment analysis of up-regulated genes and down-regulated genes were displayed respectively. DAVID 6.8 was used to perform GO and KEGG analysis of DEGs in CCRCC (Table 2). In biological processes, up-regulated DEGs are mostly involved in cell adhesion, signal transduction, immune response, especially the regulation of inflammatory response (Fig. 3A); while down-regulated DEGs are mostly involved in oxidation–reduction process proteolysis, ion transmembrane transport, response to drug, and ion transport (Fig. 3E). In terms of cellular components, up-regulated DEGs are mainly distributed in the plasma membrane, integral component of the membrane, extracellular exosome, and membrane (Fig. 3B); down-regulated DEGs are mainly distributed in integral component of membrane, plasma membrane, cellular exosomes, and integral component of plasma membrane (Fig. 3F). In terms of molecular function, up-regulated DEGs generally have protein binding, identical protein bind, ATP binding capacity, and protein homodimerization activity (Fig. 3C); down-regulated DEGs generally have identical protein homodimerization activity, calcium ion binding, sequence-specific DNA binding, and oxidoreductase activity (Fig. 3G). In the KEGG signal pathway, the up-regulation of DEGs mostly involved in the PI3K-Akt signaling pathway, pathways in cancer, focal adhesion, and HIF-1 signaling pathway (Fig. 3D); while the down-regulation of DEGs mainly involved in the metabolic pathway, Biosynthesis of antibiotics, Carbon metabolism, and Aldosterone-regulated sodium reabsorption (Fig. 3H). CCRCC is a kind of metabolic disease. Metabolic reprogramming covers different processes including aerobic glycolysis, fatty acid metabolism and the utilization of tryptophan, glutamine as well as arginine (Lucarelli et al., 2019), which has also been proved by the results of KEGG pathway enrichment analysis in the work. KEGG pathway results show that down-regulated DEGs are enriched in metabolic pathways, Glycolysis/Gluconeogenesis, Glycine, serine and threonine metabolism and so on.

In addition, GO and KEGG enrichment analyses were performed for all DEG (up-regulated and down-regulated genes) (Fig. S2). Finally, 11 common enrichment results were found from the three analyses (including up-regulated DEG, down-regulated DEG and all DEG enrichment analyses). The common enrichment results are as follows: response to drug, plasma membrane, integral component of membrane, extracellular space, extracellular exosome, cell surface, cell surface, protein homodimerization activity, identical protein binding. The results of enrichment analysis were consistent with previous studies (Wang, Yu & Chai, 2019; Tian et al., 2019). It suggests that these 11 pathways may be important in CCRCC.

Construction of PPI network and screening of key genes

To identify the key genes, the STRING online database and Cytoscape software were used to analyze all DEGs (including up-regulated genes and down-regulated genes)and construct PPI network (Fig. 4A). Based on the main role of proteins in biological functions, their interaction determines the molecular and cellular mechanisms that control the health and disease state of the organism (Safari-Alighiarloo et al., 2014). Next, using the degree algorithm in the cytohubba plug-in of Cytoscape to screen the hub genes in PPI network. The gene whose degree ≥37 could be defined hub gene, therefore the first 35 genes (Table 3) in the PPI network were chosen to be hub genes (Fig. 4D). Next, GO enrichment analysis was conducted for all DEGs and the first 35 DEGs in the PPI network. The result showed that all DEGs mainly enriched in signal transduction, oxidation–reduction process, cell adhesion, inflammatory response, plasma membrane, integral component of membrane, extracellular exosome, integral component of plasma membrane (Figs. 4B–4C). The top 35 DEGs were enriched in extracellular space, extracellular exosome, inflammatory response, an integral component of the plasma membrane, and cell surface (Fig. 4E). By comparison, it turned out that the enrichment analysis results of all DEGs in PPI network contained the results of the top 35’s. In the end, three new key genes with high degree which were ITGAX, LAPTM5 and SERPINE1, could be screened by using GEPIA and Kaplan–Meier plotter database.

The expression of key genes

Among the 35 genes, we focused on ITGAX, LAPTM5, and SERPINE1, which have not been reported to be related to the occurrence and development of CCRCC. Firstly, the expression levels of ITGAX, LAPTM5, and SERPINE1 in CCRCC tumor tissues are significantly higher than those in normal tissues adjacent to cancer according to GEPIA database (Figs. 5A–5C). Furthermore, the data of TCGA KIRC showed that compared with adjacent normal tissues, the mRNA expressions of ITGAX, LAPTM5, and SERPINE1 in 72 pairs of CCRCC tissues were significantly increased (Figs. 5D–5F). In the Oncomine Gumz renal database, the mRNA levels of ITGAX, LAPTM5 and SERPINE1 were also upregulated in CCRCC tissues when compared with adjacent normal kidney tissues (Figs. 5H–5G). According to the UALCAN database, the promoter methylation levels of ITGAX, LAPTM5 and SERPINE were decreased in CCRCC (Figs. 5J–5L). To sum up, according to the GEPIA database, Oncomine Gumz renal database and TCGA database, the expression levels of these three genes in CCRCC tumor tissues are significantly higher than those in normal tissues adjacent to cancer. It could be further speculated that the high expression of ITGAX, LAPTM5, and SERPINE1 in CCRCC tumor tissue might be related to the decrease of promoter methylation.

The association of key genes expression with clinical pathology in CCRCC

Furthermore, the relationship between the mRNA expression of ITGAX, LAPTM5 and SERPINE1 and different clinical pathology grades were measured. The results showed that their mRNA expression was significantly related to pathological grades (Figs. 6A–6C). And the expression of ITGAX, LAPTM5, and SERPINE1 mRNA in CCRCC samples are also significantly correlated with severe clinical staging (Figs. 6D–6F). Among them, the expression levels of ITGAX, LAPTM5, and SERPINE1 were higher in stage 4 and grade 4. In conclusion, ITGAX, LAPTM5, and SERPINE1 are significantly associated with clinical pathology.

Survival and diagnostic value of ITGAX, LAPTM5, and SERPINE1 in CCRCC

According to the Kaplan–Meier plotter database, the overall survival of ITGAX, LAPTM5, and SERPINE1 genes was tested (Figs. 7A–7C). The results showed that the high expression of three key genes in CCRCC was negatively correlated with prognosis. Then, the ROC curve was used to evaluate the difference between CCRCC and the normal tissues in the TCGA KIRC data. The ROC curves (Figs. 7D–7E) of these three genes showed that their area under the curve (AUC) and the 95% confidence intervals (CI) are as follows: ITGAX (AUC = 96.315, CI = 0.9423–0.9840), LAPTM5, (AUC = 94.808 CI = 0.9243–0.9719), and SERPINE1(AUC = 76.272, CI [0.7068–0.8186]). Since ROC curves had good specificity and sensitivity, ITGAX, LAPTM5, and SERPINE1 had excellent diagnostic efficiency for distinguishing tumors and normal tissues.