The RNA-sequencing transcriptome data of 80 patients with UM and corresponding clinical data were downloaded from The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/). GTF files were downloaded from Ensembl (https://asia.ensembl.org) to distinguish between lncRNAs and mRNAs for subsequent analyses. The list of immune genes was downloaded from the ImmPort database (https://www.immport.org).

The ESTIMATE algorithm in the R “estimate” package was used to calculate the TME scores of 80 UM patients. The Kaplan-Meier survival analysis was conducted to compare the difference in survival time using R “survMiner” and “survival” packages. The fraction of 22 immune cell types in each sample was estimated using CIBERSORT. The association between TME scores and tumor-infiltrating immune cells (TICs) was established using correlation analysis.

The m⁶A regulators were identified from the published literature (9–12). The m⁶A regulators contain 8 writers (METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15, and RBM15B), 13 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGF2BP1, IGF2BP2, IGF2BP3, and RBMX), and 2 erasers (FTO and ALKBH5). The expression data of m⁶A regulators were extracted from the mRNA expression data of TCGA. The m⁶A regulators, which were previously identified using the univariate Cox regression analysis, were further subjected to the LASSO Cox regression analysis using the “glmnet” package. The minimum 10-fold cross-validation was used to select the best penalty parameter λ. Then, the risk score of each patient was calculated using a linear combination of m⁶A regulators expression weighted by the multivariate Cox regression analysis. According to the median risk score, the samples were divided into high-risk and low-risk groups. Subsequently, the Kaplan-Meier survival analysis was performed to compare the survival difference. Receiver operating characteristic (ROC) analysis was performed to evaluate the prognostic value of the signature using the “timeROC” package.

We screened m⁶A regulators-related lncRNAs by Pearson’s correlation analysis. The process used the criteria of |correlation coefficient| > 0.4 and p < 0.001. The expression data of m⁶A regulators-related lncRNAs were extracted from the lncRNA expression data of TCGA. To clarify the biological characteristics of m⁶A regulators-related lncRNAs, the R “ConsensusClusterPlus” package was used to divide the samples into different clustering subtypes according to the expression of lncRNAs. Principal component analysis (PCA), Kaplan-Meier survival analysis, TME scores, and TICs profiles were performed for different clustering subtypes.

The Pearson’s correlation analysis was performed to screen common immune genes shared by 5 lncRNAs. The Kaplan-Meier survival analysis was used to compare the survival difference between high- and low-expression of immune genes. The differential expression analysis of immune genes was performed using the R “limma” package in different risk groups and clustering subtypes. A p-value < 0.05 indicated statistical significance. At present, single-cell sequencing technology has been widely used to explore the heterogeneity of TME. To characterize immune genes expression distribution in TME of UM, we search for single-cell sequencing data of UM from Tumor Immune Single-cell Hub (TISCH) (22). The GSE139829 dataset including 59,915 cell sequencing data from 11 samples was collected to perform gene expression distribution (23).

To establish the correlation of ImmuneScore and StromalScore with the survival time, we performed the Kaplan-Meier survival analysis. A high score of immune and stromal cells signified large numbers of these cells in the TME. As shown in Figure 1A , the overall survival (OS) in the low ImmuneScore group was longer than that in the high ImmuneScore group. Similarly, StromalScore and ESTIMATEScore showed a negative correlation with the OS ( Figures 1B, C ). To confirm exact changes in the genetic profiles in the TME about immune and stromal cell components, variance analysis of high and low scores was performed. As shown in Figures 1D, E , 700 common differentially expressed genes (DEGs) were upregulated, and 74 common DEGs were downregulated. To further explore the interaction between TME scores and the 22 immune cell types, we first estimated the 22 types of TICs with abundance distribution in all the tumor samples and then calculated the correlation index between the TME scores and TICs. The results showed that ImmuneScore is correlated with CD8⁺ T cells, resting memory CD4⁺ T cells, activated memory CD4⁺ T cells, helper T cells (follicular), Tregs, resting NK cells, monocytes, M0 macrophages, M1 macrophages, resting mast cells, and eosinophils ( Figure 1F ).

To clarify the biological function of m⁶A regulators in the prognosis of patients with UM, we comprehensively investigated the prognostic value of m⁶A regulators based on the expression and clinical data ( Table S1 ). Of these m⁶A regulators, seven exhibited a prognostic value based on OS ( Figure 2A ), while nine displayed a prognostic value based on progression-free survival (PFS) ( Figure 2B ). The Venn plot results indicated that 6 m⁶A regulators (RBM15B, IGF2BP2, YTHDF1, METTL16, VIRMA, and YTHDF3) were identified based on OS and PFS ( Figure 2C ). To avoid overfitting, we performed the LASSO Cox analysis and selected 3 of the 6 m⁶A regulators were to establish a risk signature ( Figures 2D, E ). Thus, we established a predictive model: risk score = (RBM15B * −0.14284) + (YTHDF3 * 0.02121) + (IGF2BP2 * −0.11533). The distribution of risk score (Upper), patients’ survival time (Middle), and heat map analysis (Bottom) of the 3 prognostic m⁶A regulators were shown based on the OS ( Figures 2F ) and PFS ( Figure 2G ). Results of the heat maps ( Figures 2F, G ) was survival curves ( Figures 3A–C ) suggested that YTHDF3 was likely to be a high-risk factor because it was upregulated in the high-risk group. However, the highly expressed RBM15B and IGF2BP2 in the low-risk group might be protective factors. As shown in Figures 3D, E , the OS and PFS of patients in the low-risk group were longer than in the high-risk group. To evaluate the prognostic accuracy of the 3 m⁶A regulators-based signature, we performed the ROC analysis based on OS and PFS. Areas under the ROC curves of 1, 2, and 3 years were 0.774, 0.811, and 0.843, respectively ( Figure 4A ). The PFS prediction of 3 m⁶A regulators-based signature was also accurate ( Figure 4B ). These results suggested that the 3 m⁶A regulators-based signature might serve as a promising parameter for prognostic prediction of UM.

