A twelve-gene signature for survival prediction in malignant melanoma patients.

BACKGROUND: Melanoma is defined as a highly mutational heterogeneous disease containing numerous alternations of the molecule. However, due to the phenotypically and genetically heterogeneity of malignant melanoma, conventional clinical characteristics remain restricted or limited in the ability to accurately predict individual outcomes and survival. This study aimed to establish an accurate gene expression signature to predict melanoma prognosis.
METHODS: In this study, we established an RNA sequencing-based 12-gene signature data of melanoma patients obtained from 2 independent databases: the Cancer Genome Atlas (TCGA) database and the Gene Expression Omnibus (GEO) database. We evaluated the quality of each gene to predict survival conditions in each database by employing univariate and multivariate regression models. A prognostic risk score based on a prognostic signature was determined. This prognostic gene signature further classified patients into low-risk and high-risk groups.
RESULTS: Based on a prognostic signature, a prognostic risk score was determined. This prognostic gene signature further divided the patients into low-risk and high-risk groups. In the chemotherapy and radiotherapy groups of the TCGA cohort and V-raf murine sarcoma viral oncogene homolog B1 (BRAF) expression group in the GEO cohort, patients in the low-risk group had a longer survival duration compared to patients in the high-risk group. Nevertheless, the immunotherapy group in the TCGA database and neuroblastoma RAS viral oncogene homolog (NRAS) expression group in the GEO database had no significant differences in statistics. Moreover, this gene signature was associated with patient prognosis regardless of whether the Breslow depth was greater than or less than 3.75 mm. Stratified gene set enrichment analysis (GSEA) revealed that certain immune-related pathways, such as the T-cell signaling pathway, chemokine signaling pathway, and primary immunodeficiency, were significantly enriched in the low-risk group of both TCGA and GEO cohorts. This information implied the immune-related properties of the 12-gene signature.
CONCLUSIONS: Our study emphasizes the significance of the gene expression signature in that it may be an indispensable prognostic predictor in melanoma and has great potential for application in personalized treatment. 2020 Annals of Translational Medicine. All rights reserved.

Introduction

Melanoma is the fifth and sixth most common and aggressive cancer and represents 5% and 4% of cases of cancer in males and females (1). Almost 232,100 cases of primary malignant cancers that were newly diagnosed are cases of cutaneous melanoma (2). Currently, the mainstay treatment of melanoma is surgical resection, comprising of local excision, wide local excision, lymphadenectomy or sentinel lymph node biopsy, followed by radiotherapy, chemotherapy and immunotherapy (3). Although the survival rate of melanoma patients has dramatically improved due to the development in a modern variety of therapeutic methods, cutaneous melanoma still has nearly 55,500 cancer deaths annually (0.7% of all cancer deaths) (2). A systematic review revealed high-risk characteristics of melanoma included multiple melanocytic naevi, family history and types of skin (4). Nowadays, a growing body of research asserts that gene alternations are integrally involved in the disease progression (5), and to some extent might monitor the effective individualized therapies of melanoma patients (6,7).

Melanoma is defined as a highly mutational heterogeneous disease containing numerous alternations of molecule (8). Generally, risk stratification and prognosis of melanoma patients are primarily informed by clinical and pathological features, including Breslow thickness, ulceration, lymphovascular invasion and perineural invasion (9). However, due to the phenotypically and genetically heterogeneity of malignant melanoma, conventional clinical characteristics remain restricted or limited in the ability to predict individual outcomes and survival accurately (10). Therefore, to improve the prognostic and therapeutic effect, different malignant melanoma-specific genomic signatures need to be identified. In many cancers, including melanoma, advances in sequencing and bioinformatics technologies have allowed for genome-wide sequence analysis. The open-access databases, such as the Cancer Genome Atlas (TCGA, https:/portal.gdc.cancer.gov/) and the Gene Expression Omnibus database (GEO, https:/www.ncbi.nlm.gov/geo/), have stored sequencing data from a multitude of melanoma samples. The objective of this study was to identify a melanoma prognostic signature using the TCGA and GEO database RNA sequencing (RNA-seq) patient data. A reliable genomic prognostic signature can supplement conventional clinical prognostic factors and further enable custom therapy.

Methods

Data collection

The raw microarray data and clinical information were acquired from GEO according to the accession number provided (GSE59455). Normalization and log2 transformation were implemented with the manufacture-provided R packages. Level 3 mRNA expression and clinical data were obtained from the TCGA cohort.

