A novel immune-related genes prognosis biomarker for melanoma: associated with tumor microenvironment

Rongzhi Huang; Min Mao; Yunxin Lu; Qingliang Yu; Liang Liao

doi:10.18632/aging.103054

Featured Nobel Articles

Elizabeth Blackburn, a member of the Editorial Board of Aging, won the Nobel Prize in Physiology or Medicine 2009, while being a member of the board. Elizabeth Blackburn co-authored a paper published in the first (inaugural) issue of Aging.
Andrew V. Schally, Nobel Prize Laureate, published his paper in Aging.
Shinya Yamanaka won the Nobel Prize in Physiology and Medicine 2012. Shinya Yamanaka co-authored a paper published in Aging.

Learn about our FREE Post-Publication Promotion Services

Research Paper Volume 12, Issue 8 pp 6966—6980

A novel immune-related genes prognosis biomarker for melanoma: associated with tumor microenvironment

Rongzhi Huang^{1,
,} , Min Mao^{1,
,} , Yunxin Lu^{1,
,} , Qingliang Yu^{1,
,} , Liang Liao^1,2, ,

¹ The First Affiliated Hospital of Guangxi Medical University, Nanning 530021, The Guangxi Zhuang Autonomous Region, China
² Department of Traumatic Orthopedics and Hand Surgery, The First Affiliated Hospital of Guangxi Medical University, Nanning 530021, The Guangxi Zhuang Autonomous Region, China

* Equal contribution

Received: August 7, 2019 Accepted: March 29, 2020 Published: April 20, 2020

https://doi.org/10.18632/aging.103054
How to Cite

Copyright © 2020 Huang 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.

Abstract

Background: Melanoma is a cancer of the skin with potential to spread to other organs and is responsible for most deaths due to skin cancer. It is imperative to identify immune biomarkers for early melanoma diagnosis and treatment.

Results: 63 immune-related genes of the total 1039 unique IRGs retrieved were associated with overall survival of melanoma. A multi-IRGs classifier constructed using eight IRGs showed a powerful predictive ability. The classifier had better predictive power compared with the current clinical data. GSEA analysis showed multiple signaling differences between high and low risk score group. Furthermore, biomarker was associated with multiple immune cells and immune infiltration in tumor microenvironment.

Conclusions: The immune-related genes prognosis biomarker is an effective potential prognostic classifier in the immunotherapies and surveillance of melanoma.

Methods: Melanoma samples of genes were retrieved from TCGA and GEO databases while the immune-related genes (IRGs) were retrieved from the ImmPort database. WGCNA, Cox regression analysis and LASSO analysis were used to classify melanoma prognosis. ESTIMATE and CIBERSORT algorithms were used to explore the relationship between risk score and tumor immune microenvironment. GSEA analysis was performed to explore the biological signaling pathway.

Introduction

Melanoma is a life-threatening malignancy with high metastasis and mortality rates [1, 2]. Approximately 232,000 new melanoma patients were diagnosed in 2011 and with 55,000 deaths recorded in the same year [3]. High mortality rates result from poor prognosis leading to late diagnosis. Therefore, there is need to come up with approaches for early diagnosis [4–6].

The TNM stage is an effective approach for detection of the cancer stage, is invaluable in cancer prognosis and informs on the right therapy approaches [7]. However, differences in the overall survival associated with TNM stage method are observed [8]. Current studies on tumors have revealed the clinical limitations of TNM stage method [9, 10]. Therefore, there is a need to explore new melanoma markers to guide the clinical treatment and improve melanoma prognosis. Gene-based biomarkers have become more popular with the advances in human gene sequencing [11, 12].

Most immune system components are implicated in the initiation and progression of melanoma [13, 14]. In tumor immunity, tumor cells act as antigens while immune cells and leukocytes infiltrates the tumor tissue function through chemotaxis for immune defense [13]. Immune escape also is an important factor in tumorigenesis [15, 16]. Currently, a myriad of new immunotherapy are used in melanoma and including PD-1, PD-L1 and CTLA-4 inhibitors [17, 18]. However, these approaches are effective only on a few patients while the majority of the patients have limited or no response to the therapy especially during melanoma progression [19, 20]. Therefore, comprehensive analyses of the correlation between immune genes and overall survival in melanoma are important in exploring the potential prognostic value of immune genes and new biomarkers.

In this study, our aim was to construct a novel immune-related genes biomarker for use in immunotherapies and melanoma prognosis. Comprehensive bioinformatics analyses were performed to explore underlying mechanisms of the biomarker. This study provides information for subsequent personalized diagnosis and treatment of melanoma.

Results

Identification of survival-related modules by WGCNA

WGCNA analysis was carried out on 950 overlapping IRGs (Figure 1). The soft-thresholding power in WGCNA was determined based on a scale-free R2 (R2 = 0.95). Six modules were identified based on the average linkage hierarchical clustering and the soft-thresholding power. The red module showed the highest correlation with the overall survival of melanoma. Additionally, the blue module was highly correlated with the overall survival of melanoma. The red module contained 22 IRGs while the blue module contained 138 IRGs (Figure 2). Data for these two modules were selected for further analysis.

Venn diagram and Histogram was used to visualize common IRGs shared between GEO dataset, TCGA dataset and IRGS. 950 IRGs overlapped in the three datasets. The value used represented the number of gene symbol covered from the ensemble IDs and probe IDs. The number of genes annotated are presented on the y-axis.

