Identification of prognostic biomarkers for malignant melanoma using microarray datasets.

Malignant melanoma is one of the most common types of cancer worldwide. Efforts have been made to elucidate the pathology of malignant melanoma. However, its molecular mechanisms remain unclear. Therefore, the microarray datasets GSE3189, GSE4570 and GSE4587 from the Gene Expression Omnibus database were used for the elucidation of candidate genes involved in the initiation and progression of melanoma. Assessment of the microarray datasets led to the identification of differentially expressed genes (DEGs), which were subsequently used for function enrichment analysis. These data were utilized in the construction of the protein-protein interaction network and module analysis was conducted using STRING and Cytoscape software. The results of these analyses led to the identification of a total of 182 DEGs, including 52 downregulated and 130 upregulated genes. The functions and pathways found to be enriched in the DEGs were GTPase activity, transcription from RNA polymerase II promoter, apoptotic processes, cell adhesion, membrane related pathways, calcium signaling cascade and the PI3K-Akt signaling pathway. The identified genes were demonstrated to belong to a set of 10 hub genes biologically involved in proliferation, apoptosis, cytokinesis, adhesion and migration. Survival analysis and Oncomine database analysis revealed that the calmodulin gene family, BAX and VEGFA genes, may be associated with the initiation, invasion or recurrence of melanoma. In conclusion, the DEGs and hub genes identified in the present study may be used to understand the molecular pathways involved in the initiation and progression of malignant melanoma. Furthermore, the present study may aid in the identification of possible targets for the diagnosis and treatment of melanoma. Copyright: © Lin et al.

Introduction

Melanoma is a malignant skin tumor derived from melanocytes, which has a high degree of malignancy and leads to high mortality. During the past 40 years, its incidence, as well as the proportion of mid- to late-stage tumors and infeasibility of surgery have increased worldwide (1).

Previous studies have demonstrated that abnormal expression and mutations of genes and proteins, including cell division cycle associated 8 (encoded by CDCA8), telomerase reverse transcriptase (TERT), B-Raf proto-oncogene (BRAF) and various tumor suppressor genes are involved in the initiation and progression of melanoma. For example, it has been reported that the CDCA8 gene is capable of promoting the malignant progression of cutaneous melanoma and is associated with poor prognosis (2). TERT promoter mutations have also been identified in up to 50% of cutaneous melanoma cases in the global population; however, their incidence in Asian populations remains unclear (3). BRAF mutations, particularly those located at codon 600, have been observed in 50% of malignant melanoma cases worldwide (4). Furthermore, overexpression of BRAF and hypermethylation of Ras binding proteins have been revealed to be associated with poor prognosis in patients with malignant melanoma (5,6). Melanoma mortality remains high due to the absence of efficient diagnostic techniques at the initial stages of the disease. Therefore, understanding the mechanisms involved in the initiation, proliferation and recurrence of this type of cancer at the molecular level is essential for the development of more effective diagnostic and treatment strategies.

Over the past few decades, microarray analyses as well as bioinformatics studies have been increasingly favored for the screening of genetic changes at the genomic level. These identification methods may be employed for the determination of differentially expressed genes (DEGs), as well as functions that may be involved in melanoma initiation and progression (7). However, the reliability of independent microarray analyses may not be high owing to the rate of false positives. Therefore, in the present study, three different mRNA microarray datasets were obtained from the Gene Expression Omnibus (GEO) database. These datasets were analyzed for the determination of DEGs between malignant melanoma and normal nevi tissues. Thereafter, Gene Ontology (GO) and Kyoto Gene and Genomic Encyclopedia (KEGG) analyses, as well as protein-protein interaction (PPI) network analysis were conducted to identify molecular processes associated with melanoma development and progression. Altogether, 182 DEGs and 10 hub genes were indicated as potential biomarkers of melanoma.

Materials and methods

Data from the microarray analyses

