2.1. Data Collection

The transcriptome profiling data for Uveal Melanoma (UM) was obtained from The Cancer Genome Atlas (TCGA-UVM project), comprising 80 UM samples. To ensure accurate comparison, transcriptome data from 187 normal tissue samples were also obtained from the Genotype-Tissue Expression (GTEx) dataset, accessible through the UCSC Xena browser (https://xenabrowser.net/, accessed on 5 March 2023).

Given the difference in the sources of data (TCGA and GTEx), batch effects were a concern, as they could obscure real biological differences between the tumor and control samples. Therefore, the Fragments Per Kilobase of exon per Million mapped fragments (FPKM) normalization method was employed. FPKM corrects for differences in sequencing depth and gene length, providing normalized expression values that are comparable between samples.

In addition to this, “limma” package in R was used to correct for batch effects using an empirical Bayes framework. This method models technical noise (such as sequencing artifacts or platform-specific effects) while preserving biological variability. Finally, to further standardize the data for analysis and comparison, log-transformation was implemented to all gene expression values to reduce the effect of outliers and make the distributions of gene expression across samples more comparable. This combined approach ensured that differences in gene expression between UM and normal tissues were primarily biological and not artifacts introduced by the data collection processes.

2.2. Bioinformatic Analysis

The “survival” package in R (version 4.2.2) was used to calculate the hazard ratio of RNAs, which can identify significant genes related to prognosis. This package was also employed in estimating the relationship between survival time and different groups, as well as survival time and significant RNAs in the tumor tissue. Moreover, “ConsensusClusterPlus” package was employed with 1000 iterations and a resampling rate of 80% to categorize these remarkable RNAs into distinct subtypes [11]. Enrichment analyses were conducted for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)using the “clusterProfiler” package [12]. Additionally, (GSEA) version 4.1.0 was utilized to discover pathways in which genes participate [13].

