Development and Validation of a Novel Metabolic Signature-Based Prognostic Model for Uveal Melanoma

Abstract

Purpose: Uveal melanoma (UM) is the most common primary malignant tumor with poor prognosis. The role of metabolism-related genes in the prognosis of UM remains unrevealed. This study aimed to establish and validate a prognostic prediction model for UM based on metabolism-related genes.

Methods: Gene expression profiles and clinicopathological information were downloaded from The Cancer Genome Atlas, and the Gene Expression Omnibus database. Univariable Cox regression, least absolute shrinkage and selection operator Cox regression, and stepwise regression were performed to establish the model. Kaplan-Meier survival analysis, receiver operating characteristic (ROC) curve analysis, and calibration and discrimination analyses were used to evaluate the prognostic model.

Results: Three metabolism-related genes, carbonic anhydrase 12, acyl-CoA synthetase long-chain family member 3, and synaptojanin 2, and three clinicopathological parameters (i.e., age, gender, and metastasis staging) were identified to establish the model. The risk score was found to be an independent prognostic factor for UM survival. High-risk patients demonstrated significantly poorer prognosis than low-risk patients. ROC analysis suggested the promising prognostic efficiency of the model. The calibration curve manifested satisfactory agreement between the predicted and observed risk. A nomogram and online survival calculator were developed to predict the survival probability.

Conclusions: The novel metabolism-based prognostic model could accurately predict the prognosis of UM patients, which facilitates the prediction of the survival probability by both ophthalmologists and patients with the online dynamic nomogram.

Translational Relevance: The dynamic nomogram links gene expression profiles to clinical prognosis of UM and is useful to evaluate the survival probability.

Introduction

Uveal melanoma (UM) is one of the most prominent primary intraocular tumors and mainly occurs between the ages of 50 and 70 years.¹ More than 90% of UM originates in the choroid; the ciliary body is a less-frequent site of origin, with an incidence of 6%, and the iris is the least-involved site, accounting for 4%.^2,3 Aronow et al.⁴ reported that the age-adjusted incidence of UM is 5.2 cases per million in the United States, which is close to that reported in European regions. It has been reported that metastasis occurred in 33% of 492 eyes with ciliary body melanoma, 25% of 7256 eyes with choroidal melanoma, and 7% of 285 eyes with iris melanoma in a 10-year follow-up study.⁵ Once metastasized, the median survival declined to only 3.7 months, and the one-year survival was only 13%.⁶ Although several risk factors for UM have been found, these common traditional clinical characteristics, such as age, gender, and staging classification, cannot accurately predict the prognosis and survival of UM patients. Therefore the identification of novel prognostic biomarkers is urgent and of utmost importance. 

Metabolic reprogramming, which nutrifies cell proliferation through regulation of energy metabolism, has been widely acknowledged as an important hallmark of cancer.⁷ Activation of oncogenes and depletion of tumor suppressor genes facilitate metabolic reprogramming in cancer cells to meet their high energy demands and accommodate changes in the tumor microenvironment.⁸ Dysregulated metabolic status is often associated with clinical outcomes. Metabolic expression subtypes with upregulated carbohydrate and nucleotide metabolism are highly correlated with poor prognosis in most human cancer types.⁹ Nevertheless, the roles of metabolism-related genes in the prognosis of UM patients remain largely unknown. 

With the development and widespread application of high-throughput screening techniques, it became achievable and accessible to obtain whole mRNA expression profiles of research subjects. The emergence of public databases, such as the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA), deposits massive microarray profiles, which facilitates researchers to investigate the molecular characterization of tumors. In this study, UM patients’ original data were retrieved from the GEO and TCGA databases. By using multiple regression analysis, we developed and validated a metabolism-related prognostic model that could accurately predict the risk of UM patients, as well as offer promising performance in assessing survival rate. 

Methods

Data Sources

