Development and validation of the prognostic value of the immune-related genes in clear cell renal cell carcinoma.

BACKGROUND: Clear cell renal cell carcinoma (ccRCC) is a highly heterogeneous tumor, resulting a challenge of developing target therapeutics. Not long ago, immune checkpoint blockade regimens combine with tyrosin kinase inhibitors have evolved frontline options in metastatic RCC, which implies arrival of the era of tumor immunotherapy. Studies have demonstrated immune-related genes (IRGs) could characterize tumor milieu and related to patient survival. Nevertheless, the clinical significance of classifier depending on IRGs in ccRCC has not been well established.
METHODS: The R package limma, univariate and LASSO cox regression analysis were used to screen the prognostic related IRGs from TCGA database. Multivariate cox regression was utilized to establish a risk prediction model for candidate genes. Quantitative real-time PCR was used to confirm the expression of candidates in clinical samples from our institution. CIBERSORT algorithm and correlation analysis were applied to explore tumor-infiltrating immune cells signature between different risk groups. A clinical nomogram was also developed to predict OS by using the rms R package based on the risk prediction model and other independent risk factors. The ICGC data was used for external validation of either gene risk model or nomogram.
RESULTS: We identified 382 differentially expressed immune related genes. Four unique prognostic IRGs (CRABP2, LTB4R, PTGER1 and TEK) were finally affirmed to associate with tumor survival independently and utilized to establish the risk score model. All candidates' expression was successfully laboratory confirmed by q-PCR. CIBERSORT analysis implied patients in unfavorable-risk group with high CD8 T cell, regulatory T cell and NK cell infiltration, as well as high expression of PD-1, CTLA4, TNFRSF9, TIGIT and LAG3. A nomogram combined IRGs risk score with age, gender, TNM stage, Fuhrman grade, necrosis was further generated to predict of 3- and 5-year OS, which exhibited superior discriminative power (AUCs were 0.811 and 0.795).
CONCLUSIONS: Our study established and validated a survival prognostic model system based on 4 unique immune related genes in ccRCC, which expands knowledge in tumor immune status and provide a potent prediction tool in future. 2021 Translational Andrology and Urology. All rights reserved.

Introduction

Clear cell renal cell carcinoma (ccRCC) representing approximately 75% of RCC cases with more than 175000 deaths per year (1). Although surgical resection is effective for localized RCC, about one-third cases suffered recurrences and metastases with worse prognosis (2). Various molecular signatures of ccRCC implies that distinct survival advantages exist in the certain subtypes (3,4). Owing to heterogeneity, discovering reliable molecular biomarkers can help to improve prognostic determination and guide clinical decision. Actually, RCC is believed to be an immunogenic tumor for long time. Interleukin 2 (IL-2) and interferon alpha (IFN-a) were used for therapeutic regimens for advanced RCC in the 1990s to early 2000s, and the incidence of complete remission was about 3–5% (5,6). Recently, trail of KEYNOTE-426 and JAVELIN Renal 101 demonstrates using PD-1 immune checkpoint inhibitor-based combination regimens as the first line setting can significantly improve advanced RCC survival, which have been approved by FDA (7,8). Therefore, further exploration of immune related molecular network in RCC definitely helps to develop comprehensive understanding of immune evasion and provide insights into making therapeutic strategy. Immune molecular regulation is the key mechanism for host innate immunity and immune surveillance. It is necessary to explore clinical significance of immune-related biomarkers, especially immune-related genes (IRGs) which could predict prognosis of patients, and potentially portrait tumor microenvironment (TME) (9,10).

