Survival associated AIRGs
Totally, 222 autophagy related genes and 1793 immunity related genes were selected for univariate cox analysis in TCGA-UVM dataset. According to the selection standard, 44 autophagy related genes were significantly associated with survival of UM patients, which contained 27 risk genes and 17 protect genes (Figure 2A). In addition, 365 immunity related genes are survival associated, which consist of 301 risk genes and 64 protect genes (Figure 2B). Apart from the overlapped genes, a total of 409 AIRGs was identified for subsequent analysis.
Identification and validation of prognostic AIRGs biomarker
Firstly, 1,000 iterations of LASSO modeling were performed to evaluate the qualified variables from the 409 prognostic AIRGs, and then these qualified AIRGs were used to construct prognostic biomarker. Totally, 12 qualified prognostic AIRGs were selected to subsequently build biomarker due to appear more than 990 times in LASSO modeling (Figure 2C). Moreover, based on the 5-year AUCs of different gene combinations, we found that a 4-AIRGs related biomarker (PRKCD, MPL, EREG and JAG2) was arrived the max value of 0.829 (Figure 2D,E). These four AIRGs were closely correlated with each other (Figure 2F). Finally, the 4-AIRGs biomarker was used to calculate risk score and the formula listed as follows: the risk score = . By applying this risk formula, a risk score for each sample will be generated and scaled to range from 0 to 1. Next, by applying the median cut-off value of the risk scores, UM patients were categorized into high-risk group (n = 40) and low-risk group (n = 40). Kaplan-Meier (KM) curves showed that patients in high-risk group have a shorter survival time than low-risk with a log-rank test p<0.001 (Figure 3A). The ROC curve was drawn and 3,5 years of AUCs were 0.916 and 0.829 respectively (Figure 3B). The distribution of risk score, OS, vital status, and expression of corresponding AIRGs in TCGA-UVM were shown in Figure 3C-E. In order to verify the robustness of the result, validated analyses were applied in another two independent datasets including GSE22138 and GSE84976. Patients in GSE22138 and GSE84976 datasets were accordingly classified into high-risk and low-risk groups. KM analyses manifested a similar result that patients in low risk group had significantly longer survival time than those in high-risk group (GSE22138 log-rank p=0.0019 and GSE84976 log-rank p<0.001 respectively) (Figure 4A and Figure 5A). The 3,5 years of AUCs in GSE22138 were 0.746, 0.713 (Figure 4B). Higher AUCs (3-year:0.836; 5-year:0.872) also observed in GSE84976 (Figure 5B). The distribution of risk score, MFS, vital status, and expression of 4 AIRGs in GSE22138 was shown in Figure 4C-E. Meanwhile, the same risk score heatmap of 4 AIRGs in GSE84976 was shown in Figure 5C-E.
Prognostic value of biomarker
To explore prognostic factors for OS or MFS in multiple datasets, the risk score of AIRGs biomarker and clinical variables were conducted by the univariate and multivariate cox regression analyses (Table 1). The results of univariate cox regression revealed that stage, age, metastasis, histological type and risk sore were significantly associated with MFS or OS. According to the multivariate cox regression analyses, only the risk score was stably remained independent prediction for OS or MFS in TCGA-VUM dataset (HR= 34.951, 95% CI = 4.891–249.759, P < 0.001), GSE22138 dataset (HR = 45.623, 95%CI = 10.025–564.298, P < 0.001) and GSE84976 dataset (HR = 3.532, 95%CI = 1.783–9.276, P =0.042). The 5 years AUCs of age, stage, metastasis, histological type and risk sore in TCGA-VUM set were 0.741, 0.778, 0.658, 0.424 and 0.829 respectively (Figure 6A). In GSE22138, the 5 years AUCs of Chromosome 3 status and risk score were 0.762 and 0.713 (Figure 6B). In addition, the 5 years AUCs of Chromosome 3 status, metastasis and risk sore in GSE84976 set were 0.912, 0.858, and 0.872 respectively (Figure 6C).
Associations between risk model with clinical characteristics and immune infiltration
The relationships between the risk score signature and various clinical characteristics were investigated in TCGA-VUM dataset. The box plot manifested that UM patients in stage IV, metastasis, epithelioid histological type and dead status have higher risk score levels than other subtypes (Figure 7A). Other clinical characteristics (age and gender) have no relationships with risk score. Thus, the risk model was intimately associated with tumor stage, metastasis, histological type and vital status. Through CIBERSORT analysis, the proportion of 22 different infiltrating immune cells was estimated in TCGA-VUM dataset. The stratifying analysis showed that the fractions of T cell CD8, T cells CD4, T cells gamma delta, B cells naïve, NK cells activated, Monocytes, Macrophages M1 and Neutrophils in high risk subgroup significantly different from those in low risk group (Figure 7B).
Gene set enrichment analysis (GSEA)
In order to explore the different hallmark pathways enriched in high and low risk group, GSEA analyses were performed. According to selected criterion, we confirmed that six positive correlated pathways were enriched in high risk group, which including epithelial mesenchymal transition, estrogen response late, myogenesis, estrogen response early, apical junction and KRAS signaling. moreover, six cancer hallmark pathways contained interferon gamma response, protein secretion, interferon alpha response, MTORC1 signaling, JAK STAT3 signaling and inflammatory response were actively associated with low risk group (Figure 7C).
High-risk subgroup more sensitive to chemotherapies
At present, chemotherapy is a common treatment for UM, therefore GDSC database was also applied to predict sensitivity to chemo drugs. Assessment by selecting criteria, 45 kinds of drugs were screened out in TCGA-UVM set. 34 and 14 kinds of drugs were sensitive to high risk group in GSE22138 and GSE84976 respectively (Table 2). The Venn plot of three datasets suggested that AMG.706 and JNK.Inhibitor.VIII were the common drugs (Figure 8A). We could observe a significant difference in the estimated IC50 between high and low risk group for these chemo drugs where high-risk group could be more sensitive to chemotherapies (Figure 8B-D).