GEO (http://www.ncbi.nlm.nih.gov/geo) (8) is a publicly available function-related genomics repository containing high-throughput gene expression data, ChIP-seq data, as well as microarrays. Three gene datasets, GSE3189 (9), GSE4570 (10) and GSE4587 (11), were downloaded from the GEO database. GSE3189 and GSE4570 were based on the ArrayGPL96 platform (Affymetrix Human Genome U133A Array), whereas GSE4587 was based on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0). The probes were later converted to their analogous gene symbols using platform information. The GSE3189 dataset comprised 45 melanoma tissue samples and 18 normal nevi tissue samples, GSE4570 contained 6 melanoma samples and 2 nevi samples, and GSE4587 contained 7 melanoma samples and 8 nevi samples.

Identification of DEGs

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r) was used for the screening of DEGs between melanoma and normal nevi tissue samples. GEO2R is a web tool used in interaction studies for the comparison of various datasets in a GEO series to identify DEGs. As mentioned above, the most common limitation of microarray analyses is the presence of false positives, which may be limited by adjusting the P-values and calculating the Benjamini-Hochberg false discovery rate (12). These adjustments led to the identification of statistically significant genes. Probe sets that had either no corresponding gene symbols or genes with numerous probe sets were eliminated or averaged, respectively. DEGs of logarithmic fold change value >1 were selected in the present study. P<0.05 was considered to indicate a statistically significant difference.

KEGG and GO enrichment analysis of the DEGs

