Using exosome-related RNA levels to predict immunotherapy response in clear cell kidney cancer

RNA‐seq and miRNA‐seq transcriptome data from 539 ccRCC tissues and 72 adjacent normal tissues and the corresponding clinical information were downloaded from the Cancer Genome Atlas (TCGA) (https://cancergenome.nih.gov/↗) database. GSO datasets (GSE22541↗ and GSE167573↗) were downloaded from the GEO (https://www.ncbi.nlm.nih.gov/↗) database. GSE22541↗ included 44 ccRCC lung metastatic tissues and 24 primary ccRCC tissues. GSE167573↗ contained 63 primary ccRCC tissues and 14 adjacent normal tissues. A total of 120 exosome‐related genes were obtained from a comprehensive database of exosomes (https://www.exocarta.org/↗).37

First, we performed a differential expression analysis to screen for exosome‐related genes that are characteristic for the ccRCC using the R package “limma” (|log FC| >0.1 and FDR <0.05). Univariate cox regression analysis was performed to identify the prognostic exosome‐related genes using R package “survival.” The intersection genes of differentially expressed genes (DEGs) and prognostic exosome‐related genes were then incorporated into the protein–protein interaction (PPI) analysis. Ten genes, with the highest number of node connections in the protein interaction network, were identified as hub exosome‐related genes in ccRCC. Moreover, hub exosome‐related gene‐to‐lncRNA correlations were generated using Spearman correlation analysis. The lncRNAs with R > 0.5 and FDR < 0.001 were identified as ERLRs.

The entire TCGA ccRCC cohort was randomly divided into the training cohort (n = 265) and the validation cohort (n = 265) with the ratio of 1:1. The training cohort was adopted for the construction of the prognostic ERLR‐based score, while the validation cohort was used to verify accuracy. In the training cohort, a univariate Cox regression analysis was performed to identify prognostic ERLRs. Next, based on prognostic ERLRs, LASSO regression analysis was used to choose the most significant genes by R package “glmnet” in the training cohort. Finally, these genes were incorporated into the multivariate Cox analysis to determine vital ERLRs and construct a score. The riskScore of each ccRCC patient was calculated using the followed formula:RiskScore=e^∑CoefiExpi,

“e” is the natural logarithm, Coef(i) is the coefficient, and Exp(i) represents the expression of ERLRs. The median of the riskScore in the training cohort was used as the cutoff value to classify all ccRCC patients into the high‐ and low‐risk groups.

A Kaplan–Meier survival curve was drawn to display the differences in prognosis between the high‐ and low‐risk groups with the log‐rank test. Next, we implemented time‐dependent receiver operating characteristic (ROC) curves at 1‐, 3‐, and 5‐year survival, using the R package “SurvivalROC” to verify the performance of the score. Moreover, stratification analyses were generated among gender (female, male), age (≤65 or >65), grade (I‐II, III‐IV), and stage (I‐II, III‐IV) to confirm that the score was consistent across these subgroups. Also, the riskScore of the different subgroups was compared using the Wilcoxon signed‐rank test to perform an internal verification.

Age, gender, grade, stage, and riskScore were included in the univariate and multivariate Cox regression analyses to screen for independent prognostic factors. A nomogram was developed using the R package “regplot” to predict the 1‐, 2‐, and 3‐year overall survival (OS) of ccRCC patients. Calibration curves were plotted using the R package “RMS” and “Survival” to evaluate the predictive power of this nomogram. Additionally, the time‐dependent ROC curves for all factors were drawn, and the area under the curve (AUC) for all factors was compared with the AUC of the score and nomogram.

Gene set enrichment analysis (ssGSEA) was used to calculate the infiltrating scores of 23 types of immune cells, immune‐related pathways, and stromal‐related pathways using the R package “GSVA.” CIBERSOFT, as another method to evaluate the infiltration scores of immune cells, was performed to calculate the infiltrating scores of 22 types of immune cells via R package “CIBERSOFT.” In addition, the “Maftools” R package was used to evaluate mutations in different risk score groups, and two TME‐related scores (Immune score and Stromal score) were calculated by the R package “ESTIMATE.”

The tumor immune dysfunction and exclusion (TIDE) algorithm (https://tide.dfci.harvard.edu↗) was applied to obtain seven immunotherapy related scores, including ips_ctla4_pos_pd1_pos score, CD8 score, CD274 score, Dysfunction score, TIDE score, merck18 score, and IFNG score.

Importance of the identified ERLRs in prognosis was analyzed through mechanical learning with R package “randomForest,” and EMX2OS was ank1 in 4ERLRs. Downstream miRNA of EMX2OS was screened based on the mircode database (https://www.mircode.org/↗), and the intersection of differentially expressed miRNA in the TCGA database was obtained to identify the target miRNA. Based on the miRwalk database (https://mirwalk.umm.uni‐heidelberg.de/↗), criterion was set with a score = 1, and downstream mRNAs of relevant miRNAs were screened in TargetScan (https://www.targetscan.org/vert_80/↗) and miRDB (https://mirdb.org/↗). This was intersected with differentially expressed mRNAs in the TCGA database to identify downstream mRNAs. Finally, correlations with the ceRNA interaction network were established. The ceRNA network was visualized by introducing edge and nodal gene information into the interaction network with Cytoscape v3.8.0. The importance of target‐miRNAs and target‐mRNAs was analyzed through mechanical learning.

To better understand the underlying function of aberrantly expressed genes, gene ontologies (GO) were collected. We used the Database for Annotation Visualization and Integrated Discovery 6.8 (DAVID) (https://david.ncifcrf.gov/↗) to perform the functional analysis. Then, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were constructed using KOBAS 3.0 (https://kobas.cbi.pku.edu.cn/anno_iden.php↗).

Based on the expression level of key RNAs, all patients were divided into low expression and high expression groups according to the median expression level. Survival curves were drawn to compare the survival differences between the two groups.

In this study, IC50 values of drugs were obtained in Genomics of Drug Sensitivity in Cancer (GDSC) using the R package “oncoPredict.” Correlations between drug IC50 values and risk scores were analyzed by Spearman's analysis to screen for potential drugs. We then compared the differences in IC50 values between the high‐ and the low‐expression groups for drugs with absolute values of correlation greater than 0.4. The results were then visualized using the R language ggplot2. Additionally, based on the expression levels of EMX2OS and TLN2, a Kaplan–Meier survival curve was drawn to display the differences of OS between the high‐ and low‐expression groups after anti‐PD‐1, anti‐PD‐L1, or anti‐CTLA‐4 treatment‐ using the log‐rank test. Next, we implemented time‐dependent ROC curves at anti‐PD‐1, anti‐PD‐L1, or anti‐CTLA‐4 treatment using R package “SurvivalROC” to verify score performance.

CcRCC tissues and normal control tissues were obtained from ccRCC patients undergoing radical nephrectomy in the First Affiliated Hospital of Wenzhou Medical University. This study was approved by the Human Research Ethics Committee of the First Affiliated Hospital of Wenzhou Medical University (KY2021‐160).

All statistical analyses were performed in R software (version 4.1.0), GraphPad Prism 8 (GraphPad Software, United States), and PASS 27.0 software. Each chapter introduces particular software packages utilized by R studio. The data were expressed as the mean ± standard deviation (SD) of at least three independent experiments, and the differences between the two groups were compared using the Student t‐test. A one‐way ANOVA with Bonferroni's multiple comparisons test was used for multiple group comparisons. Spearman's correlation coefficient test was used to evaluate the rank correlations among different variables. Kaplan–Meier survival curve was used to analyze survival differences, and the Log‐rank test was used to assess the significance of the difference in survival time between two groups. Statistical p‐values were subjected to two tailed tests, and p < 0.05 was considered statistically significant.

We identified DEGs in ccRCC by comparing gene expression levels in adjacent normal renal tissues using the TCGA database. A total of 99 DEGs were identified between ccRCC and normal samples (Figure S1A). Of these, 46 differentially expressed genes were associated with ccRCC exosomes and had prognostic values (Figure S1B,C). Then, we performed a protein–protein interaction network (PPI) analysis on these 46 genes and selected the TOP 10 genes for the following analysis according to differences in expression levels (Figure 1A). The differential expression of the TOP 10 exosome‐related DEGs in ccRCC was shown by heatmaps in Figure 1B. In total, 669 ERLRs were identified to be associated with top10 prognostic exosome‐related DEGs (Figure 1C). To construct the ERLRs prognostic score for forecasting the overall survival (OS) of ccRCC patients, we performed a LASSO regression and multivariate cox regression analysis with the 669 ERLRs, the four ERLRs (EMX2OS, AC026401.3↗, AC018690.1↗, and AL161935.1↗) were screened for the greatest prognostic value in ccRCC patients (Figure 1D,E). Univariate Cox regression analysis confirmed the prognostic value of the four ERLRs (p < 0.05) (Figure 1F). The whole ccRCC patients were randomly divided into two cohorts (training: validation = 1:1). Then, based on the correlation coefficient of the expression of lncRNAs in the multivariate Cox regression analysis, we construct a risk factor prediction linear model. According to the median ERLR‐based score, the patients in each cohort were divided into low‐risk group and high‐risk group (Figure 1H,K). The results of the OS Kaplan–Meier curve in the two cohorts were consistent (Figure 1G,J), and the survival rate of ccRCC patients in the low‐risk group was significantly better than that in the high‐risk group. A time‐dependent receiver operating characteristic was used to investigate the prognostic value of the characteristic (Figure 1I,L). The 1‐, 3‐, and 5‐year AUCs for the training group were 0.777, 0.731, and 0.778, respectively, while the 1‐, 3‐, and 5‐year AUCs for the validation group were 0.729, 0.651, and 0.694, respectively. These results showed that ERLR‐based score could be used as a predictor of ccRCC patients' outcomes.

To further explore the potential mechanism, we first used the R package “limma” to mine DEGs in low‐ and high‐groups, where a total of 99 DEGs were identified (Figure). Subsequently, gene name and log2FC of the DEGs in the groups were used for enrichment analysis of gene ontology (GO) and KEGG pathways (Figure). KEGG analysis also demonstrated that these genes were mainly involved in PPAR signaling pathways and metabolic pathways. GO analysis revealed that most DEGs were enriched in immune‐related pathways, such as B‐cell immunity, humoral immune response mediated by circulating immunoglobulin, and complement activation. The above results suggest that a better prognosis of ccRCC patients in the low‐risk group may be caused by the anticancer effects of immune‐related pathways. S2A S2B,C

For further analysis, we identified several correlations between 4‐ERLR‐based score and clinical features, including gender, age, tumor stage, and grade. Results indicated that ERLR‐based score was higher in male ccRCC patients, those with a tumor histology grade of 3–4 and stage III‐IV (Figure 2B–D), regardless of age (Figure 2A). Based on four clinical variables, we divided all patients into several subgroups to display survival using Kaplan–Meier survival analysis and found that high‐risk groups had worse survival across different subpopulations. We confirmed that the 4‐ERLR‐based score had accurate predictive ability for ccRCC patients with different ages (≤65 or >65), different genders (female or male), different tumor histology grades (1–2 or 3–4), and different stages (I‐II or III‐IV) (Figure 2E–L). These results implied that the ERLR‐based score may be intimately associated with tumor progression and could predict the outcome of multiple subgroups of ccRCC patients.

Then, univariate Cox regression analysis and multivariate Cox regression analysis showed that age, tumor stage grade, and a 4‐ERLRs score were significantly correlated with the prognosis of ccRCC patients (Figure 2M,N). Therefore, 4‐ERLR‐based score is a promising independent prognostic scoring system in ccRCC. Meanwhile, to develop a clinical utility for predicting the survival probability in ccRCC patients, we constructed a nomogram integrating 4‐ERLRs score and clinicopathological features, including age, gender, grade, and stage (Figure S3A,C). The calibration plots displayed a good performance in 1‐, 2‐, and 3‐year overall survival, indicating that the nomogram had good accuracy as an ideal model (Figure S3B,D). The nomogram predicted the probability of 1, 3, and 5 years in ccRCC patients, combined with the ERLR‐based score, providing a quantitative method for predicting the OS in ccRCC patients, and acting as a valuable clue for clinicians to make medical decisions and follow‐up plans.

To investigate the roles of the 4‐ERLR‐based score in the immune microenvironment of ccRCC patients, ESTIMATE algorithm was used to further investigate the associations between the 4‐ERLR‐based score and immune infiltration cells. Proportions of various immune cells were shown in Figure 3A, which intuitively reflected immune cells pattern differences in the low‐ and high‐risk groups. T cells and macrophages accounted for the largest components. Correlation between immune infiltration status and riskScore was assessed by seven different types of software (Figure 3B).

Then, we assessed immune infiltration status in the low‐ and high‐risk groups. We estimated the proportions of 22 immune cell types in ccRCC by running the CIBERSORT script and found that CD8 + T cells, T follicular helper cells (Tfh), and regulatory T cells (Tregs) were significantly enriched in the high‐risk group. Monocytes, macrophages M0, macrophages M1, resting dendritic cells, activated dendritic cells, and resting mast cells were mainly enriched in the low‐risk group (Figure S4A). Meanwhile, similar outcomes were found by running the ESTIMATE script (Figure S4B). At the same time, the immune score, tumor microenvironment interstitial score, and estimate score of the high‐risk group were statistically higher, indicating that the high‐risk group may be more sensitive to immunotherapy than low‐risk group (Figure 3C,D) (Figure S4E). Moreover, infiltration levels of all 10 significant immune‐associated gene sets were higher in high‐risk ccRCC patients, including APC‐co‐stimulation, CCR, checkpoint, cytolytic activity, inflammation‐promoting, MHC class I, parainflammation, T‐cell co‐inhibition, T‐cell co‐stimulation, Type I IFN Response, except for Type II IFN Response (Figure 3F). These findings suggest that the immune function of high‐risk ccRCC individuals was more active than that of low‐risk ccRCC individuals. Correlations between stromal‐related pathways and the 4‐ERLR‐based score were also analyzed. Among them, the high‐risk group was correlated with Pan‐F‐TBRS and EMT2, while the low‐risk group was associated with FGFR‐related genes, angiogenesis pathways, EMT1 and WNT targets (Figure 3E). Meanwhile, immune activation genes were upregulated in the high‐risk group (Figure S4C). Differences in TGFβ and EMT pathway‐related gene expression levels of the two groups are shown in Figure S4D, suggesting that the expression levels of PDGFRA, TWIST1, and VIM were increased in the high‐risk group, indicating a poorer prognosis in ccRCC patients.

In conclusion, the high‐risk group had higher levels of immunosuppressive cell, suggesting that the immune function of these cells was possibly impaired. At the same time, the high‐risk group had a potential state of increasing vascularization and a worse prognosis. This result provides potential immunotherapy targets and therapeutic directions.

Since the tumor mutation burden (TMB) has been recognized as a predictive biomarker for immunotherapy response in cancer, we further explored the correlation between the TMB and 4ERLRs‐based score. TMB was higher in the high‐risk group, indicating that tumors with more mutations were enriched in the high‐risk group, leading to a worse prognosis (Figure 4A,B). Combined with TMB and riskScore, high‐TMB + low‐RiskScore is slightly better than low‐TMB + high‐RiskScore (Figure 4C), so it can be concluded that riskScore plays a more significant role than TMB, and is more valuable when predicting survival and is related to the amount of tumor mutations.

We then analyzed distributional differences of somatic mutations between the low‐ and high‐risk groups. We found that VHL and PBRM1 had similar mutation frequencies and the two highest mutation frequencies in patients with ccRCC were in SETD2 and BAP1, more than three times higher in the high‐risk group than in the low‐risk group (Figure 4D,E). These would be valuable to reveal underlying mechanisms by which mutations in BAP1 and SETD2 exert diverse effects on cancer progression in ccRCC. The 4‐ERLR‐based score might have implications for predicting immunotherapeutic outcomes in ccRCC.

Immune checkpoint inhibitor (ICI), as a kind of tumor immunotherapy, is a promising strategy for ccRCC. Immune checkpoint molecules include PD‐1, PD‐L1, cytotoxic T lymphocyte‐associated protein 4 (CTLA‐4), T‐cell immunoglobulin mucin‐3 (TIM‐3), and others. The expression levels of PD‐1 and CTLA4 were upregulated in the high‐risk group, indicating that patients in the high‐risk group were promising candidates for anti‐PD1 and anti‐CTLA4 immunotherapy (Figure 5A–C). At the same time, ips_ctla4_pos_pd1_pos, CD8 score, and CD274 score were higher in the high‐risk group than in the low‐risk group, suggesting that patients in the high‐risk group were more likely to benefit from immunotherapy (Figure 5D–F). Moreover, dysfunction score, TIDE, and Merck 18 are higher in the high‐risk group, indicating that the potential immune escape is stronger and more cytotoxic T lymphocyte infiltration in the high‐risk group (Figure 5D–J). Notably, IFNG scores also are higher in the high‐risk group (Figure 5K). In conclusion, the poor prognosis of high‐risk ccRCC patients may be due to tumor immune evasion which could contribute to tumor invasion and metastasis.

The competitive endogenous RNA (ceRNA) is a well‐known regulatory mechanism of lncRNA. Among the 4ERLRs, we found EXM2OS was the first rank of importance‐through mechanical learning (Figure 6A,B). Kaplan–Meier analysis showed that ccRCC patients with high EMX2OS expression had better survival (Figure 6C). In addition, we found that EMX2OS was significantly downregulated in the ccRCC validation set GSE167573↗ and TCGA_ccRCC (Figure 6D), which was also validated in clinical samples (Figure 6E). Then, we constructed a ceRNA network composed of EMX2OS, differentially expressed miRNAs (DEmiRNAs), and differentially expressed mRNAs (DEmiRNAs) based on multiple databases. The network diagram was visualized by Cytoscape, in which EMX2OS, 6 DEmiRNAs, and 26 DEmRNAs were shown as nodes and 32 interactions were shown as edges (Figure 6F). Meanwhile, analysis of the KEGG signaling pathway and GO analysis revealed that the DEGs were mainly enriched in peroxisome proliferator‐activated receptor (PPAR) signaling pathway and complement activation classical pathway, with a close relationship to regulated amino acid metabolism, regulation of transmembrane transport, and regulation of receptor‐mediated immunity. Furthermore, five significantly correlated GO and KEGG items were determined and visualized according to genes (Figure 6G,H). Then, we found that hsa‐miR‐183‐5p, hsa‐miR‐210‐3p, and hsa‐miR‐31‐5p are more important in DEmiRNAs using mechanical learning (Figure 6I,J). Kaplan–Meier analysis showed that ccRCC patients with low expression of hsa‐miR‐183‐5p or hsa‐miR‐31‐5p had better survival, while the high expression group of hsa‐miR‐210‐3p had better survival (Figure 6K–M). The expression of hsa‐miR‐31‐5p was significantly increased in clinical samples (Figure 6N). Next, we found that TLN2 played a key role as the downstream mRNA of hsa‐miR‐31‐5p by machine learning (Figure 6O,P). Kaplan–Meier analysis also showed that ccRCC patients with high TLN2 or HEYL expression had better survival. We found that TLN2 was significantly downregulated in the ccRCC validation set GSE167573↗ and TCGA_ccRCC (Figure S6A–C), which was also validated in clinical samples (Figure 6Q,R). Therefore, our study found that the downregulation of EMX2OS and TLN2 may be related to the poor prognosis of ccRCC. The EMX2OS/hsa‐miR‐31‐5P/TLN2 axis might yield functional impact on tumor microenvironment in ccRCC and may represent a novel therapeutic target for ccRCC.

For EMX2OS, it is found that the anti‐PD‐1 AUC is 0.733 and 0.688 in the Homet cohort 2019 and Amato cohort 2020, respectively. The Anti‐CTLA‐4 AUC is 0.668 in the Nathanson cohort 2017 (Figure 7A–C). For TLN2, it is found that the anti‐PD‐1 AUC is 0.700 and 0.679 in the Homet cohort 2019 and Ascierto cohort 2016, respectively. The anti‐CTLA‐4 AUC is 0.705 in the VanAllen cohort 2015. In addition, the anti‐MAGE‐A3 AUC is 0.672 in the Dizier cohort 2013 (Figure 7G,I). Therefore, we found that the response to anti‐PD‐1 and anti‐CTLA4 treatments was closely correlated with the expression levels of EMX2OS and TLN2 in ccRCC patients (Figure 7D,J). Meanwhile, lower expression of EMX2OS or TLN2 was significantly associated with longer disease‐free survival, shorter disease‐specific survival, and shorter progression‐free survival (Figure S5A–C, G–I). Then, according to the results of drug susceptibility test (Figure 7E,F,K,L), we found that sensitivity to 6 drugs (A‐83‐01, SB52334, SB505124, Ibrutinib, Sinularin, OF‐1) for ccRCC was negatively associated with expression of EMX2OS. Moreover, we found sensitivity to six drugs (AZD1332, MK‐2206, SB52334, Linsitinib, LDN‐193189, BMS‐754807, and Cediranib) for ccRCC was negatively associated with expression of TLN2. Altogether, it appears that these drugs would have higher sensitivity in ccRCC and improve prognostic status and survival in ccRCC patients.