Research Paper Volume 12, Issue 2 pp 1446—1464
Immunological analyses reveal an immune subtype of uveal melanoma with a poor prognosis
- 1 Department of Ophthalmology, Ninth People’s Hospital, Shanghai JiaoTong University School of Medicine, Shanghai, China
- 2 Shanghai Key Laboratory of Orbital Diseases and Ocular Oncology, Shanghai, China
- 3 Department of Pathology, Ninth People’s Hospital, Shanghai JiaoTong University School of Medicine, Shanghai, China
received: October 20, 2019 ; accepted: December 25, 2019 ; published: January 18, 2020 ;https://doi.org/10.18632/aging.102693
How to Cite
Copyright © 2020 Pan et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Uveal melanoma is an aggressive intraocular malignancy that often exhibits low immunogenicity. Metastatic uveal melanoma samples frequently exhibit monosomy 3 or BAP1 deficiency. In this study, we used bioinformatic methods to investigate the immune infiltration of uveal melanoma samples in public datasets. We first performed Gene Set Enrichment/Variation Analyses to detect immunological pathways that are altered in tumors with monosomy 3 or BAP1 deficiency. We then conducted an unsupervised clustering analysis to identify distinct immunologic molecular subtypes of uveal melanoma. We used CIBERSORT and ESTIMATE with RNA-seq data from The Cancer Genome Atlas and the GSE22138 microarray dataset to determine the sample-level immune subpopulations and immune scores of uveal melanoma samples. The Kaplan-Meier method and log-rank test were used to assess the prognostic value of particular immune cells and genes in uveal melanoma samples. Through these approaches, we discovered uveal melanoma-specific immunologic features, which may provide new insights into the tumor microenvironment and enhance the development of immunotherapies in the future.
Immune heterogeneity within the tumor microenvironment has been linked to the drug sensitivity and prognosis of patients with various cancer types [1, 2]. Thus, the profiling of immune signatures might uncover biomarkers for targeted therapy and clinical outcome assessment. Recently, datasets from The Cancer Genome Atlas (TCGA) have been used to depict the immune landscapes of multiple tumor types . Researchers have used integrated approaches and multidimensional datasets to determine the infiltration levels and co-infiltration networks of various immune cell populations in tumors [3, 4]. For instance, genomic data and hematoxylin & eosin image data were used to assess the total lymphocyte infiltration and immune cell fractions of the tumor microenvironment in different cancer types in TCGA . This analysis revealed common immune subtypes, immune gene expression signatures and tumor-extrinsic features, which could be used to identify transcriptional regulatory networks in the tumor microenvironment. Moreover, an extensive immunogenomic analysis of PanCancer TCGA data from 33 diverse cancer types revealed six distinct immune subtypes and various tumor-immune cell interactions .
Uveal melanoma (UM) is the most common aggressive intraocular malignancy in adults, and originates from the uveal tract . UM is characterized by different cytogenetic alterations than cutaneous melanoma (CM), and has the potential for hepatic metastasis . Nevertheless, UM and CM share a common lineage that is determined by melanoma-specific neural crest genes . Recently, a comprehensive analysis of 80 UM cases identified four molecularly distinct subtypes. Monosomy 3 (M3) tumors were found to be enriched for genes in immune pathways such as interferon signaling, T cell invasion and cytotoxicity . Several reports have demonstrated that CM is highly infiltrated by immune cells such as CD4 and CD8 cells [9–11]. However, only a subset of UM patients exhibit similar lymphocyte infiltration of their tumors, suggesting that UM immune infiltration is heterogeneous . In liver metastases, the tumor-infiltrating lymphocyte activity is lower in UM patients than in CM patients [5, 7, 12].
The inflammatory phenotype of UM is characterized by high infiltration of lymphocytes and macrophages, and by the expression of human leukocyte antigen (HLA) Class I and II antigens . Mutation of GNA11/GNAQ was not found to significantly alter the immune infiltration and HLA Class I expression of primary UM . However, the nuclear factor kappa-light-chain-enhancer of activated B cells (NF-kB) pathway was found to be associated with the inflammatory phenotype and high HLA Class I expression, and was upregulated upon the inactivation of BAP1 in UM. The levels of the NF-kB pathway molecules NF-kB1, NF-kB2 and RELB were reported to correlate positively with the expression of HLA Class I and with the infiltration of T cells and macrophages in UM . Moreover, one of the most significant UM studies revealed that lymphocyte infiltration and tumor-associated M2 macrophage levels were associated with a poor prognosis in primary UM after adjustment for other risk factors . In another study, the authors used a digital PCR-based T cell quantification method to characterize the prognostic value of the T cell count and activated macrophage level in the microenvironment of UM . Thus, characterizing the immunological features of UM may provide novel immune biomarkers for prognostic assessment and immunotherapy.
Immunotherapy through immune checkpoint blockades has displayed promising clinical efficacy in multiple tumor types [18–20]. Enhancing the cytolytic functions of infiltrating lymphocytes can significantly improve antitumor immunity . However, due to the low levels of cytotoxic cells in the tumor microenvironment, non-responsiveness to immunotherapy remains a clinical challenge. Thus, a deep comprehension of the interplay among the immune cell subsets in the tumor microenvironment is needed for the development of effective antitumor immune therapies. Here, we sought to identify the molecular immune subtypes and infiltrating immune cells of UM, in order to discover possible candidates for immunotherapy.
Identification of M3/BAP1-specific immunological pathways using GSEA
BAP1 is frequently mutated in metastatic UM, and is associated with chromosome 3 loss (M3) [16, 21]. Thus, we investigated whether BAP1-deficient UM samples exhibited distinct immune infiltration patterns from BAP1-intact samples in TCGA. Gene Set Enrichment Analysis (GSEA) was used to develop M3/BAP1null aberrant gene signatures and reveal novel immune pathways. The immune-related gene signatures were then used to estimate the level of immune cell infiltration. Among the most significantly altered pathways in M3/BAP1null tumors, those involving the CD8/T-cell-receptor, adaptive immune system, pre-B1 lymphocytes and differentiating T lymphocytes were significantly enriched in the M3 group (Figure 1A), consistent with a previous study . In curating the differential gene expression data, we found significant gene expression changes in a variety of immune-associated processes, suggesting that M3/BAP1null tumors highly express immune pathway genes.
Figure 1. BAP1-mutant/M3 UM is enriched in immune signatures. (A) GSEA plots of gene ontology categories, including the immune response, immune system process, adaptive immune response, immune effector process, and regulation of immune system process. (B) GSVA analysis of differing immune pathways between BAP1wild type and BAP1mutant tumors. (C) Differential proportions of immune cells between D3 and M3 tumors.
Next, we performed Gene Set Variation Analysis (GSVA) to compare the sample-level infiltration of M3/BAP1null tumors and disomy 3 (D3)/BAP1intact tumors. The top 10 significantly enriched immune pathway gene sets were selected. As shown in the heatmap (Figure 1B), the “B-cell receptor signaling” pathway, the “IRF4 targets in plasma cells vs. mature B lymphocytes” pathway and the “peripheral blood mononuclear cell response to ionizing radiation (IR)” pathway were significantly enriched in M3/BAP1null tumors, suggesting that IRF4 may enhance immunity in BAP1 deficient UMs.
To further assess the immune cell subpopulations of M3/BAP1null and D3/BAP1intact tumors, we compared the relative abundance of immune cells between these two UM subtypes. As shown in Figure 1C, the levels of infiltrating CD8 T cells and T follicular helper cells were significantly higher in M3-subtype tumors (n = 42) than in D3 tumors (n = 38) (P<0.01, Mann-Whitney test). Of note, D3 tumors were previously reported to have a better prognosis than M3 tumors. In contrast, monocytes and CD4 memory resting cell levels were higher in D3 tumors than in M3 tumors (P<0.01, Mann-Whitney test).
Hierarchical clustering of immune cell-associated gene expression in UM
We then performed an unsupervised clustering analysis of 730 immune-related genes in the UM dataset of TCGA, as described in a previous study . The sample clustering revealed three clear groups of samples that separated predominantly according to the gene expression of infiltrating immune cells, here termed the Immune Low (Immune L, n=54, 67.5%), Immune Medium (Immune M, n=16, 20%) and Immune High (Immune H, n=10, 12.5%) groups (Figure 2A). As shown in the heatmap, the Immune H group expressed high levels of the majority of the immune-related genes, in contrast to the Immune L group. The Immune H group highly expressed genes associated with CD8 T cells, B cells and natural killer cells (Figure 2B). The inhibitory checkpoint molecules PD-1 and CTLA-4 and the genes directly associated with MHC CLASS I/II, Cytolytic Activity and co stimulatory Molecules were also highly expressed in the Immune H group and modestly expressed in the Immune M group (Supplementary Figure 2).
Figure 2. Immune-related gene expression in the UM dataset of TCGA. (A) Hierarchical clustering of 80 tumors based on 730 immune-related genes. Genes were median-centered. Each colored square represents the relative mean transcript abundance (log2 FPKM+1) for each sample, with the lowest levels shown in green, the median levels in black and the highest levels in red. The genetic mutation type, SCNA type, immune score, leukocyte fraction and BAP1 mutation status are shown below the array tree. (B) The expression of selected gene signatures or genes is demonstrated below the heatmap.
Next, we sought to determine the clinical and molecular features underlying the immune clustering. The age, American Joint Committee on Cancer stage (AJCC), sex and GNAQ/11 mutation status of the patients did not differ significantly between the Immune H group and the Immune L and M groups (Immune H vs. M vs. L). Interestingly, the BAP1 mutation status differed significantly among the groups (P<0.0001, Fisher’s exact test, Immune H vs. M vs. L). Similarly, a significant difference in somatic copy number alteration (SCNA1/2 vs. ¾) was found among the groups (P<0.0001, Fisher’s exact test, Immune H vs. M vs. L) (Table 1). Our analysis suggested that BAP1 mutations and chromosome alterations may determine the tumor immune state and immune infiltration of UM.
Table 1. Clinical–pathologic characteristics of the TCGA Dataset in this study.
|Age, median (range)||60 (22–86)|
|AJCC Clinical Stage|
|Epithelioid cell dominant||23(28.8%)|
|Spindle Cell dominant||57(71.2%)|
Immune subtypes correlated with immune infiltration and clinical outcomes
ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) is a tool for predicting tumor purity and assessing stromal/immune cell infiltration into tumor tissues based on gene expression data . We used the ESTIMATE algorithm to calculate the immune score and stromal score for each UM sample from TCGA, and we compared the scores of the different immune subgroups. Since the immune score reflects the overall immune infiltration based on lymphocyte gene expression, we used this score to estimate the total immune cell infiltration of each sample. The Immune L group demonstrated the lowest immune and stromal scores, whereas the Immune M and Immune H groups displayed modest and high stromal and immune scores, respectively (P<0.0001, analysis of variance [ANOVA], Immune H vs. M vs. L) (Figure 3A). The tumor purity calculated by the ESTIMATE algorithm revealed that the Immune L group was the purest subtype (mean 0.9373) and the Immune H group was the least pure subtype (mean 0.7649) (P<0.0001, ANOVA) (Figure 3A).
Figure 3. Differential infiltration of the immune subgroups. (A) Box plots comparing the distribution of immune infiltration measures in the three immune subtypes. Each box spans the interquartile range, with the lines representing the median for each group. Whiskers represent the absolute range. All outliers are included in the plot. (B) Differential proportions of immune cells in the immune subtypes. (C) Kaplan-Meier survival curves demonstrate that the Immune M and H groups of UM patients consistently exhibited worse OS and PFIs than the Immune L group (log-rank test, P<0.05).
To further examine the tumor microenvironment, we assessed the infiltration of distinct immune subpopulations in the three tumor subtypes. The Immune L group had the lowest levels of CD8 T cells and total lymphocytes, whereas the Immune M and Immune H subtypes displayed relatively high levels of these cell types (P<0.0001, Immune H vs. M vs. L). In addition, the Immune H group had lower levels of M2 macrophages and mast cells than the other two groups. The Immune L group also exhibited low infiltration of T follicular helper cells and regulatory T cells (P<0.0001, Immune H vs. M vs. L) (Figure 3B), suggesting that there are distinct subtypes of UM with low or high T cell infiltration.
Next, we investigated the prognostic impact of immune clustering on patient survival. Interesting, we observed significantly worse overall survival (OS) and progression-free intervals (PFIs) in the Immune M and Immune H groups than in the Immune L group in the cohort from TCGA, suggesting that distinct immunological features correlate with the patient prognosis in UM (Figure 3C). However, the OS and PFI did not differ significantly between the Immune M and H groups (log-rank test, P>0.05, Immune H vs. M).
Correlations of immune cells and immune scores in TCGA and Laurent datasets
Next, we performed a hierarchical clustering analysis of the pairwise correlations between the different immune subpopulations in UM samples from TCGA (Figure 4A). The heatmap revealed high correlations among several types of immune cells associated with cytotoxic T cell infiltration (i.e., lymphocytes, regulatory T cells, T follicular helper cells, M1 macrophages and CD8 T cells). We obtained similar results using Laurent microarray data (Figure 4B). Therefore, we concluded that the identified immune cells were significantly associated with immune infiltration in the tumor microenvironment.
Figure 4. The correlation between the immune score and immune cell infiltration in UM. Pairwise correlation heatmap among immune cell-type scores in the datasets from TCGA (A) and GSE22138 (B).
We then examined the correlations between other immune parameters for UM patients in TCGA and Laurent UM dataset, and found a significant correlation between the immune score and the stromal score (Spearman’s r=0.806, 0.7789, respectively; Supplementary Figure 4A and 4B). We also observed a strong correlation between the immune score and T cell infiltration (Spearman’s r=0.7633, 0.4993, respectively; Supplementary Figure 4A and 4B). Additionally, a significant inverse correlation was found between the immune score and tumor purity (Spearman’s r =-0.9758, r=-0.9651 Supplementary Figure 4A and 4B). In the Laurent UM dataset, the association between the immune score and CD8 T cell infiltration remained significant (Supplementary Figure 4A and 4B), supporting the ability of the immune score to predict immune infiltration. Therefore, the immune score based on immune marker gene expression reflected immune cell infiltration in both the UM TCGA RNA-seq dataset and an independent microarray dataset.
Prognostic value of immune cells and immune scores in TCGA and Laurent datasets
To assess the prognostic value of immune cells in UM, we first used a Univariate Cox regression model to investigate the associations of different immune subpopulations with the OS and PFI of patients in TCGA. The infiltration of M1 macrophages, activated natural killer cells and CD8 T cells was associated with a worse prognosis in the UM datasets (Figure 5A). On the contrary, the infiltration of monocytes was positively associated with patient survival (Figure 5A, 5B). In general, the infiltration of immune subpopulations involved in adaptive immunity is more likely to be an unfavorable prognostic factor than the infiltration of subpopulations involved in innate immunity.
Figure 5. The prognostic value of the immune score and immune cell infiltration in UM. (A) HRs of OS and the PFI based on the infiltration of various immune cells (as continuous variables) in all patients (left); the horizontal bars represent the 95% confidence intervals of the HRs. Statistically significant variables are shown. Each cell type was evaluated individually and rank-ordered based on the estimated HR. (B) Kaplan-Meier survival analysis based on immune score and stroma score. Patients were divided into the high and low groups based on the level of immune score and stroma score. (C) Kaplan-Meier survival analysis based on selected immune cells. Patients were divided into the high and low groups based on their expression of each cell.
To determine whether the immune score could predict the prognosis of UM patients in TCGA, we used a Kaplan-Meier curve and log-rank test to estimate the hazard ratios (HRs) of OS and the PFI. High immune scores and stromal scores were associated with worse OS in the datasets of TCGA (HR=6.721, P<0.0001, Figure 5B) and GSE22138 (HR=2.508, P=0.02, Supplementary Figure 1A). Next, we determined the associations of different immune cells with the patient prognosis. Consistent with the Univariate Cox regression analysis, our survival analysis revealed that particular immune subpopulations (M1 and M0 macrophages, CD8 T cells and T follicular helper cells) were associated with a worse prognosis, whereas total CD4 cells were associated with a better prognosis in the UM dataset of TCGA. Similar results were obtained with the GSE22138 dataset (Supplementary Figure 1B). These results suggest that UM patients with greater infiltration of immune cells and effector molecules exhibit poorer survival and may benefit from immunotherapies.
Immune marker genes predict the prognosis of UM patients
Next, we examined the prognostic value of individual immune genes in predicting patient survival. Higher mRNA levels of CD8A (P<0.0001), HLA-A (P<0.0001), HLA-B (P=0.0001), HLA-C (P=0.0003) and HLA-DRA (P=0.0001) in TCGA samples were associated with significantly shorter OS and PFIs (Figure 6A, 6B). Similarly, higher levels of HLA-A (P=0.02), HLA-B (P=0.002), HLA-C (P=0.007) and HLA-DRA (P=0.03) were associated with significantly shorter metastasis-free survival in the 63 patients in the Laurent UM dataset (Figure 6C), and CD8A expression (P=0.06) displayed a similar but non-significant trend of association with metastasis-free survival.
Figure 6. The prognostic value of individual immune genes in UM. (A) Kaplan-Meier survival curves demonstrate that elevated levels of CD8A (P<0.0001), HLA-A (P<0.0001), HLA-B (P=0.0001), HLA-C (P=0.0003) and HLA-DRA (P=0.0001) were consistently associated with worse OS in the UM dataset of TCGA. (B) Kaplan-Meier survival curves demonstrate that elevated levels of CD8A, HLA-A, HLA-B, HLA-C and HLA-DRA were consistently associated with a worse PFI in the UM dataset of TCGA (log-rank test, P<0.05). (C) Kaplan-Meier survival curves demonstrate that elevated levels of CD8A (P=0.06), HLA-A (P=0.02), HLA-B (P=0.002), HLA-C (P=0.007) and HLA-DRA (P=0.03) were consistently associated with worse metastasis-free survival in the Laurent UM dataset (n=63; log-rank test, P<0.05).
TCGA has illuminated the genomic data from bulk tumor samples, and has provided detailed information about the tumor immune microenvironment [1, 24–26]. In previous studies, low immune cell infiltration has been associated with poor clinical outcomes for patients with different cancers [1, 27]. Through gene expression profiling, researchers can identify prognostic gene signatures and detect candidate genes for targeted therapies . For example, an immune score based on gene expression data was found to correlate significantly with recurrence-free survival in thyroid cancer patients, regardless of their BRAF(V600E) status . The fraction of immune cells in clinical tumor samples can be evaluated by multiple algorithms; indeed, aside from curating samples and performing basic pathologic characterization, investigators can analyze digitized hematoxylin & eosin-stained images of TCGA samples for tumor-infiltrating lymphocytes . Using TCGA data, we identified three distinct immune subtypes of UM, with prognostic implications for immunological cancer management. The immune score was strongly associated with immune infiltration and poor outcomes, regardless of the tumor genome ploidy of the UM tumor samples. The underlying mechanism needs to be explored.
Genomic and transcriptomic data have been used to detect immune infiltration and to determine the molecular subtypes of ovarian cancer, melanoma and pancreatic cancer [29–31]. For instance, DNA sequencing data have been used to connect the neoantigen load to the T cell response and to link somatic mutations to immune infiltration [1, 32]. More recently, deconvoluted expression data have been used to measure the cytolytic activity in the tumor microenvironment and to quantify the infiltration of individual immune cell subsets [33–35]. A common theme across these studies is the integration of several types of genomic and clinical data, allowing for associations to be made among immune activity, gene expression, the mutation burden and patient survival. In this study, we examined the immune infiltration of UM samples through a single-sample GSEA and deconvolution method based on publicly available data in TCGA and the Gene Expression Omnibus. Our results suggested that the infiltration of immune cells differs markedly among immune subtypes. This analysis may ultimately reveal prognostic gene signatures and provide candidate genes for targeted therapies.
Despite the great success of immunotherapies against metastatic and late-stage melanoma [36, 37], immunotherapy has had limited success in UM . UM is considered to be an immunotherapy-resistant subtype of melanoma, and UM patients are frequently excluded from clinical trials of immunotherapies for metastatic melanoma [20, 39, 40]. Although higher cytotoxic expression patterns are associated with better anti-tumor response and better patient survival in many solid tumors [34, 41, 42], the prognostic effects of immune infiltration depend on the type of tumor, the location of the cells and the state of activation. In UM, high levels of immune cells are associated with poor prognostic factors, such as M3 and BAP1 mutation [15, 17, 43]. Immune cell infiltration occurs more frequently in epithelioid-cell-type UM, which also has a poor prognosis . Crosstalk in the tumor microenvironment can promote the inflammatory response in cancer cells. Cancer cells may also promote the type 2 differentiation of macrophages and neutrophils, and may attract myeloid-derived suppressor cells and regulatory T cells to tumor sites. Thus, we speculate that UM cells may utilize immune cells for their survival and protection from immunological attack. The immunomodulatory microenvironment in the liver could further protect escaped UM cells from systemic immune surveillance [5, 44, 45].
De Lange et al have used unsupervised clustering to investigate the gene expression profiles of 64 enucleated eyes from UM patients, and divided them into class I tumors with a good prognosis and class IIa and IIb tumors with a poor prognosis . Their study revealed an immune phenotype with a different prognosis. High expression of immune-related genes in class IIb UM suggested that the tumors were inflamed. Furthermore, study from TCGA of UM showed that the genes encoding chemotactic signals (e.g., CXCL9 and CXCL13), MHC class I (A, B, C) and MHC class II (DP, DM, DOA, DOB, DQ and DR) were upregulated in M3 patients . Consistent with the previous studies, we demonstrated that BAP1 inactivation was associated with immune infiltration and immune marker gene set expression, indicating the BAP1 may regulate tumor immunology. GSVA results suggested that IRF4 targets and BCR pathways may be induced in BAP1-deficient tumors. Loss of BAP1 expression is also associated with an increased infiltration of T cell follicular helper, Treg and CD8+ T cells, suggesting an inflammatory tumor microenvironment. Our data demonstrated that the immune cell subpopulations were differentially distributed between M3/BAP1null and D3/BAP1intact tumors, suggesting that BAP1 null tumors might be prioritized for immune checkpoint blockade therapies in UM.
At the time of this work, large, publicly available gene expression profiling datasets of UM patients treated with immune checkpoint blockers were not available. In addition, clinical trials in immunotherapy are being deployed earlier in the course of the disease, whereas the cohort in TCGA is more representative of the clinical population. We hope that large sequencing data from UM patients undergoing immune checkpoint blocker treatment will emerge in the future. Nevertheless, our analysis of the available datasets has advanced the application of genomic data to tumor immunology. The immune features reported herein should be considered for integration into prognostic models, or explored as predictors of adjuvant immune therapy responsiveness in patients with BAP1-deficient UM.
Materials and Methods
UM gene expression datasets
RNA-seq data for UM samples were generated by TCGA and downloaded from the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/). This dataset includes normalized gene expression profiles for 80 tumor samples in FPKM (Fragments Per Kilobase of transcript per Million fragments mapped). We downloaded one additional microarray gene expression dataset from the Gene Expression Omnibus database (accession number: GSE22138 ; n=63) to verify the immune scores and gene signatures and to predict the patients’ prognoses. GSE22138 included information on metastasis-free survival, while the UM dataset from TCGA included data on patients’ OS and PFI from the supplemental file of a previous study . The clinical information and molecular data for the UM samples from TCGA were also downloaded from the supplemental file of a previous study . The clinical information included the age, sex, metastatic status, histology cell type, American Joint Committee on Cancer clinical stage, mutation data, SCNA data and vital status of the patients.
We used the complete linkage method for hierarchical clustering analysis of the tumor samples, immune cell types and genes. The hierarchical clustering algorithm is agglomerative in that it joins samples based on a measure of multivariate distance, and prevents the joined samples from clustering independently again. Pair-wise joins in samples are represented as combined branches of a tree in a Dendrogram Plot. Pearson’s correlation distance method was used to determine whether samples clustered together .
ESTIMATE is a tool that predicts tumor purity and detects infiltrating stromal/immune cells in tumor tissues based on gene expression data . The ESTIMATE algorithm calculates the stromal and immune scores by performing single-sample GSEA for each sample. Data on the leukocyte fractions of the 80 tumor aliquots from TCGA were obtained from a previous study .
CIBERSORT (immune cellular fraction estimates)
The relative fractions of 22 immune cell types within the leukocyte compartment were estimated with CIBERSORT (Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts) . CIBERSORT requires a specialized knowledgebase of gene expression signatures, termed a “signature matrix,” to deconvolute cell types of interest. CIBERSORT uses a set of 22 immune cell reference profiles (LM22)  to derive a base (signature) matrix that can be applied to mixed samples to determine the relative proportions of immune cells. LM22 is a signature matrix file consisting of 547 genes that can precisely distinguish 22 mature human hematopoietic populations (from peripheral blood or in vitro culture), including seven T cell types, naïve and memory B cells, plasma cells, natural killer cells, myeloid subsets, etc. LM22 can be applied to RNA-seq data as well as to microarray data. We used CIBERSORT to re-quantify several key immune gene signatures from GSE22138, and to identify the immune cells in the dataset from TCGA, as in a previous PanCancer immune study .
Predicting patient survival with the Cox proportional hazards model
A Cox proportional hazards regression model was used to investigate the effectiveness of immune cells in predicting patients’ survival (PFI or OS). The optimal cut-off values of the proportions of immune cells in the cohort from TCGA were calculated based on their prognostic effects in X-Tile software . The UM samples were then divided into high and low groups based on the cut-off point for the fraction level of each type of immune cell. In the subsequent scoring formula, the immune cell fraction level was given a value of 0 or 1. In STATA15, a univariate Cox proportional hazards model was used to evaluate the effects of significant prognostic immune cell fractions on survival outcomes. Kaplan-Meier plots were used to visualize the differences in survival between groups. The survival proportions of the different groups were compared through a log-rank test.
GSVA was performed as previously described  with the curated gene set collection of the Molecular Signatures Database (Broad Institute, c2. CP CPG). We searched for significant pathway gene sets that were differentially enriched between the BAP1null and BAP1intact groups (FDR≤0.01). The top 10 immune-related pathways were selected and shown in the heatmap.
Gene set collections for canonical pathways (C2.CP canonical pathways) were downloaded from the Molecular Signatures Database version 6.2 (www.broad institute.org/gsea/msigdb/collections.jsp). Gene set enrichment scores were calculated  with GSEA package version 1.32.0 with RNA-seq parameters. Differential gene set enrichment was determined using the limma package. The thresholds for statistical significance are noted in the Results.
All statistical analyses in this study were performed with R version 3.2.0 (R Foundation for Statistical Computing, Vienna, Austria) and SPSS Statistics 22.0. Survival curves were generated by the Kaplan-Meier method, and the log-rank test was used to compare the survival curves. The thresholds for each Kaplan-Meier plot were determined based on the unique score distribution of each dataset in X-Tile software . We used Spearman's rank correlation coefficients to evaluate the correlation between the immune score and immune cell infiltration. The immune scores and stromal scores stratified by subgroup were compared by one-way ANOVA. A two-tailed P-value <0.05 was considered statistically significant.
AJCC: American Joint Committee on Cancer stage; ANOVA: analysis of variance; CIBERSORT: Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts; CM: cutaneous melanoma; D3: Disomy3; ESTIMATE: Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data; FPKM: Fragments Per Kilobase of transcript per Million fragments mapped; GSEA/GSVA: Gene Set Enrichment/Variation Analysis; HR: hazard ratio; IR: Ionizing radiation; M3: Monosomy3; PFI: Progression-free interval; SCNA: somatic copy number alteration; TCGA: The Cancer Genome Atlas; UM: Uveal melanoma.
Pan H designed and drafted the manuscript; Lu L analyzed the data; Cui J and Yang Y reviewed and edited the figures; Wang Z and Fan X discussed and approved the manuscript. All authors approved this manuscript.
Conflicts of Interest
The author(s) declare that they have no conflicts of interest.
The work was supported by the National Key R&D Program of China (2018YFC1106100, 2018YFC11 06101), The Science and Technology Commission of Shanghai (17DZ2260100, 19JC1410200), the National Natural Science Foundation of China (Grant U1932135, 81770934, 81702781), Shanghai Municipal Education Commission- Gaofeng Clinical Medicine Grant Support (20181810).
- 1. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, Porta-Pardo E, Gao GF, Plaisier CL, Eddy JA, Ziv E, Culhane AC, Paull EO, et al. The Immune Landscape of Cancer. Immunity. 2018; 48:812–830.e14. https://doi.org/10.1016/j.immuni.2018.03.023 [PubMed]
- 2. Cesano A, Warren S. Bringing the next Generation of Immuno-Oncology Biomarkers to the Clinic. Biomedicines. 2018; 6:E14. https://doi.org/10.3390/biomedicines6010014 [PubMed]
- 3. Saltz J, Gupta R, Hou L, Kurc T, Singh P, Nguyen V, Samaras D, Shroyer KR, Zhao T, Batiste R, Van Arnam J; Cancer Genome Atlas Research Network, Shmulevich I, et al. Spatial Organization and Molecular Correlation of Tumor-Infiltrating Lymphocytes Using Deep Learning on Pathology Images. Cell Rep. 2018; 23:181–193.e7. https://doi.org/10.1016/j.celrep.2018.03.086 [PubMed]
- 4. Varn FS, Wang Y, Mullins DW, Fiering S, Cheng C. Systematic Pan-Cancer Analysis Reveals Immune Cell Interactions in the Tumor Microenvironment. Cancer Res. 2017; 77:1271–82. https://doi.org/10.1158/0008-5472.CAN-16-2490 [PubMed]
- 5. Oliva M, Rullan AJ, Piulats JM. Uveal melanoma as a target for immune-therapy. Ann Transl Med. 2016; 4:172. https://doi.org/10.21037/atm.2016.05.04 [PubMed]
- 6. Laurent C, Valet F, Planque N, Silveri L, Maacha S, Anezo O, Hupe P, Plancher C, Reyes C, Albaud B, Rapinat A, Gentien D, Couturier J, et al. High PTP4A3 phosphatase expression correlates with metastatic risk in uveal melanoma patients. Cancer Res. 2011; 71:666–74. https://doi.org/10.1158/0008-5472.CAN-10-0605 [PubMed]
- 7. Rothermel LD, Sabesan AC, Stephens DJ, Chandran SS, Paria BC, Srivastava AK, Somerville R, Wunderlich JR, Lee CC, Xi L, Pham TH, Raffeld M, Jailwala P, et al. Identification of an Immunogenic Subset of Metastatic Uveal Melanoma. Clin Cancer Res. 2016; 22:2237–49. https://doi.org/10.1158/1078-0432.CCR-15-2294 [PubMed]
- 8. Robertson AG, Shih J, Yau C, Gibb EA, Oba J, Mungall KL, Hess JM, Uzunangelov V, Walter V, Danilova L, Lichtenberg TM, Kucherlapati M, Kimes PK, et al. Integrative Analysis Identifies Four Molecular and Clinical Subsets in Uveal Melanoma. Cancer Cell. 2017; 32:204–220.e15. https://doi.org/10.1016/j.ccell.2017.07.003 [PubMed]
- 9. Bertrand F, Montfort A, Marcheteau E, Imbert C, Gilhodes J, Filleron T, Rochaix P, Andrieu-Abadie N, Levade T, Meyer N, Colacios C, Ségui B. TNFα blockade overcomes resistance to anti-PD-1 in experimental melanoma. Nat Commun. 2017; 8:2256. https://doi.org/10.1038/s41467-017-02358-7 [PubMed]
- 10. Iglesia MD, Parker JS, Hoadley KA, Serody JS, Perou CM, Vincent BG. Genomic Analysis of Immune Cell Infiltrates Across 11 Tumor Types. J Natl Cancer Inst. 2016; 108:djw144. https://doi.org/10.1093/jnci/djw144 [PubMed]
- 11. Prat A, Navarro A, Paré L, Reguart N, Galván P, Pascual T, Martínez A, Nuciforo P, Comerma L, Alos L, Pardo N, Cedrés S, Fan C, et al. Immune-Related Gene Expression Profiling After PD-1 Blockade in Non-Small Cell Lung Carcinoma, Head and Neck Squamous Cell Carcinoma, and Melanoma. Cancer Res. 2017; 77:3540–50. https://doi.org/10.1158/0008-5472.CAN-16-3556 [PubMed]
- 12. Qin Y, Petaccia de Macedo M, Reuben A, Forget MA, Haymaker C, Bernatchez C, Spencer CN, Gopalakrishnan V, Reddy S, Cooper ZA, Fulbright OJ, Ramachandran R, Wahl A, et al. Parallel profiling of immune infiltrate subsets in uveal melanoma versus cutaneous melanoma unveils similarities and differences: A pilot study. Oncoimmunology. 2017; 6:e1321187. https://doi.org/10.1080/2162402X.2017.1321187 [PubMed]
- 13. Souri Z, Wierenga AP, Mulder A, Jochemsen AG, Jager MJ. HLA Expression in Uveal Melanoma: An Indicator of Malignancy and a Modifiable Immunological Target. Cancers (Basel). 2019; 11:E1132. https://doi.org/10.3390/cancers11081132 [PubMed]
- 14. van Weeghel C, Wierenga AP, Versluis M, van Hall T, van der Velden PA, Kroes WG, Pfeffer U, Luyten GP, Jager MJ. Do GNAQ and GNA11 Differentially Affect Inflammation and HLA Expression in Uveal Melanoma? Cancers (Basel). 2019; 11:E1127. https://doi.org/10.3390/cancers11081127 [PubMed]
- 15. Souri Z, Wierenga AP, van Weeghel C, van der Velden PA, Kroes WG, Luyten GP, van der Burg SH, Jochemsen AG, Jager MJ. Loss of BAP1 Is Associated with Upregulation of the NFkB Pathway and Increased HLA Class I Expression in Uveal Melanoma. Cancers (Basel). 2019; 11:E1102. https://doi.org/10.3390/cancers11081102 [PubMed]
- 16. Bronkhorst IH, Vu TH, Jordanova ES, Luyten GP, Burg SH, Jager MJ. Different subsets of tumor-infiltrating lymphocytes correlate with macrophage influx and monosomy 3 in uveal melanoma. Invest Ophthalmol Vis Sci. 2012; 53:5370–78. https://doi.org/10.1167/iovs.11-9280 [PubMed]
- 17. de Lange MJ, Nell RJ, Lalai RN, Versluis M, Jordanova ES, Luyten GP, Jager MJ, van der Burg SH, Zoutman WH, van Hall T, van der Velden PA. Digital PCR-Based T-cell Quantification-Assisted Deconvolution of the Microenvironment Reveals that Activated Macrophages Drive Tumor Inflammation in Uveal Melanoma. Mol Cancer Res. 2018; 16:1902–11. https://doi.org/10.1158/1541-7786.MCR-18-0114 [PubMed]
- 18. Combe P, de Guillebon E, Thibault C, Granier C, Tartour E, Oudard S. Trial Watch: therapeutic vaccines in metastatic renal cell carcinoma. Oncoimmunology. 2015; 4:e1001236. https://doi.org/10.1080/2162402X.2014.1001236 [PubMed]
- 19. Hammers H. Immunotherapy in kidney cancer: the past, present, and future. Curr Opin Urol. 2016; 26:543–47. https://doi.org/10.1097/MOU.0000000000000338 [PubMed]
- 20. Wierenga AP, Cao J, Luyten GP, Jager MJ. Immune Checkpoint Inhibitors in Uveal and Conjunctival Melanoma. Int Ophthalmol Clin. 2019; 59:53–63. https://doi.org/10.1097/IIO.0000000000000263 [PubMed]
- 21. Harbour JW, Onken MD, Roberson ED, Duan S, Cao L, Worley LA, Council ML, Matatall KA, Helms C, Bowcock AM. Frequent mutation of BAP1 in metastasizing uveal melanomas. Science. 2010; 330:1410–13. https://doi.org/10.1126/science.1194472 [PubMed]
- 22. Robertson AG, Shih J, Yau C, Gibb EA, Oba J, Mungall KL, Hess JM, Uzunangelov V, Walter V, Danilova L, Lichtenberg TM, Kucherlapati M, Kimes PK, et al, and TCGA Research Network. Integrative Analysis Identifies Four Molecular and Clinical Subsets in Uveal Melanoma. Cancer Cell. 2018; 33:151. https://doi.org/10.1016/j.ccell.2017.12.013 [PubMed]
- 23. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
- 24. Bailey MH, Tokheim C, Porta-Pardo E, Sengupta S, Bertrand D, Weerasinghe A, Colaprico A, Wendl MC, Kim J, Reardon B, Kwok-Shing Ng P, Jeong KJ, Cao S, et al, and MC3 Working Group, and Cancer Genome Atlas Research Network. Comprehensive Characterization of Cancer Driver Genes and Mutations. Cell. 2018; 174:1034–35. https://doi.org/10.1016/j.cell.2018.07.034 [PubMed]
- 25. Liu J, Lichtenberg T, Hoadley KA, Poisson LM, Lazar AJ, Cherniack AD, Kovatich AJ, Benz CC, Levine DA, Lee AV, Omberg L, Wolf DM, Shriver CD, et al. An Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome Analytics. Cell. 2018; 173:400–416.e11. https://doi.org/10.1016/j.cell.2018.02.052 [PubMed]
- 26. Li B, Cui Y, Nambiar DK, Sunwoo JB, Li R. The Immune Subtypes and Landscape of Squamous Cell Carcinoma. Clin Cancer Res. 2019; 25:3528–37. https://doi.org/10.1158/1078-0432.CCR-18-4085 [PubMed]
- 27. He B, Gao R, Lv D, Wen Y, Song L, Wang X, Lin S, Huang Q, Deng Z, Wang Z, Yan M, Zheng F, Lam EW, et al. The prognostic landscape of interactive biological processes presents treatment responses in cancer. EBioMedicine. 2019; 41:120–33. https://doi.org/10.1016/j.ebiom.2019.01.064 [PubMed]
- 28. Geraldo MV, Kimura ET. Integrated Analysis of Thyroid Cancer Public Datasets Reveals Role of Post-Transcriptional Regulation on Tumor Progression by Targeting of Immune System Mediators. PLoS One. 2015; 10:e0141726. https://doi.org/10.1371/journal.pone.0141726 [PubMed]
- 29. Zhao Y, Schaafsma E, Gorlov IP, Hernando E, Thomas NE, Shen R, Turk MJ, Berwick M, Amos CI, Cheng C. A Leukocyte Infiltration Score Defined by a Gene Signature Predicts Melanoma Patient Prognosis. Mol Cancer Res. 2019; 17:109–19. https://doi.org/10.1158/1541-7786.MCR-18-0173 [PubMed]
- 30. Yoshihara K, Tsunoda T, Shigemizu D, Fujiwara H, Hatae M, Fujiwara H, Masuzaki H, Katabuchi H, Kawakami Y, Okamoto A, Nogawa T, Matsumura N, Udagawa Y, et al, and Japanese Serous Ovarian Cancer Study Group. High-risk ovarian cancer based on 126-gene expression signature is uniquely characterized by downregulation of antigen presentation pathway. Clin Cancer Res. 2012; 18:1374–85. https://doi.org/10.1158/1078-0432.CCR-11-2725 [PubMed]
- 31. Danilova L, Ho WJ, Zhu Q, Vithayathil T, De Jesus-Acosta A, Azad NS, Laheru DA, Fertig EJ, Anders R, Jaffee EM, Yarchoan M. Programmed Cell Death Ligand-1 (PD-L1) and CD8 Expression Profiling Identify an Immunologic Subtype of Pancreatic Ductal Adenocarcinomas with Favorable Survival. Cancer Immunol Res. 2019; 7:886–95. https://doi.org/10.1158/2326-6066.CIR-18-0822 [PubMed]
- 32. Matsushita H, Sato Y, Karasaki T, Nakagawa T, Kume H, Ogawa S, Homma Y, Kakimi K. Neoantigen Load, Antigen Presentation Machinery, and Immune Signatures Determine Prognosis in Clear Cell Renal Cell Carcinoma. Cancer Immunol Res. 2016; 4:463–71. https://doi.org/10.1158/2326-6066.CIR-15-0225 [PubMed]
- 33. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol. 2018; 1711:243–59. https://doi.org/10.1007/978-1-4939-7493-1_12 [PubMed]
- 34. Narayanan S, Kawaguchi T, Yan L, Peng X, Qi Q, Takabe K. Cytolytic Activity Score to Assess Anticancer Immunity in Colorectal Cancer. Ann Surg Oncol. 2018; 25:2323–31. https://doi.org/10.1245/s10434-018-6506-6 [PubMed]
- 35. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12:453–57. https://doi.org/10.1038/nmeth.3337 [PubMed]
- 36. Eggermont AM, Blank CU, Mandala M, Long GV, Atkinson VG, Dalle S, Haydon A, Lichinitser M, Khattak A, Carlino MS, Sandhu S, Larkin J, Puig S, et al. Prognostic and predictive value of AJCC-8 staging in the phase III EORTC1325/KEYNOTE-054 trial of pembrolizumab vs placebo in resected high-risk stage III melanoma. Eur J Cancer. 2019; 116:148–57. https://doi.org/10.1016/j.ejca.2019.05.020 [PubMed]
- 37. Robert C, Ribas A, Schachter J, Arance A, Grob JJ, Mortier L, Daud A, Carlino MS, McNeil CM, Lotem M, Larkin JM, Lorigan P, Neyns B, et al. Pembrolizumab versus ipilimumab in advanced melanoma (KEYNOTE-006): post-hoc 5-year results from an open-label, multicentre, randomised, controlled, phase 3 study. Lancet Oncol. 2019; 20:1239–51. https://doi.org/10.1016/S1470-2045(19)30388-2 [PubMed]
- 38. Heppt MV, Steeb T, Schlager JG, Rosumeck S, Dressler C, Ruzicka T, Nast A, Berking C. Immune checkpoint blockade for unresectable or metastatic uveal melanoma: A systematic review. Cancer Treat Rev. 2017; 60:44–52. https://doi.org/10.1016/j.ctrv.2017.08.009 [PubMed]
- 39. Hou F, Huang QM, Hu DN, Jonas JB, Wei WB. Immune oppression array elucidating immune escape and survival mechanisms in uveal melanoma. Int J Ophthalmol. 2016; 9:1701–12. https://doi.org/10.18240/ijo.2016.12.01 [PubMed]
- 40. Achberger S, Aldrich W, Tubbs R, Crabb JW, Singh AD, Triozzi PL. Circulating immune cell and microRNA in patients with uveal melanoma developing metastatic disease. Mol Immunol. 2014; 58:182–86. https://doi.org/10.1016/j.molimm.2013.11.018 [PubMed]
- 41. Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. 2015; 160:48–61. https://doi.org/10.1016/j.cell.2014.12.033 [PubMed]
- 42. Wakiyama H, Masuda T, Motomura Y, Hu Q, Tobo T, Eguchi H, Sakamoto K, Hirakawa M, Honda H, Mimori K. Cytolytic Activity (CYT) Score Is a Prognostic Biomarker Reflecting Host Immune Status in Hepatocellular Carcinoma (HCC). Anticancer Res. 2018; 38:6631–38. https://doi.org/10.21873/anticanres.13030 [PubMed]
- 43. Bakhoum MF, Esmaeli B. Molecular Characteristics of Uveal Melanoma: Insights from the Cancer Genome Atlas (TCGA) Project. Cancers (Basel). 2019; 11:E1061. https://doi.org/10.3390/cancers11081061 [PubMed]
- 44. de Waard-Siebinga I, Hilders CG, Hansen BE, van Delft JL, Jager MJ. HLA expression and tumor-infiltrating immune cells in uveal melanoma. Graefes Arch Clin Exp Ophthalmol. 1996; 234:34–42. https://doi.org/10.1007/BF00186516 [PubMed]
- 45. Barnes TA, Amir E. HYPE or HOPE: the prognostic value of infiltrating immune cells in cancer. Br J Cancer. 2018; 118:e5. https://doi.org/10.1038/bjc.2017.417 [PubMed]
- 46. de Lange MJ, van Pelt SI, Versluis M, Jordanova ES, Kroes WG, Ruivenkamp C, van der Burg SH, Luyten GP, van Hall T, Jager MJ, van der Velden PA. Heterogeneity revealed by integrated genomic analysis uncovers a molecular switch in malignant uveal melanoma. Oncotarget. 2015; 6:37824–35. https://doi.org/10.18632/oncotarget.5637 [PubMed]
- 47. Liu X, Zhu XH, Qiu P, Chen W. A correlation-matrix-based hierarchical clustering method for functional connectivity analysis. J Neurosci Methods. 2012; 211:94–102. https://doi.org/10.1016/j.jneumeth.2012.08.016 [PubMed]
- 48. Camp RL, Dolled-Filhart M, Rimm DL. X-tile: a new bio-informatics tool for biomarker assessment and outcome-based cut-point optimization. Clin Cancer Res. 2004; 10:7252–59. https://doi.org/10.1158/1078-0432.CCR-04-0713 [PubMed]
- 49. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
- 50. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005; 102:15545–50. https://doi.org/10.1073/pnas.0506580102 [PubMed]