Identification of Prognostic Signatures for Predicting the Overall Survival of Uveal Melanoma Patients

Uveal melanoma (UM) is an aggressive cancer which has a high percentage of metastasis and with a poor prognosis. Identifying the potential prognostic markers of uveal melanoma may provide information for early detection of metastasis and treatment. In this work, we analyzed 80 uveal melanoma samples from The Cancer Genome Atlas (TCGA). We developed an 18-gene signature which can significantly predict the prognosis of UM patients. Firstly, we performed a univariate Cox regression analysis to identify significantly prognostic genes in uveal melanoma (P<0.01). Then the glmnet Cox analysis was used to generate a powerful prognostic gene model. Further, we established a risk score formula for every patient based on the 18-gene prognostic model with multivariate Cox regression. We stratified patients into high- and low-risk subtypes with median risk score and found that patients in high-risk group had worse prognosis than patients in low-risk group. Multivariate Cox regression analysis demonstrated that 18-gene model risk score was independent of clinical prognostic factors. We identified four genes whose mutations were closely to UM patients' prognosis or risk score. We also explored the relationship between copy number variation and risk score and found that high risk group showed more chromosome aberrations than low risk group. Gene Set Enrichment Analysis (GSEA) analysis showed that the different biological pathways and functions between low and high risk group. In summary, our findings constructed an 18-gene signature for estimating overall survival (OS) of UM. Patients were categorized into two subtypes based on the risk score and we found that high risk group showed more chromosome aberrations than low risk group.


Introduction
Uveal melanoma, the most common primary intraocular malignancy, mainly has an increased risk of liver metastatic involvement. Though effective therapy, surgical enucleation and radiotherapy have improved local control, however up to 50% patients will eventually die of their disease [1,2]. Uveal melanoma has been considered with the worse prognosis and poorer response to chemotherapy compared with cutaneous melanoma. Therefore, identifying the potential prognostic markers of uveal melanoma may provide information for early detection of recurrence and treatment. Although many studies have identified some important genes and pathways in the diagnosis and treatment of uveal melanoma, the prognosis of uveal melanoma patients was still very poor [3,4]. Hence, it's urgent to unveil new markers to evaluate their prognosis.
Nowadays, there are lots of emerging high-throughput sequencing technologies and databases for development of diagnostic and prognostic signatures of cancer. Although a 15-gene expression panel of UM patients was explored to determine the risk of metastasis, the signature for prognostic results remained to be further elucidated Ivyspring International Publisher [5]. Gene mutations and chromosome copy number variants strongly correlated with UM outcome. UM patients also display a pattern of mutations. GNAQ and GNA11 mutations are reported to promote cell proliferation and metastasis [6]. BAP1 mutations are related to metastasis and EIF1AX mutations are associated with favorable metastatic-free survival [7][8][9][10]. Uveal melanoma patients with chromosome 3 copy number loss are associated with a high risk of metastasis and a poor outcome [11,12]. It was reported that UM tumor progression typically displays chromosomal abnormalities. UM patients with chromosome 3 loss, 8q gain and 1p loss have reduced overall survival [12]. Therefore, exploring the gene mutations and copy number variations may provide insights into uveal melanoma prognosis.
In the present paper, we conducted uveal melanoma data by using the publicized data from the TCGA database. By using glmnet Cox model and Cox regression analysis, we developed an 18-gene prognostic model to demonstrate the association between 18-gene model and prognostic power of uveal melanoma. Meantime, we also found the specific gene mutations and chromosome copy number variations involved in uveal melanoma associated with overall survival. By utilizing 18-gene prognostic model, gene mutations and copy number variations of UM may provide evidence for the selection and determination of an individualized and targeted therapeutic treatment for each patient.

