Using a Machine Learning Approach to Identify Key Biomarkers for Renal Clear Cell Carcinoma

Introduction

The incidence of kidney cancer is increasing year by year, and it has reached 0.149% (95% CI: 0.144–0.153%) in 2018. The 5-year conditional relative survival rates were only 75.8% (95% CI: 75.5–76%), indicating a high mortality rate from this disease (SEER http://seer.cancer.gov/statfacts/html/kidrp.html). The most common and deadly subtype of renal carcinoma is KIRC, which accounts for approximately 75% of renal carcinoma.^1,2 The main cause of death in patients with KIRC is tumor metastasis.³ There are no obvious clinical features in the early stage of kidney cancer, and 25–30% of patients have already metastasized when they are first diagnosed.^4–6 KIRC patients whose local tumors have been removed by nephrectomy are still at high risk of metastasis and recurrence and are not sensitive to chemotherapy and radiotherapy, leading to poor prognosis.^7–9 Therefore, early diagnosis and treatment of this disease are very important. Besides, understanding key genetic drivers for progression can help to develop new treatments.¹⁰

The rapid development of bioinformatics helps to screen out the characteristic genes of diseases using data from public databases.^11–13 However, few studies have used two well-known machine learning algorithms at the same time, LASSO model and SVM-RFE algorithms, to identify biomarkers of KIRC. LASSO analysis, as regression analysis method, was to perform variable selection and regularization while fitting a generalized linear model, the variable selection was to selectively put variables into the model to get better performance parameters and selectively put variables into the model to get better performance parameters, regularization refers to controlling the complexity of the model through a series of parameters λ to avoid overfitting.^14,15 SVM was a generalized linear classifier (generalized linear classifier) that performed binary classification of data in a supervised learning manner, which maximized the margin between different classes in a high-dimensional space.^13,16 Here, the two machine learning algorithms were used to recognize diagnostic biomarkers data of KIRC and predicted the accurately and interpretable models.

Herein, four open access KIRC-related patient datasets from the GEO database and one dataset from TCGA dataset were downloaded. DEG screening and GO, KEGG, DO and GSEA enrichment analysis was firstly conducted and then LASSO and SVM-RFE algorithms were adopted to identify KIRC-associated key genes in training sets and validate them in the test set. The clinical prognostic analysis including the association between the expression of key genes and the overall survival, stage, grade across KIRC, the immune infiltration difference between normal samples and cancer samples, the correlation between the key genes and immune cells, immunomodulator, immune subtypes of KIRC were investigated in this research, which could provide new guidance for cancer diagnosis.

Materials and Methods

Data Collection and Download

The data sets were collected using the search term “clear cell renal cell carcinoma” as the keyword, with the search scope limited to “Homo sapiens” in GEO (Gene Expression Omnibus) database (http://www.ncbi.nlm.nih.gov/geo/),¹⁷ including KIRC samples and normal samples. After exclusion and screening, GSE6344, GSE40435, GSE781, GSE53757 datasets were included in the study. Besides, for screening biomarkers of early stages of patients with KIRC, another dataset of KIRC was obtained from TCGA datasets (The Cancer Genome Atlas, https://portal.gdc.cancer.gov/) called TCGA_KIRC. Among them, samples of patient with KIRC in early stage (including stage I and II) have been preserved and both normal and cancer samples of patients with KIRC in advanced stages (including stage III and IV) were removed. The basic clinical information has been collected. Log2 in the R software was used to correct and normalize the original data of the datasets.^18–21

Data Preprocessing and DEG Screening

The dataset format from TCGA was adjusted from fpkm to tpm which was more similar to the format of GEO chip data. The original datasets were corrected, including probe matrix and human genome array. The expression matrix was extracted using perl software and the expression matrix of GSE6344, GSE40435, GSE781, TCGA_KIRC datasets was merged using Limma and sva software package.²² And the Limma software package was also used to eliminate the differences between batches. The differences of gene expression profiles among GSE6344, GSE40435, GSE781, TCGA_KIRC datasets were analyzed by using dplyr, ggplot2, ggrepel.²¹ The genes with too low expression value were deleted and the expression profile value was converted to log2-counts per million (logCPM), and the comparison matrix was constructed by performing linear regression. Where each row represented the gene names and each column represented the sample name of this study. The eligible differential genes were screened for |log2(FC)| > 2 and P value < 0.05 based on Bayesian calculation of T-values, F-values and log-odds, and volcano plots was plotted by using the ggplot2 program package to visualize the data.^22,23

GO, KEGG and GSEA Enrichment Analysis

Gene Ontology (GO) analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, Disease Oncology (DO) analysis and Gene Set Enrichment Analysis (GSEA) were calculated by using the R package Clusterprofiler program.²⁴ P value < 0.05 was considered significant, and gene counting the identified significant analysis were classified.

Screening and Identification of Gene Prediction Model for Early Diagnosis

The data of GSE6344, GSE40435, GSE781, TCGA_KIRC were selected as the training set, GSE53757 data was used as the test set. The glmnet package was used as the binomial LASSO model of the training set.²⁵ The operating characteristic curve of the recipient was plotted and its AUC was calculated. The pROC software package in R was used to evaluate the diagnostic value of hub genes.²⁶ SVM is a generalized linear classifier (generalized linear classifier) that performs binary classification of data in a supervised learning manner, which maximized the margin between different classes in a high-dimensional space.²⁷ The e1071, kernlab and caret package were used to eliminate the recursive features of the difference genes obtained, and data calculation, finally the best gene signature was obtained. The two machine learning algorithms were used to screen the hub genes of KIRC at the same time and the same key genes were obtained by using venn package of R. Then difference analysis of key genes was performed in the test set using limma package of R. Besides, a receiver operating characteristic (ROC) curve was constructed, and AUC value was calculated to estimate the predictive value of the model both in the training set and test set.

Clinical Prognostic Analysis

The TISIDB website (http://cis.hku.hk/TISIDB/) was used to analyze associations between key genes expression overall survival, stage, grade across KIRC.

The Analysis of Immune Cell Infiltration

We conducted the correlation analysis of immune cells, assessed the proportion of immune cells in renal clear cell carcinoma tissues and normal tissues by Cibersort analysis.^22,28–32 This revealed the characteristics of renal clear cell carcinoma immune cell infiltration by heat mapping using ggplot2. Among them, we conducted a correlation analysis between the hub genes that were screened out and immune cells.

Immune Subtype Analysis

Tumor microenvironment (TME) had importantly therapeutic and prognostic significance in antitumor therapy. Based on five representative immune signals, the researchers identified six immune subtypes of TCGA tumors, which provided resources for the analysis of TME in some specific tumors. Six stable and repeatable immune subtypes were classified according to five immune expression characteristics:³³ macrophages/monocytes,³⁴ total lymphocyte infiltration (mainly T and B cells),³⁵ TGF-β response,³⁶ IFN-γ response³⁷ and wound healing.³⁸ And the six subtypes may play a key role in predicting disease outcomes, rather than relying solely on the characteristics of individual cancer types.³³ In order to detect the mRNA expression levels of key genes of six different immune subtypes in TCGA tumor types, we used Kruskal test for differential expression analysis.

Immune Checkpoint Gene Analysis

The TISIDB website (http://cis.hku.hk/TISIDB/) was used to analyze associations between key genes expression and the gene expression of the immune checkpoint.

Result

Data Preprocessing and the Screening of Differential Genes

The flowchart of this study was shown in Figure 1. After gene annotation and standardization of the data, GSE40435 dataset contains 138 samples (31,280 genes), GSE6344 dataset contains 40 samples (GSE6344_GPL96 contains 12,548 genes, GSE6344_GPL97 contains 10,603 genes), GSE781 data set contains 29 samples (GSE781_GPL96 contains 12,548 genes, GSE781_GPL97 contains 10,603 genes), GSE53757 dataset contains 86 samples (20,109 genes). TCGA_KIRC dataset contains 583 samples (56,536 genes). The basic clinical information was shown in Table 1. People diagnosed with age over 60 accounted for 52.08% of the total, people diagnosed with age over 60 accounted for 52.08% of the total. A whole-gene expression profile analysis of the combined data was shown in Figure 2A which was the differentially expressed genes (DEGs) expression heat map. Figure 2B showed the subsequent differential analysis and volcano mapping of counts genes in patient tissue and normal tissue samples. A total of 25 DEGs were screened out, including 1 upregulated gene and 24 down-regulated genes.

Table 1 The Basic Clinical Information of Included Studies

Figure 1 The flowchart of this study.

Figure 2 The gene differential expression analysis of GSE6344, GSE40435, GSE781, TCGA_KIRC data set. (A) Whole gene expression heat map of normal tissues and cancer tissues with high expression in red and low expression in blue, with clinical information (B) the DEG Volcano map shows upregulated genes in red and down-regulated genes in green.

GO, KEGG and GSEA Enrichment Analysis of Differential Genes

GO, KEGG, DO and GSEA enrichment analysis were carried out based on the obtained DEGs using the clusterProfiler package. The results showed that the biological processes (BP) of DEGs were mainly related to organic anion transport, dicarboxylic acid transport, sodium ion transport, vascular process in circulatory system and other pathways. In terms of cell composition (CC), it was mainly related to basolateral plasma membrane, basal plasma membrane, basal part of cell, collagen−containing extracellular matrix and so on. In terms of molecular function (MF), it was mainly related to active ion transmembrane transporter activity, secondary active transmembrane transporter activity, anion transmembrane transporter activity, sodium ion transmembrane transporter activity and so on (Figure 3A and B). The DEGs were mainly concentrated in neuroactive ligand−receptor interaction, complement and coagulation cascades, bile secretion, influenza A in the KEGG pathway rich concentration (Figure 4A and B). The GSEA enrichment analysis showed that organic acid catabolic process, organic acid metabolic process, small molecule metabolic process, transmembrane transport, ion transmembrane transporter activity, cell activation, lymphocyte activation, response to virus, T cell activation were mainly active in cancer samples (Figure 5A and B).

Figure 3 Barplot (A) and bubble plot (B) figures of GO enrichment analysis of 25 differential genes gotten from DEG screening.

Figure 4 Barplot (A) and bubble plot (B) figures of KEGG enrichment analysis of 25 differential genes gotten from DEG screening.

Figure 5 Figures of GSEA enrichment analysis in normal tissues (A) and cancer tissues (B).

Screening and Identification of Gene Prediction Model for Early Diagnosis

The data of GSE6344, GSE40435, GSE781, TCGA_KIRC were selected as the training set, GSE53757 data was used as the test set. The hub genes that can predict early KIRC accurately was obtained by constructing the LASSO model in the training set and selecting the min valve to screen (Figure 6A). At the same time, we used SVM-RFE algorithm to screen 25 different genes and obtained 8 hub genes (Figure 6B). Then we got four key genes by taking the intersection of the genes obtained by these two algorithms, including ACPP, SCNN1G, ANGPTL4, SLC22A7 (Figure 7A). We performed difference analysis of key genes in the test set, it suggested that the expression of ANGPTL was higher in cancer samples compared to normal samples, the expression of ACPP, SCNN1G, SLC22A7 was lower in cancer samples compared to normal samples, which was consistent with the expression pattern of the training set (Figure 7B–E). A receiver operating characteristic (ROC) curve was constructed, and AUC value was calculated to evaluate the predictive value of the model both in the training set and test set. The results showed that the AUC of ACPP, ANGPTL4, SCNN1G, SLC22A7 in the training set respectively are 0.937, 0.947, 0.923, 0.839 (Figure 8A–D), and the AUC of ACPP, ANGPTL4, SCNN1G, SLC22A7 in the test set is 0.965, 0.983, 0.961, 0.957, indicating that the built model owned excellent verification performance (Figure 8E–H).

Figure 6 The potential key genes of KIRC were screened by LASSO regression model and SVM-RFE using training set. (A) LASSO regression model screened the potential key genes of KIRC. The X-axis represented the λ value, the Y-axis represented the cross-validation error. N=10 represented that 10 key genes were screened. (B) SVM-RFE algorithm screened the potential key genes. The X-axis represented the λ value, the Y-axis represented the cross-validation error. N=8 represented that 8 key genes were screened.

Figure 7 The intersection hub genes of the two algorithm (A) and each gene differential analysis in normal and cancer tissues, including ACPP (B), ANGPTL4 (C), SCNN1G (D), SLC22A7 (E).

Figure 8 ROC curves of training set and test set for gene signature. ACPP (A), ANGPTL4 (B), SCNN1G (C), SLC22A7 (D) in training set, ACPP (E), ANGPTL4 (F), SCNN1G (G), SLC22A7 (H) in test set.

Clinical Prognostic Analysis

The clinical prognostic analysis showed that patient with low SCNN1G expression had higher overall survival and the higher the tumor stage and grade, the higher SCNN1G expression level (p<0.05). Patient with high SLC22A7 expression had higher overall survival and the higher the tumor stage, the lower SCNN1G expression level (p<0.05). The expression of other genes has no obvious relationship with the prognosis (Figure 9).

Figure 9 Clinical prognostic analysis including the association between the expression of key genes including ACPP, ANGPTL4, SCNN1G, SLC22A7 and the overall survival, stage, grade across KIRC.

The Analysis of Immune Cell Infiltration

The correlation analysis of immune cells showed that NK cells resting had a negative correlation with NK cells activated (R=−0.34), Monocyte had a negative correlation with macrophage M2 (R=−0.63), neutrophils (R=−0.34), plasma cells had a negative correlation with dendritic cells resting (R=−0.3), macrophage M2 (R=−0.3), T cells CD4 memory activated had a negative correlation with T cells regulatory (Treg) (R=−0.43), NK cells activated (R=−0.4), T cells CD4 memory resting (R=−0.33), had a positive correlation with dendritic cells activated (R=0.33), macrophage M0 had a negative correlation with dendritic cells resting (R=−0.44) (Figure 10A). The relative percentage of 22 immune cell subsets result showed that macrophages M2, NK cells, T cells accounted for the most whether in tumor samples nor normal samples (Figure 10B). The result of immune infiltration difference between normal samples and cancer samples showed that normal samples had more neutrophils, T cells CD4 memory resting (P<0.05) (Figure 10C). The correlation analysis result between the key genes and immune cells showed that the expression of ACPP had a positive correlation with eosinophils (R=0.33, P<0.05), macrophage M2 (R=0.26, P<0.05), had a negative correlation with NK cells resting (R=−0.25, P<0.05), the expression of ANGPTL4 had a positive correlation with monocytes (R=0.28, P<0.05), plasma cells (R=0.2, P<0.05), had a negative correlation with macrophages M2 (R=−0.31, P<0.05), neutrophils (R=−0.2, P<0.05), the expression of SCNN1G had a positive correlation with neutrophils (R=0.2, P<0.05), the expression of SLC22A7 had a positive correlation with T cells CD4 memory resting (R=0.2, P<0.05) (Figures 11 and 12).

Figure 10 (A) The correlation among immune cells of the samples in the training set. (B) Relative percentage of 22 immune cell subsets in normal samples and cancer samples. (C) Vioplot diagram of immune infiltration difference between normal samples and cancer samples, green as normal samples, red as cancer samples.

Figure 11 Lollipop figures of the correlation between the key genes including ACPP (A), ANGPTL4 (B), SCNN1G (C), SLC22A7 (D) and immune cells using the training set. The P value in Y-axis <0.05 was marked in red.

Figure 12 Scatter plot of the correlation between key gene expressions and immune cells using the training set.

The Immune Subtype Analysis

We used Kruskal–Wallis test to analyze the mRNA expression of key genes in 6 immune subtypes of KIRC. The expression levels in C1-C6 of ACPP (p<0.05), ANGPTL4 (p<0.05), SLC22A7 (p<0.05), SCNN1G (p<0.05) were with significant statistical difference (Figure 13). Obviously, the expression of ANGPTL4 ranked the first among the five groups and was the lowest in the C5 subtype among the six subtypes (p<0.05). The expression of ACPP was the highest in the C5 subtype among the six subtypes (p<0.05). The expression of SLC22A7 and SCNN1G were the lowest in the C5 subtype among the six subtypes (p<0.05).

Figure 13 Kruskal test analysis of key gene expression levels in immune subtypes of KIRC using TISIDB website.

The Immune Checkpoint Related Genes Analysis

The results of the immune checkpoint related genes analysis were shown in Figure 14 and 15. The gene expression of ACPP was positively correlated with the expression of most immune checkpoint genes such as BTLA (R=0.333, p<0.05), IL10 (R=0.373, p<0.05), CSF1R (R=0.535, p<0.05), PDCD1LG2 (R=0.355, p<0.05), TGFBR1 (R=0.33, p<0.05), TNFRSF8 (R=0.322, p<0.05), TNFSF13B (R=0.473, p<0.05), CD80 (R=0.369, p<0.05), CD86 (R=0.502, p<0.05), CD48 (R=0.311, p<0.05), CD28 (R=0.35, p<0.05), HLA_DPA1 (R=0.368, p<0.05), HLA_DPB1 (R=0.355, p<0.05), HLA_DQA1 (R=0.311, p<0.05), HLA_DRA (R=0.401, p<0.05). The gene expression of ANGPTL4 was positively correlated with the expression of most immune checkpoint genes such as PTNFRSF4 (R=0.383, p<0.05), TNFRSF14 (R=0.315, p<0.05), CD40 (R=0.327, p<0.05), CD70 (R=0.354, p<0.05), negatively correlated with the expression of some immune checkpoint genes such as TGFBR1 (R=−0.307, p<0.05), TNFSF15 (R=−0.453, p<0.05). The gene expression of SLC22A7 was positively correlated with the expression of most immune checkpoint genes such as HHLA2 (R=0.599, p<0.05). The gene expression of SCNN1G was positively correlated with the expression of most immune checkpoint genes.

Figure 14 The correlation analysis between key genes expression and immune checkpoint genes’ expression using TISIDB website.

Figure 15 The scatter plots of key genes expression and immune checkpoint genes’ expression using TISIDB website.

Discussion

In this study, DEG screening and GO, KEGG and GSEA enrichment analysis was firstly conducted and then the LASSO and SVM-RFE algorithms were further employed to screen out KIRC-associated hub genes in the training set and validate them in the test set. As we know, this was the first study combining LASSO and SVM-RFE algorithm to identify key biomarkers of KIRC. A total of 25 DEGs were screened out, including 1 upregulated gene and 24 downregulated genes. KEGG pathway analysis showed that 25 DEGs were significantly associated with complement and coagulation cascades. Complement and coagulation cascades were vicious cycle which were mutually induced and they can prevent immune cells from attacking cancer cells, promote immune escape, thus inducing tumor progression and metastasis.^39,40 GSEA enrichment analysis showed that lymphocyte activation and T cell activation was active in cancer group, suggesting that immune treatment might be effective.

In this study, we finally screened out 4 key genes, including ACPP, ANGPTL4, SCNN1G, SLC22A7. The AUC values of the 4 key genes in the test set were >0.95 but could only reach about 0.9 in the training set, indicating that the built model had excellent verification performance but generalization performance need to be perfected. ACPP, the prostatic acid phosphatase, was capable of dephosphorylating substrates including alkyl, aryl, and acyl orthophosphate monoesters and phosphorylated proteins when the pH of the conditions was 4–6. It dephosphorylated ERBB2 and deactivated MAPK-mediated signaling to act as the tumor suppressor of prostate cancer (http://cis.hku.hk/TISIDB/browse.php?gene=ACPP). Its expression was downloaded in the tumor samples compared to the normal samples in our study. And its expression was positively correlated with expressions of most immune checkpoint gene expression. Novel immune checkpoints, including CTLA-4, PD-1, and their ligands, programmed death-ligand 1 (PD-L1 (CD274)) and PD-L2 (CD273), provided more therapeutic options for cancer treatment,⁴¹ and their inhibitors significantly prolonged the overall survival (OS) of cancer patients.⁴² However, the objective response rate of ICIs therapy was only ~17%, indicating that a majority of patients fail to benefit from ICIs treatment.^43,44 Thus, there was a large unmet clinical need to identify the appropriate patients who may respond to ICIs therapy. Therefore, patients with high expression of ACPP might be more sensitivity to ICIs. The expression of ACPP was positively correlated the abundance of tumor-infiltrating cells such as eosinophils. The malignancy of tumor was related to several factors, including genomic instability, heterogeneity (cells types such as fibroblasts, endothelial cells, pericytes, and immune cells), and composition of the microenvironment.⁴⁵ Therefore, tumors’ heterogeneity could be changed by regulating genes’ expression, to further reduce the malignant degree of the tumor. The angiopoietin-like 4 (ANGPTL4) protein, as a member of secreted proteins superfamily, whose structure was related to factors modulating angiogenesis known as angiopoietins,⁴⁵ was involved in regulation of angiogenesis, lipid metabolism and glucose and energy homeostasis.^46–51 The overexpression of ANGPTL4 could promote tumorigenesis, angiogenesis, tumor invasion, and metastasis,^52–56 were also capable of inhibiting tumor metastasis via inhibiting vascular permeability, motility, and invasiveness⁵⁷ or via attenuating endothelial cell adhesion, migration, and sprouting.^58,59 ANGPTL4 had been reported as a potential biomarker in renal cell carcinoma,^56,59 the expression level of which would increase in KIRC. The gene expression of ANGPTL4 was negatively correlated with the expression of some immune checkpoint genes such as TGFBR1, CD274, TNFSF15, positively correlated with the expression of some immune checkpoint genes such as PDCD1, TNFRSF4, TNFRSF14, CD40, CD70, which might provide thoughts on immune checkpoint inhibitors. SCNN1G encoded one of subunits, γ subunit, which combined with α subunit and β subunit to form a nonvoltage-gated, amiloride-sensitive, sodium channels controlling fluid and electrolyte transport across epithelia in many organs. The expression of SCNN1G upregulated in Ewing’s sarcoma and fibromatosis samples compared to normal samples.⁶⁰ In our study, although the expression of SCNN1G was low in cancer samples of KIRC compared to normal samples, the survival analysis showed that the high expression group had lower survival rate. As the tumor progressed, the expression increased, indicating that elevated expression meant tumor progression and poor prognosis. Thus, high expression of SCNN1G of patients meant poor prognosis. SLC22A7 could mediate sodium-independent multispecific organic anion transport, including prostaglandin E2, prostaglandin F2, tetracycline, bumetanide, estrone sulfate, glutarate,⁶¹ dehydroepiandrosterone sulfate, allopurinol, 5-fluorouracil, paclitaxel, L-ascorbic acid, salicylate, ethotrexate, and alpha-ketoglutarate.^62,63 It was also reported that its expression was low in KIRC compared to normal samples using different dataset, which was consistent with our analysis result.⁶³ The survival analysis showed that high expression of SLC22A7 had higher survival analysis, and the more severe the disease, the lower the expression level. It indicated that high expression level was a sign of early stage of disease and better prognosis of patients.

The clinical prognostic analysis including the association between the expression of key genes and the overall survival, stage, grade across KIRC, the immune infiltration difference between normal samples and cancer samples, the correlation between the key genes and immune cells, immunemodulator, immune subtypes of KIRC were investigated in this research, which could provide new guidance for cancer diagnosis. The study had limitations. Firstly, the data in our study may be limited to the country, region, race of included research objects. Secondly, due to the data limitations, the molecular subtypes of KIRC have not been further classified. Last, our analysis needs further clinical validation.

Conclusion

This study screened the 4 key genes, including ACPP, ANGPTL4, SCNN1G, SLC22A7, whose expression showed difference among normal samples and tumor samples. Patients with higher expression of SCNN1G have a lower survival rate, the more severe the disease, the higher the expression level. Patients with higher expression of SLC22A7 have a higher survival rate, the more severe the disease, the lower the expression level. Meanwhile, the expression of key genes was related with the abundance of tumor infiltration immune cells and the expression of immune checkpoint, which could provide guidance for immunotherapy.