Identification of an immune-related risk signature for predicting prognosis in clear cell renal cell carcinoma

Identification of differently expressed IRGs in ccRCC

We analyzed 109 paired tumor and adjacent normal tissue samples from ccRCC patients available in the Gene Expression Omnibus (GEO) database and identified 2,128 differently expressed genes (DEGs), of which 1,036 were up-regulated and 1,092 were down-regulated. From 72 paired specimens in The Cancer Genome Atlas (TCGA), we further identified 4,629 DEGs, including 2,174 up-regulated and 2,455 down-regulated transcripts. Overlapping DEGs between GEO and TCGA datasets (756 up-regulated and 746 down-regulated transcripts) were next extracted (Figure 1B) and intersected with an IRG set composed of lists combined from ImmPort and InnateDB databases. A total of 326 differentially expressed IRGs was thus obtained (Figure 1C and Supplementary Table 1). Principal component analysis (PCA) of these IRGs showed varying distribution patterns, indicating distinct immune phenotypes in normal tissue and tumor samples both in GEO (Supplementary Figure 1A) and TCGA datasets (Supplementary Figure 1B). As expected, Gene Ontology (GO) analysis results revealed that inflammatory and immune pathways were most significantly enriched with these IRGs. On GO analysis, “inflammatory response”, “extracellular space”, and “receptor activity” were the most significant terms among biological processes, cellular components, and molecular functions, respectively (Supplementary Table 2). The top 20 GO terms are presented in Supplementary Figure 2A. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that “cytokine-cytokine receptor interaction” was the pathway most significantly enriched with the differentially expressed IRGs (Supplementary Figure 2B). Protein Analysis Through Evolutionary Relationships (PANTHER) pathway analysis revealed that immune-related pathways, such as interleukin signaling pathway, inflammation mediated by chemokine and cytokine signaling pathway, B cell activation, and T cell activation, were associated with these differentially expressed IRGs (Supplementary Table 3). Additional results from disease ontology (DO) analysis confirmed that this IRG set was associated with immune status and showed also involvement in kidney disease (Supplementary Table 4). The top 15 most significantly enriched diseases are presented in Supplementary Figure 2C.

Figure 1. Flow chart of the study design and identification of differentially expressed genes (DEGs). (A) Flow chart of our research project. (B) Using gene expression data of 109 pairs of normal and tumor tissues from Gene Expression Omnibus and 72 paired samples from The Cancer Genome Atlas (TCGA) database, 756 up-regulated and 746 down-regulated DEGs were identified and extracted from the two databases. (C) These DEGs were intersected with a gene set including 2103 immune-related genes (IRGs), and 326 differently expressed IRGs were identified. (D) Harrell's concordance index of the five-gene signature in the TCGA discovery, validation, and entire sets, and in the E-MTAB-3267 dataset.

Identification of hub IRGs

We next imported the identified IRGs into the STRING database to create a protein-protein interaction (PPI) network, thus obtaining 47 hub genes defined by a connectivity degree > 20 and an interaction score > 0.7. In line with the above GO analysis results, GO analysis of these hub genes revealed that “inflammatory response”, “extracellular space”, and “chemokine activity” were the terms most significantly enriched (Supplementary Table 5). Figure 2A shows the top 20 GO terms related to the hub IRGs. On KEGG pathway analysis “chemokine signaling pathway” was the most frequently enriched process (Figure 2B), and highly concordant results were obtained for diverse immune-related pathways after data analysis using the PANTHER algorithm (Supplementary Table 6). Meanwhile, DO analysis results (Supplementary Table 7) confirmed the association of these hub IRGs with kidney disease (Figure 2C).

Figure 2. Functional enrichment of 47 immune-related hub genes. (A) Top 20 most significant Gene Ontology (GO) terms identified by GO analysis. (B) Kyoto Encyclopedia of Genes and Genomes pathway analysis. FDR < 0.01 indicated significant enrichment. (C) Top 15 most significant disease ontology (DO) terms identified by DO analysis.

