1.
Introduction
Uveal melanoma(UVM)is the most common intraocular cancer in adults and the second most common melanoma. UVM has the characteristics of easy metastasis and poor prognosis. Primary UVM can be effectively controlled by surgery and radiotherapy, but the five-year survival rate is 50-70% [1]. About 50% of patients had metastasis, the prognosis after metastasis was very poor, and the median survival time was less than one year [2,3,4]. At present, the treatment of metastatic UVM is limited and ineffective [5,6]. Unlike in patients with cutaneous melanoma, immune checkpoint inhibitors are usually ineffective in patients with metastatic UVM. Due to the high metastasis rate of UVM and the high mortality rate of patients with metastasis, it is necessary to diagnose and treat UVM early, improve the survival rate of patients, and explore biological markers with prognostic value.
The tumor microenvironment(TME)is a complex ecosystem composed of cellular and non-cellular components, including tumor cells, immune cells, the extracellular matrix, and cytokines. TME acted as a pivotal part of the development of UVM [7,8]. In older mice, macrophages associated with tumors were directly related to the development of UVM [9]. Studies have found that the infiltration of immune cells in UVM is poor [10], whereas highly infiltrated macrophages and T lymphocytes can lead to poor prognosis [11]. Inflammatory mediators are critical factors in UVM metastasis [12]. Genes in activated macrophage clusters have been associated with low survival in metastatic patients [13]. Therefore, it is crucial to explore the TME of UVM to identify effective diagnostic and therapeutic markers.
Bioinformatics analysis refers to the use of computer technology and information technology to analyze experimental biological data. These data are stored in databases, such as The Cancer Genome Atlas (TCGA), which provides a convenient platform for researchers to carry out basic research. In 2013, Yoshihara et al. developed the ESTIMATE algorithm for predicting TME with the immune stromal index [14]. The researchers used this algorithm to discuss gene variation and prognostic value in a variety of tumors [15,16,17]. However, the value of the immune stromal index in UVM remains unclear.
In this study, we calculated the differentially expressed genes (DEGs) related to the immune index and the stromal index of UVM in the Gene Expression Omnibus (GEO) chip and TCGA, constructed the prognosis model using univariate and multivariate Cox regression, and analyzed the clinical significance of the immune/stromal index in UVM and its potential molecular mechanisms. This study demonstrates that the risk signature based on the immune index can be used as an indicator to evaluate the prognosis of patients with UVM.
2.
Materials and methods
2.1. mRNA sequence data processing and clinical data of the patients
We downloaded the level three Fragments Per Kilobase Million (FPKM) data from 80 patients with UVM from the TCGA database (https://gdc.cancer.gov/) and standardized FPKM data with log2 (Transcripts Per Million (TPM) + 0.001). The clinical data for UVM were downloaded which included the unique number of the patients, age, gender, eye color, clinical stage, pathology stage, T stage, metastasis, tumor site, extrascleral extension, tumor shape and prognostic index.
2.2. Calculation of immune stromal index and DEGs
We used the ESTIMATE package of R software to calculate the immune index scores and the stromal index scores based on the standardized data of mRNA-seq and GSE22138. Samples were grouped into high and low immune or stromal score groups according to the median values of immune and stromal scores. We compared the expression value of genes between high and low immune or stromal groups to calculate immune or stromal score-related DEGs. The DEGs of mRNA-seq were calculated from the count data downloaded from TCGA by LIMMA voom package. For microarray data, the normalized data were used to calculate the DEGs (log FC > 1 or < -1, adj. P < 0.05) using the LIMMA package [18]. The up-regulated and down-regulated DEGs of the immune index were intersected separately to obtain two Wayn diagrams. The high expression and low expression differential genes of the stromal index were intersected separately to obtain two Wayn diagrams.
2.3. Correlation between the immune stromal index and clinical parameters
We used an independent sample t-test and ANOVA in Graphpad software to analyze the correlation between the immune stromal index calculated by TCGA and the clinical parameters of UVM. If P < 0.05, the difference was considered statistically significant.
2.4. Functional analysis of the immune index intersection genes
In order to understand the enrichment of the intersection genes in biological processes, cellular components, and molecular functions and pathways, GO and KEGG analyses were performed by using the cluster Profiler package of R software for the intersection genes of the immune index. Protein ANalysis THrough Evolutionary Relationships (PANTHER) database was also searched for pathway annotation of the intersected genes.
P < 0.05 was considered to be a statistically significant enrichment entry. We also use the STRING website to construct a protein-protein interaction (PPI) network map of these intersecting genes to analyze the interactions between these genes.
2.5. Establishment of an immune-related prognosis model
Among 80 patients with UVM, 58 patients with survival days greater than 90 were selected for a follow-up survival analysis. We excluded patients with survival time less than 90 days because of the possibility of insufficient follow-up time or non-cancer-related deaths [19]. For the intersection genes of the immune index, firstly, the survival package was used in R software to conduct a univariate Cox regression analysis, and the genes with significant prognostic significance were selected (P < 0.01). Next, we performed a multivariate Cox regression analysis on these selected genes using the SPSS software, then selected the genes with P < 0.05 to incorporate into the prognostic model as characteristic genes. The selected characteristic genes were used to construct the prognosis model with the following formula:
where log2 (TPM + 0.001) and beta represent the standardized expression value of characteristic genes and the multivariate Cox regression coefficient, respectively [20].
2.6. Evaluation and verification of prognostic models
In the constructed prognosis model, the patients were divided into high- and low-risk groups with the median risk score as the cutoff value. The survival package, survminer package, and survival receiver's operating characteristics curves (ROC) package of R software were used to draw the Kaplan-Meier (KM) curve and time-dependent receiver's operating characteristics (tROC) curve to evaluate the three-year and five-year prognostic prediction capabilities of the model. If P < 0.05 for the KM curve and area under curve (AUC) > 0.8 for the tROC curve, this model appears to distinguish survival outcomes and predict prognosis well. To evaluate the independent prognostic value of this model compared with other clinical features, we included the clinical parameters of UVM and this model using SPSS for univariate and multivariate COX regression analyses. P < 0.05 was considered to indicate independent prognostic ability. To validate this model, we calculated the expression values of the model's signature genes in the cancer and non-cancer groups in GSE62075.
3.
Results
3.2. Immune stromal index and DEGs
The immune and stromal scores of each sample from GSE22138 and RNA-seq dataset were listed in Supplementary File 1. As shown in Figure 1, there were 1828 and 199 genes of mRNA-seq and GSE22138 with an up-regulated immune index, respectively, and 164 genes in the intersection, and there were 662 and 80 genes of mRNA-seq and GSE22138 with a down-regulated immune index, respectively, and 36 genes in the intersection (Supplementary Tables 1 and 2). There were 2152 and 136 genes with mRNA-seq and GSE22138 stromal index up-regulated genes, respectively, with 117 genes in the intersection, and there were 317 and 12 genes with mRNA-seq and GSE22138 stromal index down-regulated genes, respectively, and three genes with an intersection. All intersected immune or stromal-related DEGs were exhibited in Supplementary File 2. There were 98 common genes in the intersection of the immune index and stromal index up-regulated genes of mRNA-seq and GSE22138 and three common genes in the intersection of the immune index and stromal index down-regulated genes of mRNA-seq and GSE22138. Because the mRNA-seq and GSE22138 immune index up-regulated genes overlapped with most of the stromal index up-regulated genes, 164 genes in the intersection of the up-regulated genes in the immune index of mRNA-seq and GSE22138 were selected for subsequent prognosis analysis. The heat map of these 164 genes in the high- and low-immune index groups of mRNA-seq and GSE22138 is shown in Figure 2.
3.2. Correlation between the immune stromal index and clinical parameters
There was no correlation between the mRNA-seq immune/stromal index and clinical parameters (P < 0.05; the results are not shown) (Supplementary File 3).
3.3. Functional analysis of the immune index intersection genes
According to Table 1, these 164 genes were mainly enriched in biological processes, such as response to interferon-gamma and cellular response to interferon-gamma and the interferon-gamma-mediated signaling pathway; molecular components, such as the MHC protein complex and the side of membrane and the integral component of the lumenal side of the endoplasmic reticulum membrane; and molecular functions, such as peptide antigen binding and peptide and amide binding(P < 0.05).There was no significant enrichment of the KEGG pathway (Figure 3). Annotations from PANTHER database indicated that the 164 genes were assembled in 5-Hydroxytryptamine biosynthesis and nicotinic acetylcholine receptor signaling pathway (Table 2). The PPI diagram in Figure 3 showed the protein-level interactions of these genes.
3.4. Establishment of the immune-related prognosis model
A total of 58 patients with survival days over 90 in the mRNA sequence data met the requirements for survival analysis. After univariate Cox regression analysis of these 164 genes, there were 30 genes with significant prognostic value (P < 0.01) (Figure 4). Multivariate Cox regression of these 30 genes revealed that four genes had considerable prognostic significance (P < 0.05) (Figure 5). Therefore, using these four genes to construct a prognostic model, the formula was the expression of HLA-J * 0.764 + the expression of MMP12 * 0.552 + the expression of HES6 * 2.755 + the expression of ADAMDEC1 * (-0.505). HLA-J, MMP12 and HES6 were high risk genes and ADAMDEC1 was low risk gene. The coefficients for the univariate and multivariate Cox regression of these four genes are shown in Table 3.
3.5. Evaluation and verification of the prognostic model
The KM curve showed that the prognosis model had a considerable ability to predict differences in survival status between the high- and low-risk groups (P = 0.002) (Figure 6). The AUC values of the tROC curves of the three-year survival curve and the five-year survival curve reached 0.893 and 0.957, respectively, indicating the excellent survival prediction ability of this model (Figure 7). Because there were no uveal melanoma chips containing survival data, we validated the expression values of several characteristic genes in GSE44295 and GSE62075. The KM curve showed that the prognosis model had a significant ability to distinguish between high- and low-risk groups in terms of survival status (P = 0.002).
The samples in GSE62075 were divided into two groups, the cancer group and the non-cancer group. In GSE62075, the cancer group was the cells from humal uveal melanoma derivative cell lines; and the non-cancer group was the primary human uveal melanocytes. The results showed that the expression value of ADAMDEC1 in the non-cancer group of GSE62075 was significantly higher than that in the cancer group (P = 0.0004), consistent with its protective prognostic significance in UVM (Figure 8).
4.
Discussion
Current research on UVM focuses on the prognosis, detection, and treatment of metastatic UVM. As microarray sequencing technology continues to evolve, researchers are beginning to explore new markers and targets for clinical trials. There is no doubt that the TME affects the response to treatment. Changing its composition or function can improve the treatment effect. Therefore, better understanding of TME of UVM can provide new ideas for treatment. In this study, we attempted to analyze the TME of UVM, to find DEGs with crucial prognostic value, and to construct a prognostic model to explore the potential molecular mechanism of the immune index in UVM.
This study is the first to use the immune index and the stromal index for UVM. A total of 164 intersection genes of the immune index up-regulated genes obtained by GSE22138 and TCGA were analyzed. GO analysis revealed that these genes mainly enriched in the biological processes of response to interferon-gamma, cellular response to interferon-gamma and the interferon-gamma-mediated signaling pathway. Interferon-gamma can regulate immune response in many ways. Increased serum IFN-γ levels are associated with the spread of metastasis and poor prognosis in patients with UVM [21]. The study found that high concentrations of IFN-γ (1000 units / mL) can induce UVM cell death, and 1 micromol / L decitabine works synergistically with 10 to 1000 units / mL IFN-γ to induce massive UVM cell death [22].
By analyzing these 164 intersection genes, a prognosis model including HLA-J, MMP12, HES6, and ADAMDEC1 was constructed via univariate and multivariate COX regression. The remarkable change of HR values for HES6 and ADAMDEC1 might be explained by Simpson's Paradox, where two groups of data under certain conditions will meet certain properties when discussed separately, but once considered together, it may lead to the opposite conclusion. In univariate Cox regression analysis, "Simpson's paradox" will occur under certain circumstances because the influence of other factors is not considered. However, in multivariate Cox regression analysis, the mystery of "Simpson's paradox" can be solved by adjusting and controlling the influence of other factors. That's why we performed multivariate Cox regression analysis in addition to univariate Cox regression analysis. The large interval of HES6 of the hazard ratio of HES6 might be caused by the imbalance between binary variables, which might be addressed by increasing the number of samples in future work. In addition, the model demonstrated excellent prognostic and predictive power in this study. Before the work of the present study, some histopathological and molecular features have been proposed as prognostic indicators for UVM. Histopathological features including cell types of spindle cell nevus, spindle cell malignant melanoma, mixed cell melanoma, epithelioid cell melanoma, location of anterior place, diffuse growth pattern, mitotic figures and pigmentation were associated with poorer prognosis [23]. Molecular signatures such as miRNA signature, alternative splicing events and 15 gene expression profiling (GEP) assay have also been invented for risk stratification of UVM patients [24,25,26]. The prognostic model in the current work was constructed from the novel perspective of immune and stromal index and could compete with the above mentioned histopathological and molecular data with an AUC of 0.957. The combined detection of known histopathological features, molecular data and immune or stromal events discovered in the present work was anticipated to improve the prognostic prediction of UVM. The four genes in the prognostic model constructed in this study have been reported to be involved in vital mechanisms of human disease.
Genome-wide association studies have revealed that HLA-J is adifferential gene methylation site in patients with polycystic ovary syndrome and an allele associated with insulin resistance in childhood obesity [27,28].
MMP12 is involved in the degradation of the extracellular matrix. It inhibits the angiogenesis and metastasis of lung tumors [29]. MMP12 has been related to tumor development and invasion. High levels of MMP12 expression lead to poor prognosis in a variety of cancers, including colon cancer, hepatocellular carcinoma, and esophageal squamous cell carcinoma [30,31,32]. MMP12 has also been reported as a potential therapeutic target for lung adenocarcinoma and esophageal squamous cell carcinoma [31,33].
HES6 is linked to cell proliferation and differentiation. It can promote proliferation and cell movement. HES6 has also been relevant to the malignant progression of tumors [34], and its high expression causes poor prognosis in hepatocellular carcinoma and colorectal cancer [35,36].
ADAMDEC1 takes part in the metabolism of the extracellular matrix, such as protein and the process of fiber formation. Studies have confirmed that ADAMDEC1 is associated with various diseases, such as Crohn's disease, myocardial infarction and osteoarthritis.
Although this study developed a prognostic model for UVM based on the immune index, there are still some limitations. First, because there was no UVM chips containing survival data, we had to verify the differential expression of some characteristic genes in the cancer and non-cancer groups in only two GSE chips. In future studies, we will collect enough samples of UVM with sufficient prognosis data to verify the prognosis model of this study. Second, experimental studies on the function of prognostic features in UVM have been insufficient. Hereafter, we must further prove the mechanism of the characteristic model in the genesis and development of UVM through in vivo and in vitro experiments. Third, the expression of markers in prognostic model has not been validated by RT-qPCR, immunohistochemistry or Western blot experiments, which should be addressed in future work.
5.
Conclusions
The risk signature based on the immune index by the ESTIMATE algorithm can be used as an effective marker to evaluate the prognosis of patients with UVM. This study also provides valuable clues for the mechanism of the immune index in UVM.
Acknowledgments
The authors thank The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) for providing data.
Conflict of interest
The authors declare there is no conflict of interest.
DownLoad:
