A 14 immune-related gene signature predicts clinical outcomes of kidney renal clear cell carcinoma

Introduction

Kidney cancer is one of the most common urological tumors worldwide, with approximately 403,262 new cases and 175,098 deaths associated with this form of cancer in 2018 (Bray et al., 2018). Kidney renal clear cell carcinoma (KIRC) is the leading cause of kidney cancer mortality and the most common type, accounting for 85% of kidney cancers (Siegel, Miller & Jemal, 2017; Tseng, 2016). Considering the aging population of KIRC patients and the rising expense of treatment, KIRC is gradually becoming the focus of geriatric cancer (Tan, Filson & Litwin, 2015). Despite the rapid development of cancer treatments, mortality rates in KIRC remain stagnant. With the progression of next generation sequencing and data mining techniques, it is urgent that we explore prognostic biomarkers for KIRC, using molecular characteristics and tumor immune environments to guide patient therapy.

Over the past decade our understanding of immune components, including the impact of tumor microenvironments on patient survival and therapy response (Chen & Mellman, 2017; Grivennikov, Greten & Karin, 2010), has increased. Some studies have found that tumor-infiltrating immune cells are able to serve as either tumor suppressors or promoters in microenvironments. For example, CD8⁺ T cells have been associated with improved survival in cancer patients (Gajewski, Schreiber & Fu, 2013; Governa et al., 2017), while tumor associated macrophages and regulatory T cells demonstrate the ability to promote tumor development (Nishikawa & Sakaguchi, 2014; Noy & Pollard, 2014). Considering the complexity and significance of tumor immune microenvironments, it is imperative that we investigate immune-related biomarkers for KIRC patients. Recent studies have provided insight into the KIRC immune signature (Geissler et al., 2015; Şenbabaoğlu et al., 2016). Şenbabaoğlu et al. (2016) found mRNA signatures with the potential to be immunotherapeutic biomarkers in KIRC. However, their investigations do not include immune-related genes for analysis nor do they establish a systematic immune-related gene-risk signature for KIRC patients. Khadirnaikar et al. (2019) utilized immune associated lncRNA (long non-coding RNA) to construct prognostic subtypes in KIRC patients. Our immune clusters were more robust and independent. Additionally, they concentrated on lncRNA, not the overall immune-related genes. Smith et al. (2019) constructed endogenous retroviral signatures for KIRC patients, but they did not investigate the prognostic ability of the signature in different subtypes of patients. Therefore, it is essential that we explore a systematic prognostic signature based on tumor immune environments in KIRC.

In our study, we used RNA-seq data from The Cancer Genome Atlas (TCGA) to find immune-related genes with prognostic ability and to establish an immune-related risk signature for KIRC. To assess the clinical potential of the signature, we investigated the association between the signature, clinical parameters, and patient survival. Gene set enrichment analysis (GSEA) was performed to explore the molecular characteristics of the signature.

Materials & Methods

Statistical analysis

