Clinical characteristics of all UVM patients
The TCGA UVM dataset (N=80) was utilized to construct the prognostic immune-gene signature. GSE84976 (N=28) and GSE22138 (N=63) datasets were used as the validation cohorts. All clinical characteristics are summarized in Table 2.
Identification of the prognostic six-immune-gene signature in TCGA UVM dataset
One hundred thirty overlapping prognostic immune genes were obtained from the univariate Cox regression analysis and K-M survival analysis (Fig 1A). Then, 14 prognostic immune genes were selected by LASSO Cox regression analysis (Fig 1B-C). Finally, a six-immune-related-gene signature was identified (Fig 1D), whose risk score of each patient was generated using the following risk score formula: Risk score = Exp_JAG2 × 0.157 + Exp_ CCL18 × 0.046+Exp_ PRELID1 × 0.284 + Exp_ CXCL8 × 1.154 + Exp_ PTGER4 × 0.092- Exp_ GTPBP1 × 0.576 (Table 3). In addition, the K-M survival results of these six immune genes were presented in Fig 1E-J.
Evaluation of the six-immune-gene signature in TCGA UVM dataset
The risk score, survival time and survival status of each patient are plotted in Fig 2A-B. The heat map is displayed in Fig 2C. The overall survival time of the low-risk group patients was significantly higher than that of the high-risk group patients (Fig 2D, p < 0.001). The AUCs of the time-dependent ROC at 1-, 2- and 3 year were 0.962, 0.943, and 0.962, respectively (Fig 2E). The AUC of multi-variate ROC for the 6-immune-gene signature risk score was 0.97, which was much better than that of other clinical variables (Fig 2F). The PCA analysis showed that this six-immune-gene signature could help distinguish the high-risk patients from the low-risk patients ideally (Fig 2G).
Moreover, the results of univariate and multivariate Cox regression analyses of the age, gender, histological type, TNM stage, new tumor event, tumor basal diameter and this signature risk score, are shown in Fig 2H-I. Risk score of this newly identified signature (HR = 2.117, 95%CI = 1.525-2.938, p < 0.001) was an independent clinical prognostic risk factor (Fig 2I).
Validation of the six-immune-gene signature in the two GEO datasets
To further confirm the predictive diagnostic power and stability of this six-immune-gene signature in predicting the overall survivals of UVM patients, we validated it in two GEO datasets, including GSE84976 (N=28) (Fig 3) and GSE22138 (N=63) (Fig 4). Risk scores were also generated using the same risk score formula constructed in the TCGA UVM dataset.
The risk score, survival time and survival status of each patient in the two GEO datasets were displayed (Fig 3A-B, 4A-B). The heat map of the expression of this six-immune-gene signature was plotted (Fig 3C, 4C). The overall survival time of the low-risk group patients was significantly higher than that of the high-risk group patients (Fig 3D, 4D, p < 0.001). In addition, in GSE84976, the AUCs of the time-dependent ROC at 1-, 2- and 3 year were 0.813, 0.859, and 0.782, respectively (Fig 2E). In GSE22138, the AUCs of the time-dependent ROC at 1-, 2- and 3 year were 0.551, 0.652, and 0.629, respectively (Fig 3E). The PCA analysis indicated that this six-immune-gene signature could largely help distinguish the high-risk patients from the low-risk patients (Fig 3F, 4F).
Clinical correlations in the TCGA UVM dataset
Consistent with our expectation, patients with single subtype, higher TNM stages (stage IV), higher T stage (T4), higher M stage (M1), new tumor event (YES), and tumor basal diameter (≥15), had significant higher risk scores than those with mixed subtype, lower TNM stages (stage II+ III), lower T stage (T2+3), M0, new tumor event (NO), and tumor basal diameter (<15) (all p < 0.05) (Fig 5C-G). The risk score in age (Fig 5A) and gender (Fig 5B) did not show significant differences. These findings suggested that high risk score of this signature might be involved with the disease progression of UVM patients.
Stratification survival analysis in the TCGA UVM dataset
Compared with the patients in the low-risk group, those in the high-risk group had a worse outcome in many different subgroups, including age (Fig 6A-B), gender (Fig 6C-D), new tumor event (NO) (Fig 6E), tumor basal diameter (Fig 6G-H), TNM stage (Fig 6I-J), T3+4 (Fig 6K), M0 (Fig 6L), N0 (Fig 6M), mixed subtype (Fig 6O), and spindle cell subtype (Fig 6P) (all p < 0.05), but not in the group with new tumor event (YES) (Fig 6F), and epithelioid cell subtype (Fig 6N) (all p > 0.05). However, they all showed the same tendency.
Functional enrichment analysis in the TCGA UVM dataset
The functional enrichment analysis of GO enrichment categories, including biological process (BP), cell component (CC) and molecular function (MF), displayed the enrichment of some known immune-related pathways, including response to interferon gamma, T cell activation, interferon-gamma-mediated signaling pathway, cellular response to interferon-gamma, antigen processing and presentation of peptide antigen, and so on (Fig 7A-B). Similar result was also obtained from KEGG pathway enrichment analysis (Fig 7C-D).
ssGSEA analysis in the TCGA UVM dataset
The ssGSEA indicated that the high risk patients were enriched of many immune cells, including B cells, CD8+ T cells, DCs, macrophages, pDCs, Tfh, Th2 cells,TIL, and Treg, while the low risk patients were only enriched in aDCs (Fig 8A, all p < 0.05) . High risk patients are enriched in all immune functions (all p <0.05), except for APC-co-inhibition and Type-II IFN response (Fig 8B, all p > 0.05).
Construction and evaluation of the predictive nomogram model in the TCGA UVM dataset
This newly identified predictive nomogram model was successfully constructed by combining all the details of age, gender, histological type, TNM stage, new tumor event, tumor basal diameter and this signature risk score in TCGA UVM dataset (Fig 9A). The calibration plots suggested that no significant deviations between the observed and predicted curves were found for both 1-year and 3-year survivals (Fig 9B-C). The UVM patients in high-nomogram-score group had a worse outcome compared with those in the low nomogram-score group (p < 0.001) (Fig 9D). The AUCs of time-dependent ROC curves for 1-, 2- and 3 years were 0.977, 0.980, and 0.968, respectively (Fig 9E).
Knocking-down of CCL18 expression inhibits uveal melanoma cells proliferation, migration and invasion.
To verify our computational findings, we tested the mRNA expression level of CCL18, CXCL8, GTPBP1, JAG2, PRELID1 and PTGER4 in uveal melanoma cell lines: M17, M23 and SP6.5 and normal uveal epithelial cell: Um95. We found that compared with Um95, the expression patterns of these six genes in uveal melanoma cell lines: are consistent with our computational findings (Fig 10A). In addition, we found that CCL18 has the highest expression in the M17 cell line, so we selected the CCL18 gene to perform related functional verification in M17 cells. As shown in Figure 10B&C: siRNA-2 successfully knocked down the expression of CCL18 in M17 cells. The results of cell proliferation experiments show that low CCL18 expression can significantly inhibit the proliferation ability of M17 cells (Fig 10D). Inhibition of CCL18 expression will significantly inhibit the migration and invasion of M17 cells（Fig 11 A&B）。Furthermore, knocking down the expression of CCL18 significantly induced M17 cell apoptosis (Fig 11 C). In summary: Our results demonstrated that our computational findings were consistent with the trend of in vitro cell experiments.