Figure 1. Venn diagram and Histogram was used to visualize common IRGs shared between GEO dataset, TCGA dataset and IRGS. 950 IRGs overlapped in the three datasets. The value used represented the number of gene symbol covered from the ensemble IDs and probe IDs. The number of genes annotated are presented on the y-axis.

Weighted melanoma gene co-expression network. (A) The scale-free fit index for soft-thresholding powers. The soft-thresholding power in the WGCNA was determined based on a scale-free R2 (R2 = 0.95). The left panel presents the relationship between the soft-threshold and scale-free R2. The right panel presents the relationship between the soft-threshold and mean connectivity. (B) A dendrogram of the differentially expressed genes clustered based on different metrics. Each branch in the figure represents one gene, and every color below represents one co-expression module. (C) Distribution of average gene significance and errors in the modules associated with overall survival of melanoma patients. Based on the average linkage hierarchical clustering and the soft-thresholding power, six modules were identified. To determine the significance of each module, gene significance (GS) was calculated to measure the correlation between genes and sample traits. GS was defined as the log10 conversion of the p-value in the linear regression between gene expression and clinical data (GS = lg P). The red and blue module showed high correlation with the survival of melanoma patients. (D) A heatmap showing the correlation between the gene module and clinical traits. The red module contained 22 IRGs while the blue module contained 138 IRGs. The correlation coefficient in each cell represented the correlation between gene module and the clinical traits, which decreased in size from red to blue. The blue module showed the highest positive correlation with the survival while the red module showed the highest negative correlation with the survival.

Figure 2. Weighted melanoma gene co-expression network. (A) The scale-free fit index for soft-thresholding powers. The soft-thresholding power in the WGCNA was determined based on a scale-free R2 (R2 = 0.95). The left panel presents the relationship between the soft-threshold and scale-free R2. The right panel presents the relationship between the soft-threshold and mean connectivity. (B) A dendrogram of the differentially expressed genes clustered based on different metrics. Each branch in the figure represents one gene, and every color below represents one co-expression module. (C) Distribution of average gene significance and errors in the modules associated with overall survival of melanoma patients. Based on the average linkage hierarchical clustering and the soft-thresholding power, six modules were identified. To determine the significance of each module, gene significance (GS) was calculated to measure the correlation between genes and sample traits. GS was defined as the log10 conversion of the p-value in the linear regression between gene expression and clinical data (GS = lg P). The red and blue module showed high correlation with the survival of melanoma patients. (D) A heatmap showing the correlation between the gene module and clinical traits. The red module contained 22 IRGs while the blue module contained 138 IRGs. The correlation coefficient in each cell represented the correlation between gene module and the clinical traits, which decreased in size from red to blue. The blue module showed the highest positive correlation with the survival while the red module showed the highest negative correlation with the survival.

Construction of prognostic classifier based on IRGs

63 IRGs of the red and blue modules were identified as survival related IRGs of melanoma with the criterion of P < 0.01 (Supplementary File 1). LASSO analysis identified eight IRGs (PSME1, CDC42, CMTM6, HLA-DQB1, HLA-C, CXCR6, CD8B, TNFSF13) which were included in the classifier (Figure 3). The coefficients of the eight IRGs are shown in Table 1 and the expression levels are shown in Figure 4. The high-RS group showed a poor overall survival rate compared with low-RS group based the Kaplan-Meier analysis (Figure 5B). Time-dependent ROC curves showed that the classifier had a strong predictive ability in GSE dataset (Figure 5A). In the training cohort, the AUC was 0.679 in 1 year, 0.743 in 3 years and 0.740 in 5 years (Figure 5A).

Construction of the IRGs prognostic classifier. (A, B) Determination of the number of factors by the LASSO analysis. (C) The distribution of RS. (D) The survival duration and status of patients. (E) A heatmap of IRGs in the classifier.

Figure 3. Construction of the IRGs prognostic classifier. (A, B) Determination of the number of factors by the LASSO analysis. (C) The distribution of RS. (D) The survival duration and status of patients. (E) A heatmap of IRGs in the classifier.

Table 1. The IRGs in the prognostic classifier associated with OS in the GSE dataset.

Symbol	Univariate Cox regression analysis			LASSO coefficient
Symbol	HR	95%CI	P Value	LASSO coefficient
PSME1	0.416	0.285-0.608	5.854205e-06	-0.30396287
CDC42	0.428	0.248-0.74	0.00236537	-0.24399092
CMTM6	0.364	0.218-0.608	0.0001131757	-0.23548175
HLA-DQB1	0.692	0.592-0.809	3.711835e-06	-0.07311844
HLA-C	0.595	0.466-0.759	2.920363e-05	-0.10691953
CXCR6	0.509	0.363-0.713	8.635839e-05	-0.03143482
CD8B	0.248	0.108-0.566	0.0009273984	-0.05032655
TNFSF13	0.172	0.055-0.54	0.002576346	-0.25872281

Figure 4. Expression profile of 8 genes. (A) GSE dataset (B) TCGA dataset.

The distribution of time-dependent ROC curves and Kaplan-Meier survival based on the integrated classifier in the training and independent validation sets. ROC, receiver operator characteristic. AUC, the area under the curve. (A) ROC curve for the GSE cohort. (B) KM curve of the GSE cohort. (C) ROC curve of the TCGA cohort. (D) KM curve of the TCGA cohort. (E) 3-years correlation ROC curve in the TCGA cohort for the comparison of the classifier prognostic accuracy and clinical characteristics. (F) 5-years correlation ROC curve in the TCGA cohort for the comparison of the classifier prognostic accuracy and clinical characteristics.