Identification of a five-IRG signature in the TCGA discovery set

We used the least absolute shrinkage and selection operator (LASSO) Cox regression model to identify the most valuable prognostic IRGs based on their association with patient survival in our TCGA discovery set (Supplementary Figure 3). Five IRGs, namely SHC1, IRF7, KDR, JAK3, and CXCL5 were thus selected (see detailed information in Supplementary Table 8). Univariate Cox analysis results for these IRGs are shown in Table 1. All five IRGs were up-regulated and correlated negatively with patient survival (P < 0.05) except for KDR, which showed a positive relationship. The same results were obtained on Kaplan-Meier (K-M) analyses (Supplementary Figure 4A). The prognostic value of these 5 genes was further validated in the Gene Expression Profiling Interactive Analysis (GEPIA) database. Here, consistent results were obtained, except for SHC1 (Supplementary Figure 4B). The relative expression of each gene in our IRG signature was further verified by qRT-PCR in 35 paired ccRCC/normal specimens collected at our institution. Results corroborated significantly higher expression in tumor samples (Supplementary Figure 5), in agreement with the database validation approach.

Based on the respective expression levels, we then established the following risk score formula:

Risk score = (0.08995984 * SHC1) + (0.05754872 * IRF7) + (-0.13910054 * KDR) + (0.01889022 * JAK3) + (0.02791388 * CXCL5)

The C-index of the risk score in the TCGA discovery set was 0.666 (95% confidence interval [CI], 0.603-0.729; Figure 1D). Using X-tile software, the optimum cut-off value to distinguish high-risk (n = 130) from low-risk (n = 127) patients was 0.135 (Supplementary Figure 6A). High-risk patients showed poorer prognosis, while low-risk patients had better overall survival (Figure 3A, left panel). Time-dependent receiver operating characteristic (ROC) analysis was performed to evaluate the accuracy of the five-gene signature in predicting patient survival. The area under the curve (AUC) was 0.753 at 1 year, 0.686 at 3 years, and 0.637 at 5 years (Figure 3A, middle panel). K-M analysis corroborated that high-risk patients had a significantly lower survival rate than low-risk cases (P < 0.001, Figure 3A, right panel).

Figure 3. Validation of the prognostic risk signature. Left panel: Distribution of the risk signature based on survival status. High-risk and low-risk patients were distributed above and below the x-axis, respectively. Pink and green colors indicate dead and alive patients, respectively. Middle panel: Time-dependent ROC curves were performed to evaluate the accuracy of the risk signature. Right panel: Kaplan-Meier survival curves were performed to assess patients’ prognosis. (A) TCGA discovery set. (B) Validation set. (C) Entire set. (D) E-MTAB-3267.

Validation of the five-IRG signature

The prognostic value of the five-IRG signature was further evaluated in three validation sets (TCGA internal validation set, TCGA entire set, and the E-MTAB-3267 dataset). We calculated the risk score for each patient using the same formula. Patients in TCGA internal and entire sets were classified into high-risk and low-risk groups with the cut-off value (0.135) used in TCGA discovery set. Patients in the E-MTAB-3267 dataset were divided into two groups with a cut-off value of 0.448 (Supplementary Figure 6B) determined using X-tile software. The C-indexes of risk scores in TCGA internal, TCGA entire, and E-MTAB-3267 sets were 0.686 (95% CI, 0.623-0.749), 0.677 (95% CI, 0.634-0.720), and 0.591 (95% CI, 0.481-0.701), respectively (Figure 1D). Consistent with the results from the TCGA discovery set, high-risk patients in the three validation sets had poorer prognosis than those in the low-risk group (Figure 3B–3D, left panel). The ROC curves of the five-gene signature in the three validation sets showed good performance (Figures 3B–3D, middle panel). Survival analysis in the three validation sets confirmed lower survival rate in the high-risk groups (Figures 3B–3D, right panel).