In this study, we identified 4 immune-related genes (CRABP2, LTB4R, PTGER1 and TEK) through integrated analyses of mRNA expression data from TCGA database and independently assessed. Multivariate Cox proportional hazards models were constructed by 4 genes and validated the accuracy in an external ICGC dataset. Moreover, we investigated a high proportion of CD8 T cell, regulatory T cell and NK cell in the unfavorable-risk group. High levels of immune response suppressors (PD-1, CTLA4, TNFRSF9, TIGIT and LAG3) were observed in unfavorable-risk group and positively correlated with risk score. According to aforementioned data, a nomogram was well established for clinical use and also externally validated its superior power by ICGC data.

We present the following article in accordance with the TRIPOD reporting checklist (available at http://dx.doi.org/10.21037/tau-20-1348).

Methods

Study design and dataset information

The work flow of our study is shown in . The expression profile and clinical data of 539 ccRCC patients in TCGA-KIRC dataset were downloaded from TCGA portal (online URL: https://cdn.amegroups.cn/static/public/tau-20-1348-1.pdf) and 91 RCC patients from ICGC database were downloaded from ICGC portal (online URL: https://dcc.icgc.org/). Immune-related genes list was downloaded from the ImmPort database (online URL: https://www.immport.org/home/). All data were preprocessed in R software (online URL: https://www.r-project.org/; version 3.6.0;). 518 patients in TCGA-KIRC cohort and 91 patients in ICGC cohort with clinical information (Table S1 and S2) were screened for subsequent analyses. All procedures performed in this study were in accordance with the Declaration of Helsinki (as revised in 2013).

Figure 1

The workflow of our research project. Differential analysis was conducted in TCGA-KIRC dataset to obtain differentially expressed with the following cut-off: adjusted P value <0.05 and absolute log2FC >1.5; 382 immune-related genes were identified after taking intersection of the lists of DEGs and IRGs; 4 IRDEGs (CRABP2, LTB4R, PTGER1 and TEK) were finally identified and used to constructed a prognostic model after univariate cox and LASSO cox analyses. Further analyses were conducted to validate the robustness of model and explored the potential mechanism.

Identification of IRDEGs in TCGA-KIRC dataset

Differential analysis was conducted in TCGA-KIRC dataset through limma package (11,12), with the following cut-off: adjusted P value <0.05 and absolute log2FC >1.5. The differentially expressed genes list and immune-related genes list from the ImmPort database were uploaded into the Venn diagram online software (online URL: http://bioinformatics.psb.ugent.be/webtools/Venn/) to obtain the IRDEGs. The heatmap of IRDEGs expression was performed by pheatmap package.

Construction and validation of the risk model

518 patients in TCGA-KIRC dataset were included as a training set while 91 patients in the ICGC database were assigned as a validation set. Univariate Cox proportional hazards regression analysis was applied to identify the significant prognostic factors associated with OS, and Lasso regression was used to exclude overfitting genes. The candidate genes were analyzed in a multivariate Cox proportional hazards regression analysis to estimate their relative contributions to survival prediction. Subsequently, a prognostic model was constructed: risk score = expression of gene1 × β1 + expression of gene2 × β2 +…… + expression of genen × βn (13,14). According to the median risk score, patients were divided into two groups (favorable-risk group and unfavorable-risk group), and we applied the Kaplan-Meier and log-rank methods to test whether the survival distribution of different groups was equal. Receiver operating characteristic (ROC) curves were used to assess the predictive value of the risk model according to the areas under the respective ROC curves (AUCs). Time-dependent ROC curve analysis was conducted by using the survival ROC package (15).

RNA extraction and qRT-PCR analysis

Total RNA of 35 pairs of ccRCC and normal tissues RNA were extracted using a Trizol reagent, and 500 ng of RNA was used to synthesize cDNA, and qRT-PCR was performed on ABI system. The primer sequences are listed in Table S3.

Estimation of TME infiltration between groups

CIBERSORT, a deconvolution algorithm to characterize different cell compositions of the samples based on the immune gene signature sets, including 547 genes and 22 immune cell subtypes (16). We downloaded the result of 518 patients in TCGA cohort calculated by CIBERSORT algorithm from TIMER 2.0(Online URL: http://timer.cistrome.org/) (17) to estimate the infiltration of 22 different immune cell subtypes in the TME for further investigation of the composition and difference between favorable-risk and unfavorable-risk group. Each sample had been calculated a proportion in each cell subtype to estimate the relative abundance of TME immune infiltrating cells. Wilcox test was applied to compare the infiltration proportion of the 22 cell types between the unfavorable-risk group and the favorable-risk group.

Construction and validation of the clinical nomogram

We constructed a nomogram, which was widely used to predict the survival probability of patients in clinical (18), with the incorporation of age, gender, TNM stage, Fuhrman grade, necrosis and risk score through R rms package. We also used ROC curves to evaluate the predictive performance of nomogram at 1-, 3- and 5-year. In addition, calibration curves were used to evaluate the accuracy of the predicted survival time for 3- and 5-year OS, and decision curve analysis (DCA) was performed to evaluate the clinical application benefit between different variables.

Statistical analysis

All statistical analyses were performed by R software with the cut-off of P<0.05. Univariate and multivariate Cox proportional hazards regression analyses were used for identifying of prognosis-related IRDEGs and independent prognostic factors. Spearman correlation test was applied to analyze the correlation between the risk score and the expression of immune checkpoint genes. Survival data were calculated using the Kaplan-Meier method and the log-rank test. The relative expression level of four IRDEGs was analyzed by paired t test.

Results

Identification of IRDEGs in ccRCC

After differential analysis in TCGA-KIRC cohort, 3,169 DEGs were detected (adjusted P value <0.05 and absolute log2FC >1.5), among which 1,635 genes were upregulated, and 1534 genes were downregulated (). Taking the intersection of DEGs and the immune-related genes list, 382 IRDEGs were identified (), with 253 IRDEGs upregulated and 129 IRDEGs downregulated. The heatmap to visualize the expression of 382 IRDEGs in normal samples and tumor samples is shown in , and the result of differential analysis are shown in https://cdn.amegroups.cn/static/public/tau-20-1348-2.xlsx and Table S4.

Figure 2

Differential analysis of TCGA-KIRC cohort and LASSO regression. (A) The volcano plot of TCGA-KIRC cohort. Red plots represent upregulated genes while blue plots represent downregulated genes both with adj.P <0.05. (B) Venn diagram of DEGs and IRGs. (C) The heatmap of 382 IRDEGs in ccRCC and normal samples. Each column represents one sample and each row represents one gene. The gradual color ranging from blue to red represents the changing process from down to up regulation. (D) Plot of LASSO coefficient profiles. (E) Plot of partial likelihood deviance for the 382 IRDEGs.

Construction of prognostic model in TCGA cohort

175 IRDEGs was calculated to be significantly associated with OS after univariate Cox regression analysis (P<0.05) (Table S5). Lasso regression was used to filter genes to obtain 8 candidate genes (), which were subsequently included in the multivariate Cox regression analysis. A prognostic gene signature consisting of 4 genes was ultimately constructed with the P value <0.05 in multivariate Cox regression analysis ( and Table S6). Among these 4 genes, TEK was identified as a protective gene because of its hazard ratios (HR) value of <1, while CRABP2, LTB4R and PTGER1 were considered to be predictive genes of poor prognosis. Based on the analysis result, we constructed a computational formula: risk score = (0.074× expression level of CRABP2) + (0.165× expression level of LTB4R) + (0.052× expression level of PTGER1) + (-0.203× expression level of TEK), and the expression level was obtained by the log2-transformed FPKM+1 of each gene. Subsequently, a total of 518 patients in TCGA cohort were divided into two groups (unfavorable-risk group and favorable-risk group) according to the median risk score. shows the distribution of risk scores, patient survival status and the four gene expression levels in the 518 patients, which were sorted by the risk score of the four-gene signature. As the risk score increasing, the expression of the 4 IRDEGs also changed accordingly, and the prognosis of patients also became worse. Besides, as and 3C show, there were obvious differences in both OS and PFS between the two groups (P<0.0001). On the other hand, a time-dependent ROC was used to assess the prognostic value of the four-gene signature in the training set. The AUCs of the signature were respectively 0.744, 0.734, and 0.753 for the 1-, 3- and 5-year OS () while for the 1-, 3- and 5-year PFS were 0.696, 0.711, and 0.734(), indicating our risk model had a good performance on predicting prognosis. Besides, further analysis showed risk score was an independent prognostic factor (HR: 3.137; 95% CI: 2.383-4.131; P<0.001) ().

Table 1

Univariate analysis and Multivariate analysis of the 4 IRGs and signature

	Univariate analysis			Multivariate analysis
	HR (95% CI)	P value		HR (95% CI)	P value	Coefficient
Gene symbol
CRABP2	1.077 (1.001–1.158)	0.0508		0.672	0.0210	−0.3979
LTB4R	1.180 (1.001–1.390)	0.0029		1.217	0.0300	0.1961
PTGER1	1.054 (1.008–1.102)	P<0.001		1.199	0.0200	0.1816
TEK	0.639 (0.711–0.936)	P<0.001		0.754	0.0150	−0.2817
Factors
Age	1.030 (1.016–1.043)	P<0.001		1.031(1.016–1.045)	P<0.001
Stage		P<0.001
I	Reference			Reference
II	1.273 (0.683–2.371)			1.230 (0.657–2.304)	0.518
III	2.656 (1.755–4.019			1.937 (1.272–2.952)	0.002
IV	6.685 (4.537–9.851)			5.239 (3.524–7.787)	P<0.001
Four-gene signature				3.137 (2.383–4.131)	P<0.001

Figure 3

The 4-IRG prognostic signature in the ccRCC patients (TCGA cohort). (A) From top to bottom are the risk score distribution, the patients’ survival status distribution, and the heatmap of the 4 genes for low and unfavorable-risk groups, in which each column represents one sample and each row represents one gene. (B, C) The Kaplan-Meier curves of OS and PFS for low and unfavorable-risk groups. (D, E) The ROC curves for predicting OS and PFS in training set by the risk score and the AUC of 1, 3, 5-year.

External validation of prognostic model in ICGC cohort

We used the data from ICGC database as an external validation (). The distribution of risk scores and survival status of patients, as well as the expression level of the 4 IRDEGs in ICGC cohort, were shown in and it was observed that the expression had an obvious difference between unfavorable-risk group and favorable-risk group. Besides, Kaplan-Meier analysis indicated that unfavorable-risk group was significantly associated with a poor prognosis (P<0.05), consistently with the above results (). And the AUCs of ROC analysis were respectively 0.635, 0.638 and 0.635 at 1-, 3- and 5-year OS value, indicating the stability of risk model in different cohorts ().

Figure 4

The 4-IRG prognostic signature in the ICGC cohort. (A) From top to bottom are the risk score distribution, the patients’ survival status distribution, and the heatmap of the 4 genes for low and unfavorable-risk groups (B) The Kaplan-Meier curves for low and unfavorable-risk groups. (C) The ROC curves for predicting OS in validation set by the risk score and the AUC of 1, 3, 5-year.

Validation in clinical samples

We used 35 pairs of ccRCC (Table S7) and normal tissues to detect the expression of the four IRDEGs. The results showed LTB4R was high expressed in ccRCC compared with normal tissues, and CRABP2, PTGER1, TEK were low expressed, consistent with the expression data of TCGA database ().

Figure 5

Relative expression level of the 4 IRGs in paired ccRCC and normal tissues.

TME immune cell infiltration analysis and immune checkpoints analysis

We summarized the result of 518 ccRCC patients calculated by CIBERSORT algorithm (Figure S1) and compared all the immune cell subtypes in two groups (Table S8). Infiltration proportion of partial cell subtypes have an obvious difference between two groups, among which mainly CD8 T cell, follicular helper T cell, regulatory T cell, activated NK cell, M0 Macrophage have a higher infiltration proportion in the unfavorable-risk group, while M1 Macrophage, M2 Macrophage and other cell types have a lower proportion (). We further explored the expression of the T cell exhaustion-related markers and immunomodulators (PD-1, CTLA4, TNFRSF9, TIGIT, LAG3) in two groups and found all markers in the unfavorable-risk group were upregulated, indicating an immunosuppressive and exhausted phenotype in the unfavorable-risk group (). Subsequent correlation analysis also showed a positive correlation between risk score and the above markers (). Based on the above analyses, we found two groups had a significant distinct pattern of immune infiltration, which may lead to different survival benefits.

Figure 6

TME immune cell infiltration characteristics of 22 immune cell subtypes in unfavorable-risk and favorable-risk groups. (A) The violin plot of the abundance of immune cell subtypes in two groups. The asterisks on the top represented the P value of Wilcoxon test (ns P>0.05; *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001). (B) The violin plot of the expression of T cell exhaustion-related markers and common immune checkpoint in two groups. (C) Spearman correlation analysis of risk score and above markers.

Construction and validation of the nomogram

We constructed a nomogram containing age, gender, TNM stage, Fuhrman grade, necrosis and risk score to predict 3- and 5-year survival probability. Each variable had a corresponding score (Table S9), and an overall score could be finally calculated to predict the survival probability at the corresponding time (). To validate the performance of nomogram, we conducted the ROC analysis, and the result showed respective AUCs were 0.811 and 0.795 in the TCGA cohort (). The calibration curves showed good consistency between the actual and predicted outcomes of 3- and 5-year OS (). Decision curve analysis (DCA) was also conducted, and all variable curves were above the two solid curves. The curve of the nomogram was above the curve of risk score at 3- and 5-year (), indicating nomogram had a better clinical net benefit. We also constructed another nomogram containing age, gender, TNM stage, and risk score. And the AUCs were 0.811, 0.786 while in the validation set were 0.728, 0.713 at 3-, 5-year, suggesting its stability and effectiveness (Figures S2, S3). Calibration curves and DCA also showed the robustness of our nomogram.

Figure 7

Construction and validation of the nomogram in TCGA-KIRC cohort. (A) The clinical nomogram of the ccRCC patients (TCGA cohort). (B) The ROC curves of nomogram for predicting OS in training set and the AUC of 1, 3, 5-year. (C) The 3-year and 5-year calibration curves of the nomogram in training set. (D) The 3-year and 5-year DCA plots of the nomogram. The grey, red dotted lines respectively represent net benefits of nomogram and risk score at different threshold probabilities.

Discussion

Immune related genes (IRGs) play an important role in tumor immune infiltration as well as tumor progression in ccRCC (19,20) and strongly influence complicate soluble factors secretion, which correlate with therapeutic response and clinical outcome (21). IRGs based prognostic model have been successfully developed for hepatocellular, colorectal, lung and bladder cancers (22-25). In this study, we screened and validated CRABP2, LTB4R, PTGER1 and TEK from TCGA database as potent IRGs to predict the survival risk in ccRCC patients. Multivariate Cox proportional hazards models were constructed to stratify patients based on the 4 genes signature. Laboratory q-PCR confirmed LTB4R upregulated and CRABP2, PTGER1, TEK downregulated in 35 pairs of ccRCC and normal tissues.

LTB4R is a receptor of Leukotriene B4 and is expressed mainly in leukocytes like granulocytes, macrophages and eosinophils (26). Several studies implied it was involved in CD8 T cells recruiting (27,28). The neutrophilic influx induced by LTB4 increases the pro-tumorigenic activity of tumor-associated neutrophils through releasing reactive oxygen species, inflammatory cytokines and injuring innate immune response (29,30). CRABP2 was responding in retinoic acid (RA) transduction as a tumor suppressive pathway (31). However, artificially overexpressing CRABP2 in Caki-2 cells did not exhibit a significant change in RA sensitivity. Our data indicated CRABP2 was lowly expressed in ccRCC samples, which was consistent with previous study (32). Although the exactly role of CRABP2 in RCC is not clear yet, our data showed high CARBP2 expression was an independent predictor factor for worse prognosis. Further investigations are necessary to define other molecules involved in CARBP2 mediated RA signaling and metabolism in RCC. PTGER1 is one of the receptors of prostaglandin and it couples with G-proteins to activate protein kinase C (33). Previous study implied blocking PTGER1 could suppress immunosuppressive function of Treg and subsequently inhibit tumor growth in colon cancer (34). In ccRCC, our result also exhibited high PTGER1 expression was correlated with worse prognosis. TEK encodes Tie2, which cooperate with VEGFs as critical regulators of vascular development (35). Actually, mechanisms of angiogenesis are extremely complex depending on different tumor stage and content, therefore, the prediction role of Tie2 is inconsistent among different tumor types (36). Increasing Tie2 expression correlated with high metastasis risk and poor survival among breast cancer and glioblastoma patients (37,38). However, in ccRCC, our data implied downregulation of TEK associated with a poor prognosis which was also demonstrated in previous studies (39,40). Low expression of TEK has been noted in aggressive ccRCC for years (41). Recently, when compared gene prevalence between nonmetastatic and metastatic ccRCC by target next generation sequencing from 106 sporadic cases, higher frequencies of TEK mutations involved in metastatic cohort (42). Since Tie2 signaling influences vascular permeability, low expression of Tie2 may potentiate inflammatory cells migration into tumor microenvironment (43). Inflammatory cytokines such as TNF-α, IL-6, CXCL8 induces a more aggressive tumor phenotype via immune surveillance and form premetastatic niche. Overall, we assume TEK is a tumor suppressor in ccRCC, further studies are needed in future.

Solid tumors usually disrupt tumor target immune response by subvert immune surveillance. Tumor immune signature is highly correlated with tumor prognosis and response to immunotherapy. When using a 34-gene expression signature, ccRCC can be characterized into high angiogenesis tumor with improved prognosis or high immunocytes tumor with worse survival (44). Tumor infiltrating lymphocytes in ccRCC were analyzed by gene expression and cytometry phenotyping, the result implied more poorly cytotoxic CD8 T cell, Treg infiltrate in unfavorable risk group (45). In this current study, we also investigated high proportion of CD8 T cell, Treg, NK cell in unfavorable risk group through CIBERSORT algorithm, which indicated that our novel prediction model could properly distinct patient into different immunological features. Besides, a group of immunomodulators (PD-1, CTLA4, TNFRSF9, TIGIT and LAG3) was significantly correlated with our risk score, and confirmed our risk model stably stratified patients from immune evasion perspective.

In order to increase prediction accuracy, we developed a clinical nomogram with age, gender, TNM stage, Fuhrman grade, necrosis and risk score. This nomogram obtained an AUC of 0.846, 0.811 and 0.795 in predicting the possibility of survival at 1-, 3- and 5-year respectively. As there is lack of Fuhrman grade, necrosis information in ICGC database, we removed these two factors from original nomogram for validation. The modified nomogram consistently obtained a relatively high AUC of 0.755, 0.728 and 0.713 in survival prediction at 1-, 3- and 5-year separately in ICGC data. Hence, based on the 4 immune related genes CRABP2, LTB4R, PTGER1 and TEK, we successfully constructed a prognostic risk model for ccRCC and externally validated its accuracy. Defective T-cells and aberrant expression suppressive immunomodulators lead tumor be more aggressive. Owing to data we obtained from public database, further independent validation in prospective studies is needed.

The article’s supplementary files as