Construction of exosome-related genes risk model in kidney cell carcinoma predicts prognosis and immune therapy response

Data acquisition

Exosome-related data containing 15 KIRC and 118 healthy samples were downloaded from exoRBase database (http://www.exorbase.org/). The mRNA expression profile and clinical data of the TCGA-KIRC dataset, including 535 tumour tissues and 72 normal tissues, and TCGA-KIRP, including 285 tumour tissues, were downloaded from the public website (https://xena.ucsc.edu/public). The TCGA-KIRC and KIRP datasets’ mRNA expression values were transformed into log2(TPM+1) units. E-MTAB-1980 dataset containing 101 RCC samples was downloaded from the Array Express database (https://www.ebi.ac.uk/arrayexpress). E-MTAB-1980 is an external RCC data that contains clinical data and normalized mRNA expression data. Four immune therapy-related datasets, GSE135222 (n = 27), GSE78820 (n = 28), GSE79671 (n = 36), and GSE42664 (n = 45), were downloaded from the Gene Expression Omnibus (GEO) database, respectively (https://www.ncbi.nlm.nih.gov/geo/).

Differential expression analysis

Differential expression analysis was performed using the limma R package to screen for differentially expressed genes (DEGs). |log2(Fold change)| >1 and P-value < 0.05 were used as filtering thresholds for DEGs.

Survival analyses, gDS, and target development levels (TDLs)

Based on the survival R package, the mRNA expression data matched with clinical survival data was used to conduct uni-variable Cox regression analysis. The risk regression model assessed the predictive effect to identify the prognosis-related genes. The ideal cut-off value of gene expression was calculated by the surv_cutpoint R function from the survminer R package. Then, the data was divided into the high- and low-expression groups. Kaplan-Meier survival curves were depicted to show the prognostic difference between the high- and low-expression groups. gDS score of each gene in kidney cancer cells was obtained from the CCLE database (https://sites.broadinstitute.org/ccle/) and was used to demonstrate the importance of each gene for cells. Each gene’s PubTator score and target development levels were obtained from a public website (https://pharos.nih.gov/targets).

Pathway correlation analysis

The 50 hallmark pathways’ activities were figured using the Gene Set Variation Analysis (GSVA) by subjecting the gene sets acquired from the public website - MsigDB database (https://www.gsea-msigdb.org/gsea/index.jsp) to the GSVA R package. Then, we calculated the correlation between activity score and gene expression using the Pearson method.

Construction of PRERGs risk signature

Three machine learning algorithms, including GBM, Coxboost, and Boruta, were performed to select potentially essential genes. Then, the dataset was separated into training and testing cohorts in terms of a ratio of 6:4. The training cohort was used to conduct Lasso-penalised Cox regression analyses based on 10-fold cross-validation from glmnet R-package. Genes were finally identified, and risk scores were calculated using the following equation:

$$\text{Risk score}=\text{sum}{\text{(coef}}_{\text{gene}}\times {\text{Expression}}_{\text{gene}}\text{)}$$

The receiver operating characteristic (ROC) curve was applied to estimate the predictive efficacy of the risk score. Then, the training cohort was divided into high- and low-risk groups based on the median risk score value. Kaplan-Meier survival curves were depicted to determine the prognostic difference between high- and low-risk groups. The same analyses were conducted in the testing and two external validation cohorts (KIRP and E-MTAB-1980 datasets).

Clinical nomogram model

In multivariate Cox regression analysis, the risk score was adjusted by age, gender, stage, and grade. Then, based on the stepwise methods, the final variables were screened and utilized to construct a nomogram model to predict 3-/5-/7-year overall survival probability. Time-dependent ROC curve, calibration curve, and decision curve analysis were utilized to estimate the predicted efficacy and accuracy of the nomogram model.

PPI network analysis and co-expression analysis

The protein-protein interaction analysis was performed by the public website (https://string-db.org/). The clinical-actionable genes were downloaded from previous research [22]. The correlation between clinical-actionable genes and PRERGs was estimated by Pearson methods in light of their mRNA expression, and the correlation was then shown through a network diagram.

Pathway enrichment analysis

Gene ontology (GO), including cell component (CC), biological progression (BP), and molecular function (MF), were analyzed using the clusterProfiler R package. Gene set enrichment analysis (GSEA) of hallmark pathways was performed using GSEA software.

Tumour immune infiltration analysis

Based on RNA-seq data from TCGA specimens, the CIBERSORT algorithm was used to quantify the proportion and distribution of tumour-infiltrating immune cells (TIICs) by presenting expression matrix and reference gene sets. We used the xCell R package to assess each sample’s immune, stromal, and tumour microenvironment scores. The gene sets of immune cells were downloaded from previous research [23]. Then, the single sample gene set enrichment analysis (ssGSEA) method was used to calculate the activities of immune cells.

Statistical analysis

The unpaired t-test was used to compare the differences between the two groups. The association between variables and overall survival (OS) was assessed using univariate and multivariate Cox regression analyses. The survival probability differences of the two groups in Kaplan-Meier were assessed using the log-rank test. Pearson correlation coefficient (PCC) was used to quantify the correlation between the two groups. |PCC| >0.3 and p-value < 0.05 were regarded as significant.

Data availability

All datasets generated for our research are introduced in the article.

Exosome-related genes identification

The Flow diagram is shown in Figure 1. We first conducted the differential expression analysis between KIRC and normal samples from TCGA and exoRBase. 3617 differential expression genes (DEGs) in TCGA-KIRC and 723 DEGs in exoRBase were eventually identified (Figure 2A). We further determined 18 common up-regulated genes and 44 common down-regulated genes in both datasets based on the Venn diagrams; these genes were named exosome-related genes (ERGs) (Figure 2B). The heatmaps and t-SNE analysis revealed that the mRNA expression of these ERGs could significantly distinguish normal and tumour samples (Figure 2C–2F). Therefore, we included these genes for further exploration.

Figure 1. Flow chart of this research.

Figure 2. Identification of exosome-related genes in KIRC. (A) The volcano figure shows the differentially expressed genes between KIRC and normal control in exoRBase and TCGA-KIRC datasets. (B) The Venn diagram shows the common up-regulated genes and down-regulated genes in exoRBase and TCGA-KIRC datasets. (C, D) Heatmaps show the expression levels of ERGs in tumours and normal tissues. (E, F) T-SNE analyses were performed to show the distribution of tumour tissues and normal tissues based on ERG expression.

Clinical importance evaluation of ERGs

Next, we performed a univariate Cox regression analysis to verify the relationship between prognosis and ERGs. 26 prognosis-related ERGs (PRERGs) were finally determined. Among these, 9 PRERGs were recognized as risk factors (TNFSF14, CD7, CDC6, BATF, CD82, XK, MEST, XCL2, and CDT1) while others as protective factors (LRBA, CTDSPL, OPHN1, SNRPN, APP, MFAP3L, ASAP2, FAXDC2, CA2, DNAJC6, VEPH1, BMP6, SELP, SUCNR1, DLL4, SPX, and ABO) (Figure 3A). The patients with high expression of risk genes showed a worse survival probability than those with low expression. In the same way, patients with low expression of protective genes showed a poorer prognosis than those with high expression (Supplementary Figure 1). Then, we used the gene dependency score (gDS) to assess the “importance” of 26 PRERGs. Several genes showed a low gDS in most kidney cancer cell lines, such as CD82, CD7, MEST, SELP, BMP6, CA2, DNAJC6, and VEPH1. CDC6 and CDT1, especially, represented the lowest gDS in kidney cancer cells, indicating that deletion of these two genes severely affects tumour cell viability (Figure 3B). In addition, we reviewed the target development levels (TDLs) and PubTator scores of these genes. We found APP, CA2, SELP, and ABO possess a higher PubTator score which is over 2000. PubTator scores of DLL4, CD7, LRBA, CD82, BMP6, and SNRPN are 200-2000, while others are lower than 200. Next, we analyzed the drug development status of PRERGs and their target development levels (TDLs). The TDLs were evaluated by Illuminating the Druggable Genome Knowledge Management Center, NIH. Three genes, APP, CA2, and SELP, have been developed targeted drugs, and SUCNR1 has been developed Tchem. Most genes do not have known approved drugs, and their small molecule activities do not satisfy the cutoffs of Tclin and Tchem. However, they are annotated with experimental evidence or gene ontology and are defined as Tbio. MFAP3L and FAXDC2 have very few studies and are defined as Tdark (Figure 3C). We further explored the relationship between these genes and cancer-related pathways by calculating the Pearson correlation coefficient (PCC) between the mRNA expression value of PRERGs and the activity score of cancer-related pathways. These PRERGs showed distinct associations with various pathways. CDC6 and CDT1, for instance, positively correlated with various cycle- and proliferation-related pathways, such as E2F targets and G2M checkpoint, indicating that they regulate cancer development by activating these pathways (Figure 3D). CD7, MEST, and XCL2 are highly correlated with IL6/JAK/STAT3, inflammatory response, and interferon-gamma/alpha response, revealing their peculiar regulation mechanism. Notably, CDC6 and CDT1 showed lower gDS and PubTator scores, are referred to as Tbio in Target development levels, and positively correlated with cancer-related pathways, which is worthy of further experimental exploration in RCC.

Figure 3. Importance evaluation of PRERGs in KIRC. (A) The forest diagram shows prognosis-related ERGs (PRERGs). (B) gDS score map shows the importance of a single gene to kidney cancer cells. (C) PubTator score and target development level (TDL) of PRERGs. The size of the circle represents different PubTator scores; the colour of the circle represents different TDLs. (D) The heatmap reveals a correlation between 50 cancer-related hallmark pathways activities and expression of PRERGs in KIRC.

PRERGs risk signature construction based on machine learning

To explore the relationship between PRERGs and the prognosis of RCC patients, we conducted further screening by utilizing three machine learning methods, including GBM, Boruta, and Cox-Boost. 9 PRERGs were finally determined (Figure 4A). Then, we divided the TCGA-KIRC dataset into training and testing cohorts and conducted a Lasso-cox regression analysis based on the mRNA expression of these nine genes in the training cohort (Supplementary Figure 2A, 2B). A 5-gene risk signature (PRERGs risk signature) was finally determined, including TNFSF14, CA2, CDC6, LRBA, and MFAP43L. Based on their lasso coefficients and mRNA expression values, the risk score per patient was calculated by the formula as follows:

$$\begin{array}{l}\text{Risk score}=\text{TNFSF14}\times 0.1508474-\text{CA2}\times \\ 0.01339085+\text{CDC6}\times 0.18011894-\text{LRBA}\times \\ 0.31228396-\text{MFAP3L}\times 0.17567860\end{array}$$

High-risk patients showed a poorer prognosis than low-risk patients (Figure 4B). The risk score showed a high accuracy in predicting the prognosis of KIRC patients (5 years ROC = 0.7229, 7 years ROC = 0.7274). In the high-risk group, the number of patients proved to be dead was far more than the number of patients in the low-risk group. In addition, we observed that CDC6 and TNFSF14 are highly expressed in the high-risk group, while CA2, LRBA, and MFAP3L are highly expressed in the low-risk group. Further validation was performed in the testing cohort, E-MTAB-1980, and KIRP datasets, and the results bear resemblance to the training cohort (For testing cohort: 5 years ROC = 0.7341, 7 years ROC = 0.7295; for E-MTAB-1980: 5 years ROC = 0.7529, 7 years ROC = 0.6735; for TCGA-KIRP: 5 years ROC = 0.6758, 7 years ROC = 0.6645) (Figure 4C–4E). In addition, we found that these genes interacted with different proteins and positively or negatively correlated with different clinically actionable genes, indicating that they might regulate cancer by different mechanisms and interact with different substrates (Supplementary Figure 3A, 3B).

Figure 4. Construction of PRERGs signature. (A) The flow chart shows the machine-learning process. (B–E) Kaplan-Meier survival curve shows the survival differences in high- and low-risk groups. ROC curve shows the predicted accuracy of the risk score on the prognosis of the patients. The risk dot plot shows the number of dead and alive patients with gradually increasing risk scores. Heatmap shows the expression of 5 PRERGs that constructed the risk signature in the high- or low-risk group.

PRERGs risk signature correlated with clinical information

Next, we excavated the relationship between clinical status and risk score in three RCC datasets. We found that PRERG’s risk score was associated with pathological stage and grade of TCGA-KIRC. Specifically, patients with advanced stages and grades are mainly distributed in the high-risk group (Figure 5A). The risk score was then adjusted by other clinical variables such as age, gender, grade, and stage during multivariate Cox regression analysis. The results showed that the PRERG risk score was an independent prognostic factor (Figure 5B). What is more, the risk scores were practical in different clinical subgroups. Patients with high risk always showed worse outcomes than patients in low-risk groups (Figure 5C). The same analyses were performed in another two datasets and revealed similar results (Supplementary Figures 4 and 5). Notably, the risk score showed no correlation with stage and grade in the E-MTAB-1980 dataset, which might be due to small sums of samples. However, we still observed that more patients with advanced stages and grades were distributed in the high-risk group than in the low-risk group. Despite different datasets, PRERGs still acted as an independent prognostic factor, and the high-risk group invariably showed a lower survival rate than the low-risk group, indicating that PRERGs were steady and robust in predicting the prognosis of RCC patients.

Figure 5. Clinical relevance and robustness of PRERGs risk score in TCGA-KIRC dataset. (A) The heatmap shows the distribution of clinical variables in the high- or low-risk groups. (B) The results of multiple Cox-regression analysis in the TCGA-KIRC dataset. (C) Kaplan-Meier shows the prognostic differences between high- and low-risk groups in classified clinical variables. ^*P-value < 0.05; ^**P-value < 0.01; ^***P-value < 0.001; ^****P-value < 0.0001.

Clinical model construction

Using multivariate Cox regression analysis and stepwise regression method in the TCGA-KIRC dataset, a nomogram model was subsequently established to estimate 3-/5-/7-year survival probability (Figure 6A). Time-dependent ROC and calibration curves revealed that the model possesses a discriminative accuracy (3 years AUC = 0.8079; 5 years AUC = 0.7733; 7 years AUC = 0.7510). Decision curve analysis (DCA) illustrates that this nomogram model could receive more clinical benefit than single variables (Supplementary Figure 6A). Based on the nomogram model, a risk score was finally determined, and the patients in different risk groups showed significant survival differences (Figure 6A). The same analyses were performed in another two external validation cohorts, and the results proved to be the same trends (Figure 6B, 6C; Supplementary Figure 6B, 6C). We then compared our PRERG signature with nine other risk signatures [24–32]. PRERGs risk signature showed the top 3 of AUC for 3-/5-/7-year prognosis of the KIRC patients (Supplementary Figure 7A–7C), indicating that our risk signature has a higher predicted accuracy.

Figure 6. Nomogram construction in three RCC datasets. (A–C) A nomogram was constructed based on PRERGs and clinical variables to predict the 3-/5-/7-year survival probability of the patients in the TCGA-KIRC, TCGA-KIRP, and E-MTAB-1980 datasets. ROC and calibration analysis were performed to estimate the predicted accuracy. The Kaplan-Meier survival curve was depicted to show the prognosis difference of high- or low-risk groups which was calculated by nomogram.

Potential molecular mechanism

To demonstrate the potential molecular mechanism between the high PRERGs risk group and the low PRERGs risk group, we performed enrichment analysis by submitting the differentially expressed genes screened from the high-/low-risk group to the clusterProfiler R package. We found cell components, such as blood microparticles, immunoglobulin complex, immunoglobulin complex circulating, plasma lipoprotein particle, and protein-lipid complex, highly enriched in the high-risk group. In the case of biological progression, antibacterial humoral response, antimicrobial humoral response, positive regulation of T cell migration, regulation of CD4-positive alpha-beta T cell differentiation, and regulation of T cell migration are highly enriched in the high-risk group. Meanwhile, molecular function, antigen binding, chemokine activity, chemokine receptor binding, fatty acid binding, and immunoglobulin receptor binding are highly enriched in high-risk groups (Supplementary Figure 8A–8C). By presenting the mRNA expression profile ordered by risk score per sample to GSEA software, we found no hallmark pathways significantly enriched in a high-risk group. In contrast, several pathways showed a correlation with the low-risk group. The top 5 hallmark pathways that were enriched in the low-risk group are adipogenesis, heme metabolism, fatty acid metabolism, bile acid metabolism, and peroxisome pathways, indicating that low enrichment of these genes in the high-risk group might be the reason causing poor outcomes in the patients (Supplementary Figure 8D).

Immune infiltration and therapeutic response

We further explored the infiltration levels of different immune cells in these two groups. As expected, we found that CD4+ T cells and CD8+ T cells show higher fraction and infiltration levels in high-risk groups. However, we observed that Tregs are highly enriched in the high-risk group (Figure 7A, 7B). Furthermore, the immune and microenvironment scores in the high-risk group show a high enrichment, while the stroma score shows no significant difference in these two risk groups (Figure 7C). The immune checkpoints PDCD1 (also known as PD1) and CTLA4 show a higher expression in the high-risk group, while CD274 (also known as PDL1) shows a higher expression in the low-risk group (Figure 7C). To explore whether the PRERGs risk signature could predict immune therapy response, we further carried out Tumor Immune Dysfunction and Exclusion (TIDE) (http://tide.dfci.harvard.edu/) analysis by presenting a normalized gene expression profile of TCGA-KIRC [33]. The results reveal that most responders (R) harboured in low-risk groups, though the distribution did not show statistical significance (Supplementary Figure 9). We then validated the speculation using four other external immune therapy cohorts (including GSE135222, GSE78820, GSE79671, and GSE42664). We extracted the mRNA expression of PRERGs from normalized datasets and then calculated the risk score. Then, we matched the risk score with the clinical data of each dataset; among these, GSE135222 and GSE78820 contain survival data that enable us to perform survival analysis. The high- and low-risk groups were determined by ideal cutoff value, and we also observed that high-risk group patients showed low survival probability in these two cohorts, even though the difference was not significant statistically, that might be due to the small sums of samples (Figure 7D, 7E). We further analyzed another three immune therapeutic cohorts (GSE78820, GSE79671, and GSE42664); conformably, no responders in these datasets showed a higher risk score, while those proved to be responders represented a lower risk score (Figure 7F, 7H and 7J). We calculated the proportion of responders and non-responders harboured in the high- or low-risk group. The results revealed that a higher proportion of no-responders were distributed in the high-risk group. In contrast, only a slight collection of non-responders was located in the low-risk group (Figure 7F–7K). These results suggested that the PRERG risk signature might act as a predictor for predicting the reaction of patients who received immune therapy.

Figure 7. Relevance between tumour immune infiltration and PRERGs risk signature. (A) Fraction of immune-related cells in the high- or low-risk group. (B) Activity score of immune-related cells in the high- or low-risk group. (C) Immune score, stromal score, microenvironment score, and three immune checkpoint expressions in the high- or low-risk group. (D, E) Kaplan-Meier survival curves were portrayed to show the survival difference in the high- or low-risk group in the GSE135222 and GSE78820 datasets. (F, G) Boxplot shows the risk scores of CR/PR and SD/PD groups in the GSE78820 dataset. The barplot shows the distribution of high- or low-risk groups in CR/PR and SD/PD groups. (H, I) Boxplot shows the risk score of the response and no-response groups in the GSE79671 dataset. The barplot shows the distribution of high- or low-risk groups in the response and no-response groups. (J, K) Boxplot shows the risk score of CR/PR and SD/PD groups in the GSE42664 dataset. The barplot shows the distribution of high- or low-risk groups in CR/PR and SD/PD groups. Abbreviations: CR: complete response; PR: partial response; SD: stable disease; PD: progressive disease. ^*P-value < 0.05; ^**P-value < 0.01; ^***P-value < 0.001; ^****P-value < 0.0001.