Association with clinicopathological factors and sub-group analysis

Results of the correlation analysis between clinical factors and our IRG signature’s risk score is shown in Table 2. The risk score of the signature was correlated with gender (P < 0.001), stage (P < 0.001), and grade (P < 0.001), but not with age (P = 0.347). Indeed, we found that male gender, advanced stage, high grade, and high-risk patients tended to have higher risk score (P < 0.05, Supplementary Figure 7). Subgroup analysis was performed to further assess whether the five-IRG signature had prognostic value for survival within specific clinical parameters. The results showed that the our risk signature was still a powerful tool for predicting survival in younger (age < 65) and older (age ≥ 65) patients, male or female patients, low (Fuhrman grade 1 & 2) or high (Fuhrman grade 3 & 4) grade patients, and advanced disease stage (stage III & IV) patients (P < 0.05, Supplementary Figure 8). In the subgroup of early stage (stage I & II), the prognosis of the high-risk group also tended be poorer than the low-risk group, although this difference did not reach significance (P = 0.056). Overall, the risk score based on our IRG signature was able to predict overall survival independently of specific clinical factors in ccRCC patients.

Distinct immune phenotypes characterize high-risk and low-risk ccRCC patients

Considering the distinct prognosis in high-risk and low-risk ccRCC patients, we sought to explore whether phenotypic differences in immune cell populations are found upon segregation of cases based on our IRG signature’s specific risk score. PCA based on the 47 hub IRGs (Figure 4A) and the 326 differentially expressed IRGs (Figure 4B) showed different distribution patterns, indicating remarkably different immune phenotypes for high- and low-risk groups. Gene Set Enrichment Analysis (GSEA) was also performed to explore whether immunological pathways differ between risk groups. As expected, significant alterations in immunological pathways were detected in high-risk patients, compared with low-risk ones (Figure 4C). Enriched pathways and processes in tumors from high-risk ccRCC patients included primary immunodeficiency, intestinal immune network for IgA production, cytokine-cytokine receptor interaction, complement and coagulation cascades, Nod-like receptor signaling pathway, and cell cycle.

Figure 4. Principal component analysis (PCA) and gene set enrichment analysis (GSEA). (A) PCA based on 47 immune-related hub genes (IRGs) showing distinct immune phenotypes in high- and low-risk patient groups. (B) Distribution patterns for the two risk groups based on 326 differentially expressed IRGs. (C) GSEA results showing significant enrichment of immune-related phenotype in high-risk patients.

Immune landscapes differ between low-risk and high-risk ccRCC patients

To assess whether the five-IRG signature could accurately reflect the status of the tumor immune microenvironment, we estimated infiltration rates for 22 immune cell types in ccRCC patients through CIBERSOFT algorithm, and investigated potential differences between the low- and high-risk groups. As shown in Figure 5A, the distribution of immune cell types varied within and between groups, suggesting that the immune landscape in ccRCC cases might be an inherent feature of each individual.

Figure 5. Immune landscapes in high-risk and low-risk patients. (A) Relative proportions of 22 immune cell types in high- and low-risk patients. (B) Correlation matrix of relative proportions of the 22 immune cell types. (C) Box plots showing differential immune cell infiltration status between high- and low-risk patients. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

The relative proportions of the 22 immune cell types were weakly to moderately correlated with risk classification (Figure 5B). Patients in the high-risk group had significantly higher representation of CD8 T⁺ cells, T follicular helper cells, regulatory T cells (Tregs), and M0 macrophages, and significantly lower abundance of natural killer (NK) cells, plasma cells, dendritic cells, mast cells, and eosinophils (P < 0.05, Figure 5C), indicating a weakened immune phenotype in high-risk patients. These data suggest that the five-IRG signature could serve as an indicator of immune status in ccRCC.