Figure 5. The distribution of time-dependent ROC curves and Kaplan-Meier survival based on the integrated classifier in the training and independent validation sets. ROC, receiver operator characteristic. AUC, the area under the curve. (A) ROC curve for the GSE cohort. (B) KM curve of the GSE cohort. (C) ROC curve of the TCGA cohort. (D) KM curve of the TCGA cohort. (E) 3-years correlation ROC curve in the TCGA cohort for the comparison of the classifier prognostic accuracy and clinical characteristics. (F) 5-years correlation ROC curve in the TCGA cohort for the comparison of the classifier prognostic accuracy and clinical characteristics.

Verification of the prognostic classifier in TCGA cohort

We used the TGCA cohort to validate the predictive ability of the classifier. Kaplan-Meier analysis showed that the high-RS group had a poor overall survival (P<0.0001, Figure 5D). Time-dependent ROC curves showed that the classifier had a good accuracy with 0.642 in 1 year, 0.636 in 3 years and 0.645 in 5 years (Figure 5C). Moreover, the classifier had better predictive power and accuracy compared with other clinical features (Figure 5E, 5F). In Addition, the classifier was an independent factor in multivariate Cox analysis. Results of univariate and multivariate analyses in prognostic factors and overall survival were showed in Table 2.

Table 2. Univariate and multivariate analyses of prognostic factors and overall survival of melanoma patients in TCGA cohort.

Characteristics	Univariate Cox regression analysis			multivariate Cox regression analysis
Characteristics	HR	95%CI	P Value	HR	95%CI	P Value
Age	1.025	1.015-1.035	3.63e-07	1.021	0.01-4.06	4.83e-05
Gender	0.877	0.655-1.175	3.79e-01	1.088	0.16-0.54	0.592
Local invasion	0.988	0.955-1.021	4.65e-01	0.987	0.02--0.66	0.511
Lymph node metastasis	1.087	1.032-1.145	1.74e-03	1.092	0.03-3.3	0.00096
Distant metastasis	1.161	0.887-1.52	2.78e-01	1.429	0.14-2.55	0.0107
TNM stage	1.000	0.964-1.038	9.80e-01	0.982	0.02--0.75	0.455
Multi-IRGs Classify	1.588	1.315-1.919	1.61e-06	1.704	0.1-5.13	2.94e-07

Immune infiltration score between high and low RS group

Kaplan-Meier analysis showed that different immune scores had differential overall survival in melanoma samples (Figure 6A, 6B). The immune score showed a significant difference between high and low-RS group (Figure 6C, 6D).

(A) Impact of immune score on overall survival in melanoma based on KM analysis. (A) GSE cohort. (B) TCGA cohort. (C, D) Association with immune score, stromal score and risk score. The high-RS group showed lower immune score and stromal score comparing with low-RS group. (C) GSE cohort. (D) TCGA cohort.

Figure 6. (A) Impact of immune score on overall survival in melanoma based on KM analysis. (A) GSE cohort. (B) TCGA cohort. (C, D) Association with immune score, stromal score and risk score. The high-RS group showed lower immune score and stromal score comparing with low-RS group. (C) GSE cohort. (D) TCGA cohort.

Immune cell subtypes between high and low RS group

The 22 immune cell proportions of melanoma are shown in Figure 7A, 7B. Macrophages M0, Macrophages M2 and T cells CD8 accounted for a large proportion of melanoma immune cell infiltration. High and low RS groups showed differential immune cells expression (Figure 7C, 7D).

(A, B) The mean proportion of 22 immune cells in GSE cohort. Macrophages M0, Macrophages M2 and T cells CD8 account for a large proportion of melanoma immune cell infiltration. (A) GSE cohort. (B) TCGA cohort. (C, D) Violin plot showing the relationship between risk score with immune score and stromal score. Red color represents high-RS group while blue color represents low-RS group. Differential immune cell type expression was observed between the high and low-RS groups. (C) GSE cohort. (D) TCGA cohort.

Figure 7. (A, B) The mean proportion of 22 immune cells in GSE cohort. Macrophages M0, Macrophages M2 and T cells CD8 account for a large proportion of melanoma immune cell infiltration. (A) GSE cohort. (B) TCGA cohort. (C, D) Violin plot showing the relationship between risk score with immune score and stromal score. Red color represents high-RS group while blue color represents low-RS group. Differential immune cell type expression was observed between the high and low-RS groups. (C) GSE cohort. (D) TCGA cohort.

GSEA analysis

GSEA analysis showed 14 significant KEGG pathways associated with risk score, including Rap1 signaling pathway, Ras signaling pathway, Herpes simplex virus 1 infection, Regulation of actin cytoskeleton, MAPK signaling pathway, Neuroactive ligand-receptor interaction, Human cytomegalovirus infection, Human T-cell leukemia virus 1 infection, Human immunodeficiency virus 1 infection, Kaposi sarcoma-associated herpesvirus infection, Chemokine signaling pathway, Epstein-Barr virus infection, Tuberculosis and Cytokine-cytokine receptor interaction (Figure 8).

Figure 8. GSEA analysis.

Discussion

