Research Paper Volume 12, Issue 9 pp 8484—8505
Identification of transforming growth factor beta induced (TGFBI) as an immune-related prognostic factor in clear cell renal cell carcinoma (ccRCC)
- 1 Department of Urology, Zhongnan Hospital of Wuhan University, Wuhan 430071, China
Received: October 29, 2019 Accepted: April 16, 2020 Published: May 14, 2020https://doi.org/10.18632/aging.103153
How to Cite
Copyright © 2020 Du et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Clear cell renal cell carcinoma (ccRCC) is the most common subtype among kidney cancer, which has poor prognosis. The aim of this study was to screen out novel prognostic biomarkers and therapeutic targets for immunotherapy, and some novel molecule drugs for ccRCC treatment. Immune scores ranged from -1109.36 to 2920.81 and stromal scores ranged from -1530.11 to 1955.39 were firstly calculated by applying ESTIMATE algorithm. Then 17 DEGs associated with immune score and stromal score were further identified. 6 candidate hub genes were screened out by performing overall survival (OS) and disease-free survival analyses based on TCGA-KIRC data, one of which including TGFBI was further regarded as hub gene associated with prognosis by calculating the R2 (R2 = 0.011, P = 0.018) and AUC (AUC = 0.874). The prognostic value of TGFBI was validated by performing OS, CSS, and PFS analyses based on GSE29609 and E-MTAB-3267. CMap analysis suggested that 3 molecule drugs might be novel choice for ccRCC treatment. Further analysis demonstrated that CNVs of TGFBI was associated with OS of patients with ccRCC. TGFBI expression was also correlated with histologic grade, pathologic stage, and immune infiltration level, significantly. TGFBI was the most relevant gene with OS among the candidate hub genes, which might be novel DNA methylation biomarkers for ccRCC. In conclusion, our findings indicated that TGFBI was correlated with prognosis of patients with ccRCC, which might be novel prognostic biomarkers, and targets for immunotherapy in ccRCC. Three small molecule drugs were also identified, which showed strong potential for ccRCC treatment.
Immunotherapy is one of the treatment methods for malignant tumors at present, which mainly uses the immune effects of autoimmune or alloimmune cells in patients to improve the symptoms, prolong the survival and improve the prognosis [1, 2]. In recent years, immunotherapy has become a novel treatment for cancers, whose effectiveness and safety have been gradually conformed . With the development of precision medicine and immunotherapy for cancers, nowadays, more and more researchers focus on finding out more accurate therapeutic targets for immune treatment [4, 5].
Kidney cancer is one of the most common malignances around the world, which is with poor prognosis . According to recent statistics from the International Agency for Research on Cancer (IARC), part of the World Health Organization (WHO), there were 403,262 new cases of kidney cancer and 175,098 associated deaths worldwide in 2018. Renal cell carcinoma (RCC) accounted for approximately 90% of kidney cancers, the commonest histological subtype of which was clear cell renal cell carcinoma (ccRCC) [7, 8]. Although surgical treatment was the most effective therapy for localized ccRCC, there was a lack of drugs for adjuvant treatment . What was worse, there was no effective treatment method for advanced ccRCC . Thus, in this study, we tried to find out some novel prognostic biomarkers, which might be novel targets for immunotherapy.
For the first time, in this study, we firstly calculated immune score and stromal score of each case from TCGA-KIRC data, by applying Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) algorithm (a method provided by Yoshihara et al.) . Then we screened out 17 differentially expressed genes (DEGs) associated with immune score and stromal score. Based on these DEGs, 3 small molecule drugs were obtained, which showed strong potential for ccRCC treatment. Finally, transforming growth factor beta induced (TGFBI) was screened out by using four kinds of survival analyses and two independent datasets, which were significantly associated with prognosis of patients with ccRCC. 9 different kinds analyses were further performed to explore the potential value of TGFBI in various aspects.
In conclusion, our finding indicated that TGFBI had great effects for assessing prognosis of patients with ccRCC, which might be a novel prognostic biomarker and target for immunotherapy. Moreover, three molecule drugs were screened out, which might be novel choice for clinicians for ccRCC treatment.
Immune scores and stromal scores were correlated with clinical features of patients with ccRCC
Among all the 530 ccRCC downloaded from TCGA database, 35.1% (n = 186) samples were female, 64.9% (n = 344) samples were male. As for the neoplasm histologic grade, patients with Gx, G1, G2, G3, G4 grade accounted for 0.9% (n = 5), 2.7% (n = 14), 43.1% (n = 227), 39.1% (n = 206), 14.2% (n = 75), respectively. Tumors on the left accounted for 47.1% (n = 249) and the right side accounted for 52.9% (n = 280). Pathologic stage included 265 (50.0%) patients of stage I, 57 (10.8%) of stage II, 123 (23.2%) patients of stage III, and 83 (15.7%) patients of stage IV. When talking about person neoplasm cancer status, patients with tumor free accounted for 68.9% (n = 354) of the total number, and 31.1% (n = 160) cases were with tumor. After calculating immune score and stromal score of each ccRCC from TCGA-KIRC, immune scores ranged from -1109.36 to 2920.81 meanwhile stromal scores ranged from -1530.11 to 1955.39 as the result suggested (Supplementary Table 1).
Further analysis demonstrated that immune score was associated with gender (t = -2.220, P = 0.027, Figure 1A), neoplasm histologic grade (F = 9.470, P < 0.001, Figure 1B), pathologic stage (F = 7.390, P = 0.009, Figure 1C), significantly. Moreover, stromal score was significantly related to neoplasm histologic grade (F = 3.020, P = 0.038, Figure 1G), and pathologic stage (F = 226.080, P < 0.001, Figure 1I).
Figure 1. Distribution of immune scores of gender (A), neoplasm histologic grade (B), laterality (C), pathologic stage (D), and person neoplasm cancer status (E). Distribution of stromal scores of gender (F), neoplasm histologic grade (G), laterality (H), pathologic stage (I), and person neoplasm cancer status (J).
Immune scores and stromal scores were correlated with OS
To explore the relationship between immune score (or stromal score) and survival, we performed survival analysis in this part. As Figure 2A showed, 1179.98 was set as the optimal cutoff for grouping. ccRCC patients were divided into high immune score group (n = 234) and low immune score group (n = 296). The result suggested that ccRCC patients with high immune score had worse OS compared with these with low immune score (P = 0.0024, Figure 2B). Meanwhile, an illustration of optimal cutoff identification for stromal score is shown in Figure 2D (stromal score cutoff = 794.10). ccRCC patients were divided into two groups (high stromal score group: n = 208; low stromal group: n = 322). Similarly, low stromal score of cases was associated with better OS (P = 0.0250, Figure 2C). We also explored the correlation between scores and DFS, unfortunately, there was no significant relationship between them, as Figure 2C, 2F suggested.
Figure 2. (A) An illustration of optimal cutoff identification for immune score. Survival analysis of the association between immune score and overall survival (B), disease-free survival (C) time in ccRCC. (D) An illustration of optimal cutoff identification for stromal score. Survival analysis of the association between stromal score and overall survival (E), disease-free survival (F) time in ccRCC.
17 DEG associated with immune score and stromal score screening
Based on “limma” in R software, 337 DEGs associated with immune score were screened out, including 334 up-regulated DEGs and 3 down-regulated DEGs (Figure 3A, 3C). Adjust P value and log2FC of each immune-related DEG were showed in Supplementary Table 2, in detail. Furthermore, 218 DEGs (204 up-regulated and 14 down-regulated) related to stromal score were picked out, accurately. Supplementary Table 3 showed the detailed information of each stromal-related DEG. Finally, 17 DEGs overlapped in immune-related DEGs and stromal-related DEGs were identified for further analysis (Figure 3E).
Figure 3. Differentially expressed genes (DEGs) analysis in ccRCC. (A) Volcano plot visualizing the immune-related DEGs. (B) Volcano plot visualizing the stromal-related DEGs. (C) Heatmap of immune scores of high score vs low score (P < 0.05, fold change > 1). (D) Heatmap of stromal scores of high score vs low score (P < 0.05, fold change > 1). (E) Identification of common DEGs between immune-related DEGs and stromal-related DEGs.
Function and pathway enrichment analysis
The results of GO analysis demonstrated that 17 DEGs were significantly enriched in 114 BPs (Supplementary Table 4), 4 CCs (Figure 4B), and 8 MFs (Figure 4C). The top 10 enriched BPs were regulation of cell-cell adhesion, regulation of immune effector process, leukocyte apoptotic process, positive regulation of cell-cell adhesion, negative regulation of cell adhesion, regulation of leukocyte cell-cell adhesion, regulation of T cell activation, leukocyte cell-cell adhesion, negative regulation of leukocyte apoptotic process, and protein activation cascade (Figure 4A). Moreover, DEGs were significantly correlated with four KEGG pathways including complement and coagulation cascades, cytokine-cytokine receptor interaction, staphylococcus aureus infection, and viral protein interaction with cytokine and cytokine receptor suggested by Figure 4D.
3 small molecule drugs might be novel choices for ccRCC treatment
The highly associated molecule drugs were identified by using CMap. Totally 6 molecule drugs were screened out (Table 1). Among them, three small molecule drugs including vincamine (mean = -0.493, n = 6, P < 0.001), clenbuterol (mean = 0.556, n = 5, P = 0.004), betazole (mean = 0.422, n = 5, P = 0.039) showed strong potential to treat ccRCC.
Table 1. Results of CMap analysis based on DEGs in ccRCC.
|cmap name||mean||n||enrichment||p||specificity||% non-null|
Identification of 6 candidate hub genes
The result of OS analysis demonstrated that 7 DEG including C1R (HR = 1.900, P = 2.1E-05, Figure 5A), C1S (HR = 2.100, P = 1.5E-06, Figure 5B), IGLL5 (HR = 1.500, P = 0.016, Figure 5H), MMP7 (HR = 1.400, P = 0.031, Figure 5L), SERPINF1 (HR = 1.500, P = 0.006, Figure 5N), SLC38A5 (HR = 2.100, P = 1.5E-06, Figure 5O), and TGFBI (HR = 1.600, P = 0.002, Figure 5P) were associated with OS of patients with ccRCC. In addition, high expressions of C1R (HR = 2.300, P = 7.9E-06, Figure 6A), C1S (HR = 2.300, P = 6E-06, Figure 6B), CP (HR = 1.500, P = 0.037, Figure 6F), MMP7 (HR = 1.700, P = 0.003, Figure 6L), PRIMA1 (HR = 1.600, P = 0.007, Figure 6M), SERPINF1 (HR = 1.800, P = 0.002, Figure 6N), SLC38A5 (HR = 1.800, P = 0.002, Figure 6O), and TGFBI (HR = 1.800, P = 0.001, Figure 6P) were significantly correlated with DFS of ccRCC patients. Finally, 6 genes including C1R (complement C1r), C1S (complement C1s), MMP7 (matrix metallopeptidase 7), SERPINF1 (serpin family F member 1), SLC38A5 (solute carrier family 38 member 5), and TGFBI (transforming growth factor beta induced) were selected as candidate hub genes for subsequent analysis.
Figure 5. Overall survival analyses on DEGs based on the TCGA-KIRC data. (A) C1R. (B) C1S. (C) CCL19. (D) CCL21. (E) CD163. (F) CP. (G) FGG. (H) IGLL5. (I) IL2RA. (J) IL7R. (K) JCHAIN. (L) MMP7. (M) PRIMA1. (N) SERPINF1. (O) SLC38A5. (P) TGFBI. (Q) VSIG4.
TGFBI was identified as hub gene
We firstly calculated the R2 to evaluate the relationship between candidate hub genes and overall survival (OS) days. The result demonstrated that all the six genes (including C1R (R2 = 0.014, P = 0.007, Figure 7A), C1S (R2 = 0.013, P = 0.010, Figure 7B), SERPINF1 (R2 = 0.012, P = 0.012, Figure 7D), SLC38A5 (R2 = 0.019, P = 0.001, Figure 7E), and TGFBI (R2 = 0.011, P = 0.018, Figure 7F)) except MMP7 (R2 = 0.003, P = 0.229, Figure 7C) were negatively associated with OS days. Then we calculated AUC for each candidate hub genes. Only TGFBI reached the standard of AUC ≥ 0.80 (AUC = 0.874, Figure 7F), which was regarded as hub gene in this study.
Patient living days gradually decreased with an increasing TGFBI expression
We firstly compared the mRNA expression levels of TGFBI between tumors and normal tissues, the result suggested that expressions of TGFBI in tumors were significantly higher than these in normal tissues (Figure 8A). Moreover, Expression of TGFBI (F = 6.34, P = 3.17E-04, Figure 8B) was significantly associated with tumor stage. High expression of TGFBI always related to higher tumor stage, as Figure 8B showed.
Figure 8. (A) Expression comparison of hub gene in ccRCCs and normal tissues by GEPIA. (B) Stage plot of TGFBI across different pathological stages in the TCGA-KIRC data. (C) Survival analysis of the association between TGFBI expression and progression free survival time in ccRCC (E-MTAB-3267). Survival analysis of the association between TGFBI expression and overall survival (D), cancer specific survival (E), progression free survival (F) time in ccRCC (GSE29609).
To validate the prognostic value of TGFBI, we performed three different survival analyses (OS, CSS, PFS) for all the candidate hub genes based on GSE29609. As shown in Figure 8D, high expressions of TGFBI (P = 0.030) were associated with worse OS. Meanwhile, there was a trend that patients with high expression of TGFBI had worse CSS (P = 0.066, Figure 8E) and PFS (P = 0.062 suggested by E-MTAB-3267 (Figure 8C), P = 0.054 suggested by GSE29609 (Figure 8F)). Summary above, we thought higher TGFBI expression was significantly associated with a shorter survival time, which showed great prognostic value in ccRCC.
Hub gene validation
Based on ccRCC data from Oncomine database, TGFBI expression comparison between tumors and normal tissues across 7 analyses were identified. The result demonstrated that the mRNA expression of TGFBI (P = 9.89E-05, Figure 9A, 9B) was lower in normal tissues compared with ccRCCs. Furthermore, we explored the translational-level expression of hub genes. As Figure 9C, 9D showed, the translational-level expressions of TGFBI were higher in ccRCCs compared with normal tissues.
mRNA expressions of TGFBI in cancer cell lines
Based on GEPIA and CCLE, we explored hub genes expressions in all types of cancer cell lines. The expressions of TGFBI were different in these types of cancers. The expression of TGFBI was higher in GBM, HNSC, KIRC, KIRP, LGG, PAAD, SARC, STAD compared with corresponding normal tissues (Figure 10A). Moreover, we found that not only the mRNA expression (Figure 10B) but also the CNV level (Figure 10C) of TGFBI ranked as the first highest of all types of cancer cell lines. All the results made the hub genes reliable.
Genetical alteration of TGFBI
According to the result, TGFBI altered in 75 (14%) of 533 ccRCC patients (Figure 11A). And the main type was amplification (Figure 11B). Meanwhile, there was no differences of TGFBI expression between cases with CNVs and cases without CNVs (Figure 11C). When talking about the association between CNVs of TGFBI and survival of patients with ccRCC, there was a trend that amplification of TGFBI caused better OS (Figure 11D) and DFS (Figure 11E) of patients with ccRCC.
Figure 11. A summary of mutations and CNVs of hub genes. (A) Genetic alterations associated with hub genes and expression heatmap of hub genes based on the data from TCGA. (B) The total alteration frequency of TGFBI in TCGA-ACC is illustrated. Correlation between different CNV patterns and mRNA expression levels of TGFBI (C), respectively. (D, E) Survival analysis of ccRCC patients with CNVs of TGFBI based on TCGA-ccRCC data.
Associations between TGFBI expression and clinical factors
Cases with complete clinicopathological data from TCGA database were included for this analysis. The result demonstrated that TGFBI expression was significantly associated with stromal score (P < 0.001), immune score (P < 0.001), gender (P < 0.001), neoplasm histologic grade (P < 0.001), pathologic stage (P = 0.011), and person neoplasm cancer status (P < 0.001) (Table 2).
Table 2. Associations between TGFBI expression and clinicopathological factors of patients with ccRCC (based on TCGA-KIRC).
Prognostic value of TGFBI
According to the result (Figure 12A), immune score (Hazard Ratio = 1.000, 95%CI of ratio: 1.000-1.000, P = 0.041), TGFBI (Hazard Ratio = 1.124, 95%CI of ratio: 1.054-1.198, P < 0.001), age (Hazard Ratio = 1.629, 95%CI of ratio: 1.209-2.196, P = 0.001), laterality (Hazard Ratio = 0.706, 95%CI of ratio: 0.523-0.953, P = 0.023), neoplasm histologic grade (Hazard Ratio = 2.661, 95%CI of ratio: 1.889-3.750 P < 0.001), and pathologic stage (Hazard Ratio = 3.841, 95%CI of ratio: 2.795-5.278, P < 0.001) were influence features of OS as suggested by univariate Cox analysis. Even being adjusted by other features, TGFBI (Hazard Ratio = 1.071, 95%CI of ratio: 0.998-1.148, P = 0.049) was still relevant to OS among patients with ccRCC suggested by multivariate Cox analysis (Figure 12B).
TGFBI expression was correlated with immune infiltration level in ccRCC
Immune infiltration level played important role in survival in tumors. Therefore, we explored the relationship between hub genes and immune infiltration level. We found that TGFBI expression was negatively relevant to tumor purity (cor = -0.250, P = 5.14E-08) and positively related to macrophages (cor = 0.232, P = 7.08E-07), Figure 13A). Summary above we found TGFBI expression was significantly correlated with tumor purity in ccRCC, which suggested that TGFBI played an important role in immune infiltration in ccRCC.
Figure 13. Correlation of TGFBI expression with immune infiltration level in ACC. (A) TGFBI expression was negatively relevant to tumor purity and positively related to macrophages. (B) Survival analyses across the six different tumor-infiltrating immune cells including B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells.
TGFBI might be potential DNA methylation biomarkers for ccRCC
We also explored the relationship between hub genes and methylation in this study by using MEXPRESS. The result demonstrated that the promoter region of TGFBI (Figure 14) showed higher methylation levels in ccRCC compared with normal tissues.
Figure 14. Visualization of the TCGA data for TGFBI in ccRCC using MEXPRESS. Samples were ordered by their expression value. This figure showed the correlation between hub gene expression and promoter methylation, with the Pearson correlation coefficients on the right (* p < 0.05, ** p < 0.01, *** p < 0.001).
Hub genes related KEGG signaling pathways
The result of GSEA (Table 3) suggested that TGFBI was significantly associated with two KEGG signaling pathways including cytokine cytokine receptor interaction (nominal P = 0.006, |ES| = 0.623, gene size = 257 and FDR = 14.763%), and chemokine signaling pathway (nominal P = 0.018, |ES| = 0.604, gene size = 185 and FDR = 14.794%).
Table 3. Genet set enrichment analysis (GSEA) of C1S and TGFBI.
As the most common subtype in kidney cancer, ccRCC was difficult to be cured. Patients with ccRCC always had poor prognosis . Because of the lack of effective therapy target and molecule drugs for ccRCC treatment, the objective of this study was to screen out novel prognostic biomarkers and small molecule drugs in ccRCC.
Tumor microenvironment (TME) was the internal environment of tumor cells producing and survival, which described the non-cancerous cells present in the tumor . These cells included but not limited to fibroblasts, immune cells and cells that comprise the blood vessels . Some proteins (such as these produced by all of the cells present in the tumor which support the growth of the cancer cells) were also included in TME . Among all the non-cancerous cells, immune cells and stromal cells were two major types, which had been proved to be correlated with prognosis of cancers . Thus, in this study, we firstly used ESTIMATE algorithm for immune score and stromal score calculation. Then we divided cases from TCGA database into two groups (high-immune score group, and low-immune score group/high-stromal score group, and low-stromal score group). DEGs between the two groups were selected. 337 immune-related DEGs including 334 up-regulated and 3 down-regulated and 218 stromal-related DEGs including 204 up-regulated and 14 down-regulated were identified. 17 DEGs overlapped in the two kinds of DEGs were identified for further analysis. 6 genes including C1R (complement C1r), C1S (complement C1s), MMP7 (matrix metallopeptidase 7), SERPINF1 (serpin family F member 1), SLC38A5 (solute carrier family 38 member 5), and TGFBI (transforming growth factor beta induced) were picked out as candidate hub genes by performing OS and DFS analyses based on TCGA-KIRC data. Then we picked out one hub gene (TGFBI) by calculating the R2 and AUC for each candidate hub gene, which were validated in mRNA level and translational-level. The prognostic value of TGFBI was validated by applying OS, CSS, and PFS analyses based on GSE29609 and E-MTAB-3267. In addition, high expression of TGFBI was correlated to higher tumor stage. The result of TGFBI genetical alteration demonstrated that hub gene altered in 75 of 533 ccRCC patients, and the main type was amplification. Further analysis suggested that TGFBI expression was associated with gender, neoplasm histologic grade, and pathologic stage. Cox regression analysis also indicated that TGFBI was influence feature of OS of patients with ccRCC.
Considering about that immune infiltration level were significantly correlated with survival in tumors, we explored the relationship between hub gene expression level and immune infiltration level in ccRCC. The result demonstrated that TGFBI expression was correlated with tumor purity in ccRCC, which indicated that TGFBI played an important role in immune infiltration in ccRCC. We also explored the correlation between TGFBI and methylation around in the promoter region. The methylation levels of the promoter region of TGFBI were higher in ccRCC compared with normal tissues.
According to the function analyses, 17 DEGs were significantly enriched in 114 BPs, 4 CCs (blood microparticle, external side of plasma membrane, side of membrane, and extracellular matrix) and 8 MFs (main type: serine hydrolase activity). The top 10 BPs were regulation of cell-cell adhesion, regulation of immune effector process, leukocyte apoptotic process, positive regulation of cell-cell adhesion, negative regulation of cell adhesion, regulation of leukocyte cell-cell adhesion, regulation of T cell activation, leukocyte cell-cell adhesion, negative regulation of leukocyte apoptotic process, and protein activation cascade. Summary above we found that these DEGs showed strong correlation with immune response and tumor immune microenvironment.
To provide clinicians some novel choice for ccRCC treatment, CMap analysis was performed and the result indicated that 3 small molecule drugs including vincamine, clenbuterol, betazole showed strong potential for ccRCC treatment.
A literature review for hub genes was carried out in order to understand hub genes deeper and better. TGFBI encoded an RGD-containing protein that binds to type I, II and IV collagens, which played an important role in the regulation of a variety of BPs . Ozawa et al. found that TGFBI expression in esophageal squamous cell carcinoma was associated with poor prognosis and hematogenous metastasis recurrence . Chen et al. demonstrated that TGFBI was an important factor in epithelial-mesenchymal transformation (EMT) and malignant progression of prostate cancer . Combining the literature review with this study, we believed that TGFBI played an important role in many events during immune response and tumor immune microenvironment.
However, this study also had certain limitations. Firstly, CD8+ cells played key role in tumor elimination. But in this study, the result (showed in Figure 13B, panel 2) demonstrated that CD8+ T cell infiltration had no or little significant effect on cumulative survival. Perhaps because of the unreasonable grouping and the TCGA data itself. Therefore, we will use novel dataset for this analysis in the short future. Secondly, CSS and PFS analyses of hub genes not showed significant results based on GSE29609 as we expected. Perhaps because of the few numbers of cases in GSE29609, which is related to too few data sets and the data itself. Thirdly, though we designed this research well and used strict thresholds for each part in this study, we did not conduct experiments to verify the results. Therefore, we will explore the related mechanisms of TGFBI in ccRCC through in vivo and in vitro experiments in further analysis.
In summary, we identified 17 DEGs associated with immune score and stromal score by using TCGA-KIRC data. TGFBI was screened out as hub gene calculating the R2 and AUC for each candidate hub gene, which was validated in mRNA level and translational-level. The prognostic value of TGFBI was validated by applying OS, CSS, and PFS analyses based on GSE29609 and E-MTAB-3267. Further analysis indicated that TGFBI might be a novel prognostic biomarker and therapy target for immunotherapy. In addition, three molecule drugs including vincamine, clenbuterol, and betazole were identified, which showed strong potential to treat ccRCC.
Materials and Methods
Data collection and data preprocessing
Microarray data of ccRCC (TCGA-KIRC data) was downloaded from The Cancer Genome Atlas (TCGA) database (https://genomecancer.ucsc.edu/). After excluding unqualified samples, totally 530 ccRCC samples from ccRCC patients were used in this study, which contained complete clinical information. Moreover, GSE29609  performed on GPL1708 was retrieved from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). This dataset including 39 ccRCCs with complete clinical information. Furthermore, E-MTAB-3267 including 53 ccRCCs with complete survival information was retrieved from ArrayExpress database (https://www.ebi.ac.uk/arrayexpress/). Both the two datasets were used for validating our findings.
Figure 15 showed the flow diagram of this study. For the TCGA-KIRC data displayed as counts number, normalized and log2 transformation were performed based on R package “DEseq.2” . For GSE29069, the raw expression data was normalization and transformation by using R package “affy” . For E-MTAB-3267, we directly downloaded the normalized expression matrix for ArrayExpress database.
Immune score and stromal score calculation and their correlation with clinical phenotype
Based on TCGA-KIRC data, we firstly calculated immune score and stromal score of each sample by using ESTIMATE algorithm, which was applied by “estimate”  in R software. Then we compared the score difference between different gender (female vs male), neoplasm histologic grade (Gx, G1, G2, G3, G4), laterality (left vs right), pathologic stage (I, II, III, IV), and person neoplasm cancer status (tumor free vs with tumor), in order to explore the relationship between immune score or stromal score and important clinical phenotype. Moreover, survival analysis was performed to explore the correlation between immune score (or stromal score) and survival by using R package “survival” .
Differentially expressed gene (DEG) screening
To pick out biomarkers associated with immune score and stromal score, we identified DEGs in this part. We firstly used R package “maxstat” for identifying the best score cutoff for grouping samples most significantly by applying a method called maximally selected rank statistics. To screen out DEGs related to immune score, we firstly divided 530 ccRCC samples into high-immune group and low-immune group based on the optimal immune score cutoff calculated by “maxstat”. Then “limma”  package in R was used for DEG screening. As for DEGs correlated to stromal score, we also divided 530 samples into two groups (high-stromal, and low-stromal) based on the optimal stromal score cutoff and screened out DEGs based on “limma” mentioned above. In this study, genes with Adjust P value < 0.05 and |log2 fold change (FC) | ≥ 1.0 were regarded as DEGs. The DEGs which exhibited the same expression trends in immune related DEGs and stromal related DEGs were picked out for further analysis.
Function enrichment analysis
To explore potential functions of DEGs, Gene Ontology (GO)  enrichment analysis and Kyoto encyclopedia of Genes and Genomes (KEGG)  pathway analysis were performed through “clusterProfiler”  in R software. Gene sets were regarded as significantly enriched gene sets when P < 0.05.
Candidate molecule drug identification
Connectivity map (CMap)  could help researchers to quickly identify molecule drugs with high correlation with diseases (https://portals.broadinstitute.org/cmap/). Thus, based on the DEGs we screened out before, we performed Connectivity map (CMap) analysis to explore molecule drugs (which had high correlation with ccRCC). In this study, drugs were thought to be highly correlated to ccRCC when number of instances (n) > 5 and P value < 0.05. Moreover, small molecule drugs with |mean| ≥ 0.40 were thought to have great potential for ccRCC treatment in this study.
Candidate hub gene identification
In this study, we aimed to screen out some stromal-immune related prognostic biomarkers in ccRCC. Therefore, we performed two kinds of survival analyses (overall survival (OS) and disease-free survival (DFS)) for each DEG to screen out candidate hub genes based on Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/) . GEPIA collected 516 ccRCC samples with complete survival information, we divided 516 ccRCCs into high-expression group (n = 258) and low-expression group (n = 258) by evaluating the median value of each DEG. DEGs with P value < 0.05 in both the two survival analyses were thought as candidate hub genes.
Hub gene identification
After screening out candidate hub genes, by using R package “basicTrendline” , we firstly calculated the R2 to see if candidate hub genes were associated with OS. P value < 0.05 was significant in this analysis. Furthermore, in order to evaluate the potential of distinguishing normal tissues and ccRCCs, we plotted receiver operating characteristic (ROC) curve and calculated AUC (area under curve) for each candidate hub gene (based on R package “pROC” ). In this study, we thought genes with P value < 0.05 in both the two analyses and AUC ≥ 0.80 were hub genes.
Hub gene expressions comparison in tumors and normal tissues and in different pathologic stages
In this part, we firstly explore the expression levels of hub genes in tumors and normal tissues by using GEPIA. Unpaired t test was performed to measure the statistical significance. Based on TCGA-KIRC data with complete stage information, tumor stage (I, II, III and IV) boxplots were also performed. We used One-way Analysis of Variance (ANOVA) test to measure the statistical significance.
Prognostic value of hub genes validation
Based on E-MTAB-3267 and GSE29609, three different survival analyses including OS, cancer specific survival (CSS), and progression free survival (PFS) were performed to validate the prognostic value of hub genes. P value < 0.05 was thought significantly.
Oncomine analysis and translational-level expression validation
After screening out the hub genes, we assessed the mRNA expression levels of hub genes in ccRCC and normal tissue based on Oncomine database (https://www.oncomine.org/) . Student t test was used to measure the statistical significance. Moreover, we validated the translation-level expression levels of hub genes by using The Human Protein Atlas database (https://www.proteinatlas.org/) .
Cancer cell line encyclopedia (CCLE) analysis
To understand the hub genes better, we explored the mRNA expression levels and copy number variation (CNV) levels of hub genes in 40 types of tumors based on CCLE database (https://portals.broadinstitute.org/ccle/).
Hub gene genetical alteration
To explore mutations and CNVs of hub genes, based on CBio Cancer Genomics Portal (http://www.cbioportal.org/) [32, 33], we first explored the genetic alterations of hub genes. Combined with relative mRNA expression levels of hub genes, the correlation between CNVs and mRNA expression levels of hub genes was also explored. Furthermore, we identified the relationship between CNVs and survival (OS, and DFS) of ccRCC patients by performing survival analysis.
Exploring relationship between hub genes and clinical features of patients with ccRCC
Based on TCGA-KIRC data, we divided 530 ccRCCs into high-expression group (n = 265) and low-expression group (n = 265) by evaluating the median value of each hub gene expression. Then χ2 test or ANOVA was performed to analyze the correlations between hub gene expressions and clinical features (age, gender, laterality, neoplasm histologic grade, pathologic stage, and person neoplasm cancer status).
Cox proportional hazards regression analysis
In this part, expression levels of hub genes and some important features (age, gender, laterality, neoplasm histologic grade, pathologic stage, and person neoplasm cancer status) were included for univariable Cox analysis of overall survival (OS) by using TCGA-KIRC data. Factors were included for multivariate Cox analysis when P value < 0.05. To do this, we could make sure if these hub genes were irrelevant to other clinical features for predicting OS of ccRCC.
Exploring correlation between hub genes and immune microenvironment
By using TIMER (https://cistrome.shinyapps.io/timer/) , we explored the correlation between hub genes expressions and the abundance of immune infiltrates. Six kinds of tumor-infiltrating immune cells (B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells) were included for this analysis . Hub genes were thought to be significant associated with infiltrating level of an immunocyte when |correlation coefficient (cor) | ≥ 0.2 and P value < 0.05. Survival analysis was also performed between high infiltrating levels and low infiltrating levels of immunocyte.
Exploring relationship between hub gene expression and methylation around in the promoter region
Considering about that MEXPRESS [35, 36] was a webtool for DNA methylation, expression, and clinical data visualization (https://mexpress.be/), we used this web tool to explore correlation between hub gene expression and methylation level. We calculated the Pearson correlation for difference (between expression value and methylation level) evaluation. For clinical factors contained two levels, P value was calculated to measure the difference. Furthermore, false discovery rate was calculated to compare the difference for clinical parameters which contained more than two levels.
Gene set enrichment analysis (GSEA)
To explore the potential functions of hub genes, GSEA  was conducted by using TCGA-KIRC data. 530 ccRCCs were firstly divided into two groups (high-expression, and low-expression) based on the median of hub genes expression levels. We set “c2.cp.kegg.v7.0.symbols.gmt” as the reference gene sets. KEGG signaling pathways with nominal P < 0.05, |ES| > 0.6, gene size ≥ 100 and FDR < 25% were considered to be significantly enriched in this study.
ANOVA: Analysis of Variance; AUC: area under curve; BP: biological process; C1R: complement C1r; C1S: complement C1s; CC: cellular component; CCLE: cancer cell line encyclopedia; ccRCC/KIRC: Clear cell renal cell carcinoma; CMap: Connectivity map; CNV: copy number variation; CSS: cancer specific survival; DAVID: Database for Annotation, Visualization, and Integrated Discovery; DEGs: Differentially expressed genes; DFS: Disease-free survival; EMT: epithelial-mesenchymal transformation; ESTIMATE: Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data; FDR: False discovery rate; GO: Gene ontology; GEO: Gene expression omnibus; GEPIA: Gene Expression Profiling Interactive Analysis; GSEA: Gene set enrichment analysis; IARC: the International Agency for Research on Cancer; KEGG: Kyoto encyclopedia of genes and genomes; MF: molecular function; OS: Overall survival; PFS: progression free survival; RCC: renal cell carcinoma; ROC: receiver operating characteristic; TCGA: The Cancer Genome Atlas; TGFBI: transforming growth factor beta induced; TME: tumor microenvironment; TOM: Topological overlap matrix; WHO: World Health Organization.
T. L., G. D., and X. Y. conceived and designed the study, G. D., X. Y., and Z. C. performed the analysis procedures, G. D., X. Y., Z. C., R. Z., K. T., J. B., H. W., and T. L. analyzed the results, T. L., and X. Y. contributed analysis tools, G. D. contributed to the writing of the manuscript. All authors reviewed the manuscript.
We would like to acknowledge the TCGA and GEO database developed by the National Institutes of Health (NIH).
Conflicts of Interest
No conflicts of interest existed in this study.
This work was supported by the National Key Research and Development Plan of China (Grant No. 2016YFC0106305).
- 1. Marabelle A, Tselikas L, de Baere T, Houot R. Intratumoral immunotherapy: using the tumor as the remedy. Ann Oncol. 2017 (suppl_12); 28:xii33–43. https://doi.org/10.1093/annonc/mdx683 [PubMed]
- 2. Yang Y. Cancer immunotherapy: harnessing the immune system to battle cancer. J Clin Invest. 2015; 125:3335–37. https://doi.org/10.1172/JCI83871 [PubMed]
- 3. Sprooten J, Ceusters J, Coosemans A, Agostinis P, De Vleeschouwer S, Zitvogel L, Kroemer G, Galluzzi L, Garg AD. Trial watch: dendritic cell vaccination for cancer immunotherapy. Oncoimmunology. 2019; 8:e1638212. https://doi.org/10.1080/2162402X.2019.1638212 [PubMed]
- 4. Xie F, Xu M, Lu J, Mao L, Wang S. The role of exosomal PD-L1 in tumor progression and immunotherapy. Mol Cancer. 2019; 18:146. https://doi.org/10.1186/s12943-019-1074-3 [PubMed]
- 5. Martini DJ, Kline MR, Liu Y, Shabto JM, Williams MA, Khan AI, Lewis C, Collins H, Akce M, Kissick HT, Carthon BC, Shaib WL, Alese OB, et al. Adiposity may predict survival in patients with advanced stage cancer treated with immunotherapy in phase 1 clinical trials. Cancer. 2020; 126:575–82. https://doi.org/10.1002/cncr.32576 [PubMed]
- 6. Ochocki JD, Khare S, Hess M, Ackerman D, Qiu B, Daisak JI, Worth AJ, Lin N, Lee P, Xie H, Li B, Wubbenhorst B, Maguire TG, et al. Arginase 2 Suppresses Renal Carcinoma Progression via Biosynthetic Cofactor Pyridoxal Phosphate Depletion and Increased Polyamine Toxicity. Cell Metab. 2018; 27:1263–1280.e6. https://doi.org/10.1016/j.cmet.2018.04.009 [PubMed]
- 7. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018; 68:394–424. https://doi.org/10.3322/caac.21492 [PubMed]
- 8. Linehan WM, Ricketts CJ. The Cancer Genome Atlas of renal cell carcinoma: findings and clinical implications. Nat Rev Urol. 2019; 16:539–52. https://doi.org/10.1038/s41585-019-0211-5 [PubMed]
- 9. Subramanian P, Haas NB. Recent advances in localized RCC: A focus on VEGF and immuno-oncology therapies. Urol Oncol. 2018; 36:23–30. https://doi.org/10.1016/j.urolonc.2017.09.008 [PubMed]
- 10. Powles T, Albiges L, Staehler M, Bensalah K, Dabestani S, Giles RH, Hofmann F, Hora M, Kuczyk MA, Lam TB, Marconi L, Merseburger AS, Fernández-Pello S, et al. Updated European Association of Urology Guidelines: Recommendations for the Treatment of First-line Metastatic Clear Cell Renal Cancer. Eur Urol. 2018; 73:311–15. https://doi.org/10.1016/j.eururo.2017.11.016 [PubMed]
- 11. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
- 12. Wu T, Dai Y. Tumor microenvironment and therapeutic response. Cancer Lett. 2017; 387:61–68. https://doi.org/10.1016/j.canlet.2016.01.043 [PubMed]
- 13. Liccardi G, Pentimalli F. Cancer, immunity and inflammation. Report from the CDD Cambridge Conferences 2018 and 2019. Cell Death Dis. 2019; 10:798. https://doi.org/10.1038/s41419-019-2032-0 [PubMed]
- 14. Hanahan D, Coussens LM. Accessories to the crime: functions of cells recruited to the tumor microenvironment. Cancer Cell. 2012; 21:309–22. https://doi.org/10.1016/j.ccr.2012.02.022 [PubMed]
- 15. Kim JE, Kim SJ, Lee BH, Park RW, Kim KS, Kim IS. Identification of motifs for cell adhesion within the repeated domains of transforming growth factor-beta-induced gene, betaig-h3. J Biol Chem. 2000; 275:30907–15. https://doi.org/10.1074/jbc.M002752200 [PubMed]
- 16. Ozawa D, Yokobori T, Sohda M, Sakai M, Hara K, Honjo H, Kato H, Miyazaki T, Kuwano H. TGFBI Expression in Cancer Stromal Cells is Associated with Poor Prognosis and Hematogenous Recurrence in Esophageal Squamous Cell Carcinoma. Ann Surg Oncol. 2016; 23:282–89. https://doi.org/10.1245/s10434-014-4259-4 [PubMed]
- 17. Chen WY, Tsai YC, Yeh HL, Suau F, Jiang KC, Shao AN, Huang J, Liu YN. Loss of SPDEF and gain of TGFBI activity after androgen deprivation therapy promote EMT and bone metastasis of prostate cancer. Sci Signal. 2017; 10:eaam6826. https://doi.org/10.1126/scisignal.aam6826 [PubMed]
- 18. Edeline J, Mottier S, Vigneau C, Jouan F, Perrin C, Zerrouki S, Fergelot P, Patard JJ, Rioux-Leclercq N. Description of 2 angiogenic phenotypes in clear cell renal cell carcinoma. Hum Pathol. 2012; 43:1982–90. https://doi.org/10.1016/j.humpath.2012.01.023 [PubMed]
- 19. Anders S. Differential gene expression analysis based on the negative binomial distribution. Journal of Marine Technology & Environment. 2009; 2.
- 20. Gautier L, Cope L, Bolstad BM, Irizarry RA. affy—analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. 2004; 20:307–15. https://doi.org/10.1093/bioinformatics/btg405 [PubMed]
- 21. Therneau TM. survival: survival Analysis. Technometrics. 2015; 46:111–12. https://doi.org/10.1017/S0954579414001333 [PubMed]
- 22. Smyth GK. limma: Linear Models for Microarray Data. Bioinformatics & Computational Biology Solutions Using R & Bioconductor. 2011:397–420. https://doi.org/10.1007/0-387-29362-0_23
- 23. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, et al, and The Gene Ontology Consortium. Gene ontology: tool for the unification of biology. Nat Genet. 2000; 25:25–29. https://doi.org/10.1038/75556 [PubMed]
- 24. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000; 28:27–30. https://doi.org/10.1093/nar/28.1.27 [PubMed]
- 25. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–87. https://doi.org/10.1089/omi.2011.0118 [PubMed]
- 26. Lamb J, Crawford ED, Peck D, Modell JW, Blat IC, Wrobel MJ, Lerner J, Brunet JP, Subramanian A, Ross KN, Reich M, Hieronymus H, Wei G, et al. The Connectivity Map: using gene-expression signatures to connect small molecules, genes, and disease. Science. 2006; 313:1929–35. https://doi.org/10.1126/science.1132939 [PubMed]
- 27. Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 2017; 45:W98–102. https://doi.org/10.1093/nar/gkx247 [PubMed]
- 28. Hothorn T, Zeileis A. Generalized maximally selected statistics. Biometrics. 2008; 64:1263–69. https://doi.org/10.1111/j.1541-0420.2008.00995.x [PubMed]
- 29. Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, Müller M. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics. 2011; 12:77. https://doi.org/10.1186/1471-2105-12-77 [PubMed]
- 30. Rhodes DR, Yu J, Shanker K, Deshpande N, Varambally R, Ghosh D, Barrette T, Pandey A, Chinnaiyan AM. ONCOMINE: a cancer microarray database and integrated data-mining platform. Neoplasia. 2004; 6:1–6. https://doi.org/10.1016/S1476-5586(04)80047-2 [PubMed]
- 31. Uhlén M, Fagerberg L, Hallström BM, Lindskog C, Oksvold P, Mardinoglu A, Sivertsson Å, Kampf C, Sjöstedt E, Asplund A, Olsson I, Edlund K, Lundberg E, et al. Proteomics. Tissue-based map of the human proteome. Science. 2015; 347:1260419. https://doi.org/10.1126/science.1260419 [PubMed]
- 32. Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, Jacobsen A, Byrne CJ, Heuer ML, Larsson E, Antipin Y, Reva B, Goldberg AP, et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2012; 2:401–04. https://doi.org/10.1158/2159-8290.CD-12-0095 [PubMed]
- 33. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, Sun Y, Jacobsen A, Sinha R, Larsson E, Cerami E, Sander C, Schultz N. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal. 2013; 6:pl1. https://doi.org/10.1126/scisignal.2004088 [PubMed]
- 34. Li B, Severson E, Pignon JC, Zhao H, Li T, Novak J, Jiang P, Shen H, Aster JC, Rodig S, Signoretti S, Liu JS, Liu XS. Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol. 2016; 17:174. https://doi.org/10.1186/s13059-016-1028-7 [PubMed]
- 35. Koch A, De Meyer T, Jeschke J, Van Criekinge W. MEXPRESS: visualizing expression, DNA methylation and clinical TCGA data. BMC Genomics. 2015; 16:636. https://doi.org/10.1186/s12864-015-1847-z [PubMed]
- 36. Koch A, Jeschke J, Van Criekinge W, van Engeland M, De Meyer T. MEXPRESS update 2019. Nucleic Acids Res. 2019; 47:W561–65. https://doi.org/10.1093/nar/gkz445 [PubMed]
- 37. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005; 102:15545–50. https://doi.org/10.1073/pnas.0506580102 [PubMed]