The bulk RNA-seq of uveal melanoma (UM) data as well as clinical information were extracted from open access resources (TCGA and GEO database). The uveal melanoma dataset (TCGA-UVM) was acquired from UCSC Xena (http://xena.ucsc.edu). The UM-associated transcriptomes (GSE84976 and GSE22138) were obtained from the GEO website (https://www.ncbi.nlm.nih.gov/geo). Then, fragments per kilobase per million (FPKM) format of gene expression in TCGA-UVM was converted into transcripts per kilobase million (TPMs). The raw count of GEO datasets was also transformed to the TPMs format. Next, these datasets were merged into a large metadata by the “ComBat” algorithm, which can decrease batch effects caused by different platforms and non-biological technical biases. A simplified workflow for the current investigation is depicted in Figure 1 .

Based on the 1,000 permutations of LM22 signature, the “CIBERSORT” method was used to calculate the percentage of 22 different types of immune cells in the UM TME. Next, the “ConsensusClusterPlus” program was applied to stratify the UM patients into three separate subgroups in order to show the biological importance of different immune cell infiltrations. This program employed the unsupervised clustering “Pam” approach and Euclidean and Ward’s linkage, which was performed 1,000 times to confirm classification consistency.

UM patients in metadata were categorized into subgroups based on infiltrated immune cells. We used the “limma” method to screen immune cell-related DEGs and set the significant cutoff criteria to p < 0.05 (Adjust) and |log2FC| ≥1.5. To further investigate the underlying biological mechanism of DEGs, GO and KEGG annotation were conducted and the enriched items with adjust p < 0.05 were considered as significant.

To begin with, the expressed values of differently expressed genes were used to divide the UM patients by using unsupervised clustering. Genes with positive and negative correlations to genomic clusters were respectively defined as TIC gene features α and β. Next, the important gene sets in features α and β were identified by using the Boruta algorithm, which subsequently employed PCA calculation. Finally, principal component 1 was respectively extracted to represent the feature score. TIC scores in UM samples were calculated by the following formula: ΣPC1 (α) - ΣPC1 (β).

UM patients in the TCGA-UVM cohort had their mutation data deposited in the UCSC Xena website. We downloaded the Mutation Annotation Format file and analyzed it using the “Maftools” package. OncoPrint plots were used to display the mutation landscape of high- and low-TIC score subgroups. The different alteration frequencies between high- and low-TIC score subgroups were explored. Besides, Maftools was used to compute tumor mutational burden (TMB) in the TCGA-UVM cohort; the relationship between TMB and TIC score was then investigated.

To test the prognostic value and immune checkpoint therapy response of TIC scores, three independent melanoma cohorts, namely, TCGA-SKCM, GSE35640, and GSE78220, were downloaded and analyzed. The transcriptional profile of TCGA-SKCM and related information were acquired from the UCSC Xena website. The gene expression profiles of GSE35640 and GSE78220 as well as related information were obtained from the GEO database. In all, a number of 169 melanoma patients who accepted immune checkpoint therapy were analyzed to measure the TIC score.

To determine the possible therapeutic agents, three different drug response databases, namely, CTRP (https://portals.broadinstitute.org/ctrp.v2.1), PRISM (https://www.theprismlab.org/), and GDSC (https://www.cancerrxgene.org/), were used to investigate the associations between drug response and TIC scores. Firstly, Spearman correlation between TIC scores and AUC values was used to identify potential therapeutic agents (CTRP: r < −0.20; PRISM: r < −0.20; GDSC: r < −0.40). Next, differential drug response between the high-TIC score group (upper decile) and the low-TIC score group (lower decile) was determined to identify compounds with a higher AUC in the low-TIC score group.

All statistical tests were conducted using R software version 3.6.0 with packages. The “ConsensusClusterPlus” package was used to conduct clustering. The “Boruta” package was performed to Boruta algorithm. The “CIBERSORT” package was used to estimate the immune cells. The “Survival” package made it possible to perform Kaplan–Meier (KM) survival analysis. The correlation analyses were assessed by Spearman test. Subgroup analyses was conducted by Kruskal−Wallis or Wilcoxon test. The chi-square test was used to analyze the different mutations across groups. p < 0.05 or Adjust p < 0.05 was considered to be statistically significant.

A total of 171 uveal melanoma (UM) samples were pooled into a big cohort that was studied in this research. We firstly removed the batch effect caused by the various platforms via the “ComBat” method. Before the removal of batch effect, the clusters of different platforms were clustered more closer than after removal ( Figure 2A ). Next, the CIBERSORT method was performed to determine the percentage of 22 immune cells in the immune microenvironment of UM. The UM patients were further categorized into subgroups based on the similar proportions of infiltrated immune cells. A steadily increasing trend in the cumulative distribution function (CDF) value was seen as a sign of stable clustering. We finally discovered three stable subtypes by applying an unsupervised clustering (k = 3), which contained Cluster 1 (51 UMs), Cluster 2 (66 UMs), and Cluster 3 (54 UMs). The relationship between subtypes and clinical characteristics was investigated and displayed in a complete heatmap ( Figure 2B ). The chi-square test showed significant variations in vital status, metastasis, chromosome 3 status, stage, and histological type among subtypes. Moreover, KM curves revealed that UM patients in Cluster 1 have a worst overall survival rate among subtypes with log rank test p < 0.0001 ( Figure 2C ). We further analyzed the immune cell composition of the UM microenvironment to better understand the fundamental biological distinctions that led to different clinical presentations. The association heatmap visualized the comprehensive interaction of immune cells in the UM microenvironment ( Figure 2D ). Among these subtypes, Cluster 1 had a high number of plasma cells, CD8 T cells, activated CD4 T cells, follicular helper T cells, gamma delta T cells, M1 macrophages, and resting dendritic cells, while Cluster 2 and Cluster 3 had a high number of naïve B cells, memory B cells, resting CD4 memory T cells, activated NK cells, monocytes, M0/M2 macrophages, activated dendritic cells, resting mast cells, eosinophils, and neutrophils ( Figure 2E ). In addition, the expressions of three essential immune checkpoint molecules (PD-1, PD-L1, and CTLA-4) were investigated in each subgroup. The Kruskal–Wallis test indicated that Cluster 1 had a higher expression of PD-1, PD-L1, and CTLA-4 than Cluster 2 and Cluster 3 ( Figure 2F ).

We used “limma” differential analysis to evaluate the DEGs among these subgroups to uncover the underlying biological properties of diverse immune phenotypes. A total of 1,054 DEGs were identified for unsupervised clustering analysis. UM samples were subsequently divided into three immune gene-related clusters (A–C). Cluster B was linked to a better prognosis with a log rank test p = 0.04 ( Figure 3B ). The feature α gene set was made up of 108 DEGs that were positively associated with the gene cluster, whereas feature β was made up of 946 DEGs that were negatively associated with the gene cluster. The relationship between the immune gene-related cluster and clinical characteristics was depicted in a heatmap of DEG expression ( Figure 3A ). The significantly enriched GO terms in feature α and feature β gene sets are illustrated in Figures 3C, D . Besides, the expression levels of PD-1, PD-L1, and CTLA-4 were investigated in each cluster. The Kruskal–Wallis test indicated that the expressions of PD-1, PD-L1, and CTLA-4 were generally higher in Cluster A than those in Cluster B and Cluster C ( Figure 3E ). Among the three gene clusters, Cluster A had a high infiltration of CD8 T cells, follicular helper T cells, gamma delta T cells, and M1 macrophages, whereas Cluster B and Cluster C had a high number of naïve B cells, resting memory CD4 T cells, activated NK cells, monocytes, resting dendritic cells, and neutrophils ( Figure 3F ).

To find the qualified gene indicator for the computation of TIC scores, we used the Boruta method to respectively select the important gene sets in feature α and feature β. A total of 109 DEGs ( Table S1 ) consisting of 13 DEGs in feature α and 96 DEGs in feature β were selected for PCA calculation. Furthermore, we acquired the TIC scores for each UM patient based on the earlier formula. After that, the optimal cutoff value of TIC score was used to divide UM patients into two subgroups: high score and low score. Log-rank test in KM curve suggested that high-TIC score patients had a good overall survival time than those in the low-TIC score group ( Figure 4A ). To assess the robustness of TIC scores, stratification analyses were also conducted. TCGA-UVM, GSE22138, and GSE84976 were accordingly separated into high-score and low-score groups. KM curves prove that high-TIC score patients had a good prognosis than low-TIC score patients regardless of TCGA-UVM ( Figure 4B ), GSE22138 ( Figure 4C ) and GSE84976 ( Figure 4D ).

An increasing number of results indicated that high mutation load closely associates with increased neoantigen expression and gives more chances for the immune system to recognize the tumor. Therefore, we used the “maftools” method to investigate the possible relationships between mutation load and TIC score. To begin, participants in the TCGA-UVM cohort were divided into two groups based on their TIC scores: low and high. The oncoPrint plots summarized the top 20 mutated genes in the low- and high-score groups ( Figure 4E ). GNAQ, GNA11, BAP1, SF3B1, and EIF1AX had 52%, 39%, 19%, 19%, and 19% mutation, respectively, in the high-score samples, while GNA11, BAP1, GNAQ, SF3B1, and DEPDC5 had 54%, 46%, 46%, 31%, and 8%, mutation, respectively, in the low-score samples. Despite tumor burden mutation (TMB) being an indicator for responsiveness to immunotherapy in various tumors, Spearman’s test revealed that TIC scores were not correlated with TMBs ( Figure 4F ). Furthermore, a forest plot indicated that BAP1 highly mutated in the low-score group and EIF1AX highly mutated in the high-score group ( Figure 4G ).

The possible associations of TIC scores with clinical parameters, molecular indicators, and biological signal pathways were further investigated. We firstly classified UM samples into low- and high-score groups, and then Alluvial diagram of TIC scores subgroup illustrated that high TIC score patients who had a large proportion of disomy in chromosome 3 status with no metastasis and alive status ( Figure 5A ). The different distributions of TIC scores in clinical subtypes were also estimated by Kruskal−Wallis or Wilcoxon test, which indicated that immune gene cluster ( Figure 5B ), infiltrated immune cell cluster ( Figure 5C ), metastatic status ( Figure 5D ), vital status ( Figure 5E ), and chromosome 3 status ( Figure 5F ) were intimately associated with TIC scores. Besides, to figure out the immunological activity and tolerance levels of each group, we selected PDCD1, CD274, CTLA4, LAG3, IDO1, HAVCR2, CD40, and CD40LG as immune checkpoint-related genes and picked BTLA, CXCL9, GZMA, CD8A, TBX2, PRF1, TNF, and TIGIT as immunological activity-related genes. Wilcoxon test suggested that the immune checkpoint-related and immunological activity-related genes were all significantly increased in the high-score group ( Figure 5G ). Finally, GSEA was used to discover the various signal pathways that were enriched in each TIC score group. The top five pathways were respectively depicted in the high-score group ( Figure 5H ) and the low-score group ( Figure 5I ) based on the selection standard and the ranking pathways enriched in each phenotype.

The relationship between TIC scores and the expression patterns of current immune checkpoint inhibitor-associated genes (PD-1, PD-L1, and CTLA-4) was examined further to better understand the possible response to immunotherapy. According to Spearman’s correlation tests, the TIC scores were significantly positively linked to the expression levels of PD-1 ( Figure 6A ), PD-L1 ( Figure 6B ), and CTLA-4 ( Figure 6C ). To investigate the effects of interaction between TIC scores and immune checkpoint genes on patients’ survival, these UM patients were categorized into four subtypes based on TIC scores and immune checkpoint genes, and KM curve analyses were conducted to determine the different survival times among the four subtypes. The log-rank tests indicated that TIC scores can effectively separate patients’ prognosis with contradictory expression levels of PD-1 ( Figure 6D ), PD-L1 ( Figure 6E ), and CTLA-4 ( Figure 6F ). UM patients with the worst prognosis had a low TIC score and a high level of immune checkpoint genes, whereas UM patients with high TIC scores and low levels of immune checkpoint genes had the highest survival rates among the four subtypes.

Because of the intimate connections between TIC scores and immune checkpoint inhibitor-associated genes, we hypothesized that it may be used to predict UM immunotherapy response. Therefore, we firstly used the TIDE module in an online website (http://tide.dfci.harvard.edu/) to estimate TIDE scores for all UM patients ( Figure 6G ). The chi-square test demonstrated that the high TIC score group has a higher response rate (40% vs. 21.57%) than the low-score group ( Figure 6H ). Furthermore, we used subclass mapping analysis to contrast the expression profiles of the high-/low-score subgroups with a previously published dataset of 47 melanoma patients who accepted the immune checkpoint inhibitor therapy (CTLA-4 and PD-1). Surprisingly, we discovered that the high-TIC score group has a better chance of responding to anti-PD-1 therapy. UM patients in the low-TIC score category, on the other hand, are unresponsive to anti-CTLA-4 or anti-PD-1 treatment ( Figure 6I ).

To support our hypothesis, melanoma patients from the TCGA-SKCM, GSE35640, and GSE78220 datasets who received immunotherapy were divided into low- and high-TIC score groups accordingly. In particular, melanoma patients in the high-score category had a longer survival time than those in the low-score group whether in TCGA-SKCM cohort ( Figure 7A ) or the GSE78220 cohort ( Figure 7E ). Besides, the high-score groups had a larger proportion of response rate than low-score groups whether in the TCGA-SKCM cohort ( Figure 7B ), the GSE78220 cohort ( Figure 7F ), and the GSE35640 cohort ( Figure 7I ). The box plots of TIC score distribution in TCGA-SKCM ( Figure 7C ), GSE78220 ( Figure 7G ), and GSE35640 ( Figure 7J ) prove that immunotherapy response groups had higher TIC scores than nonresponse groups. Eventually, ROC curves of TIC scores in TCGA-SKCM ( Figure 7D ), GSE78220 ( Figure 7H ), and GSE35640 ( Figure 7K ) suggested that the TIC scores had a relatively high accuracy to predict the response of immunotherapy.

Recently, numerous studies are searching for new promising therapeutic agents for advanced UMs due to resistance to standard chemotherapeutics. Therefore, we used multiple drug response databases to identify possible chemotherapy drugs for low-TIC score patients with poor prognosis. Eventually, we discovered that 6 compounds were negatively associated with TIC scores, which contained selumetinib and paclitaxel in the CTRP database ( Figure 8A ), cobimetinib and TAK-733 in the PRISM database ( Figure 8B ), and dasatinib and staurosporine in the GDSC database ( Figure 8C ). In addition, Wilcoxon test suggested that all of the six compounds had higher AUC values in the low-TIC score group than in the high-TIC score group.