The Database for Annotation, Visualization and Integrated Discovery (DAVID; version 6.8; http://david.ncifcrf.gov) (13) is an online bioinformatics database that integrates biological data with analytic tools, and offers substantial gene- and protein-related information, thus contributing to the extraction of biological data. KEGG database is a online tool used to study advanced functions and biological processes of genes via high-throughput sequencing (14). GO, a significant web-based tool used for annotating genes and analyzing the biological processes that these genes are involved in, was also employed for DEG enrichment (15). The biological functions of the DEGs were analyzed using DAVID. P<0.05 was considered to indicate a statistically significant difference.

Construction of a PPI network and module analysis

The Search Tool for the Retrieval of Interacting Genes (STRING; version 10.5; http://string-db.org) database was used to analyze the PPI network of genes (16). Analysis of functionally relevant interactions among the proteins encoded by DEGs may provide valuable insights into the mechanisms of genesis or progression of various diseases. In the present study, the STRING database was used for the construction of the DEG PPI network. A statistically relevant interaction was defined using STRING (combined score >0.4). Cytoscape (version 3.7.0), an open access bioinformatics software, is a platform used to study networks of molecular interactions (17). Furthermore, Molecular Complex Detection (MCODE; version 1.5.1), a Cytoscape plugin, clusters a given network based on topology for the determination of compact connected portions (18). Cytoscape was used for the identification of the PPI network, while MCODE was used to identify the most significant interaction. The MCODE selection criteria were as follows: i) MCODE score >5; ii) MCODE degree cut-off level=2; iii) node score cut-off level=0.2; iv) max depth=100 and v) k-score=2.

Hub genes selection and analysis

After the construction of the PPI network, genes (MCODE degrees ≥10 using Cytoscape) were identified as hub genes. The gene network and genes that were co-expressed within this network were determined using cBioPortal (http://www.cbioportal.org) (19,20). Investigation of the biological processes associated with the hub genes was conducted using the Biological Networks Gene Oncology (version 3.0.3) Cytoscape plugin (21). The overall survival and disease-free survival analyses of the hub genes were conducted using the Kaplan-Meier function in cBioPortal. The gene expression levels of the hub genes between melanoma and normal nevi tissues were evaluated using Oncomine, an online database (http://www.oncomine.com) (22,23). Differences in expression between melanoma and normal nevi tissues were analyzed by t-test. P<0.05 was considered to indicate a statistically significant difference.

Results

DEG identification in melanoma

With logarithmic fold change value >1 and P<0.05, the DEGs (2,484 in GSE3189, 1,300 in GSE4570 and 6,759 in GSE4587) were identified. The intersection of the 3 datasets comprised of 182 genes, 52 of which were downregulated and 130 were upregulated between melanoma and normal nevi tissues, as demonstrated in the Venn diagram (Fig. 1A).

Figure 1.

Venn diagram of the three gene microarray datasets, PPI network obtained using bioinformatics algorithms and statistically significant DEG module. (A) DEGs with a fold change >2 and P<0.05 were selected among the gene expression profiling GSE4587, GSE3189 and GSE4570 datasets. A common set of 182 genes was identified at the intersection of the 3 datasets. (B) Cytoscape was used for the construction of the PPI network of the DEGs. (C) Module of maximum relevance, consisting of 10 nodes and 29 edges, extracted from PPI network. Genes with upregulated expression are shown in pink, whereas genes with downregulated expression are shown in blue. PPI, protein-protein interaction; DEG, differentially expressed gene; DES, desmin; VEGFA, vascular endothelial growth factor A, CALM1, calmodulin 1; CALM2, calmodulin 2; CALM3, calmodulin 3; FN1, fibronectin 1; PRKCA, protein kinase C α; IGF1, insulin-like growth factor 1; RB1, retinoblastoma transcriptional corepressor 1.

Enrichment analyses of the DEGs using KEGG and GO

Biological classification of the DEGs, achieved via functional and pathway enrichment analyses, was performed using DAVID. The results of the GO analysis revealed that the DEGs in biological processes were enriched in transcription from RNA polymerase II promoter, cell adhesion, GTPase activity, apoptotic processes and transcription. At the molecular level, changes in functions were mainly enriched in actin and protein kinase binding. Regarding the ‘cellular component’, DEGs were mainly enriched in the cell membrane, cytoskeleton and extracellular region (Table I). KEGG pathway analysis demonstrated that the downregulated DEGs were mostly enhanced in the estrogen signaling cascade, melanogenesis and the calcium signaling pathway, whereas the upregulated DEGs were mostly enriched in microRNAs in focal adhesion and pathways in cancer, as well as the PI3k-Akt signaling pathway.

Table I.

Enrichment study of DEGs between malignant melanoma tissue and normal nevus tissue.

A. Downregulated DEGs
Term and description	Count in gene set	P-value
GO:0030801~Positive regulation of cyclic nucleotide metabolic process	3	3.39×10−5
GO:0005980~Glycogen catabolic process	4	4.72×10−5
GO:0051343~Positive regulation of cyclic-nucleotide phosphodiesterase activity	3	1.12×10−4
GO:0001975~Response to amphetamine	4	1.56×10−4
GO:1901841~Regulation of high voltage-gated calcium channel activity	3	2.35×10−4
GO:1901844~Regulation of cell communication by electrical coupling involved in cardiac conduction	3	3.13×10−4
GO:0007190~Activation of adenylate cyclase activity	4	3.35×10−4
GO:0060316~Positive regulation of ryanodine-sensitive calcium-release channel activity	3	4.01×10−4
GO:0006936~Muscle contraction	5	4.75×10−4
GO:0043647~Inositol phosphate metabolic process	4	5.41×10−4
GO:0021762~Substantia nigra development	4	6.12×10−4
GO:0019233~Sensory perception of pain	4	7.28×10−4
GO:0060315~Negative regulation of ryanodine-sensitive calcium-release channel activity	3	7.31×10−4
GO:0032516~Positive regulation of phosphoprotein phosphatase activity	3	8.62×10−4
GO:0005513~Detection of Ca2+	3	1.00×10−3
GO:0010801~Negative regulation of peptidyl-threonine phosphorylation	3	1.00×10−3
GO:0051412~Response to corticosterone	3	1.67×10−3
GO:0010880~Regulation of release of sequestered Ca2+ into cytosol by sarcoplasmic reticulum	3	1.67×10−3
GO:0010881~Regulation of cardiac muscle contraction by regulation of the release of sequestered Ca2+	3	1.87×10−3
GO:0060314~Regulation of ryanodine-sensitive calcium-release channel activity	3	1.87×10−3
GO:0055117~Regulation of cardiac muscle contraction	3	2.28×10−3
GO:0035307~Positive regulation of protein dephosphorylation	3	2.28×10−3
GO:0031954~Positive regulation of protein autophosphorylation	3	2.28×10−3
GO:0051000~Positive regulation of nitric-oxide synthase activity	3	2.50×10−3
GO:0032465~Regulation of cytokinesis	3	2.50×10−3
GO:0043547~Positive regulation of GTPase activity	8	3.00×10−3
GO:0001822~Kidney development	4	3.11×10−3
GO:0043065~Positive regulation of apoptotic process	6	3.45×10−3
GO:0050999~Regulation of nitric-oxide synthase activity	3	3.49×10−3
GO:0010800~Positive regulation of peptidyl-threonine phosphorylation	3	3.76×10−3
GO:0043388~Positive regulation of DNA binding	3	4.04×10−3
GO:0022400~Regulation of rhodopsin-mediated signaling pathway	3	4.33×10−3
GO:0002027~Regulation of heart rate	3	5.59×10−3
GO:0071902~Positive regulation of protein serine/threonine kinase activity	3	6.27×10−3
GO:0007223~Wnt signaling pathway, calcium-modulating pathway	3	7.74×10−3
GO:0045893~Positive regulation of transcription, DNA-templated	7	7.81×10−3
GO:0034704~Calcium channel complex	3	2.75×10−3
GO:0030017~Sarcomere	3	6.61×10−3
GO:0005876~Spindle microtubule	3	8.36×10−3
GO:0043274~Phospholipase binding	4	2.72×10−5
GO:0031997~N-terminal myristoylation domain binding	3	3.24×10−5
GO:0072542~Protein phosphatase activator activity	3	1.61×10−4
GO:0030235~Nitric-oxide synthase regulator activity	3	2.99×10−4
GO:0008440~Inositol-1,4,5-trisphosphate 3-kinase activity	3	4.78×10−4
GO:0031996~Thioesterase binding	3	9.59×10−4
GO:0031432~Titin binding	3	9.59×10−4
GO:0043539~Protein serine/threonine kinase activator activity	3	1.60×10−3
GO:0015276~Ligand-gated ion channel activity	3	6.33×10−3
GO:0001105~RNA polymerase II transcription coactivator activity	3	6.68×10−3
GO:0019901~Protein kinase binding	6	8.04×10−3
hsa05010:Alzheimer's disease	6	9.66×10−4
hsa04915:Estrogen signaling pathway	5	1.07×10−3
hsa04922:Glucagon signaling pathway	5	1.07×10−3
hsa04916:Melanogenesis	5	1.11×10−3
hsa04020:Calcium signaling pathway	6	1.28×10−3
hsa04728:Dopaminergic synapse	5	2.75×10−3
hsa05214:Glioma	4	3.29×10−3
hsa05031:Amphetamine addiction	4	3.43×10−3
hsa04720:Long-term potentiation	4	3.43×10−3
hsa04744:Phototransduction	3	6.84×10−3
B, Upregulated DEGs
Term and description	Count in gene set	P-value
GO:0007155~Cell adhesion	13	1.02 ×10−4
GO:0007010~Cytoskeleton organization	7	1.01×10−3
GO:0001501~Skeletal system development	6	2.92×10−3
GO:0045944~Positive regulation of transcription from RNA polymerase II promoter	16	3.94×10−3
GO:0008584~Male gonad development	5	4.48×10−3
GO:0002576~Platelet degranulation	5	6.18×10−3
GO:0042475~Odontogenesis of dentin-containing teeth	4	6.96×10−3
GO:0090190~Positive regulation of branching involved in ureteric bud morphogenesis	3	7.87×10−3
GO:0050679~Positive regulation of epithelial cell proliferation	4	8.85×10−3
GO:0016020~Membrane	29	8.21×10−4
GO:0005886~Plasma membrane	45	8.97×10−4
GO:0005856~Cytoskeleton	10	1.05×10−3
GO:0031012~Extracellular matrix	8	4.34×10−3
GO:0005925~Focal adhesion	9	5.57×10−3
GO:0031093~Platelet α granule lumen	4	6.37×10−3
GO:0005576~Extracellular region	21	6.61×10−3
GO:0005913~Cell-cell adhesion junction	8	6.92×10−3
GO:0043034~Costamere	3	7.39×10−3
GO:0005887~Integral component of plasma membrane	19	7.89×10−3
GO:0003779~Actin binding	11	2.55×10−5
GO:0005200~Structural constituent of cytoskeleton	5	7.42×10−3
GO:0050839~Cell adhesion molecule binding	4	9.33×10−3
hsa04510:Focal adhesion	9	9.83×10−4
hsa04512:ECM-receptor interaction	6	1.72×10−3
hsa04151:PI3K-Akt signaling pathway	11	2.13×10−3
hsa05205:Proteoglycans in cancer	8	3.71×10−3
hsa05200:Pathways in cancer	11	5.44×10−3
hsa05206:MicroRNAs in cancer	9	7.47×10−3

GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes; MM, malignant melanoma.

Construction of the PPI network and module analysis

The DEG PPI network was constructed to find novel protein interactions. A total of 122 nodes and 266 protein interaction pairs were identified (Fig. 1B). The most significant interaction of 10 nodes and 29 protein interaction pairs was identified using Cytoscape (Fig. 1C).

Selection and analysis of hub genes

Ten genes (MCODE degrees ≥10 using Cytoscape) were revealed as hub genes. The names, acronyms and roles of these hub genes are presented in Table II. The cBioPortal online platform was used to analyze the hub gene network, as well as their co-expression genes (Fig. 2). The cBioPortal network contains 60 nodes, including 10 hub genes and the 50 most frequently altered neighbor genes. The results of the hub gene biological process analysis are demonstrated in Fig. 3. The biological processes ‘cell division’, ‘cytokinesis, release of sequestered ion into cytosol’, ‘release of sequestered calcium ion into cytosol by sarcoplasmic reticulum’, ‘ryanodine-sensitive calcium-release channel activity’ and ‘ion transmembrane transporter activity’ were significantly enriched. Subsequently, survival analysis of the hub genes was performed using Kaplan-Meier analysis, as presented in Fig. 4. Melanoma patients with BAX alterations were identified to have poor disease-free as well as overall survival (Fig. 4), indicating that BAX may serve a significant role in the initiation or progression of melanoma. Oncomine analysis of melanoma and normal nevi tissues revealed that BAX was significantly overexpressed in melanoma tissues in the different datasets (Fig. 5). In the Oncomine database, the mRNA expression levels of BAX (P<0.0001), CALM1 [which encodes calmodulin (CaM1)] (P<0.0001), CALM3 (CaM3; P<0.0001), FN1 (fibronectin 1; P<0.0001), PRKCA (protein kinase C α; P=0.0320), RB1 (RB transcriptional corepressor 1; (P<0.0001) and VEGFA (vascular endothelial growth factor A; P<0.0001) genes were demonstrated to be higher in melanoma tissues than in the normal nevi tissues. By contrast, the mRNA expression level of IGF1 (insulin-like growth factor 1) was found to be lower in melanoma tissues. There was no statistically difference in the expression of DES (P=0.0628) and CALM2 (P=0.0552) between melanoma and normal nevi tissues.

Table II.

Functional roles of the 10 hub genes with a degree ≥10.

No.	Gene symbol	Full name	Function
1	VEGFA	Vascular endothelial growth factor A	Induction of proliferation and migration of vascular endothelial cells; essential for both physiological as well as pathological angiogenesis.
2	BAX	BCL2-associated X, apoptosis regulator	It acts as either an anti- or proapoptotic regulator and is involved in several cellular activities.
3	CALM1	Calmodulin 1	Expression of a member of the EF-hand calcium-binding protein family
4	CALM2	Calmodulin 2	Role in proliferation, cell cycle progression, and signaling cascades.
5	CALM3	Calmodulin 3	Calcium binding; enzymatic co-factor; role in cell cycle regulation as well as cytokinesis.
6	FN1	Fibronectin 1	Fibronectin expression; involvement in cell adhesion and migration processes.
7	PRKCA	Protein kinase C α	Essential roles in numerous cellular processes, including cell volume control, checkpoints of the cell cycle, cell adhesion, as well as cell transformation.
8	IGF1	Insulin-like growth factor 1	Expression of protein with similar structure as well as function to insulin; member of a family of proteins involved in regulating growth and development.
9	RB1	RB transcriptional corepressor 1	Expression of negative regulatory protein of the cell cycle; the first tumor suppressor gene to be identified.
10	DES	Desmin	Expression of a muscle-specific class III intermediate filament. Homopolymers of this protein form a stable intracytoplasmic filamentous network which connects myofibrils to each other as well as to the plasma membrane.

Figure 2.

Hub gene interaction network. cBioPortal was used to examine the hub genes and their co-expression genes. Blue arrows represent controlling state change of genes, green arrows represent controlling expression of genes and brown lines represent complex with other genes. Nodes with a bold black outline represent the hub genes, while nodes with a thin black outline represent co-expression genes. DES, desmin; VEGFA, vascular endothelial growth factor A; CALM1, calmodulin 1; CALM2, calmodulin 2; CALM3, calmodulin 3; FN1, fibronectin 1; PRKCA, protein kinase C α; IGF1, insulin-like growth factor 1; RB1, retinoblastoma transcriptional corepressor 1.

Figure 3.

BiNGO analysis of the biological process of genes with a degree ≥10. BiNGO was used for hub gene biological process analysis. The color depth of the nodes represents the corrected P-value of the ontologies. The magnitude of the nodes represents the genes that participate in the ontologies. BiNGO, Biological Networks Gene Oncology.

Figure 4.

Disease-free survival and overall survival analyses of the hub genes using Kaplan-Meier analysis on the cBioPortal online platform. (A) Disease-free survival and (B) overall survival. P<0.05 was considered to indicate a statistically significant difference. DES, desmin; VEGFA, vascular endothelial growth factor A; CALM1, calmodulin 1; CALM2, calmodulin 2; CALM3, calmodulin 3; FN1, fibronectin 1; PRKCA, protein kinase C α; IGF1, insulin-like growth factor 1; RB1, retinoblastoma transcriptional corepressor 1.

Figure 5.

Expression of the hub genes in the GSE3189 dataset obtained using Oncomine online analysis. 0, No value (normal skin); 1, Skin nevus; 2, Cutaneous melanoma. DES, desmin; VEGFA, vascular endothelial growth factor A; CALM1, calmodulin 1; CALM2, calmodulin 2; CALM3, calmodulin 3; FN1, fibronectin 1; PRKCA, protein kinase C α; IGF1, insulin-like growth factor 1; RB1, retinoblastoma transcriptional corepressor 1.

Discussion

Melanoma is a highly metastatic type of cancer, which exhibits strong resistance to both chemotherapy and radiotherapy (24–26). It has been observed to develop rapidly during the early phase. Melanoma cells acquire mutations during this phase. Subsequently, mutant cells invade the dermal layer and trigger angiogenesis. This ‘angiogenic switch’ is involved in invasiveness and is characterized by the regulation of genes, including VEGFA, VEGF receptor genes and related angiogenic signaling pathways (27). However, the molecular pathways underlying melanoma remain largely unknown. The development of biomarkers with improved accuracy is essential for the effective diagnosis and treatment of melanoma. Genetic changes in melanoma may be observed using microarray technology, which may also be beneficial for the identification of novel biomarkers in other diseases.

The present study used 3 mRNA microarray datasets to identify DEGs between malignant melanoma tissues and normal nevi tissues. A total of 182 DEGs were identified, which included 52 downregulated and 130 upregulated genes. DEG interactions were studied using GO and KEGG enrichment analyses. The upregulated genes were found to be associated with the cell membrane, cytoskeleton, extracellular region, actin binding, mitotic cell cycle and PI3K-Akt signaling cascades, while the downregulated genes were found to be involved in progressive regulation of transcription, protein kinase binding, melanogenesis, and the estrogen and calcium signaling pathways. Literature retrieval results indicated that the associations between malignant melanoma and these molecular mechanisms (oocyte meiosis, protein kinase binding and estrogen signaling pathway) have not been reported widely. A total of 10 DEGs with degrees ≥10 in the PPI network were selected as hub genes. The PPI network revealed that BAX directly interacts with MYCN, IGF1, RUNX1 (runt-related transcription factor 1), MSH6, PRKDC (protein kinase, DNA-activated, catalytic subunit), CALM1, CALM2, CALM3, TP63, VEGFA and RB1, indicating a key role for BAX in melanoma. Furthermore, the PPI network in DEGs is a novel gene interaction observed in malignant melanoma. It is well known that VEGFA and BAX are involved in tumor malignancy (28,29), a finding confirmed in the present study. VEGFA is related to angiogenesis, a process which is required for tumor growth and metastasis (30). Vasculogenic mimicry (VM) is an endothelial vessel supply system in cancers that VM reflects the aggressive ability of tumor cells (31). Recently, the c-Myc gene was reported to promote tumorigenesis of melanoma by promoting vasculogenic mimicry via the Bax signaling pathway (24). VEGFA overexpression has also been observed in lung, pancreatic and other cancers (32–34). In addition to its role in cell cycle progression, VEGFA is a potent angiogenic factor, which is required for oncogenesis (35,36). CALM1, CALM2 and CALM3 belong to the CaM gene family (37,38). They all encode a similar CaM protein, with differences at the nucleotide level. CaM serves an essential role in disease pathogenesis via the Ca2+ signaling pathway (39,40). Furthermore, it is involved in apoptosis by balancing the expression levels of the proapoptotic protein, BAX, with those of the antiapoptotic protein, Bcl-2 (41).

Among the identified hub genes, BAX overexpression was found to be associated with the lowest survival rate. The protein encoded by this gene is a part of the Bcl-2 protein family that consists of antiapoptotic and proapoptotic members. BAX gene expression is associated with shorter patient survival, chemoresistance and recurrence in melanoma (42,43). Thus, it is regarded as a target for anticancer agents. The expression of the hub genes in relation to both overall and disease-free survival was evaluated. BAX alteration was found to significantly affect both overall survival and disease-free survival. BAX could induce the decrease in overall survival and disease-free survival. Moreover, clinical studies have reported that a shorter survival period is significantly associated with BAX gene overexpression (44,45). However, the expression levels of the other hub genes in overall survival were not statistically significant compared to BAX. This result may have occurred due to the fact that survival analysis in cBioPortal is performed based on a relationship between gene mutation and prognosis. However, gene overexpression may arise via either mutation or amplification. Accordingly, hub gene overexpression in melanoma may occur due to gene amplification rather than mutation, thus creating the need for further research in order to confirm the association between melanoma and the hub genes. Oncomine analysis demonstrated that mRNA expression levels of BAX, CALM1, CALM3, FN1, PRKCA, RB1 and VEGFA were higher in melanoma tissue than in normal nevi tissues, whereas the mRNA expression level of IGF1 was lower in melanoma tissues. Previous studies have reported that CALM2 levels in gastric, breast and other cancer tissues are higher than those in normal tissues, and that this gene may subsequently be used as a prognostic target (46–48). Oncomine analysis in melanoma indicated that CALM2 expression in melanoma was not higher. Further investigation is therefore required to confirm CALM2 expression in cancers. The DES gene encodes desmin, one of the first muscle-specific proteins to be expressed during the early phases of skeletal and cardiac muscle differentiation (49). Certain studies have suggested that desmin expression is elevated in colorectal cancer, and that it may be used as a novel prognostic predictor (50,51). DES expression is increased in osteogenic melanoma, however the expression of DES is low in other types of melanoma (52). Thus, the reason why there was no significant difference in the Oncomine database analysis results may be that the Oncomine database did not include the classification of melanoma. Similarly, FN1, PRKCA, RB1 and IGF1 had been reported to influence tumorigenesis and initiation in other types of cancer (53,54), which was consistent with our study in melanoma. It was speculated that the Ca2+ signaling pathway is associated with malignant melanoma. The disruption of the homeostasis of this pathway during tumorigenesis leads to abnormal expression of BAX (55). Disruption of Ca2+ signaling pathway homeostasis also promotes tumor angiogenesis via the overexpression of VEGFA.

There were certain limitations in the present study. Firstly, the GES3189 dataset comprised 45 melanoma samples and 18 normal nevi samples, whereas the sample quantity of the other two databases was insufficient. However, since three databases were used to choose the overlap in Venn diagram as DEGs, this may improve the credibility of the analysis results. Secondly, certain biomarkers associated with melanoma were identified; however, further experimental studies, including immunohistochemistry, animal testing and clinical trials, are required to validate these findings. Despite these limitations, there is few report that upregulation of the CALM gene family (CALM1, CALM2 and CALM3) in malignant melanoma was associated with poor prognosis. In addition, the present study was the first to report associations between the identified DEGs and hub gene interactions in malignant melanoma.

In conclusion, the present study identified DEGs that may be associated with the initiation or progression of melanoma. The 182 DEGs and 10 hub genes that were identified may be considered as potential biomarkers of melanoma. Nonetheless, additional research is required to further understand the biological functions of these genes in melanoma.