We also evaluated the correlation between the IRG signature’s risk score and the expression of T-cell markers, including CD4 and CD8A, and immune checkpoint genes such as PD-1 and its ligands (PD-L1, PD-L2), cytotoxic T-lymphocyte associated protein 4 (CTLA-4), lymphocyte-activation gene 3 (LAG3), and CD47. We found that the risk score was positively correlated with the expression of CD4, CD8A, PD-1, CTLA-4, LAG3, and CD47, and negatively correlated with PD-L1 expression (P < 0.05, Supplementary Figure 9). Next, we contrasted these findings in low- and high-risk ccRCC patients. In strong agreement with the risk classification derived from our five-IRG signature, we found that CD4, CD8A, PD-1, CTLA-4, LAG3, and CD47 were significantly overexpressed in high-risk patients, compared to low-risk ones (P < 0.05; Figure 6). These data further suggest the presence of an immunosuppressive tumor microenvironment in high-risk patients, which might help explain their poor prognosis.

Figure 6. Expression of T-cell and immune checkpoint markers in low-risk and high-risk patients. Group differences were assessed by Wilcoxon test.

Improved prognostic prediction through integrative analysis of the IRG signature and clinicopathological factors

To assess whether inclusion of clinicopathological factors could increase the prognostic accuracy of our five-IRG signature’s risk score, patients’ age, gender, Fuhrman grade, and stage were accessed to perform univariate and multivariate Cox analyses in the discovery, validation, and entire TCGA sets. Tumor grade and stage were defined according to the Fuhrman classification system and the 7^th edition American Joint Committee on Cancer TNM classification, respectively [16]. Regression analysis results are shown in Table 3. The multivariate analysis showed that the risk score based on the IRG signature was an independent prognostic factor in the discovery set (hazard ratios (HR): 3.964; 95% CI: 1.852-8.486; P < 0.001), the validation set (HR: 3.957; 95% CI: 1.727-9.067; P = 0.001), and the entire TCGA set (HR: 3.896; 95% CI: 2.228-6.813; P < 0.001). Additionally, we found that age (HR: 1.030; 95% CI: 1.015-1.044; P < 0.001), grade (HR: 1.269; 95% CI: 0.999-1.612; P = 0.051), and stage (HR: 1.653; 95% CI: 1.124-2.432; P = 0.011) were independent prognostic factors in the entire set.

Based on the results of multivariate analysis of the entire set, C-indexes and Akaike information criterions (AICs) were calculated to evaluate the power of selected parameters, i.e. risk score, age, grade, and stage (Table 4). The combination of the five-IRG risk score with clinical factors had a higher C-index (0.772; 95% CI: 0.737-0.807) and a lower AIC than the risk score or the clinical factors alone. This indicated that combining our IRG risk score with clinical variables can improve prognostic accuracy for ccRCC. To provide a quantitative method to predict survival probability of ccRCC patients after surgery, we constructed a nomogram integrating the risk score of the IRG signature and clinical factors (Figure 7A). Calibration plots indicated that the nomogram showed good performance for predicting 1-year, 3-year, and 5-year survival probabilities (Figure 7B).

Figure 7. Nomogram and calibration plots for prediction of patients’ survival in the entire TCGA set. (A) Nomogram combining the five-IRG risk signature with clinical factors for prediction of 1-year, 3-year, and 5-year survival rates. (B) Calibration plots showing high predictive accuracy of the nomogram.

Study design and data collection