Melanoma is a fatal skin cancer that affects many people worldwide each year [21]. Currently, immunotherapy is a successful treatment option for melanoma [22]. Notably, many researchers demonstrates the role of the immune cells on tumor cells [23, 24]. Moreover, immune components in melanoma tissue can be used to evaluate therapeutic efficacy and melanoma prognosis in patients [25]. In this study, 63 IRGs were found to be associated with melanoma prognosis, of which eight IRGs were adopted to construct a classifier. The classifier showed reliable predictive value and accuracy. In addition, we explored the relationship between RS and the prognosis value in melanoma. The findings showed differences in immune cell infiltration and multiple signaling pathways between high and low-RS group.

The PSME1, CDC42, CMTM6, HLA-DQB1, HLA-C, CXCR6, CD8B and TNFSF13 RGs were used in the classifier. These IRGs were reported to be associated with tumor prognosis in previous studies. Cell division cycle 42 (CDC42) protein, a member of Rho GTPases, activates multiple cellular processes by regulating actin cytoskeleton [26]. In addition, CDC42 facilitates the invasion and migration of melanoma cells [27–29]. Therefore, CDC42 inhibitors have been effective in melanoma treatment [30, 31]. CMTM6 is a ubiquitously expressed protein encoded by two distinct gene clusters located on chromosome 16 and chromosome 3 [32]. It enhances PD-L1 expression and anti-tumor immunity. Therefore, CMTM6 is a potential biomarker and therapeutic target for melanoma patients [33–35]. Among the HLA class I antigens, HLA-C locus recognizes the inhibitory killer cells and suppresses the functions of NK cells in melanoma patients [36–38]. Furthermore, the frequency of HLA-DQB1*0301 and HLA-DQB1*0303 alleles are highly expressed in melanoma patients [39]. Moreover, melanoma patients with DQBI*0301 allele have thicker primary tumor and are more likely to have local or distant metastatic disease [40]. Besides, the chemokine co-receptor CXCR6 was identified as a new biomarker associated with asymmetric self-renewal of tissue-specific stem cells. CXCR6 + cells cause rapid increase in tumor mass compared with CXCR6- cells [41]. TNFSF13, a member of the TNF superfamily, was reported to indicate the proliferative or survival state in tumor cells [42]. The multi-IRGs classifier established in this study showed high predictive value and accuracy through various analyses.

The degree of immune infiltration significantly affected melanoma survival. Previous studies demonstrate that immune cells in the tumor microenvironment can be used in the prognostic assessment of multiple tumors, such as glioblastoma, breast cancer, and melanoma [43–45]. In this study, the expression of eight genes affected immune infiltration scores. Patients with higher immune scores had better prognosis. This finding implies that prognosis value of risk score is associated with melanoma immune system.

To further explore the immune and risk score, we used the CIBERSORT algorithms to calculate the immune cell subtype in R platform. Our result showed that the two risk score groups expressed differential immune cell subtypes. Ali et al. demonstrated that imbalance in immune cell component ratio is highly correlated with poor prognosis and low survival in cancer patients [46, 47]. A previous study reported that CD8⁺ T cells produces granulocyte and perforin to kill tumor cells [48]. In our study, the immune cells found in melanoma mainly comprised macrophages M0, macrophages M2 and T cells CD8. In this study, T cell CD8 levels were low whereas M0 and M2 macrophages levels were high in the high-risk group. This implies that imbalance of T cell CD8 and M0, M2 macrophage ratio may reduce the survival rate of patients in the high-risk group. High expression of CD8+T cells may improve the prognosis of melanoma patients as well as reduce the risk factors.

GSEA analysis showed differences in 14 important signaling pathways between high and low RS groups. Inhibition of MAPK signaling pathway improved melanoma immune microenvironment by enhancing the melanoma antigen expression and down-regulating immunosuppressive cytokines [49, 50]. Additionally, chemokine signaling pathway participates in tumor growth. Some chemokines, such as CCR10 and CXCR3, have been shown to play an important role in the proliferation and metastasis of melanoma cells [51, 52].

In this study, LASSO regression analysis was used to establish a novel classifier using multiple IRGs and the classifier was verified using an independent cohort. Currently, few studies have used ESTIMATE and CIBERSORT algorithms to explore immune infiltration in melanoma. In this study, we use these algorithms to explore immune infiltration in melanoma using the R software. These preliminary results could provide a perspective for exploring the role of immune infiltration in melanoma. However, this study has the following limitations. First, the reliability of our molecular mechanism analysis results is limited due to lack of vitro or vivo experiments. Second, this study was a retrospective study, therefore, prospective study should be carried out to validate the findings of our study.

In conclusion, we successfully constructed a multi-IRGs classifier with the powerful predictive function. Differences in the overall survival of high and low risk groups are implicated in immune infiltration, tumor microenvironment and the interaction of multiple signaling pathways. This study provides additional information on the analysis of melanoma pathogenesis and clinical treatment.

Materials and Methods

Data Procession