The ‘estimate’ package was utilized to compute the immunoScore using the ESTIMATE method (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) algorithm [14]. Moreover, cell type identification was conducted to estimate the TIME infiltration levels (CIBERSORT, https://cibersort.stanford.edu/) [15].

2.3. Statistical Analysis

Statistical analysis was performed using R version 4.2.2. The t-test was initially applied to identify the significant RNAs between normal and tumor samples. Following that, Pearson correlation test was employed to examine the association between the various subtypes and clinical features, such as stage and age. Subsequently, the LASSO (Least Absolute Shrinkage and Selection Operator) regulation was utilized to construct the risk model: risk score = sum of coefficients ×expression levels of RNAs [16]. The Survival curves were generated via the Kaplan–Meier method. Additionally, the receiver operating characteristic (ROC) curves was employed to reflect the model’s predictive ability. And finally, the univariable cox regression was utilized to identify the significant RNAs associated with prognosis [17].

3.1. Expression Levels of RNAs in UM and Normal Tissue

The detailed workflow of the present study is shown in Supplementary Figure S1. After gene expression analysis, 32 genes with significant expression levels between normal and tumor tissue were found (p < 0.01, Figure 1A). Of these, the expression levels of 20 RNAs were found to be higher in uveal melanoma tissue as compared to their normal counterparts, including HLA-DRB1, PHLDA2, PROSER2, TOLLIP-AS1, FLYWCH2, OGDHL, TUBB4A, AC004925.1, RMND5A, OTUD1, SCFD2, RBM15B, WWC2-AS2, ZMAT3, AC112907.2, PTPRH, MFSD3, AC011481.1, KCNE4, MT1A. The remaining genes, like UPP1, RAB4B, CAMK1, IL32, FOXP1, PHYKPL, HSPE1-MOB4, ACOX2, APIP, NUCB2, MITD1 and SPCS2, were lower in uveal melanoma than normal tissue.

3.2. Identification of Significant RNAs Related to Prognosis

Univariate Cox regression revealed that 25 RNAs were found to be significantly correlated with the overall survival (OS) of individuals with uveal melanoma (p < 0.05, Supplementary Figure S2). In addition, the forest plot indicated that 15 genes showed low hazard ratios, implying that patients with higher expression levels of these genes had shorter OS (p < 0.001, Figure 1B). For example, FLYWCH2, TUBB4A and RBM15B all had lower expression levels, resulting in longer OS (p < 0.001, Figure 1B). On the other hand, the remaining RNAs displayed high hazard ratios, with higher expression levels leading to shorter OS, such as CAMK1, PHLDA2, and HLA-DRB1 (p < 0.001, Figure 1B). Among all the high-risk genes, HSPE1-MOB4 was found to have the highest hazard ratio, suggesting a possible crucial role in survival rate of uveal melanoma patients. Furthermore, there were some genes with significantly positive correlations, especially between TUBB4A and RBM15B, as well as HLA-DRB1 and IL32. However, numerous genes exhibited negative correlations, such as PD-L1 and RMND5A, and RBM15B and NUCB2 (Figure 1C).

3.3. Relationship between the Survival Time and RNAs Subtype

Similarity analysis of RNAs expression levels, as well as ambiguous clustering measure percentages, indicated that k = 2 was the optimal clustering stability (Figure 2A and Supplementary Figure S3). This resulted in two different subtypes, namely cluster 1 (n = 53) and cluster 2 (n = 27). Interestingly, the survival rate of the former was found to be significantly greater than that of the latter (p < 0.001, Figure 2B). Following a comparison of clinical characteristics such as age and stage, notable differences were observed between the two clusters. To begin with, the expression levels of genes including ZMAT3, OGDHL, FOXP1, TUBB4A, RBM15B, AC112907.2, FLYWCH2, TOLLIP−AS1, WWC2−AS2, APIP, RMND5A, OTUD1, ACOX2, AC004925.1 and KCNE4 were significantly lower in cluster 2 compared to cluster 1. Conversely, the remaining genes exhibited a contrasting result. Furthermore, the T stage values varied remarkably between both clusters (Figure 2C).

3.4. Enrichment Analysis about Prognostic RNAs

Subsequently, the prognostic RNAs of two clusters were further evaluated for enrichment analysis. Results of the Gene Ontology (GO) analysis showed that RNAs were preferentially implicated in Biological Process (BP) and Molecular Function (MF) (Figure 3A,B). Enriched pathways related to immune regulation, such as the regulation of mRNA splicing via spliceosome, positive regulation of B cell receptor signaling, and negative regulation of lymphocyte proliferation, suggest an immune-active microenvironment. These processes contribute to a less aggressive tumor profile with better prognosis (Figure 3C). Key genes involved in cluster 1 UVM include FOXP1 and TUBB4A, which are linked to enhanced immune surveillance and anti-tumor immune responses. The cluster 2 phenotype displayed enrichment in pathways such as mRNA processing, regulation of monocyte differentiation, and dendritic cell antigen processing and presentation, suggesting a more immune-suppressive tumor microenvironment. Cellular components like the intermediate filament and endoplasmic reticulum-plasma membrane contact site are more prominent in cluster 2 UVM, highlighting a more aggressive and metastatic potential, aligned with poor prognosis. Genes such as HSPE1, CAMK1, and HLA-DRB1 contribute to the immune evasion mechanisms and enhance the tumor’s ability to progress and resist immune responses. Additionally, the GSEA analysis highlighted that in cluster 2, most genes were involved in pathways such as mTOR (mammalian target of rapamycin) (normalized enrichment score = 1.647, normalized p-value < 0.001) and P53 pathway (normalized enrichment score = 1.688, normalized p-value = 0.002) (Figure 3D,E).

3.5. Construction and Prediction of Risk Model of UM

Using LASSO regression, five RNAs were identified to build a risk model for predicting the overall survival (OS) of patients with uveal melanoma, along with the expression levels of significant RNAs in the UM training cohort (Supplementary Figure S4). This risk model was expressed as follows:

Based on the median overall risk score, patients with UM were categorized into two groups: high-risk and low-risk groups. The survival rate was significantly greater in low-risk group as compared to high-risk group in both the training and test cohorts (p < 0.001, Figure 4A). In test group, we obtained the similar result (Figure 4B).

Additionally, ROC curve analysis was performed to assess the 5 risk features. In training cohort, the Area Under the Curve (AUC) value was 0.92. On the other hand, the AUC value was 0.71 in the test set (Figure 4C,D). This finding supports the effectiveness of the model in accurately predicting the survival rate of patients based on the transcriptome profiling data.

The expression of HSPE1-MOB4 and CAMK1 were observed to be significantly increased in the high-risk group, while the expression of FOXP1, WWC2-AS2, and ACOX2 were significantly decreased in low-risk group. This was determined by using the risk score to assess the expressions of each gene in both high and low-risk groups (as shown in Figure 4E,F).

3.6. Correction between Subtypes and Level of Immune Cell Infiltration in the TIME

The relationships between subtypes and levels of immune cell infiltration in TIME were studied and it was found that there were six significant immune cells in both clusters 1 and 2, including T cells CD8, T cells CD4 memory activated, Monocytes, Macrophages M1, Dendritic cells resting, and Mast cells resting (p < 0.05, Supplementary Figure S5). However, only the infiltration of Monocytes was higher in cluster 1 than cluster 2 (p < 0.001). In contrast, the remaining five cell types showed lower levels in cluster 1 than cluster 2 (Figure 5A). Moreover, some immune markers such as ImmuneScore of cell fraction, ESTIMATEScore and StromalScore were significantly increased in cluster 2 (p < 0.001, Figure 5B–D).

3.7. Relation between Clinical Characteries and Risk Score

Moreover, when performing the survival analysis, differences between high and low-risk group were observed in clinical features such as age (p < 0.001, Figure 6A,B), gender (p < 0.002, Figure 6C,D), M stage (p < 0.001, Figure 6E), N stage (p < 0.001, Figure 6F) and T stage (p < 0.001, Figure 6G). On top of that, Figure 6H illustrated that ImmunoScore and cluster had remarkable dissimilarities in both high and low-risk groups (Supplementary Figure S6).

3.8 Risk Score Associated with the Immune Cell Infiltration in the TIME

Furthermore, the relationship among the risk score and six types of immune cells in the tumor immune microenvironment (TIME) of UM was also investigated. The analysis revealed a positive correlation between the risk score and some immune cells, such as dendritic cells resting, macrophages M1, T cells CD4 memory activated, and T cells gamma delta (p < 0.05, Figure 7A–D). Conversely, some immune cells fractions had a negative relation to risk score, such as B cells naïve, Monocytes, NK cells activated, Plasma cells (p < 0.05, Figure 7E–H).