The study was performed in accordance with the tenets of the Declaration of Helsinki. The gene expression profiles and corresponding clinicopathological information of 80 UM patients were downloaded from TCGA (https://portal.gdc.cancer.gov/) and the Xena database (http://xena.ucsc.edu/)¹⁰ as the training cohort (TCGA-UVM). Simultaneously, 63 UM patients from the GSE22138¹¹ dataset were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) as the validation cohort. The metabolism-related gene signatures were derived from the molecular signatures database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/).¹² 

Precursory Screen of Gene Profiles

Raw transcriptome data were processed with the affy¹³ package and sva¹⁴ package in R software version 4.0.3 (R Foundation for Statistical Computing, Vienna, Austria) for background correction and normalization. Strawberry Perl programming language version 5.32.1.1 (The Perl Foundation, Walnut, CA, USA) was used to annotate the expression matrix of the two cohorts, extract the profiles of metabolism-related gene signatures, and merge expression profiles with clinical survival information of each patient. Then the metabolism-related genes of the TCGA-UVM and GSE22138 cohorts were identified through the annotation of MSigDB. 

Establishment of Prognostic Model

To analyze the relationship between the above identified metabolism-related genes and the survival time of UM patients, univariable Cox regression and least absolute shrinkage and selection operator (LASSO) Cox regression were performed to screen the survival-related metabolic genes and build a preliminary model. Then, stepwise regression was performed to develop an optimized model that combined risk genes with clinicopathological parameters. The survival R package and glmnet R package were used to analyze and investigate the cutoff threshold of the risk score. The performance of the optimized model was assessed by the area under the receiver operating characteristic curve (AUC) using 500 iterations of bootstrap cross-validation. 

Evaluation of Prognostic Model

The Survminer R package was applied to generate the Kaplan-Meier survival curve. The pheatmap¹⁵ R package was employed to generate the risk plot. Time-dependent prognostic value was illustrated with a receiver operating characteristic (ROC) curve and assessed by AUC using the timeROC R package. Moreover, the calibration plot, net reclassification index, integrated discrimination improvement, and concordance index were used to evaluate the prognostic model. 

Construction of Nomogram and on-Line Calculator

A nomogram was built to display the probability of the clinical outcomes through a simple diagram of the statistical prediction model and to facilitate individualized estimations of prognosis for specific patients. We adopted the factors in the optimized model to develop a nomogram with the rms R package. Furthermore, an online UM survival prediction calculator was constructed using the DynNom, Shiny, and rsconnect R packages. 

Gene Set Enrichment Analysis (GSEA)

To investigate the concomitant signaling pathways involved in UM patients with high-risk scores, we performed GSEA with GSEA software version 4.1 (Broad Institute, Inc., Cambridge, MA, USA). The Kyoto Encyclopedia of Genes and Genomes (KEGG) gene set (c2.cp.kegg.v7.4.symbols) of MSigDB was used as the reference. Normalized P value < 0.05 and normalized enrichment score (|NES| >1) were considered to be significantly enriched. 

Statistical Analysis

The statistical work of the present study was conducted using R software and EmpowerStats version 2.0 (X&Y Solutions Inc., Boston, MA, USA). Spearman correlation was applied to assess the correlation between the risk metabolism-related genes. Univariable and multivariable Cox regression analyses were performed to determine the independent prognostic factors that were related to overall survival. P value < 0.05 was considered statistically significant. 

Results

Basic Demographics and Clinical Characteristics of UM Patients

The flow chart of this study is briefly illustrated in Figure 1. Due to missing TNM staging or a short follow-up period (<30 days), 9 patients in the TCGA-UVM cohort were excluded from this study, while the remaining 71 UM patients were enrolled to train the prognostic model. All 63 UM patients from GSE22138 were included in the validation cohort. The TNM staging of the patients were categorized according to the American Joint Committee on Cancer (eighth edition). The clinicopathological features of the two cohorts are listed in Table 1. 

Figure 1.

 

Brief flow chart of this study. NRI, net reclassification index; IDI, integrated discrimination improvement; C-index, concordance index.

Table 1.

 

Basic Demographics and Clinical Features of the Two UM Cohorts

 

Construction of a Prognostic Model Based on Three Metabolism-Related Gene Signatures

A total of 861 metabolism-related gene signatures from the training cohort were screened for further analysis. Twenty-three metabolism-related genes were previously identified to be associated with the overall survival time of TCGA-UVM patients through univariable Cox regression (Fig. 2A; Supplementary Table S1). All hazard ratios (HRs) of these 23 genes were greater than 1, suggesting that these genes were risk factors rather than protective factors. Then, LASSO Cox regression (Figs. 2B, 2C) was performed to build a preliminary prognostic model with nine metabolism-related gene signatures (Table 2). 

Figure 2.

 

Construction of the prognostic model. (A) Forest plot showing the hazard ratio of 23 metabolism-related gene signatures associated with overall survival time of the TCGA-UVM cohort through univariable Cox regression. (B) Selection of tuning parameter (lambda) in the LASSO model with 100-fold cross-validation. (C) LASSO coefficient profiles of 23 metabolism-related genes in the TCGA-UVM cohort. Each coefficient profile plot was generated versus the log (lambda) sequence. (D) ROC analysis of the preliminary model (black curve) and the optimized model (red curve) with 500 iterations of bootstrap resampling. (E, F) The risk score distribution of the patients from the TCGA-UVM cohort (E) and GSE22138 cohort (F) were plotted in ascending order and colored as high-risk (red) and low-risk (green).

Table 2.

 

Nine Metabolism-Related Gene Signatures Associated With Prognosis in UM Patients

 

Furthermore, stepwise regression was used to establish an optimized prognostic model that consisted of three risk genes and three clinical parameters (Table 3). The three genes included in the model were carbonic anhydrase 12 (CA12), acyl-CoA synthetase long chain family member 3 (ACSL3), and synaptojanin 2 (SYNJ2). The three clinicopathological parameters were age, gender, and metastasis staging (M) without multicollinearity relations. Three patients in the TCGA-UVM cohort were excluded from further analysis due to their Mx stage. The risk score of the optimized model was calculated by the following formula: risk score = 0.09780 × expression of CA12 + 0.20650 × expression of ACSL3 + 0.24480 × expression of SYNJ2 + 0.04084 × age + 3.64541 × (M staging: M0 = 0, M1 = 1) + 0.89711 × (gender: female = 0, male = 1). ROC analysis with 500 iterations of bootstrap resampling revealed that the optimized model was significantly preferable to the preliminary model due to a higher AUC and concordance index (P < 0.001, Fig. 2D; Table 4). Subsequently, the UM patients of the two cohorts were divided into high- and low-risk groups based on the cutoff value of 5.03839 for the risk score in the training cohort. As displayed in Figures 2E and 2F, the numbers of high-risk and low-risk patients in the TCGA-UVM cohort were equal due to our settings, whereas the number of high-risk patients in GSE22138 was greater than that of low-risk patients. 

Table 3.

 

Three Metabolism-Related Gene Signatures and Three Clinicopathological Parameters Constitute the Optimized Model

 

Table 4.

 

Comparison of the Optimized Model and the Preliminary Model With ROC Analysis

 

Cox Regression Analyses Revealed Independent Prognostic Value of the Risk Score

Univariate and multivariate Cox regression analyses were performed among the available covariables to investigate the independent prognostic value of the risk score. In the TCGA-UVM cohort, univariate Cox regression analysis showed that the risk score was significantly associated with overall survival time (Fig. 3A, P < 0.001, HR = 2.041, 95% confidence interval [CI] = 1.515-2.751). After adjusting confounding factors, we found that the risk score remained an independent predictor for overall survival in the multivariate Cox regression analysis (Fig. 3B, P < 0.001, HR = 3.967, 95% CI = 1.761-8.937). In the GSE22138 cohort, consistent results were observed in univariate (Fig. 3C) and multivariate (Fig. 3D) Cox regression analyses (P < 0.001, HR = 6.584, 95% CI = 2.719-15.946; P < 0.001, HR = 14.013, 95% CI = 2.874-68.332, respectively). 

Figure 3.

 

Prognostic value of the risk score and three metabolism-related gene signatures. In the TCGA-UVM cohort, both univariate (A) and multivariate (B) Cox regression analyses indicate that the risk score was a powerful independent predictor associated with overall survival (P < 0.001, HR = 2.041, 95% CI = 1.515-2.751; P < 0.001, HR = 3.967, 95% CI = 1.761-8.937, respectively). In the GSE22138 cohort, both univariate (C) and multivariate (D) Cox regression analyses indicate that the risk score was a powerful independent predictor associated with overall survival (P < 0.001, HR = 6.584, 95% CI = 2.719-15.946; P < 0.001, HR = 14.013, 95% CI = 2.874-68.332, respectively). Spearman correlation show weak correlations between the three metabolism-related gene signatures in the TCGA-UVM cohort (E) and GSE22138 cohort (F). A heatmap shows that the three metabolism-related gene signatures were upregulated in the high-risk group in the TCGA-UVM cohort (G) and GSE22138 cohort (H).

In addition, Spearman correlation showed that the correlations of the three metabolism-related gene signatures in the two cohorts were all weak (Figs. 3E, 3F), which met the requirement of establishing a predictive model. Moreover, the expression of the three metabolism-related gene signatures was upregulated in the high-risk group in both the TCGA-UVM and GSE22138 cohorts (Figs. 3G, 3H). 

The Prediction Model Based on Metabolism-Related Gene Signatures Showed Strong Power in Prognostic Evaluation

Kaplan-Meier survival analysis suggested that the patients in the TCGA-UVM cohort in the high-risk group had poorer overall survival than those in the low-risk group (Fig. 4A, P < 0.001). As shown in Figure 4B, none of the low-risk group patients in the GSE22138 cohort died during the follow-up period. Consistently, significantly poorer overall survival was detected in the high-risk group (P < 0.001). The mortality proportion in the high-risk group was higher than those of the low-risk group in the training cohort (35% vs. 6%, Fig. 4C) and validation cohort (71% vs. 0%, Fig. 4E). Likewise, patient status distribution showed that more death events occurred in the high-risk groups of the training cohort (Fig. 4D) and validation cohort (Fig. 4F). 

Figure 4.

 

Prognostic evaluation performance of the prediction model. Kaplan-Meier curve of the training cohort (A) and validation cohort (B). (C) Bar plot shows the mortality proportions in the high- and low-risk groups of the training cohort. (D) Dot plot displays patient status distributions and death event occurrences in the high- and low-risk groups of the training cohort. (E) Bar plot shows the mortality proportions in the high- and low-risk groups of the validation cohort. (F) Dot plot displays patient status distributions and death event occurrences in the high- and low-risk groups of the validation cohort. ROC analysis reveals the predictive efficiency of the optimized model in one-, three-, and five-year survival probability in the training cohort (G) and validation cohort (H). (I) Calibration curve of the optimized model showed satisfactory agreement between the predicted risk (black curve) and observed risk (red diagonal).

Furthermore, the predictive power for one-, three-, and five-year survival probability achieved AUCs of 0.959, 0.825, and 0.944 in the training cohort (Fig. 4G) and 0.899, 0.990, and 0.995 in the validation cohort (Fig. 4H), indicating the promising prognostic efficiency of the optimized model. The calibration curve manifested satisfactory agreement between the predicted risk and observed risk (Fig. 4I). As displayed in Table 5, integrated discrimination improvement suggested that the predictive power of the optimized model was 32.40% higher than that of the preliminary model (P = 0.0070), whereas the net reclassification index indicated 54.90% improvement in the power of the optimized model to correctly classify high- and low-risk patients compared with the preliminary model (P = 0.0070). 

Table 5.

 

Improvement of the Predictive Power of the Optimized Model Compared With the Preliminary Model

 

Nomogram and On-Line Calculator Developed for the Established Model Facilitate Prediction of Survival Probability

A nomogram was established to predict one-, two-, three-, four-, and five-year overall survival probability using the variables for calculating the risk score, including age, gender, M stage, and the expression of CA12, ACSL3, and SYNJ2 (Fig. 5A). Each variable was allocated points in accordance with its risk contribution. We further developed an online calculator for UM prognostic prediction which can be accessed at https://keshi2021.shinyapps.io/dynnomapp/. Ophthalmologists and patients could also scan the QR code in Figure 5B to access the website. 

Figure 5.

 

A nomogram and online calculator were developed for the prediction model. (A) The nomogram was established by using the parameters of age, gender, M stage, and the expression of CA12, ACSL3, and SYNJ2. (B) QR code of our online calculator website to predict the survival probability. (C) The user interface of the online calculator, including the input section on the left and output section on the right. The forest plot on the right section displays multiple prediction results.

The interface of our web calculator, consisting of two sections, is displayed in Figure 5C. The left input section comprises the six variables of the optimized model, and users can drag circular buttons or drop-down boxes for parameter setting. In checking the “Predicted survival at this follow-up,” users could select different follow-up time points by dragging the circular button. After clicking the “Predict” button, the right output section shows the survival plot and predicted probability, accompanied by 95% CI. The calculator supports displaying multiple predictions on the same plot. A reminder is required that after finishing the prediction, the users need to click the “Quit” button on the left lower position; otherwise, the next user might encounter an error. In the case that users encounter an error when predicting, they could click the “Quit” button and reload the website. 

GSEA Demonstrated Upregulation of Metabolism-Related Pathways in the High-Risk Groups

GSEA was performed between the high- and low-risk groups in the two cohorts. We found that the KEGG-Proteasome pathway and KEGG-ATP-binding cassette transporter pathway were upregulated in the high-risk groups of both cohorts (Figs. 6A, 6B). In addition to these two pathways, six metabolism-related pathways were significantly upregulated in the high-risk patients of the TCGA-UVM cohort, including KEGG-glutathione metabolism, KEGG-galactose metabolism, KEGG-amino sugar and nucleotide sugar metabolism, KEGG-pyrimidine metabolism, KEGG-glycerolipid metabolism, and KEGG-porphyrin and chlorophyll metabolism (Fig. 6C), whereas the KEGG-nicotinate and nicotinamide metabolism pathways were significantly upregulated in the GSE22138 cohort (Fig. 6D). 

Figure 6.

 

The significantly enriched KEGG pathways in the high-risk group of the TCGA-UVM cohort and GSE22138 cohort determined by GSEA. Common upregulated KEGG-Proteasome pathway and KEGG-ABC transporters pathway in the high-risk group of the TCGA-UVM cohort (A) and GSE22138 cohort (B). (C) Six upregulated metabolism-related KEGG pathways in the TCGA-UVM cohort. (D) One upregulated metabolism-related KEGG pathway in the GSE22138 cohort. ABC, ATP-binding cassette.

Discussion

UM is a worldwide intraocular malignancy that is especially common in elderly individuals.¹⁶ UM manifests a high mortality rate, and approximately half of cases metastasize, mainly to the liver.^1,17 Once diagnosed with metastasis, the overall survival time is as low as approximately one year.¹⁸ To date, no proven treatment can improve overall survival for UM patients. Metabolic process dysregulation is a putative hallmark in the pathogenesis and progression of tumors.¹⁹ Nevertheless, few studies have explored the association between metabolism-related gene signatures and the risk and prognosis of UM with high-throughput sequencing.^20,21 

In the present study, we systematically explored the associations between metabolism-related gene signatures in UM tissues and the overall survival of UM patients from the TCGA database. After a series of regression analyses, a novel prognostic model integrating three metabolism-related gene signatures and three clinicopathological parameters was established and validated in an external GEO cohort with promising performance. 

Three clinicopathological parameters, including age, gender, and M stage, have been recognized as risk factors for UM in previous publications.^22,23 UM is rarely observed in children and teenagers but is more commonly diagnosed in elderly people, with a gradual increase in age-specific incidence that reaches a peak at 70 years of age and a plateau after 75 years.²² A systematic review revealed that the age-adjusted incidence rate of UM is higher in male patients than in female patients.²⁴ As mentioned above, the presence or absence of metastasis is closely related to the prognosis of UM patients. Therefore the inclusion of these three clinical characteristics into our model was consistent with the findings of these previous reports. 

CA12, a member of the carbonic anhydrase family of zinc metalloenzymes, catalyzes the reversible interconversion between carbon dioxide and bicarbonate.²⁵ The role of CA12 in tumorigenesis remains controversial. Li et al.²⁶ reported that upregulated CA12 expression was related to good prognosis of breast cancer, whereas high expression of CA12 in glioblastoma predicted poor prognosis.²⁷ In our study, the HR value of CA12 was 1.1027 (Table 3; 95% CI: 1.0559-1.1517, P < 0.0001), indicating its risk factor role in UM. ACSL3, an isozyme of the long-chain fatty-acid-coenzyme A ligase family, plays a crucial role in lipid and fatty acid metabolism.²⁸ In non-small cell lung cancer, ACSL3 was highly expressed and promoted tumor cell proliferation, migration, and invasion, thus identifying it as an unfavorable prognostic biomarker.²⁹ Moreover, ACSL3 was overexpressed in hepatocellular carcinoma and hepatic gastrointestinal metastasis tissues as compared with healthy controls.³⁰ SYNJ2, a member of the inositol-polyphosphate 5-phosphatase family, regulates various critical cellular processes, including membrane trafficking, nucleation of actin filaments, and the activity of ion channels and transporters.³¹ Csolle et al.³² reported that SYNJ2 was upregulated in breast cancer and promoted cancer development via AKT-dependent and AKT-independent mechanisms. In a case-control study, SYNJ2 variants were suggested to play a crucial role in the progression of colorectal cancer.³³ However, fewer studies have reported the correlation of these three metabolism-related gene signatures with UM. In our present work, the results of univariable and LASSO Cox regression suggested the increased risk effect of these three genes. 

The risk score is a generally used approach for the development of a valid signature. A series of validation analyses, including ROC and Kaplan-Meier analyses, indicated that the risk score of our model offers the ability to accurately predict the prognosis of UM patients. The calibration plot and discrimination analysis demonstrated the excellent predictive performance of the model, especially the C index value of 0.883 (Table 4; 95% CI: 0.841-0.926), suggesting the perfect discrimination capacity of our model to distinguish high- and low-risk patients. 

Our newly established model offers some advantages as compared with the existing UM prognostic model^34–38 (e.g., visualized nomograms and online calculators facilitate informed treatment and management decisions for ophthalmologists and patients). Second, other models need to examine nine,³⁴ 10,³⁵ 22,³⁶ or 18³⁷ genes or five alternative splicing events,³⁸ whereas only three genes were required for prediction by using our model, thus enabling greater convenience and speed. Third, calibration and discrimination analyses were implemented to prove the predictive power of our model and make it more accurate for prediction. 

There are some limitations in this study. First, some researchers have suggested that the latitude of the patient's country,³⁹ skin color,⁴⁰ occupational cookery status,⁴¹ and radiation exposure⁴² would also be risk factors for UM. However, we did not include these metrics in the process of predictive model establishment because this information was missing from the available profiles. Second, our data were derived from online public databases and have limitations such as missing information of disease-free survival, as well as insufficient sample size, which limits the scope of the model. For example, the minimum age of patients for which our nomogram and web calculator were applicable was 22 years old. Third, this study requires eyeball samples that were obtained from enucleation surgery. Patients of T1 stage and clinical stage І were not included because enucleation surgery were not suitable for these patients with relative early stage of the tumor. 

In conclusion, we established and validated a novel prognostic model for the survival prediction of UM patients, facilitating ophthalmologists and patients to accurately evaluate the personalized survival probability of UM patients and distinguish high‐risk individuals who might need more frequent follow‐up strategies. These three metabolism-related genes could serve as potential biomarkers and therapeutic targets for UM. 

Acknowledgments

Supported by grants from the National Natural Science Foundation of China (81730026, and 81960158), Science and Technology Commission of Shanghai Municipality (20Z11900400), Shanghai Collaborative Innovation Center for Translational Medicine (TM201917), Shanghai Hospital Development Center (SHDC2020CR2040B), the Cultivation Program of National Natural Science Foundation for Outstanding Youth (20202ZDB01014), and the Key Program of Youth Science Foundation of Jiangxi Province (20202ACBL216009). 

Disclosure: K. Shi, None; X. Li, None; J. Zhang, None; X. Sun, None