GSE65904 gene expression profiles were retrieved from the gene expression omnibus database (GEO: https://www.ncbi.nlm.nih.gov/geo/). In this study, the samples with no follow-up information or follow-up time less than 1 day were excluded. 210 melanoma samples were retrieved for subsequent analysis. Further, the probe IDs were converted to gene symbols using the illuminaHumanv4.db R package. The probe IDs with the highest mean value were reversed when more than one probe had a matched gene symbol. The GEO expression file was converted into log2 (expression) for further analysis. Additionally, the RNA-FPKM data and clinical data of melanoma samples were retrieved for external validation analysis using the TCGA biolinks R package. Samples with no follow-up information or follow-up time less than 1 day were excluded. The expression file of patients with the highest mean value was reversed when more than one expression file had matched patients. 428 melanoma samples in TCGA were used for analysis.

Immune-related gene extraction

Immune-related genes (IRGs) data were retrieved from the ImmPort database (https://immport.niaid.nih.gov) (Supplementary File 2). Overlapping immune-related genes from the GEO dataset, TCGA dataset and IRGs were selected for further analysis.

Weighted gene co-expression network analysis

GEO expression file was used for weighted gene co-expression network analysis (WGCNA) using WGCNAR package. WGCNA was used to explore the relationship between the clinical features with expression modules. Module eigengenes (MEs) were defined as the first principal component of each gene module and adopted as the representative of all genes in each module. Gene significance (GS), as the mediator p-value (GS = lg P) for each gene, represented the degree of linear correlation between gene expression of the module and clinical features. Survival-related modules were defined according to P≤0.01 and the higher GS value was extracted for further analysis.

LASSO analysis

Univariate Cox regression analysis was performed to explore the impact of each gene on overall survival. The IRGs of survival-related modules with P<0.01 were identified as survival-related IRGs and integrated into the Least Absolute Shrinkage and Selection Operator (LASSO) regression for identification of prognostic risk signatures. The risk score (RS) of each sample was calculated using the formula: risk score = Σexpgenei* βi.

The Kaplan-Meier curve analysis was further conducted to evaluate the relationship between the risk score and overall survival. The median value was used as the cut-off. Univariate and multivariate Cox regression analysis were performed to study the relationship between the index and the clinical features. To validate the accuracy and predictive ability of the signature, it was included in the TCGA dataset. The area under the curve (AUC) of the ROC curve was calculated and compared to examine the classifier performance using time ROC R package.

Comparison of the degree of immune cell infiltration between high and low RS groups

To explore the relationship between risk score and melanoma prognosis, we analyzed the relationship between risk score and tumor microenvironment. The tumor microenvironment comprises a variety of cell types, including immune cells, mesenchymal cells, endothelial cells, inflammatory mediators, and extracellular matrix (ECM) molecules [53]. We used the ESTIMATE algorithm to determine the immune score of each sample using R software and further compared the difference in degree of immune cell infiltration between high and low-risk groups by Wilcoxon test.

Comparison of 22 immune cell subtypes between high RS and low RS groups

To explore the differences of immune cell subtypes, CIBERSORT package was used to assess the proportions of 22 immune cell subtypes based on expression file [54]. The perm was set at 1000. Samples with P < 0.05 in CIBERSORT analysis result were used in further analysis. Mann-Whitney U test was used to compare differences in immune cell subtypes in the high RS and low RS groups.

Gene Set Enrichment analysis (GSEA)

To identify signaling pathway that are differentially activated between the high RS and low RS groups, we selected an ordered list of genes through limma R package and conducted Gene Set Enrichment Analysis (GSEA) with adjusted p < 0.05 using the cluster filer R package.

Statistical analysis

All analyses were carried out by R version 3.5.2 and corresponding packages. Kaplan-Meier analysis was further conducted to evaluate the relationship between immune score and overall survival using the survimer R package. The median value was set as the cut-off. The glmnet R package was used for LASSO analysis.

Availability of data and materials

The GSE65904gene expression profiles were retrieved from GEO (https://www.ncbi.nlm.nih.gov/geo). The TCGA data were retrieved from GDC data portal (https://portal.gdc.cancer.gov/). The immune-related genes (IRGs) data were retrieved from the ImmPort database (https://immport.niaid.nih.gov). The R software (https://www.r-project.org/) was used for all statistical analyses.

Supplementary Materials

Supplementary File 1

Supplementary File 2

Abbreviations

IRGs: Immune-related genes; TME: Tumor microenvironment; GEO: Gene Expression Omnibus; TCGA: The cancer genome atlas project; LASSO: Least Absolute Shrinkage and Selection Operator; MEs: Module eigengenes; GS: Gene significance; ROC: Receiver operating characteristic curve; AUC: Area under the curve; RS: Risk score; OS: Overall survival.

Author Contributions

Data curation, Rongzhi Huang and Min Mao; Methodology, Rongzhi Huang; Software, Min Mao and Yunxin Lu; Verification, Yunxin Lu and Qingliang Yu; Visualization, Qingliang Yu and Liang Liao; Writing – original draft, Qingliang Yu and Yunxin Lu; Writing – review and editing, Liang Liao.

Acknowledgments

The authors thank TCGA and GEO for sharing the melanoma data.

Conflicts of Interest

All authors have read and approved submission of the manuscript. There authors declare no conflict of interest in relation to the submission.

Funding

This work was supported by the funding by Guangxi health commission appropriate technology promotion project (S2018064) and Guangxi university young and middle-aged teachers' basic research ability promotion project (2019KY0114).

References

1. Rastrelli M, Tropea S, Rossi CR, Alaibac M. Melanoma: epidemiology, risk factors, pathogenesis, diagnosis and classification. In Vivo. 2014; 28:1005–11. [PubMed]
2. Wei D. A multigene support vector machine predictor for metastasis of cutaneous melanoma. Mol Med Rep. 2018; 17:2907–14. https://doi.org/10.3892/mmr.2017.8219 [PubMed]
3. Ferlay J, Soerjomataram I, Dikshit R, Eser S, Mathers C, Rebelo M, Parkin DM, Forman D, Bray F. Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012. Int J Cancer. 2015; 136:E359–86. https://doi.org/10.1002/ijc.29210 [PubMed]
4. Cho YR, Chiang MP. Epidemiology, staging (new system), and prognosis of cutaneous melanoma. Clin Plast Surg. 2010; 37:47–53. https://doi.org/10.1016/j.cps.2009.07.001 [PubMed]
5. Avilés-Izquierdo JA, Molina-López I, Rodríguez-Lomba E, Marquez-Rodas I, Suarez-Fernandez R, Lazaro-Ochaita P. Who detects melanoma? Impact of detection patterns on characteristics and prognosis of patients with melanoma. J Am Acad Dermatol. 2016; 75:967–74. https://doi.org/10.1016/j.jaad.2016.07.009 [PubMed]
6. Lee CS, Thomas CM, Ng KE. An Overview of the Changing Landscape of Treatment for Advanced Melanoma. Pharmacotherapy. 2017; 37:319–33. https://doi.org/10.1002/phar.1895 [PubMed]
7. Bertero L, Massa F, Metovic J, Zanetti R, Castellano I, Ricardi U, Papotti M, Cassoni P. Eighth Edition of the UICC Classification of Malignant Tumours: an overview of the changes in the pathological TNM classification criteria-What has changed and why? Virchows Arch. 2018; 472:519–31. https://doi.org/10.1007/s00428-017-2276-y [PubMed]
8. Perakis SO, Thomas JE, Pichler M. Non-coding RNAs Enabling Prognostic Stratification and Prediction of Therapeutic Response in Colorectal Cancer Patients. Adv Exp Med Biol. 2016; 937:183–204. https://doi.org/10.1007/978-3-319-42059-2_10 [PubMed]
9. Edge SB, Compton CC. The American Joint Committee on Cancer: the 7th edition of the AJCC cancer staging manual and the future of TNM. Ann Surg Oncol. 2010; 17:1471–4. https://doi.org/10.1245/s10434-010-0985-4 [PubMed]
10. Galon J, Pagès F, Marincola FM, Thurin M, Trinchieri G, Fox BA, Gajewski TF, Ascierto PA. The immune score as a new possible approach for the classification of cancer. J Transl Med. 2012; 10:1. https://doi.org/10.1186/1479-5876-10-1 [PubMed]
11. Gu X, Li B, Jiang M, Fang M, Ji J, Wang A, Wang M, Jiang X, Gao C. RNA sequencing reveals differentially expressed genes as potential diagnostic and prognostic indicators of gallbladder carcinoma. Oncotarget. 2015; 6:20661–71. https://doi.org/10.18632/oncotarget.3861 [PubMed]
12. Liu Y, Jing R, Xu J, Liu K, Xue J, Wen Z, Li M. Comparative analysis of oncogenes identified by microarray and RNA-sequencing as biomarkers for clinical prognosis. Biomark Med. 2015; 9:1067–78. https://doi.org/10.2217/bmm.15.97 [PubMed]
13. Angell H, Galon J. From the immune contexture to the Immunoscore: the role of prognostic and predictive immune markers in cancer. Curr Opin Immunol. 2013; 25:261–67. https://doi.org/10.1016/j.coi.2013.03.004 [PubMed]
14. Buonaguro FM, Pauza CD, Tornesello ML, Hainaut P, Franco R, Tommasino M. Cancer Diagnostic and Predictive Biomarkers 2016. Biomed Res Int. 2017; 2017:7362721. https://doi.org/10.1155/2017/7362721 [PubMed]
15. Liu Y, Cao X. Immunosuppressive cells in tumor immune escape and metastasis. J Mol Med (Berl). 2016; 94:509–22. https://doi.org/10.1007/s00109-015-1376-x [PubMed]
16. McGranahan N, Rosenthal R, Hiley CT, Rowan AJ, Watkins TBK, Wilson GA, Birkbak NJ, Veeriah S, Van Loo P, Herrero J, Swanton C; TRACERx Consortium. Allele-Specific HLA Loss and Immune Escape in Lung Cancer Evolution. Cell. 2017; 171:1259–71.e11. https://doi.org/10.1016/j.cell.2017.10.001 [PubMed]
17. Becht E, Giraldo NA, Dieu-Nosjean MC, Sautès-Fridman C, Fridman WH. Cancer immune contexture and immunotherapy. Curr Opin Immunol. 2016; 39:7–13. https://doi.org/10.1016/j.coi.2015.11.009 [PubMed]
18. Silva IP, Long GV. Systemic therapy in advanced melanoma: integrating targeted therapy and immunotherapy into clinical practice. Curr Opin Oncol. 2017; 29:484–92. https://doi.org/10.1097/CCO.0000000000000405 [PubMed]
19. Braun DA, Burke KP, Van Allen EM. Genomic Approaches to Understanding Response and Resistance to Immunotherapy. Clin Cancer Res. 2016; 22:5642–50. https://doi.org/10.1158/1078-0432.CCR-16-0066 [PubMed]
20. Chen DS, Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature. 2017; 541:321–30. https://doi.org/10.1038/nature21349 [PubMed]
21. Owens B. Melanoma. Nature. 2014; 515:S109. https://doi.org/10.1038/515S109a [PubMed]
22. Valpione S, Campana LG. Immunotherapy for advanced melanoma: future directions. Immunotherapy. 2016; 8:199–209. https://doi.org/10.2217/imt.15.111 [PubMed]
23. Giavina-Bianchi MH, Giavina-Bianchi PF, Festa C. Melanoma: tumor microenvironment and new treatments. An Bras Dermatol. 2017; 92:156–66. https://doi.org/10.1590/abd1806-4841.20176183 [PubMed]
24. Fischer GM, Vashisht Gopal YN, McQuade JL, Peng W, DeBerardinis RJ, Davies MA. Metabolic strategies of melanoma cells: Mechanisms, interactions with the tumor microenvironment, and therapeutic implications. Pigment Cell Melanoma Res. 2018; 31:11–30. https://doi.org/10.1111/pcmr.12661 [PubMed]
25. Ladányi A. Prognostic and predictive significance of immune cells infiltrating cutaneous melanoma. Pigment Cell Melanoma Res. 2015; 28:490–500. https://doi.org/10.1111/pcmr.12371 [PubMed]
26. Hall A. Rho family GTPases. Biochem Soc Trans. 2012; 40:1378–82. https://doi.org/10.1042/BST20120103 [PubMed]
27. Tucci MG, Lucarini G, Brancorsini D, Zizzi A, Pugnaloni A, Giacchetti A, Ricotti G, Biagini G. Involvement of E-cadherin, beta-catenin, Cdc42 and CXCR4 in the progression and prognosis of cutaneous melanoma. Br J Dermatol. 2007; 157:1212–16. https://doi.org/10.1111/j.1365-2133.2007.08246.x [PubMed]
28. Gadea G, Sanz-Moreno V, Self A, Godi A, Marshall CJ. DOCK10-mediated Cdc42 activation is necessary for amoeboid invasion of melanoma cells. Curr Biol. 2008; 18:1456–65. https://doi.org/10.1016/j.cub.2008.08.053 [PubMed]
29. Vaysse A, Fang S, Brossard M, Wei Q, Chen WV, Mohamdi H, Vincent-Fetita L, Margaritte-Jeannin P, Lavielle N, Maubec E, Lathrop M, Avril MF, Amos CI, et al. A comprehensive genome-wide analysis of melanoma Breslow thickness identifies interaction between CDC42 and SCIN genetic variants. Int J Cancer. 2016; 139:2012–20. https://doi.org/10.1002/ijc.30245 [PubMed]
30. Lin CM, Lin YL, Ho SY, Chen PR, Tsai YH, Chung CH, Hwang CH, Tsai NM, Tzou SC, Ke CY, Chang J, Chan YL, Wang YS, et al. The inhibitory effect of 7,7″-dimethoxyagastisflavone on the metastasis of melanoma cells via the suppression of F-actin polymerization. Oncotarget. 2016; 8:60046–59. https://doi.org/10.18632/oncotarget.10960 [PubMed]
31. Maldonado MD, Dharmawardhane S. Targeting Rac and Cdc42 GTPases in Cancer. Cancer Res. 2018; 78:3101–11. https://doi.org/10.1158/0008-5472.CAN-18-0619 [PubMed]
32. Han W, Ding P, Xu M, Wang L, Rui M, Shi S, Liu Y, Zheng Y, Chen Y, Yang T, Ma D. Identification of eight genes encoding chemokine-like factor superfamily members 1-8 (CKLFSF1-8) by in silico cloning and experimental validation. Genomics. 2003; 81:609–17. https://doi.org/10.1016/S0888-7543(03)00095-8 [PubMed]
33. Burr ML, Sparbier CE, Chan YC, Williamson JC, Woods K, Beavis PA, Lam EY, Henderson MA, Bell CC, Stolzenburg S, Gilan O, Bloor S, Noori T, et al. CMTM6 maintains the expression of PD-L1 and regulates anti-tumour immunity. Nature. 2017; 549:101–05. https://doi.org/10.1038/nature23643 [PubMed]
34. Mezzadra R, Sun C, Jae LT, Gomez-Eerland R, de Vries E, Wu W, Logtenberg ME, Slagter M, Rozeman EA, Hofland I, Broeks A, Horlings HM, Wessels LF, et al. Identification of CMTM6 and CMTM4 as PD-L1 protein regulators. Nature. 2017; 549:106–10. https://doi.org/10.1038/nature23669 [PubMed]
35. Mamessier E, Birnbaum DJ, Finetti P, Birnbaum D, Bertucci F. CMTM6 stabilizes PD-L1 expression and refines its prognostic value in tumors. Ann Transl Med. 2018; 6:54. https://doi.org/10.21037/atm.2017.11.26 [PubMed]
36. Campillo JA, Martínez-Escribano JA, Muro M, Moya-Quiles R, Marín LA, Montes-Ares O, Guerra N, Sánchez-Pedreño P, Frías JF, Lozano JA, García-Alonso AM, Alvarez-López MR. HLA class I and class II frequencies in patients with cutaneous malignant melanoma from southeastern Spain: the role of HLA-C in disease prognosis. Immunogenetics. 2006; 57:926–33. https://doi.org/10.1007/s00251-005-0065-2 [PubMed]
37. Konjević G, Mirjacić Martinović K, Jurisić V, Babović N, Spuzić I. Biomarkers of suppressed natural killer (NK) cell function in metastatic melanoma: decreased NKG2D and increased CD158a receptors on CD3-CD16+ NK cells. Biomarkers. 2009; 14:258–70. https://doi.org/10.1080/13547500902814658 [PubMed]
38. Kandilarova SM, Paschen A, Mihaylova A, Ivanova M, Schadendorf D, Naumova E. The Influence of HLA and KIR Genes on Malignant Melanoma Development and Progression. Arch Immunol Ther Exp (Warsz). 2016 (Suppl 1); 64:73–81. https://doi.org/10.1007/s00005-016-0437-3 [PubMed]
39. Bateman AC, Turner SJ, Theaker JM, Howell WM. HLA-DQB1*0303 and *0301 alleles influence susceptibility to and prognosis in cutaneous malignant melanoma in the British Caucasian population. Tissue Antigens. 1998; 52:67–73. https://doi.org/10.1111/j.1399-0039.1998.tb03025.x [PubMed]
40. Lee JE, Reveille JD, Ross MI, Platsoucas CD. HLA-DQB1*0301 association with increased cutaneous melanoma risk. Int J Cancer. 1994; 59:510–13. https://doi.org/10.1002/ijc.2910590413 [PubMed]
41. Taghizadeh R, Noh M, Huh YH, Ciusani E, Sigalotti L, Maio M, Arosio B, Nicotra MR, Natali P, Sherley JL, La Porta CA. CXCR6, a newly defined biomarker of tissue-specific stem cell asymmetric self-renewal, identifies more aggressive human melanoma cancer stem cells. PLoS One. 2010; 5:e15183. https://doi.org/10.1371/journal.pone.0015183 [PubMed]
42. Mhawech-Fauceglia P, Kaya G, Sauter G, McKee T, Donze O, Schwaller J, Huard B. The source of APRIL up-regulation in human solid tumor lesions. J Leukoc Biol. 2006; 80:697–704. https://doi.org/10.1189/jlb.1105655 [PubMed]
43. Fortis SP, Mahaira LG, Anastasopoulou EA, Voutsas IF, Perez SA, Baxevanis CN. Immune profiling of melanoma tumors reflecting aggressiveness in a preclinical model. Cancer Immunol Immunother. 2017; 66:1631–42. https://doi.org/10.1007/s00262-017-2056-1 [PubMed]
44. Priedigkeit N, Watters RJ, Lucas PC, Basudan A, Bhargava R, Horne W, Kolls JK, Fang Z, Rosenzweig MQ, Brufsky AM, Weiss KR, Oesterreich S, Lee AV. Exome-capture RNA sequencing of decade-old breast cancers and matched decalcified bone metastases. JCI Insight. 2017; 2:95703. https://doi.org/10.1172/jci.insight.95703 [PubMed]
45. Jia D, Li S, Li D, Xue H, Yang D, Liu Y. Mining TCGA database for genes of prognostic value in glioblastoma microenvironment. Aging (Albany NY). 2018; 10:592–605. https://doi.org/10.18632/aging.101415 [PubMed]
46. Ali HR, Chlon L, Pharoah PD, Markowetz F, Caldas C. Patterns of Immune Infiltration in Breast Cancer and Their Clinical Implications: A Gene-Expression-Based Retrospective Study. PLoS Med. 2016; 13:e1002194. https://doi.org/10.1371/journal.pmed.1002194 [PubMed]
47. Bense RD, Sotiriou C, Piccart-Gebhart MJ, Haanen JB, van Vugt MA, de Vries EG, Schröder CP, Fehrmann RS. Relevance of Tumor-Infiltrating Immune Cell Composition and Functionality for Disease Outcome in Breast Cancer. J Natl Cancer Inst. 2016; 109:djw192. https://doi.org/10.1093/jnci/djw192 [PubMed]
48. Tsukumo SI, Yasutomo K. Regulation of CD8⁺ T Cells and Antitumor Immunity by Notch Signaling. Front Immunol. 2018; 9:101. https://doi.org/10.3389/fimmu.2018.00101 [PubMed]
49. Deken MA, Gadiot J, Jordanova ES, Lacroix R, van Gool M, Kroon P, Pineda C, Geukes Foppen MH, Scolyer R, Song JY, Verbrugge I, Hoeller C, Dummer R, et al. Targeting the MAPK and PI3K pathways in combination with PD1 blockade in melanoma. Oncoimmunology. 2016; 5:e1238557. https://doi.org/10.1080/2162402X.2016.1238557 [PubMed]
50. Kakavand H, Rawson RV, Pupo GM, Yang JY, Menzies AM, Carlino MS, Kefford RF, Howle JR, Saw RP, Thompson JF, Wilmott JS, Long GV, Scolyer RA, Rizos H. PD-L1 Expression and Immune Escape in Melanoma Resistance to MAPK Inhibitors. Clin Cancer Res. 2017; 23:6054–61. https://doi.org/10.1158/1078-0432.CCR-16-1688 [PubMed]
51. Murakami T, Cardones AR, Finkelstein SE, Restifo NP, Klaunberg BA, Nestle FO, Castillo SS, Dennis PA, Hwang ST. Immune evasion by murine melanoma mediated through CC chemokine receptor-10. J Exp Med. 2003; 198:1337–47. https://doi.org/10.1084/jem.20030593 [PubMed]
52. Jenkins MH, Brinckerhoff CE, Mullins DW. CXCR3 signaling in BRAFWT melanoma increases IL-8 expression and tumorigenicity. PLoS One. 2015; 10:e0121140. https://doi.org/10.1371/journal.pone.0121140 [PubMed]
53. Runa F, Hamalian S, Meade K, Shisgal P, Gray PC, Kelber JA. Tumor microenvironment heterogeneity: challenges and opportunities. Curr Mol Biol Rep. 2017; 3:218–29. https://doi.org/10.1007/s40610-017-0073-7 [PubMed]
54. 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]

Featured Nobel Articles

Research Paper Volume 12, Issue 8 pp 6966—6980