STON1 expression was significantly reduced in KIRC
We first identified the transcription level of STON1 in different TCGA tumors with the TIMER2.0 database. A total of 61 human cancer types were explored. We found that most cancers, including KIRC, kidney renal papillary cell carcinoma (KIRP), kidney chromophobe (KICH), bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), colon adenocarcinoma (COAD), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), stomach adenocarcinoma (STAD), thyroid carcinoma (THCA), and uterine corpus endometrial carcinoma (UCEC), showed decreased STON1 mRNA, while cholangiocarcinoma (CHOL), glioblastoma multiforme (GBM), and head and neck squamous cell carcinoma (HNSN) displayed upregulated STON1 mRNA levels (Figure 1a). Moreover, results from three GEO datasets—GSE16441, GSE16449, and GSE71963—also revealed aberrantly downregulated STON1 mRNA in KIRC (P<0.01, Figure 1b, 1c, 1d). All three series were normalized before the analysis (Figure S1). Finally, the protein level of STON1 in KIRC was also found to be reduced using the CPTAC database (P < 0.0001) (Figure 1e).
Correlation between STON1 expression and clinicopathological characteristics of KIRC patients
The KIRC GDC Expression Matrix was downloaded from the UCSC Xena, which is one of the most comprehensive clinical oncology databases [33]. A total of 491 KIRC patients were included in the cohort after excluding non-KIRC diagnoses and missing clinical information, and were classified into a high expression group (n = 246) and a low expression group (n = 245) according to the median value of STON1. Notably, STON1 mRNA levels were related to grade (P = 0.022), TNM stage (P = 0.039), distant metastasis (P = 0.019), and vital status (P = 5.46e−6). However, no correlation was found between STON1 mRNA levels and age (P = 0.136), sex (P = 0.949), invasion depth (P = 0.074), and lymph node metastasis (P = 0.385) (Table 1).
STON1 was a protective prognostic factor in KIRC
To further analyze the prognostic value of STON1, we performed a survival analysis using four oncological databases. The results indicated that STON1 expression levels were strongly associated with OS in KIRC. Patients with high expression of STON1 had a better outcome compared with those with low levels of STON1 in the OncoLnc database (log rank P = 9.78e−8) (Figure 2a). Similar results were acquired from other oncology databases, including the UALCAN database, the Kaplan–Meier database, and GEPIA2 with TCGA data (Figure 2b, 2c, 2d).
Correlations of STON1 levels with tumor mutation burden and mismatch repair genes in KIRC
Given the sensitive connection between TMB, microsatellite instability (MSI), and immunotherapeutic response, we focused on their relationship with STON1 expression levels. Marginal scatter plots revealed that alterations in STON1 expression were always accompanied by changes in mismatch repair proteins (MLH1, PMS2, MSH2, MSH6, and EpCAM). Moreover, the four mismatch proteins seem to be positively associated with STON1 (Figure 3a, P < 2.2e−16). Additionally, STON1 mRNA level was negatively associated with the TMB score (P = 0.012) (Figure 3b).
Association between STON1 expression and immune cells in KIRC
STON1 gene copy numbers were accompanied by changes in immune cell infiltration levels, including that of B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells (Figure 4a). Subsequently, we examined the relationship between STON1 expression levels and the infiltration levels of 26 immune-related cell types. Our data illustrated that the abundance of different subtypes of immune cells was significantly positively (effector memory CD8 T cells, effector memory CD4 T cells, type 1 T helper cells, type 2 T helper cells, regulatory T cells, memory B cells, natural killer cells, immature dendritic cells, eosinophils, mast cells, and neutrophils) or negatively (activated CD8 T cells, CD56 bright natural killer cells, CD56 dim natural killer cells, myeloid-derived suppressor cells, activated dendritic cells, and plasmacytoid dendritic cells) related to STON1 expression (Figure 4b, 4C). These results demonstrated that STON1 was positively linked with most immune cells. Furthermore, STON1 mRNA expression was obviously highest in the C3 subtype (inflammatory) (Figure 4d), which was reported to have the best prognosis in KIRC [34]. This result indicated that STON1 may have a protective role in KIRC by exhibiting upregulated expression in the C3 immune subtype.
STON1 may mediate immune response through immunomodulators
Thirty-one immunomodulators were obviously related to STON1 in KIRC, including 21 immunostimulators (C10orf54, CD27, CD80, CXCL12, CXCR4, ENTPD1, HHLA2, IL6, IL6R, LTA, NT5E, PVR, RAET1E, TMIGD2, TNFRSF14, TNFRSF17, TNFRSF18, TNFRSF25, TNFSF14, TNFSF15, and TNFSF9) (Figure 5a) and 10 immunoinhibitors (ADORA2A, CD274, CTLA4, IDO1, KDR, LAG3, PDCD1, PDCD1LG2, PVRL2, and TGFBR1) (Figure 5b). We subsequently explored the possible signaling pathways through which STON1 may modulate the immune response in KIRC using GO and KEGG enrichment analyses. Thirty-one immunomodulators and 57 closely related genes were uploaded onto the DAVID. Biological process analysis revealed that these genes were enriched in the tumor necrosis factor (TNF)-mediated signaling pathway, T cell co-stimulation, and immune response. Molecular functional analysis revealed that they were enriched in TNF receptor binding, TNF-activated receptor activity, and cytokine activity. Cellular component analysis indicated that STON1 was mainly enriched in the plasma membrane (Figure 6a). KEGG enrichment analysis further confirmed that these genes are tightly related to cytokine–cytokine receptor interactions (Figure 6b). The associations between 31 immunomodulators and 57 closely related genes are presented in a protein-to-protein interaction network (Figure 6c)
Prognostic value of STON1-related immunoregulators in KIRC
In total, 504 patients were randomly divided into a training set (n = 328) and a validation set (n = 176) at a ratio of 65%–35%, which was downloaded from GDC. To explore the prognostic value of STON1-related immunomodulators, univariate Cox analysis was first performed on the training set, which identified 15 of the abovementioned immunoregulators with statistical significance (Table 2). We then input these 15 immunomodulators into a stepwise multivariate Cox regression model using the training set. CTLA4 (P = 0.0177), KDR (P = 0.0135), PDCD1 (P = 0.024), and HHLA2 (P < 0.001) were independent prognostic factors, as shown in Figure 7a. Here, risk scores were obtained using the following formula:
Risk score = EXP(CTLA4) × 0.25599+ EXP(KDR) × (−0.18132) + EXP(LAG3) × 0.18505 + EXP(PDCD1) × (−0.25896) + EXP(HHLA2) × (−0.16182) + EXP(TNFRSF25) × 0.13274 + EXP(TNFSF14) × 0.13404.
The Kaplan–Meier survival curve indicated that the high-risk score group had worse OS than the low-risk score group in KIRC (Figure 7b). Through time-dependent ROC curves, the Cox model displayed better predictive ability. The AUC of 1-, 3- and 5-year OS was 0.732, 0.733, 0.825, respectively (Figure 7d). The risk factor analysis graph presents the complete landscape of risk scores, patient status, and gene expression profiles (Figure 8a). We used the validation set to verify the predictive value of the seven-gene prognostic signature. Consistent with these results, the Kaplan–Meier curves of the validation set showed that the low-risk group had a better outcome than the high-risk group (P = 7e−8) (Figure 7c). The AUC of 1-, 3- , 5-year OS of our model was 0.795, 0.724, 0.784, respectively (Figure 7e). The risk factor analysis chart is also shown in Figure 8b.
Construction of a prognostic nomogram according to risk scores and clinical parameters
Univariate Cox proportional hazards regression analysis demonstrated that risk score (HR = 1.289, 95% CI = 1.217–1.365, P = 5.81e−18), stage (HR = 1.819, 95% CI = 1.532–2.160, P = 8.62e−12), M (HR = 3.858, 95% CI = 3.573–5.784, P = 6.47e−11), T (HR = 1.892, 95% CI = 1.529–2.342, P = 4.74e−9), and age (HR = 1.039, 95% CI = 0.021–1.056, P = 9.52e−6) were distinctly associated with OS in the training set (Figure 9a). Meanwhile, multivariate Cox proportional hazards regression analysis was implemented to assist with the identification of independent prognostic factors. The independent prognostic factors were risk score (HR = 1.197, 95% CI = 1.127–1.271, P = 6.71e−9), age (HR = 1.038, 95% CI = 1.020–1.057, P = 3.52e−5), and stage (HR = 1.699, 95% CI = 1.420–2.033, P = 6.76e−9) (Figure 9b). Therefore, a comprehensive nomogram prognostic model was constructed according to all independent prognostic factors screened by multivariate Cox proportional hazards regression analysis of the training set to predict the probability of 1-, 3-, and 5-year OS (Figure 9c). In the training set, AUC of 1-, 3-, and 5-year OS were 0.889, 0.804, 0.790, respectively (Figure 9d). In the validation set, AUC of 1-, 3-, and 5-year OS were 0.882, 0.847, 0.809, respectively (Figure 9e). The AUC of the training set and validation set illustrated our model with good predictive ability for OS. Furthermore, the calibration curves, which reflected the consistency between the predicted risk and the actual risk, showed that the nomogram-predicted probability was well-matched with the ideal line for 1-, 3-,5-year OS in the training set and validation set (Figure 10a). The same methods were applied to the validation set to further validate the authenticity of the above results (Figure 10b). The concordance probability (C-index) can be used to measure and compare the discriminative power of a risk prediction model. The C-indexes of our nomogram model were 0.783 and 0.795 in the training set and validation set, respectively. Moreover, we also discovered that our nomogram prognostic model showed a better net benefit for predicting OS compared with age, stage, and risk score in the training set and validation set (Figure 10c, 10d).