Boxplots were created with the R package, ggplot2 (v3.0.0). The R package ComplexHeatmap (v1.18.1) was used to create heatmaps (Gu, Eils & Schlesner, 2016). We counted C-index with R packages “survcomp” (Harrell Jr et al., 1982; Schröder et al., 2011). The student’s t-test was performed for statistical comparison. We chose R to conduct statistical analysis (https://www.r-project.org/). P-values lower than 0.05 were considered statistically significant. The main code of analysis was pushed to github (https://github.com/huchua/KIRC).

Results

Establishment and validation of the immune-related gene signature in KIRC

Figure 1 illustrates the workflow we used to develop an immune-related gene-risk signature. The immune-related gene signature was constructed within the KIRC training data set, while we applied the testing set to validate the signature. After 1,000 iterations, seven unique gene models were selected (Fig. 2A, Table S1). The model selected as the immune-related gene-risk signature consisted of 14 genes and ranked the highest in frequency 282 times. Parameter selection in the LASSO-cox model was log lambda = − 2.641 and alpha = 1. The univariate and multivariate Cox analysis of the 14 immune-related genes are shown in Tables 1 and 2. The principle component analysis of the 14 immune-related gene signature displayed a different distribution pattern between low- and high-risk groups when comparing the training and testing (Figs. 2B–2C). This indicated that the low- and high-risk groups had different immune phenotypes. In the training and testing, the c-index was 0.7862 and 0.6534, respectively (P < 0.001; Fig. 2D). In a time-dependent receiver operating characteristic curve (ROC) created for training and testing datasets, area under the curve (AUC) values at 1, 3, and 5 years were 0.679, 0.63, 0.627 and 0.65, 0.596, 0.568, respectively (Fig. S2). The 14 immune-related genes are AR, BID, BMP8A, CCL7, CCR10, FGF17, GDF1, IL20RB, IL4, KLRC2, LHB, SEMA3A, SEMA3G, and TXLNA. K-M (Kaplan–Meier) survival curves and gene expression of the 14 immune-related genes are shown in Fig. 3 and Fig. S3. Because of the differences between the log-rank test and the univariate cox analysis, the results of the univariate cox analysis of other genes except GDF1 are very consistent with the results of the K-M (Kaplan–Meier) survival curves.

Table 1:

Univariate Cox analysis of 14 immune relate genes in all cohort.

Variable	HR	95% CI	p value
KLRC2	20.657	9.116–46.81	<0.001
CCL7	2.757	2.090–3.637	<0.001
SEMA3A	1.748	1.437–2.126	<0.001
SEMA3G	0.617	0.534–0.714	<0.001
CCR10	1.759	1.381–2.241	<0.001
BMP8A	1.813	1.298–2.532	<0.001
FGF17	7.335	3.202–16.803	<0.001
GDF1	1.483	1.062–2.070	0.021
IL4	109.106	34.328–346.781	<0.001
LHB	5.053	3.240–7.881	<0.001
TXLNA	2.329	1.500–3.618	<0.001
AR	0.548	0.463–0.649	<0.001
IL20RB	1.225	1.145–1.310	<0.001
BID	3.505	2.439–5.035	<0.001

DOI: 10.7717/peerj.10183/table-1

Table 2:

Multivariate Cox analysis of 14 immune relate genes in all cohort.

AR	0.638	0.556–0.731	<0.001		IL20RB	1.268	1.102–1.46	0.001
age	1.352	1.143–1.599	<0.001		age	1.414	1.195–1.672	<0.001
stage	3.725	2.708–5.124	<0.001		stage	3.275	2.351–4.562	<0.001
gender	0.943	0.687–1.294	0.718		gender	1.025	0.744–1.412	0.879
BID	1.411	1.211–1.645	<0.001		IL4	1.458	1.294–1.643	<0.001
age	1.407	1.188–1.666	<0.001		age	1.39	1.171–1.651	<0.001
stage	3.173	2.284–4.41	<0.001		stage	3.426	2.484–4.725	<0.001
gender	1.039	0.755–1.43	0.814		gender	0.999	0.726–1.375	0.997
BMP8A	1.213	1.057–1.391	0.006		KLRC2	1.402	1.233–1.594	<0.001
age	1.454	1.229–1.721	<0.001		age	1.383	1.169–1.637	<0.001
stage	3.607	2.619–4.969	<0.001		stage	3.429	2.483–4.737	<0.001
gender	1.052	0.765–1.447	0.754		gender	0.958	0.694–1.322	0.792
CCL7	1.329	1.202–1.471	<0.001		LHB	1.285	1.168–1.414	<0.001
age	1.455	1.226–1.728	<0.001		age	1.393	1.178–1.647	<0.001
stage	3.485	2.529–4.802	<0.001		stage	3.519	2.55–4.854	<0.001
gender	1	0.726–1.377	1		gender	1.17	0.844–1.622	0.347
CCR10	1.235	1.111–1.372	<0.001		SEMA3A	1.214	1.107–1.33	<0.001
age	1.382	1.166–1.639	<0.001		age	1.445	1.218–1.715	<0.001
stage	3.844	2.793–5.292	<0.001		stage	3.549	2.572–4.898	<0.001
gender	1.034	0.751–1.424	0.836		gender	1.108	0.802–1.529	0.535
FGF17	1.235	1.112–1.372	<0.001		SEMA3G	0.652	0.549–0.775	<0.001
age	1.399	1.181–1.657	<0.001		age	1.379	1.164–1.633	<0.001
stage	3.722	2.705–5.121	<0.001		stage	3.29	2.384–4.539	<0.001
gender	1.109	0.806–1.525	0.527		gender	0.973	0.707–1.339	0.867
GDF1	1.086	0.975–1.211	0.135		TXLNA	1.285	1.098–1.505	0.002
age	1.411	1.191–1.672	<0.001		age	1.433	1.208–1.7	<0.001
stage	3.76	2.735–5.17	<0.001		stage	3.696	2.687–5.085	<0.001
gender	1.103	0.796–1.529	0.555		gender	1.212	0.87–1.688	0.255

DOI: 10.7717/peerj.10183/table-2

Discussion

Our study used the TCGA database and immune-related genes to establish a KIRC signature consisting of 14 immune-related genes. We found that patients in the high-risk group showed a positive association with M stage, T stage, and advanced pathological staging. Additionally, the signature exhibited strong prognostic abilities and independently predicted KIRC patient prognosis. Functional analysis highlighted the significance of our signature based on its involvement in many important pathways.

In our investigation, we used a total of 14 immune-related genes to construct our signature. In the signature, CCL7 increased the peripheral blood mononuclear cell recruitment in renal cell cancer through the inhibition of let-7d (Riihimaki et al., 2014). Wyler et al. (2014) demonstrated the ability of CCL7 to recruit monocytes through CCR2, promoting renal cell cancer metastasis to the brain. This indicates that CCL7 is a major factor in the development of KIRC in tumor immune microenvironments and may be a potential immunotherapy target. Fibroblast growth factor receptor 17 (FGF17) is another immune-related gene in our signature. FGF17 demonstrates a variety of functions in cancer development. Gauglhofer et al. (2011) found that FGF17 was involved in the paracrine and autocrine signaling of hepatocellular carcinoma and promoted the neoangiogenesis of hepatocellular carcinoma. Heer et al. (2004) showed that FGF17 was overexpressed in prostate cancer and participated in prostate carcinogenesis. Our study showed that FGF17 plays an important role in KIRC based on tumor immunology. However, the underlying mechanism of FGF17 in KIRC immune microenvironments requires further investigation. In our signature, IL-4 is another well-studied immune-related biomarker for KIRC. IL-4, released by immune cells, controls the expression of B7-H1, thus altering T cell responses in KIRC (Quandt et al., 2014). The investigation conducted by Chang et al. (2015) demonstrated that IL-4 expression can predict KIRC patient recurrence and survival outcomes. Collectively, these investigations reinforce the significance of the 14 immune-related gene-risk signature in KIRC.

Our signature is associated with the survival outcome of KIRC patients and clinical parameters, including pathological staging, M stage, and T stage. No correlation was found between the signature and N-stage, possibly due to a lack of N-stage information for many patients. According to the immune-related gene-risk signature, our study found that clinical cohorts in KIRC have different immune-related risk factors. This signature also reflects differences in tumor immune microenvironments and predicts survival outcomes in KIRC patients; thus, demonstrating the clinical significance of our signature and its possible use as a survival predictor in KIRC.

Despite of the relatively low C-index of our testing set (0.6534), the C-index of the testing set falls within a similar range as those of other related studies. For instance, in the study of Bailiang Li et al., they developed an immune signature in non-small cell lung cancer with the C-index of 0.64 (Li et al., 2017). Besides, our signature achieved a similar C-index of testing set with the immune signature in ovarian cancer (0.625) (Shen et al., 2019). More importantly, the C-index of our testing set showed a similar accuracy with the C-index of clinical staging systems in renal cancer (0.62) (Qu et al., 2018).

Our study expands on the signatures association with several important pathways, especially the metabolism pathway. This may reflect the mutual interaction between tumor metabolism and tumor immunology in KIRC. Pearce et al. (2013) showed that metabolic reprogramming can influence the fate and function of T cells in tumors. Our study further indicates the importance of metabolism pathways in KIRC immune microenvironments.

We acknowledge the limitation of our evaluation scheme, including randomly dividing the full data set into training and testing data set, which result in the inherent bias for the specific study. Another limitation of our research is the lack of independent validation data sets in our evaluation scheme. Therefore, our study should be further validated through a prospective cohort data to further illustrate the robustness of the model.

Conclusions

This investigation utilized RNA-seq data from the TCGA database to construct a 14 immune-related gene-risk signature with the ability to independently predict survival outcomes in KIRC; thus, providing novel clinical applications and possible immune targets for KIRC.

Supplemental Information