A flow chart of the study design is shown in Figure 1A. Raw data from mRNA expression microarrays profiling a total of 109 ccRCC and matched normal tissue samples were analyzed. The corresponding datasets included GSE53757 (72 pairs), GSE36895 (23 pairs), and GSE66270 (14 pairs), all downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). The three GEO datasets were produced using the Affymetrix Human Genome U133 Plus 2.0 Array platform (Affymetrix, Santa Clara, CA, USA). Expression data of 72 paired ccRCC and adjacent normal tissue samples in TCGA database were obtained from UCSC Xena (https://xenabrowser.net/). The GEO and TCGA datasets were used to identify DEGs, thus increasing DEG detection reliability. IRG lists were downloaded from the ImmPort database (http://www.ncbi.nlm.nih.gov/geo/) [48] and the InnateDB database (https://www.innatedb.com/) [49], which provided a total of 1534 and 1378 human IRGs, respectively. After removal of duplicate genes, a combined gene set that included 2103 unique genes (Supplementary Table 9) was used to screen hub IRGs.

To construct an IRG signature, gene expression data of all ccRCC samples, and their associated clinical and survival information, were download from UCSC Xena. Criteria for study inclusion were: (1) Repeated tumor samples in the same patient were removed. (2) Patients with unknown survival status and follow-up information, and those who died within a follow-up period of 30 days were excluded. (3) Patients with unknown disease stage or grade were excluded. Finally, 514 patients meeting the inclusion criteria were randomly divided into two groups, the discovery set (n = 257) and the internal validation set (n = 257). The discovery set was used to establish an IRG signature to predict patients’ prognosis, while the internal validation set and the entire TCGA set were used for internal validation. The E-MTAB-3267 dataset from ArrayExpress database (https://www.ebi.ac.uk/arrayexpress/), which includes 53 ccRCC samples, was used for external validation. Detailed information is presented in Supplementary Table 10.

Data preprocessing

For the three GEO datasets, Robust Multi-array Averaging (RMA) and the ComBat algorithm were used for background correction, data normalization, and to remove batch effects. Then, we carried out PCA and generated a heat map to visualize expression patterns and correction effects (Supplementary Figure 10); normalization was sequentially applied. For E-MTAB-3267 microarray data, RMA background correction and normalization were performed for processed signals. The “affy’’ R package and Affymetrix annotation files were used, respectively, to summarize and annotate the probes. For TCGA datasets, gene expression values were presented as the log2(x+1) transformed RSEM-normalized counts.

Collection of clinical samples

A total of 35 matched tumor/normal tissue specimens were collected from ccRCC patients after surgery at The First Affiliated Hospital of Anhui Medical University (Hefei, China). All patients signed informed consent forms. Histological diagnosis of ccRCC was confirmed by two pathologists. Samples were immediately frozen in liquid nitrogen and stored at -80°C until RNA extraction. The study was approved by the Ethics Committee of Human Research of The First Affiliated Hospital of Anhui Medical University (No. PJ2019-14-22).

Reverse-transcription and qRT-PCR

Total RNA from ccRCC tissues was extracted using Trizol Reagent BD (Invitrogen, Carlsbad, CA). The reverse transcription reactions were performed using a PrimeScriptTM RT reagent kit (Takara, Kusatsu, Japan) according to the manufacturer’s instructions, and qRT-PCR was prepared at a final volume of 20 μl using a SYBR Green Mix (Takara, Kusatsu, Japan) with primers synthesized by Sangon Biotech (Sangon, Shanghai, China). The reactions were measured on an ABI7500 platform (Thermo, Massachusetts, USA). The 2−ΔΔCT method was used to determine relative gene expression levels, and GAPDH was used as an internal control to normalize the data. Each reaction was performed in triplicate. The primers used for these reactions are presented in Supplementary Table 11.

Screening of differentially expressed IRGs

Differentially expressed mRNAs between tumor and adjacent normal tissue samples were screened by the “limma” package of R software 3.4.2 (https://www.r-project.org/). |log2 fold change| > 1 and false discovery rate (FDR) < 0.05 were set as the cut-off criteria. Overlapping DEGs in paired samples from GEO and TCGA datasets, characterized by upregulation and downregulation, were extracted and intersected with an IRG set to obtain differentially expressed IRGs.

Identification of hub IRGs

The differentially expressed IRGs were imported into the STRING database (https://string-db.org/) to construct a PPI network. Genes with a connectivity degree of > 20 and interaction score > 0.7 in the PPI were defined as hub genes for further analysis.

Construction and validation of an IRG signature

After identifying hub IRGs, we performed LASSO Cox regression analysis [50] in the TCGA discovery set to select the best gene model for predicting prognosis in ccRCC patients. LASSO Cox regression analysis was performed using the R package “glmnet” and the optimal values of penalty parameters were determined by 10-fold cross-validation. Subsequently, we performed univariate Cox regression and survival analysis for the selected genes, and their prognostic values for overall survival were further validated using the GEPIA database (http://gepia.cancer-pku.cn/) based on available ccRCC data [51]. The IRG risk score model for each patient was determined by a linear combination of gene expression weighted by the regression coefficient from LASSO Cox regression analysis [52, 53]. Then, we used X-tile software version 3.6.1 (Yale University, New Haven, CT, USA) [54] to select the optimum cut-off value for the risk score of the IRG signature based on the patients’ survival information in the TCGA discovery set. This cut-off value was next used to divide patients into high- and low-risk groups.

To assess the predictive accuracy of the IRG signature, Harrell's concordance index (C-index) was calculated and time-dependent ROC analysis was performed. The AUC at different cut-off times was calculated using the “survival ROC” package in R [55]. K-M survival curves were used to assess survival differences between high- and low-risk groups using the “survminer” package in R. PCA was performed to assess gene expression patterns.

Estimation of relative abundance of immune cell types

The CIBERSORT algorithm (https://cibersort.stanford.edu) [56], an approach to quantify the relative abundance of immune cell types based on specific gene expression profiles, was used to assess the distribution (normalized to 1) of 22 immune cell types in ccRCC samples from TCGA. Moreover, the P-value, correlation coefficient and root mean squared error were also presented to evaluate the accuracy of the results in each patient. A total of 500 ccRCC Patients with P-value < 0.05 were retained, and were then used to compare immune cell composition between low-risk and high-risk patients.

Gene enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/) [57] was used to perform GO and KEGG pathway analysis with the cut-off criterion of FDR < 0.01. PANTHER (http://pantherdb.org) was also used to performed pathway analysis [58]. Moreover, we performed DO analysis using the “cluster Profiler” package in R with the cut-off criteria of P < 0.05 and Q < 0.05.

Gene set enrichment analysis (GSEA)

According to the risk score of the IRG signature, 514 ccRCC samples downloaded from TCGA were divided into high- and low-risk groups. To compare immune phenotypes between the two groups and to unmask potential functions of gene signature members, GSEA (http://software.broadinstitute.org/gsea/index.jsp) was carried out. The annotated gene set list c2.cp.kegg. v5.2.symbols.gmt was selected as the reference gene set. FDR < 0.25 and a nominal P < 0.01 were considered as cut-off criteria.

Nomogram construction and validation

A nomogram integrating the IRG signature and various clinicopathological factors was established for predicting patients’ prognosis using the “rms” package in R [59]. Calibration curves were established to assess the accuracy of the nomogram. The C-index and AIC, representing the effect of prognosis factors, were calculated and compared to assess the predictive accuracy of the model.

Statistical analysis

Boxplots and PCA plots were generated using the “ggplot2” package in R. The heat map was generated using the “pheatmap” package in R. Student’s t test was used to compare subgroups and paired data. Pearson’s chi-square test was applied to analyze differences between the discovery and validation sets and the association between the risk score and clinical parameters. Pearson’s correlation test was used to analyze the correlation between the IRG signature and the expression of immune checkpoint genes. K-M survival curves were compared using log-rank test. Univariate Cox regression was conducted to estimate the HR for different factors. Multivariate Cox regression was performed to determine independent factors. All statistical analyses were conducted on R software. P < 0.05 was considered significant.