3.1 Subgroup analysis of m6A regulators
As a result, a total of thirteen m6A RNA methylation regulators changed mRNA expression values and clinicopathological characteristics of UM were obtained from TCGA and GEO. Based on “ClassDiscovery” algorithm, 80 UM patients from TCGA and 28 UM patients from GEO can be identified two clusters of groups, respectively. (Figure 1A and B). Then, we contrasted the clinical features of these two subgroups, namely, C1 and C2. The subgroups analysis of clinical characteristics showed that Chromosome 3 status and subtype of UM have significant differences (Table 1). The others clinical characteristics like age, gender, TNM and stage have no statistical significance. To find out the potential correlation of overall survival with C1 and C2. Kaplan-Meier survival analysis was performed and the curves showed that overall survival of samples in C2 is longer than the samples in the C1 group (Figure 1C, D). Then, expression levels of thirteen m6A RNA methylation regulators in UM patients with different C1/2 groups were shown in Figure 1E, F.
3.2 Gene mutation and m6A regulators
Then, we assessed the relationship between gene mutation and m6A regulators. we firstly identified the top 20 mutated genes in TCGA database, which was calculated by ‘maftools’ R package (Figure 2A). The heatmap of m6A regulators expression and 20 highly variant mutated genes indicated that SF3B1, CYSLTR2 and ADAMTSL1 were the most significantly regulated the expression of m6A regulators. (Figure 2B). Kaplan-Meier analysis of these 3 mutant genes showed that only SF3B1 have a significant difference with overall survival. And then studied the relationship between SF3B1, CYSLTR2 and ADAMTSL1 mutation status and expression levels of each m6A RNA methylation regulator in TCGA database, respectively. The results showed that there are significant differences between with mutant-SF3B1 and wildtype-SF3B1 for the expression levels of ALKBH5, FTO, WTAP, YTHDF1, YTHDF2, YTHDC2 and KIAA1429 respectively (Figure 2C). Compared mutant-CYSLTR2 to wildtype-CYSLTR2, the expression levels of ALKBH5, FTO, METTL14, WTAP, YTHDF2, YTHDC2, ZC3H13, KIAA1429 and RBM15 are significantly different (Figure 2D). The subgroup analysis of mutant-ADAMTSL1 and wildtype-ADAMTSL1 also showed that m6A regulators of ALKBH5, METTL14, WTAP, YTHDF2, YTHDC1, YTHDC2, KIAA1429 and RBM15 are also significantly different (Figure 2E).
3.4 Clustered molecular subtype of uveal melanoma
The above results revealed that the clustered molecular subtype was intimately related to the prognosis of uveal melanoma. For better understanding of the interrelations among the thirteen m6A regulators, we also analyzed the interrelation (Figure. 3A) and correlation (Figure 3C) among these regulators. ALKBH5 seems to be the hub gene of the ‘Eraser’, and correlated or co-expressed with METTL3, WTAP, YTHDF2, M ETTL14, YTHDF1, YTHDC1, YTHDC2, RBM15, KIAA1429. The correlation analysis showed that these regulators were significantly positively correlated with each other. Principal components analysis showed that C1 samples and C2 samples in TCGA datasets could be well differentiated based on the expression of m6A regulators (Figure 3B). To investigate biologic pathways shared by the different C1/2 subtype, we performed GSEA analysis. According to the following criteria: p value<0.05 and | NES | ≥1. 49 BP terms were differentially enriched in C1 expression phenotype. The top 5 BP terms indicated that pathways are commonly enriched T cell mediated pathways, including positive regulation of T cell mediated cytotoxicity, antigen processing and presentation of endogenous antigen, regulation of T cell mediated cytotoxicity, positive regulation of T cell mediated immunity and regulation of T cell mediated immunity (Figure 3D). What’s more, The GSEA analysis of cancer malignant hallmarks of tumors showed that 9 terms including mTORC1 signaling, oxidative phosphorylation, interferon-a response and apoptosis signaling were significantly associated with the C1 subgroup expression phenotype. (Figure 3E).
3.5 Identification and confirmation of m6A regulators signature
For better predict the clinical and pathologic outcomes of UM with m6A regulators. Firstly, best survival analysis was used to evaluate associations between m6A regulators and OS in TCGA dataset. Totally, three m6A regulators were seeded out (Figure 4A). Then, we used the three m6A regulators to constructed risk system by multivariate cox regression analysis. The risk system reckons a risk score for each patient. The distributions of the risk scores, OS, vital status, and expression levels of corresponding 3 m6A regulators in TCGA dataset were shown in Figure 4B-D. Next, UM samples were divided into a high-risk group (n = 40) and a low-risk group (n = 40) by applying the median value of the risk scores. Kaplan-Meier curves revealed that low risk group have a significant longer survival time than high risk (Figure 4E). The ROC curve showed that the 5 years of AUC was 0.645 (Figure 4F). To verify the predictive ability of the three m6A regulators, validation analysis was performed in GEO dataset. The distributions of the risk scores, OS, vital status, and expression levels of corresponding 3 m6A regulators in GEO dataset were shown in Figure 4G-I. The curve of Kaplan-Meier revealed that there is a significant difference between high-risk and low-risk group with log-rank test of p=0.044 (Figure 4J). The 5 years of AUC was 0.677 (Figure 4K). The subgroups analysis of clinical characteristics between low- and high- risk groups showed that chromosome 3 status, subtype, and vital status in TCGA and GEO have significant differences (Table 2).
3.6 Associations between risk score and clinical variables
The associations between risk score of m6A regulators and clinical variables such as chromosome 3 status, mutated SF3B1, mutated BAP1 and subtype were explored. The box plots showed that monosomy 3 have higher risk scores than disomy 3 (Figure 5A), wildtype of SF3B1 own higher risk scores than mutant (Figure 5B), and subtype 4 of UM have the highest scores than other subtypes (Figure 5D). While, subgroup of mutated BAP1 manifested that there is on significant difference between mutant and wildtype (Figure 5C). To evaluate the associations between risk score and immune microenvironment, “MCPcounter” package in R was applied to calculate the immune scores of immune cells. Subgroup analysis of immune cells showed that significant differences were founded in T cells, CD8 T cells, cytotoxic lymphocytes, natural killer (NK) cells, monocytic lineage and myeloid dendritic cells between high and low risk subgroups (Figure 5E). Besides, univariate and multivariate logistic regression were used to compare the prognostic value of risk score and other clinical variables in TCGA and GEO datasets. The forest plots indicated that age, stage, histology, subtype, chromosome 3 status, metastasis and risk sore were significantly associated with OS in univariate analysis, but only the risk score were significantly correlated with OS in multivariate analysis (Figure 5F-G). The 5 years AUC of age, stage, histology, subtype, chromosome 3 status and risk score in TCGA were 0.591, 0.535, 0.351, 0.788, 0.791 and 0.645 respectively (Figure 5H). As for GEO, the 5 years AUC of chromosome 3 status and risk score were 0.698 and 0.677 (Figure 5I).