Screening for signatures of the prognostic gene

Our analysis charts are shown in . To evaluate the relationships between overall survival and the level of expression of each gene in both cohorts, the proportional hazards regression analysis of univariate Cox analysis was applied. For the correlation analysis, only genes with P values <0.01 were considered candidate genes, and those common genes were applied to both TCGA and GEO cohorts to formulate the predictive model. The relative contribution of each gene to survival prediction in both cohorts was evaluated by step-wise multivariate Cox’s regression model that the common candidate genes were suitable for. The genes that were common to both cohorts and had correlation with survival were finally covered in the prognostic signature. According to the estimated regression coefficients in multivariate Cox’s regression model, a prognostic risk score for predicting overall survival was then calculated as follows,

Figure 1

Study outline.

where n is the number of prognostic genes, expi is the expression level of prognostic gene i, and βi is the regression coefficient of gene i. Using the median risk score in TCGA cohort as the cutoff value, melanoma patients in each cohort were divided into high-risk and low-risk groups.

Gene set enrichment analysis (GSEA)

The downloaded GSEA software from the Broad Institute (http://www.broadinstitute.org/gsea/index.jsp) was used to identify pathways that were significantly enriched between low-risk and high-risk tumor samples. Gene set permutations were used 1,000 times for each analysis. The normalized enrichment score (NES) value was calculated for each gene set. The nominal P<0.05 and NES >1.5 were used to find the pathway with significant enrichment.

Statistical analysis

Kaplan-Meier analysis was used to demonstrate the survival in the high-risk and low-risk groups in TCGA and GEO cohort, and log-rank test was applied to evaluate the statistical significance. With the Cox’s proportional hazard regression model, univariate and multivariate analyzes were conducted to show whether the gene signature was independent of other clinical characteristics. The time-dependent receiver operating characteristic (ROC) curve and the area under the curve (AUC) were applied to characterize the predictive accuracy of the gene signature. The significance cut-off was fixed at a P value <0.05 for both Cox regression analysis and the log-rank test.

Results

Identification of prognostic gene signature

There were 56,863 and 38,639 genes identified in the TCGA and GEO databases, respectively. After data filtering, 19,524 and 18,379 genes were applied for univariate Cox’s analyses in order to screening candidate genes related to the overall survival rate. By using the cut-off values of P<0.01, we narrowed down 3,248 and 808 candidate genes in the TCGA and GEO databases, respectively, among which 105 genes were present in both databases. After using step-wise multivariate Cox’s analysis, 20 and 54 independent survival predicted genes were recognized respectively by TCGA and GEO databases, respectively. In conclusion, 12 related genes comprising of DOCK10, CCDC7, SLA2, TNFSF14, PWP2, PFKFB3, UBE2L6, PHF11, TRIM14, REC8, SENP6, and NLRP4 were covered in both database and were further analyzed in a subsequent article (). The twelve genes are summarized in .

Figure 2

Correlation between the 12-gene signature and overall survival of malignant melanoma patients. (A) Venn diagram showing overlap among prognostic genes from TCGA and GSE59455 cohorts in univariate and step-wise multivariate cox analysis; (B) the expression heatmap of the 12 prognostic genes; (C) Kaplan-Meier curves of overall survival of the low- and high-risk groups; (D) ROC curve for 1-year survival prediction by the 12-gene signature; (E) ROC curve for 2-year survival prediction by the 12-gene signature.

Table 1

General information of 12 genes for constructing the prognostic signature.

Ensembl ID	Gene name	Gene type	Chromosome	Gene start (bp)	Gene end (bp)	Coefficient
ENSG00000135905	DOCK10	protein coding	2	224765090	225042689	0.171
ENSG00000216937	CCDC7	protein coding	10	32446082	32882864	−0.506
ENSG00000101082	SLA2	protein coding	20	36612318	36646207	−0.398
ENSG00000125735	TNFSF14	protein coding	19	6658126	6670595	0.761
ENSG00000241945	PWP2	protein coding	21	44107323	44131181	−0.298
ENSG00000170525	PFKFB3	protein coding	10	6144918	6254648	−0.223
ENSG00000156587	UBE2L6	protein coding	11	57551655	57568330	−0.314
ENSG00000136147	PHF11	protein coding	13	49495610	49528992	0.007
ENSG00000106785	TRIM14	protein coding	9	98072014	98119420	0.416
ENSG00000100918	REC8	protein coding	14	24171931	24180408	−0.344
ENSG00000112701	SENP6	protein coding	6	75601509	75718281	−0.178
ENSG00000160505	NLRP4	protein coding	19	55836540	55881855	−0.246

Validation of the 12-gene signature in the two independent malignant melanoma cohorts

In our study, we developed a risk-score prognostic model for survival prediction following the expression degree of the 12 surviving relevant genes and their relative contributions based on multiple regression analysis. This 12-gene signature-based prognostic model was confirmed in the TCGA and GEO databases and established a risk assessment criterion for each melanoma patient. We divided patients into high- and low-risk groups in both cohorts using the median risk score in the TCGA cohort as the cutoff value. The 12 gene expression values are shown in . Kaplan-Meier survival curves with log-rank test proved that patients in the low-risk group had a significantly longer survival duration compared to the high-risk group (). The 1- and 2-year survival rates were evaluated by risk scores and are displayed in , with AUC values of 0.81 and 0.785 for the TGGA database, and 0.719 and 0.705 for the GSE59455 cohort, respectively. The results show that the 12-gene signature had high sensitivity and specificity in predicting survival conditions and was appropriate for predicting the survival of patients with melanoma.

The 12-gene signature is an independent prognostic factor of survival

We decided to apply the Cox’s regression model to determine whether the prognostic value of the 12-gene signature was an independent factor compared to other clinical factors in each database. The results are illustrated in . Univariate regression analysis of the TGGA cohort indicated that the 12-gene risk score (high/low: HR =2.955, 95% CI: 2.122–4.114; P=0.000), Stage (≥stage III/< stage III: HR =1.762, 95% CI: 1.276–2.434; P=0.001), Clinical T (≥T2b/0.05). Multivariate regression analysis revealed an indispensable relationship of the 12-gene signature and Clinical T with survival compared with other clinical factors. All in all, the 12-gene risk score (high/low: HR =2.842, 95% CI: 2.031–3.977; P=0.000) was an independent prognostic factor, while Clinical T (≥T2b/

Table 2

Univariate and multivariate Cox regression analyses of the gene signature in two independent cohorts.

Variables	Group	Patients	Univariate analysis			Multivariate analysis
			HR (95% CI)	P		HR (95% CI)	P
TCGA
Age	<58	165	1.000 (reference)			1.000 (reference)
	≥58	181	1.265 (0.923–1.734)	0.144		1.174 (0.843–1.635)	0.342
Gender	Male	214	1.000 (reference)			1.000 (reference)
	Female	132	0.961 (0.690–1.340)	0.816		0.941 (0.674–1.314)	0.721
Stage	<stage III	191	1.000 (reference)			1.000 (reference)
	≥stage III	155	1.762 (1.276–2.434)	0.001		0.625 (0.215–1.819)	0.388
Clinical_T	<T2b	122	1.000 (reference)			1.000 (reference)
	≥T2b	224	1.877 (1.350–2.608)	0.000		1.586 (1.116–2.255)	0.010
Clinical_N	N0	199	1.000 (reference)			1.000 (reference)
	≥N1	147	1.834 (1.326–2.536)	0.000		2.772 (0.968–7.946)	0.058
Clinical_M	M0	335	1.000 (reference)			1.000 (reference)
	M1	11	1.736 (0.709–4.249)	0.228		1.725 (0.662–4.493)	0.264
Risk score	Low	173	1.000 (reference)			1.000 (reference)
	High	173	2.955 (2.122–4.114)	0.000		2.842 (2.031–3.977)	0.000
GSE59455
Age	<68	61	1.000 (reference)			1.000 (reference)
	≥68	61	1.637 (1.126–2.379)	0.010		1.561 (0.896–2.721)	0.116
Gender	Male	84	1.000 (reference)			1.000 (reference)
	Female	38	0.806 (0.548–1.186)	0.274		0.866 (0.520–1.444)	0.582
Mutation	WT	64	1.000 (reference)			1.000 (reference)
	NRAS	27	0.646 (0.408–1.023)	0.063		0.559 (0.328–1.024)	0.056
	BRAF	31	0.892 (0.578–1.376)	0.606		0.666 (0.384–1.158)	0.150
Sun exposure score	<3	43	1.000 (reference)			1.000 (reference)
	≥3	52	0.931 (0.614–1.410)	0.735		1.012 (0.622–1.648)	0.962
Breslow depth	<3.75	40	1.000 (reference)			1.000 (reference)
	≥3.75	40	1.812 (1.150–2.856)	0.010		1.907 (1.095–3.322)	0.023
Risk score	Low	65	1.000 (reference)			1.000 (reference)
	High	57	2.812 (1.866–4.237)	0.000		3.646 (2.094–6.349)	0.000

In the GEO cohort, univariate analysis showed that age (≥68/<68: HR =1.637, 95% CI: 1.126–2.379; P=0.010), Breslow depth (≥3.75/<3.75: HR =1.812, 95% CI: 1.150–2.856; P=0.010) and 12-gene risk score (high/low: HR =2.572, 95% CI: 1.733–3.816; P=0.000) were significantly associated with overall survival, while gender, V-raf murine sarcoma viral oncogene homolog B1 (BRAF)/neuroblastoma RAS viral oncogene homolog (NRAS) mutation and sun exposure score had no significant correlation (P>0.05). Multivariate Cox’s regression analysis revealed that Breslow depth and 12-gene signature had a remarkable prognostic value compared with other clinical factors. In short, the 12-gene risk score (high/low: HR =3.475, 95% CI: 2.030–5.94; P=0.000) was an independent supportive factor, while Breslow depth (≥3.75/<3.75: HR =1.907, 95% CI: 1.095–3.322; P=0.023) was also an independent prognostic factor.

Stratification analysis: prognostic value of 12-gene signature stratified by clinical factors

In our study, BRAF mutation, NRAS mutation and Breslow depth were considered as prognostic-associated factors. Chemotherapy, radiotherapy and immunotherapy are not listed in as certain survival-associated clinical factors. Nonetheless, these treatments still have significant values in determining the prognosis of patients. After adequate evaluation, we decided to stratify the patients in each database into subgroups comprising chemotherapy, radiotherapy, immunotherapy, Breslow depth, BRAF and NRAS mutation, and each subgroup was further distributed into the high-risk and low-risk groups by using the 12-gene signature. In the chemotherapy and radiotherapy groups of the TGGA cohort and BRAF expression group in the GEO cohort, patients in the low-risk group had a longer survival duration compared to patients in the high-risk group (, ; P<0.05), which suggested that the 12-gene risk score is capable of predicting survival of patients who receive chemotherapy and radiotherapy, or patients with BRAF genetic mutations. Although this gene signature was associated with patient prognosis regardless of whether the Breslow depth was greater than or less than 3.75 mm, its prognostic value was more significant when Breslow depth was over 3.75 mm (; P<0.05). Nonetheless, due to the P value being greater than 0.05, the immunotherapy group in TCGA database and the NRAS expression group in the GEO database had no significant difference in statistics, indicating that the 12-gene risk score was not a prognostic factor for patients receiving immunotherapy or having NRAS mutations (, ; P>0.05).

Figure 3

Kaplan-Meier analysis of the overall survival of patients stratified by treatment in TCGA cohort. (A) Chemotherapy; (B) radiotherapy; (C) immunotherapy.

Figure 4

Kaplan-Meier analysis of the overall survival of patients stratified by gene mutations or different Breslow depth. (A) Stratification analysis in GSE59455 cohort based on gene mutations (BRAF and NRAS); (B) Stratification analysis for Breslow depth in TCGA cohort.

Identification of KEGG pathways

We contrasted the expression profile of high-risk or low-risk patients with melanoma. Stratified GSEA revealed that certain immune-related pathways, such as the T-cell signaling pathway, chemokine signaling pathway, and primary immunodeficiency, were significantly enriched by TCGA and GEO cohorts in the low-risk group. This information implied the immune-related properties of the 12-gene signature ().

Figure 5

Immune-related KEGG pathways enriched in the low-risk group from TCGA and GSE59455 cohorts.

Discussion

Numerous articles proposed the opinions that gene alterations play a remarkable role in oncogenesis, treatment, and prognosis, implying an underlying application of gene signatures in the diagnosis and prognosis of melanoma. However, almost none of the systematic articles have specifically delineated the relationships between genes and melanoma prognosis. Therefore, we used regression analyses to test the predictive ability of these accurate and robust melanoma-associated genes. By analyzing the RNA sequence data from TCGA and GEO databases, we made use of 12 genes to formulate a risk score model to predict the prognosis of melanoma patients. Based on the univariate and multivariate analyses, our study serves as a proof-of-concept that this 12-gene signature, along with Clinical T and Breslow depth, was an independent prognostic factor which competitively predicted the survival of melanoma patients.

Because of the heterogeneity of melanoma, it is challenging for us to target different melanoma patients by using the same therapeutic strategies. On the other hand, observing patients with different genetic landscape that have similar histopathological features is very common in clinical practice. Therefore, risk stratification and individualized therapies are urgently needed. Omics-based molecular techniques include genomics, transcriptomics, proteomics, metabolomics and epigenomics (11). A better understanding of the special types gene mutations would help develop the most specialized theranostic method for the different patients and more precisely predict patient’s clinical courses.

Univariate analysis in the TCGA cohort found that age was not significantly associated with overall survival, while the GEO cohort reported a substantial correlation. Age is a significant predictor factor of patient survival at diagnosis or prognosis, and distant metastases tend to more commonly occur in patients older than 70 years than in younger patients (12-14). However, the multivariate analysis confirmed in our cohort that age was not an independent survival predictor. The reason might be that the elderly are not receiving proper treatment in time.

Melanoma classification is based on tumor thickness (T stage), regional lymph node (N stage), and the presence of distant metastasis (M stage) (15,16). Univariate and multivariate analyses showed that stage (≥stage III/

Interestingly, the GSEA results suggest that the 12-gene signature is immune-related, and the specific reasons for this may be a subject for further study. We identified several gene mutations that were connected with the pathogenesis of melanoma. For instance, approximately one-half of the metastatic melanomas have mutations in BRAF, and V600E is the most common mutation in BRAF. Several articles put forth the view that BRAF mutation status is vital in deciding whether patients will benefit from BRAF inhibitor therapies (26-28). NRAS mutations which are commonly found in 15–20% of melanomas are associated with the poorer prognosis, and immune checkpoint inhibitors may be regarded as an effective option for NRAS-mutant melanoma (29,30). The GEO database includes gene mutations from a number of melanoma patients. However, univariate and multivariate analyses in our study showed that these mutations had no significant prognostic value, possibly because most of them are more common in patients with secondary melanoma while most samples in our study had primary melanoma.

Our 12 prognostic panel genes have established roles in tumorigenesis. For instance, DOCK10 could be driven by senescence-related genes and DOCK10-mediated Cdc42 activation was associated with amoeboid invasion of melanoma cells (31,32). RNA expression of TNFSF14 could regulate T-cell responses to tumor cells, prevent melanogenesis in melanocytes and be applied in the therapy of cutaneous pigment disorders (33,34). With regard to PFKFB3 inhibitors, these might influence the progression and development of a tumor and have been identified as an attractive target for therapeutic strategies (35). Moreover, UBE2L6 was reported to influence the chemosensitivity in esophageal cancer cells (36). Increased PHF11 was confirmed to have an association with increased total IgE levels and asthma (37). As for TRIM14, it was reported to promote the invasion and migration in several types of tumor including glioblastoma (38). Yu et al. (39) reported that REC8 mutation was an independent risk factor for the shortened survival of gastric cancer patients.

Furthermore, SENP6 was confirmed as having predictive strength in liver transplantation tolerance (40), while NLRP4 has demonstrated the ability to predict asthma exacerbation in adult asthmatics who smoke (41). Overall, our study suggested that these 12 gene signatures were independent factors to predict overall survival in patients of melanoma. The strength of our work was to establish a gene-associated risk model and to evaluate the survival of melanoma patients by combing the TCGA/GEO data and clinical features. It is more stable and reliable to make choices in those survival-related genes that appear in the TCGA and GEO cohorts simultaneously, rather than in those that appear in a single database. Hence, the co-occurrence of these 12 genes (DOCK10, CCDC7, SLA2, TNFSF14, PWP2, PFKFB3, UBE2L6, PHF11, TRIM14, REC8, SENP6, NLRP4) was identified as a novel factor which has great potential as an indicator for patient prognosis of melanoma.

In brief, by analyzing RNA sequence-based gene expression signatures in TCGA and GEO patients, our results demonstrated the capacity of the 12-gene signature to predict survival rates in melanoma. Furthermore, the stratified and multivariate analysis confirmed that the gene panel was an independent factor with other pathological and clinical characteristics, and as a consequence is an underlying prognostic predictor of melanoma.