The correlation between the risk score and the immune cell infiltration levels was calculated to establish the association between the 3 m⁶A regulators-based signature with the TME. ImmuneScore, StromalScore, and ESTIMATEScore in the high-risk group were significantly higher than in the low-risk group ( Figures 5A–C ). The results showed that the risk score was significantly negatively correlated with the infiltration levels of naive B cells ( Figure 5D ), eosinophils ( Figure 5F ), monocytes ( Figure 5G ), and plasma cells ( Figure 5H ). Only the infiltration level of resting dendritic cells was positively correlated with the risk score ( Figure 5E ). These results indicated that the 3 m⁶A regulators were involved in the immune cell infiltration of UM.

Using the correlation analysis, we identified 514 lncRNAs based on RNA-Seq data and constructed a network between m⁶A regulators and lncRNAs ( Figure 6A ). A total of 66 of the 514 lncRNAs had a prognostic value based on the OS ( Figure S1A ), while 70 of the 514 lncRNAs had a prognostic value based on the PFS ( Figure S1B ). The Venn plot results showed that 38 lncRNAs were common according to OS and PFS ( Figure 6B ). Based on the similarity identified by consensus clustering using the ‘ConsensusClusterPlus’ package, we found that k = 2 was the optimal clustering stability value ( Figures 6C, D ). The 80 UM samples were well-differentiated into two subtypes according to the expression of 38 lncRNAs ( Figure 6E ).

The TME scores of cluster 1 were lower than those of cluster 2 ( Figures 7A–C ), while the OS ( Figure 7D ) and PFS ( Figure 7E ) of cluster 1 were notably longer than those of cluster 2. Subsequently, the 22 immune cell levels for the two subtypes were calculated. The results showed that cluster 1 had higher immune infiltration levels of plasma cells, and monocytes while there were higher immune cell infiltration levels of activated memory CD4⁺ T cells, follicular helper T cells, M1 macrophages, and resting dendritic cells ( Figure 7F ).

First, we identified the common lncRNAs from 5 lists that were DEGs in ImmuneScore (Immune-DEGs) and StromalScore (Stromal-DEGs), m⁶A regulators-related lncRNAs (m⁶A-lncRNAs), and immune gene-specific lncRNAs based on OS (OS-immune) and PFS (PFS-immune). A total of 5 lncRNAs, namely AC008555.4, AC018529.1, AC104129.1, CYTOR, and MIR4435-2HG, were found to be common across the 5 lists ( Figure 8A ). High expression of AC018529.1 ( Figure 8B ), MIR4435-2HG ( Figure 8C ), AC104129.1 ( Figure 8D ), and CYTOR ( Figure 8F ) was related to the poor prognosis of patients, while high expression of AC008555.4 was associated with good prognosis ( Figure 8E ). Next, 13 immune genes having 5 common lncRNAs were screened using the correlation analysis ( Figure 9A ). The 13 immune genes were ADGRE5, C2, CD79B, CTSC, GEM, JAG2, LYN, MAFB, MBP, MR1, PREX1, RUNX1, and TCF12. Except for C2, the other 12 immune genes were associated with a poor prognosis ( Figures 9B–M ). To further verify whether the immune genes were differentially expressed in different groups based on m⁶A regulators and lncRNAs, we extracted the expression data from RNA-Seq data and performed the differential expression analysis. The expression levels of 13 immune genes in the high-risk group were upregulated than those in the low-risk group ( Figure 10 ). Similarly, the 13 immune genes levels in cluster 2 were upregulated compared with cluster 1 ( Figure 11 ). For further validation, we used TISCH to depict the expression distribution of 13 immune genes in the TME of UM. The overall distribution of the 3 cell types in the GSE139829 dataset was shown in Figure 12A . Through analysis, we found that the expression distribution of ADGRE5 ( Figure 12B ), CTSC ( Figure 12E ), LYN ( Figure 12H ), MAFB ( Figure 12I ), PREX1 ( Figure 12L ), and RUNX1 ( Figure 12M ) were abundant in immune cells. CD79B ( Figure 12D ), GEM ( Figure 12F ), and JAG2 ( Figure 12G ) expression distribution were concentrated in malignant cells. The gene expression of C2 ( Figure 12C ), MBP ( Figure 12J ), MR1 ( Figure 12K ), and TCF12 ( Figure 12N ) were evenly distributed in immune cells and malignant cells. These results suggested that these immune genes may be the downstream regulators of m⁶A regulators and related lncRNAs participated in TME remodeling.