Uveal melanoma data collection and analysis from TCGA
We retrospectively collected uveal melanoma gene expression data and corresponding clinical information of TCGA from UCSC Xena Public Data Hubs (https://xena.ucsc.edu/public-hubs/). All expression values were transformed with log2 (FPKM+1). Copy number profiles from Affymetrix Genome-Wide Human SNP Array 6.0 platform and segments were mapped to hg38 genome assembly (https://xena.ucsc.edu/public-hubs/). Mutation data containing somatic variants was stored in the form of Mutation Annotation Format (MAF) and was downloaded from Genomic Data Commons (GDC) (https://portal.gdc.cancer.gov/). There were 80 patients in our analysis.

Construction of a prognostic gene model
Univariable Cox regression analysis was performed to identify gene which was considered statistically significant if its p value was less than or equal to 0.01. We further used glmnet Cox model to trained these identified genes for 2000 times independently (R package, glmnet, version 2.0-16), and finally selected a best suitable prognostic model with the highest frequency [13,14]. With this prognostic gene model, the prognostic genes were strictly included. Subsequently, a risk score was established using the sum of each prognostic gene expression values weighted by their regression coefficients as previously described [15,16]. Based on this formula, the risk score for each patient was calculated. Then the uveal melanoma cohort was separated into high-and low-risk groups using the median risk score as a cut-off [17,18]. N represented the number of prognostic model genes, Coei was the regression coefficient value of each gene and Expi represented the expression of prognostic model genes.

Gene mutations and copy number variants analysis
The significant mutation genes reported in uveal melanoma were used to investigate the differences between high and low risk group [9]. The relationship between risk score and chromosomal copy number variations (CNVs) were explored. The R package 'copynumber' was used to show the CNVs according the order of risk scores [19].

Functional enrichment analysis of prognostic gene model
The Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways between high and low risk group were identified by The Gene Set Enrichment Analysis (GSEA) (https://pypi.org/project/gseapy/) [20].

Statistical analysis
The time-independent receiver operating characteristic curve (ROC) was to assess the risk prediction power of the prognostic gene model for overall survival (R package, survivalROC, version 1.0.3). Kaplan-Meier survival analysis together with a log-rank test was implemented to compare the difference of overall survival between high-and low-risk cohort using R package "survival" [21]. R package "glmnet" was used to perform Cox proportional hazards regression with the least absolute shrinkage (glmnet, version 2.0-16). Multivariate Cox regression analysis was performed to provide evidence that the prognostic gene model was independent of other clinicopathological factors [22]. All analysis was performed with R software (version 3.3.0). P<0.05 was considered statistically significant.

Construction of an 18-gene prognostic model
As a result, 4388 significantly prognostic genes in 80 uveal melanoma samples were identified by using univariate Cox regression analysis (p ≤0.01). And then a Cox proportional hazards regression model from glmnet was used to select the key genes played important roles in the prognosis of patients. Ten-fold cross-validation for the Cox model was performed to estimate the penalty parameter to minimize the risk of overfitting. After 2000 iterations, we ranked fifteen significant prognostic gene models according to their frequency ( Figure 1A, Supplementary 1). We selected the one with the highest frequency as our final prognosis gene model, including 18 genes, of which are S100A13, FAM189A2, ZBED1, RNF208, TCTN1, SIRT3, AC092821.1, AL137784.1, ZNF497, GRIN2A, AC010442.3, AC023790.2, FABP5P1, CA12, PARP8, MIR4655, MGLL, and MMP9 ( Figure 1B, Table 1). Kaplan-Meier survival curves analysis of the 18 genes for overall survival of uveal melanoma were shown in Figure 2. We further calculated a risk score for each patient in uveal melanoma samples with the risk score formula. Then patients were divided into high-and low-risk groups using the median risk score value as cut-off to investigate the prognostic role of the 18-gene signature model in overall survival. Next, principal components analysis (PCA) based on the 18 gene expression demonstrated a different separation pattern between high-and low-risk group ( Figure  1C), indicating distinct phenotypes among uveal melanoma. Further, ROC curve was performed to assess the prognostic efficiency of the model. As a prediction figure, the final model of prediction power for one year, three years and five years achieved an area of under curve of 0.803, 0.873 and 0.801, separately ( Figure 1D). Utilizing our predictive model of uveal melanoma, we could achieve early prediction of distinguishing poor and good prognosis of patients.

Prognostic role of an 18-gene model in uveal melanoma
After constructing an 18-gene prognostic model, we further performed a multivariate Cox regression analysis to confirm its power for independently predicting prognosis. The results showed that the 18-gene model risk score was independent of stage, gender, and recurrence for overall survival in UM patients (HR=9.298, 95%CI=4.381-19.734, P<0.001, Figure 3A, Table 2). The hazard ratio (HR) value of risk score is greater than one, also much greater than the HR values of stage, gender, and recurrence, which means that UM patients whose risk score is high will have a worse prognosis. After separating UM patients into high-and low-risk groups, we sought to investigate the correlation between two groups and overall survival of uveal melanoma. The Kaplan-Meier curve demonstrated that uveal melanoma patients with high risk score had poorer survival than those with low risk score patients (log-rank test, p<0.001) ( Figure 3B). Consistently with the patients observed in the high-risk group had higher recurrence rate and larger dead numbers compared with low-risk group ( Figure 1B). What's more, we also investigated the 18-gene prognostic signature in UM patients at pathological stage II and stage III-IV. As a result, the 18-gene signature model could predict UM patients with different prognosis in different subgroups including pathological stage II, stage III-IV ( Figure 3C, 3D).

Identification of gene mutations and copy number variants in uveal melanoma
Uveal melanoma displays gene mutations and copy number variations that correlate strongly with clinical outcome which are not present in cutaneous melanoma. Therefore, we investigated if there are gene mutation differences between high-and low-risk groups. Ten significantly mutated genes were detected: GNAQ, GNA11, BAP1, SF3B1, EIF1AX, COL14A1, CYSLTR2, MACF1, MUC16 and MYOF. As shown in Figure 4A, patients with SF3B1 mutations and EIF1AX mutation were mainly distributed in low-risk group, while patients with BAP1 mutations mainly emerged in high-risk group. BAP1 mutations were identified in 28% uveal melanoma patients, and the types of mutations included frame shift deletion, missense mutation, nonsense mutation, splice site and frame shift insert. SF3B1 and EIF1AX mutations are mainly missense mutations.   Patients with SF3B1 and EIF1AX wildtypes have higher risk score than patients with mutation types ( Figure 4D, 4E). The Kaplan-Meier curve showed that the overall survival in SF3B1 mutation group was significantly longer than in wildtype group (p=0.0071, Figure 4H), while the overall survival of EIF1AX in mutation cohort versus wildtype cohorts had no statistical difference (p=0.2547, Figure 4I). BAP1 risk score was much higher in mutation group than in wildtype group (p=6×10 -4 , Figure 4C), adversely compared to GNAQ, SF3B1, and EIF1AX. Notably, there was no significantly prognostic difference between BAP1 mutated group and wildtype group (p=0.0875, Figure 4G). Meanwhile, we found that patients with GNAQ mutation have more favorable prognosis than patients in GNAQ wildtype group (p=0.0426, Figure 4B, 4F). As shown in Figure 5, we found that UM patients in high-risk group have much more chromosome 3 deletion and chromosome 8 gain. These specific gene mutations and chromosome copy number variations (CNVs) can also be a prognostic biomarker for uveal melanoma patients.

Biological process identified by Functional enrichment analysis
The stratification power of the 18-gene prognostic model in predicting overall survival of UM could be ascribed to their roles in tumor development and metastasis progress. Therefore, GSEA analysis was performed to identify associated enriched genes or pathways in datasets. We selected top 10 ranked KEGG sets and GO sets results shown in Figure 6 and Figure 7. In the results, we found that the "p53 signaling pathway" were enriched in the high-risk group patients. Several studies have investigated this pathway is associated with the prognosis of uveal melanoma. Young Sun et al. demonstrated that p53 signaling pathway played a major role in mediating cellular response to UM radiation-induced DNA damage and may be defective in UM. Further, other cancer related pathways and genes such as "Toll-like receptor signaling pathway", "cytokine receptor interaction", "wide pore channel activity", "B cell differentiation", "death receptor activity", "regulation of cell adhesion mediated by integrin", "cellular response to heat" and "Notch receptor processing", were enriched in high-risk UM groups. In conclusion, functional enrichment analysis results showed that this 18-gene prognostic model may play a significant role in UM development and biological progress.

Discussion
In our study, we used univariate Cox regression analysis and glmnet Cox analysis method to deal with the 80 uveal melanoma samples from the TCGA. As a result, we identified an 18-gene prognostic model, including S100A13, FAM189A2, ZBED1, RNF208, TCTN1, SIRT3, AC092821.1, AL137784.1, ZNF497, GRIN2A, AC010442.3, AC023790.2, FABP5P1, CA12, PARP8, MIR4655, MGLL, and MMP9 ( Figure 1B, Table 1). As we know, there were several known cancer-specific biomarkers of the model that had already been translated into clinical prognostic signatures. For example, low levels or no expression of S100A13 may be one of the key predictive markers to identify melanoma patients responding to dacarbazine\temozolomide chemotherapy [23]. Daniela et al. demonstrated S100A13 as a new angiogenic indicator that facilitates human UM progression and as a promising prognostic marker [24]. Jasmine et al. showed that SIRT3 overexpression played a pro-proliferative function in melanoma and was essential for melanoma growth and survival [25]. Chen et al. revealed that MMP9 played vital roles in facilitating the metastasis of uveal melanoma cells and the involvement of MMPS in cancer progression has been reported in various cancer cell types [26][27][28]. Moreover, Harbour et al. presented a 15-gene signature that could distinguish UM patients between two groups with high risk and low risk of metastasis without regarding to copy number variations and gene mutation status [5]. A. Gordon Robertson et al. performed a comprehensive multi-platform analysis of 80 cases of primary uveal melanoma from the Cancer Genome Atlas, finally they identified and characterized four different subtypes with unique genomic aberrations, transcriptional features and clinical outcomes [9]. Their study mainly aimed at revealing distinct molecular and clinical UM profiles, emphasizing the need for stratified UM patient management. While our study built up a prognostic gene model to predict the overall survival of UM patients, aiming at patients' prognosis but not molecular subtypes. Although we also investigated the gene mutations and copy number variations of UM patients. However, importantly, we analyzed them in order to find if there were different gene mutations and copy number variations between highand low-risk cohorts based on the 18-gene prognostic model, rather than exploring their relationships with UM prognosis independently. If more gene mutations and copy number variations exist in high-risk group consistent with the previous study, which show our prognostic gene model is more accurate. In our work, we first demonstrated an 18-gene signature could act as prognostic signatures for uveal melanoma, which could provide insights into new prognostic biomarkers exploration. Multivariate Cox regression analysis showed that the 18-gene risk score can be an independent prognostic factor.
Uveal melanoma displays gene mutations and copy number variations that correlate strongly with clinical outcome. Therefore, we investigated if there were gene mutation differences between high-and low-risk groups of UM patients. In the 80 uveal melanoma samples, 50% had GNAQ and 44% had GNA11. In our study, most of uveal melanoma harbored mutations in GNAQ, GNA11 as previously reported [9,29,30]. These gene mutations were not identified in the cutaneous melanoma, which was consistent with the literature [31]. Mutations in GNAQ (50%) and GNA11 (44%) are early events that promote cell proliferation, suggesting that activated G-protein signaling plays a crucial role in early uveal melanoma development [12]. In addition to mutated genes mentioned above, uveal melanoma also harbored additional SF3B1 (22%), EIF1AX (12%), COL14A1 (4%), CYSLTR2 (4%), MACF1 (4%), MUC16 (4%) and MYOF (4%). BAP1, SF3B1 and EIF1AX gene mutations in UM patients are related to poor, median and favorable prognosis, respectively [7,8,32]. While our study identified UM patients to have a favorable prognosis harboring GNAQ and SF3B1 mutations ( Figure 4F, 4H). BAP1 inactivation being related to poor prognosis in uveal melanoma is well studied [9,33,34]. However, In our Kaplan-Meier analysis, the difference in overall survival among patients with BAP1 mutation and patients without mutation was not statistically significant (P = 0.0875), indicating that BAP1 mutation was not significantly associated with OS ( Figure 4G). This discrepancy may be contributed to the small number of UM patients of the study. Except for exploring related gene mutations with UM overall survival, we also investigated chromosome copy number variants between high-and low-risk groups. Our result showed that UM patients with chromosome 3 loss and chromosome 8 gain should be considered in high-risk group. Consistently, Scholes et al. demonstrated that half uveal melanoma patients with loss copy of chromosome 3, combining with other chromosome variations and clinical information, such as 6p and 8q gain, which will largely improve prognostic accuracy [11,35,36]. Prior studies have also shown poorer clinical outcomes to be associated with higher chromosome 8q copy number [37][38][39].
Biological pathway profiling showed that tumorigenesis and progress related processes such as p53 signaling, Notch receptor signaling, regulation of cell adhesion mediated by integrin and other pathways were enriched in high-risk UM patients. The results implied that specific strategies targeted these biological pathways may achieve therapeutic efficacy in high-risk UM group of poor prognosis. Further clinical studies on these pathways are needed.
In summary, we identified an 18-gene prognostic model associated with UM overall survival, combining with determination of chromosome copy number variants and gene mutations can be a useful tool to predict UM patient's prognosis. However, several limitations to this study exist. First, our study is retrospective, and prospective study should be further validated. Second, we just constructed our model in one cohort owing to the small number of UM samples. Validating this 18-gene prognostic model in a larger cohort of uveal melanoma patients can make the prognosis more convincing.