Prognostic genes of melanoma identified by weighted gene co-expression network analysis and drug repositioning using a network-based method.

Melanoma is one of the most malignant types of skin cancer. However, the efficacy and utility of available drug therapies for melanoma are limited. The objective of the present study was to identify potential genes associated with melanoma progression and to explore approved therapeutic drugs that target these genes. Weighted gene co-expression network analysis was used to construct a gene co-expression network, explore the associations between genes and clinical characteristics and identify potential biomarkers. Gene expression profiles of the GSE65904 dataset were obtained from the Gene Expression Omnibus database. RNA-sequencing data and clinical information associated with melanoma obtained from The Cancer Genome Atlas were used for biomarker validation. A total of 15 modules were identified through average linkage hierarchical clustering. In the two significant modules, three network hub genes associated with melanoma prognosis were identified: C-X-C motif chemokine receptor 4 (CXCR4), interleukin 7 receptor (IL7R) and phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit γ (PIK3CG). The receiver operating characteristic curve indicated that the mRNA levels of these genes exhibited excellent prognostic efficiency for primary and metastatic tumor tissues. In addition, the proximity between candidate genes associated with melanoma progression and drug targets obtained from DrugBank was calculated in the protein interaction network, and the top 15 drugs that may be suitable for treating melanoma were identified. In summary, co-expression network analysis led to the selection of CXCR4, IL7R and PIK3CG for further basic and clinical research on melanoma. Utilizing a network-based method, 15 drugs that exhibited potential for the treatment of melanoma were identified. Copyright: © Wang et al.

Introduction

Melanoma is a tumor that originates in melanocytes of the skin or other parts of the body (1). The main function of melanocytes is to produce melanin via melanogenesis, a multistep biochemical process regulated by L-tyrosine, L-DOPA and other hormones (2,3). Melanogenesis leads to the upregulation of hypoxia-inducible factor 1, which modulates the cellular metabolism of melanoma (4). A previous study has demonstrated that pigmentation level is associated with the overall and disease-free survival time of patients with stage III and IV melanoma (5). In the United States, >91,000 individuals were diagnosed with cutaneous melanoma in 2018, and >9,000 patients succumbed to the disease in the same period (6). Since melanoma tends to spread lymphogenously and hematogenously, patients with inoperable metastatic melanoma exhibit median survival times between 8 and 12 months (7). Therefore, melanoma poses a serious threat to life.

Gene mutations in melanoma may activate multiple signaling pathways that regulate proliferation, epithelial-mesenchymal transition, invasion and metastasis in an abnormal manner (8). For example, BRAF mutations, predominantly V600E, occur in 40–50% of all melanomas, whereas NRAS proto-oncogene, GTPase and neurofibromin 1 mutations occur in ~20 and 15% of melanomas, respectively (9). Targeted therapy and immunotherapy have been demonstrated to be effective treatment methods (10,11). BRAF/mitogen-activated protein kinase kinase inhibitors, as well as antibodies against cytotoxic T-lymphocyte-associated protein 4 and programmed cell death protein 1 have been used for treatment of metastatic melanoma, with patient response rates ranging between 20 and 70% (12). Although these breakthrough treatments have prolonged progression-free survival to a certain extent, drug resistance still limits their effectiveness (13). For example, immune-based therapy is subject to limitations, such as the prevention of the generation of an immunosuppressive environment (14). Therefore, there remains a need for novel markers of prognosis and novel therapeutic drugs for melanoma.

Weighted gene co-expression network analysis (WGCNA) is widely used to analyze genetic expression data, locate modules of highly correlated genes and identify potential biomarkers, as well as therapeutic targets. Thus, the present study aimed to utilize WGCNA to identify novel biomarkers associated with melanoma prognosis. Additionally, the present study aimed to determine the proximity between disease-associated proteins and drug targets in the human protein-protein interactome in order to identify potential drugs for the treatment of melanoma.

Materials and methods

Data processing

Melanoma transcriptome dataset GSE65904 (15) was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo). GSE65904 comprised 214 samples from patients with melanoma, no non-tumor tissue samples or healthy subjects were included. Illumina HumanHT-12V4.0 expression beadchip was used as the sequencing platform. Clinical information of patients, including sex, age, tumor stage, distant metastasis and survival state, was collected. The GEO query package in R v2.52.0 (https://git.bioconductor.org/packages/GEOquery) was used to process the data. If the expression of a gene was not significant compared with the background value (standard probe) in >25% of all samples (P>0.05), the probe was removed from further analysis. A total of 10,566 genes were obtained.

Weighted co-expression network construction

The top 50% most differentially expressed genes (5,283 genes) were selected for WGCNA analysis following analysis of variance using R 3.3.2 (https://www.r-project.org/) (16). These genes were used for screening and cluster analysis of all samples, as well as to identify outliers, following which one patient was removed from the study (Fig. 1). The gene expression data of the patients was used to construct the co-expression network, and the WGCNA algorithm was utilized for analysis (16). To ensure that the nodes of the constructed co-expression network conformed to the power rate distribution, appropriately soft threshold was selected (β=3), which enabled the deletion of low mutual correlation relationships. The distribution of network nodes conformed to the power rate distribution at β=3. Further investigation of the distribution of node degrees in the co-expression network revealed that the degree of nodes conformed to the power law distribution. This indicated that the constructed co-expression network was a scale-free network, conforming to the characteristics of common biological networks. The average linkage hierarchical clustering method (17) was used to cluster all genes.

Figure 1.

Cluster dendrogram of 214 melanoma samples. The GSE65904 dataset was used. The red line indicates the outlier to rule out biased samples. The black lines represent each sample in the dataset, and the numbers represent corresponding GSM of the patient.

Identification of clinically significant modules

To obtain the gene modules that were associated with clinical phenotypes, the correlation between modules and clinical phenotypes was determined. Module eigengenes (MEs) were considered as characteristics of all genes in a certain module. The association between MEs and clinical characteristics was analyzed to determine a clinically significant module for further use.

Gene Ontology (GO) and pathway enrichment analysis

The ClusterProfiler package (https://github.com/GuangchuangYu/clusterProfiler) in R v3.12.0 was used to determine the functions of the enriched genes from the two modules in Fig. 3 (black and turquoise modules) in GO (18) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (19) pathway analysis, respectively. Genes in the clinically significant module were categorized into three functional groups: Biological process (BP), cellular component (CC) and molecular function (MF).

Figure 3.

Identification of modules associated with the progression of melanoma. (A) Dendrogram of all differentially expressed genes clustered based on a dissimilarity measure. (B) Heatmap of the association between module eigengenes and the progression of melanoma. DMFS, distant metastasis-free survival; Dss, disease-specific death survival; ME, module eigengene. The red color of each box represents the positive association between the module and trait whereas the green color of each box represents the negative associations. The association of the module and trait is calculated to be between −1 and 1.

Identification and validation of hub genes

To identify genes associated with melanoma prognosis, the association between each gene and clinical characteristics was evaluated, as well as the association between each gene and core modules, such as module membership (MM) and gene significance (GS). MM is defined by the correlation between the gene expression profile and MEs, whereas GS is defined by the association between a gene and external traits. Genes with |MM+GS|=5% in the aforementioned modules (black and turquoise modules) in Fig. 3 were selected as potentially prognostic genes; all other genes were removed. To further analyze the association between these genes, the remaining candidate genes were input into STRING (https://string-db.org/) to construct a protein-protein interaction (PPI) network using Cytoscape v3.2 (20).

To verify whether the identified genes were associated with tumor progression and prognosis, the association between each gene and survival was determined using the R survival package v2.41-3 (https://cran.r-project.org/web/packages/survminer/index.html). Clinical and RNA-sequencing data from 417 patients with melanoma were downloaded from The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/) using the TCGA biolinks package in R v2.12.3 (https://git.bioconductor.org/packages/TCGAbiolinks). Overall survival was analyzed using the log-rank test. In addition, the ggpubr package v0.2.1 (http://cran.r-project.org/web/packages/ggpubr/index.html) was used to demonstrate the mRNA expression of hub genes in primary and metastatic tumor, and the two groups were compared by Student's t-test. Receiver operating characteristic (ROC) curve and area under the curve (AUC) values were obtained using the pROC package v1.15.0 (http://cran.r-project.org/web/packages/ROCR) to evaluate the efficiency of the genes in distinguishing metastatic and non-metastatic tumors.

Screening candidates for treatment

Drug-target information of Food and Drug Administration (FDA)-approved drugs was obtained from DrugBank (https://www.drugbank.ca/). The exclusion of drugs that had no known targets in the interactome resulted in a total of 1,269 unique drugs and 1,185 targets selected for further analysis. Notably, only pharmacological targets (‘Targets’ section in DrugBank), excluding enzymes, carriers and transporters typically shared among different drugs, were considered. The protein interaction information was obtained from a previously published study, which contained data from 15 databases (21). Among these, 15,969 nodes and 217,160 mutual relationships were identified in the PPI networks. The prognostic genes of melanoma were mapped to the PPI network. The Igraph package v1.2.4.1 (https://igraph.org/) was used to estimate the shortest distance between each target and a particular prognostic gene for each FDA-approved drug (21). Standardization-based approximation indicated that lower values were associated with an increased likelihood that the drug may act on melanoma and prevent its progression.

Results

Weighted co-expression network construction and key module identification

Following a cluster analysis of all samples, one sample in GSE65904 was removed from subsequent analysis due to bias (Fig. 1; Table I). To ensure a scale-free network, it must satisfy R2> 0.8 (Fig. 2A), and the mean connectivity should be conserved as much as possible (Fig. 2B). Furthermore, the degree distribution of nodes in the co-expression network was investigated further and the degree of nodes conforms to power law distribution (Fig. 2C and D). The WGCNA package in R was used to place genes with similar expression patterns into modules through average linkage clustering; a total of 15 modules were identified (Fig. 3A). The black module exhibited the strongest association with tumor metastasis-free survival and disease-specific death survival (Fig. 3), whereas the turquoise module exhibited the strongest association with tumor stage. Therefore, these two modules were considered to be clinically significant and were selected for further analysis.

Table I.

Summary of number and corresponding GSM in GSE65904.

Number	GSM
1	GSM1608593
2	GSM1608594
3	GSM1608595
4	GSM1608596
5	GSM1608597
6	GSM1608598
7	GSM1608599
8	GSM1608600
9	GSM1608601
10	GSM1608602
11	GSM1608603
12	GSM1608604
13	GSM1608605
14	GSM1608606
15	GSM1608607
16	GSM1608608
17	GSM1608609
18	GSM1608610
19	GSM1608611
20	GSM1608612
21	GSM1608613
22	GSM1608614
23	GSM1608615
24	GSM1608616
25	GSM1608617
26	GSM1608618
27	GSM1608619
28	GSM1608620
29	GSM1608621
30	GSM1608622
31	GSM1608623
32	GSM1608624
33	GSM1608625
34	GSM1608626
35	GSM1608627
36	GSM1608628
37	GSM1608629
38	GSM1608630
39	GSM1608631
40	GSM1608632
41	GSM1608633
42	GSM1608634
43	GSM1608635
44	GSM1608636
45	GSM1608637
46	GSM1608638
47	GSM1608639
48	GSM1608640
49	GSM1608641
50	GSM1608642
51	GSM1608643
52	GSM1608644
53	GSM1608645
54	GSM1608646
55	GSM1608647
56	GSM1608648
57	GSM1608649
58	GSM1608650
59	GSM1608651
60	GSM1608652
61	GSM1608653
62	GSM1608654
63	GSM1608655
64	GSM1608656
65	GSM1608657
66	GSM1608658
67	GSM1608659
68	GSM1608660
69	GSM1608661
70	GSM1608662
71	GSM1608663
72	GSM1608664
73	GSM1608665
74	GSM1608666
75	GSM1608667
76	GSM1608668
77	GSM1608669
78	GSM1608670
79	GSM1608671
80	GSM1608672
81	GSM1608673
82	GSM1608674
83	GSM1608675
84	GSM1608676
85	GSM1608677
86	GSM1608678
87	GSM1608679
88	GSM1608680
89	GSM1608681
90	GSM1608682
91	GSM1608683
92	GSM1608684
93	GSM1608685
94	GSM1608686
95	GSM1608687
96	GSM1608688
97	GSM1608689
98	GSM1608690
99	GSM1608691
100	GSM1608692
101	GSM1608693
102	GSM1608694
103	GSM1608695
104	GSM1608696
105	GSM1608697
106	GSM1608698
107	GSM1608699
108	GSM1608700
109	GSM1608701
110	GSM1608702
111	GSM1608703
112	GSM1608704
113	GSM1608705
114	GSM1608706
115	GSM1608707
116	GSM1608708
117	GSM1608709
118	GSM1608710
119	GSM1608711
120	GSM1608712
121	GSM1608713
122	GSM1608714
123	GSM1608715
124	GSM1608716
125	GSM1608717
126	GSM1608718
127	GSM1608719
128	GSM1608720
129	GSM1608721
130	GSM1608722
131	GSM1608723
132	GSM1608724
133	GSM1608725
134	GSM1608726
135	GSM1608727
136	GSM1608728
137	GSM1608729
138	GSM1608730
139	GSM1608731
140	GSM1608732
141	GSM1608733
142	GSM1608734
143	GSM1608735
144	GSM1608736
145	GSM1608737
146	GSM1608738
147	GSM1608739
148	GSM1608740
149	GSM1608741
150	GSM1608742
151	GSM1608743
152	GSM1608744
153	GSM1608745
154	GSM1608746
155	GSM1608747
156	GSM1608748
157	GSM1608749
158	GSM1608750
159	GSM1608751
160	GSM1608752
161	GSM1608753
162	GSM1608754
163	GSM1608755
164	GSM1608756
165	GSM1608757
166	GSM1608758
167	GSM1608759
168	GSM1608760
169	GSM1608761
170	GSM1608762
171	GSM1608763
172	GSM1608764
173	GSM1608765
174	GSM1608766
175	GSM1608767
176	GSM1608768
177	GSM1608769
178	GSM1608770
179	GSM1608771
180	GSM1608772
181	GSM1608773
182	GSM1608774
183	GSM1608775
184	GSM1608776
185	GSM1608777
186	GSM1608778
187	GSM1608779
188	GSM1608780
189	GSM1608781
190	GSM1608782
191	GSM1608783
192	GSM1608784
193	GSM1608785
194	GSM1608786
195	GSM1608787
196	GSM1608788
197	GSM1608789
198	GSM1608790
199	GSM1608791
200	GSM1608792
201	GSM1608793
202	GSM1608794
203	GSM1608795
204	GSM1608796
205	GSM1608797
206	GSM1608798
207	GSM1608799
208	GSM1608800
209	GSM1608801
210	GSM1608802
211	GSM1608803
212	GSM1608804
213	GSM1608805
214	GSM1608806

Figure 2.

Determination of soft-thresholding power in weighted gene co-expression network analysis. (A) Scale-free fit index of various soft-thresholding powers. (B) Mean connectivity of various soft-thresholding powers. (C) Histogram of connectivity distribution at β=3. (D) Scale-free topology at β=3. β, soft thresholding power; k, connectivity.

GO and KEGG pathway enrichment analysis

The genes in the clinically significant modules were categorized into functional groups: BP, CC and MF. The genes in the black module were mainly enriched in ‘antigen processing and presentation’, ‘antigen processing and presentation of peptide antigen’ and ‘antigen processing and presentation of exogenous peptide antigen’ in the BP group, ‘endocytic vesicle membrane’, ‘Golgi-associated vesicle’ and ‘COPII-coated ER to Golgi transport vesicle’ in the CC group, and ‘amide binding’, ‘peptide binding’ and ‘antigen binding’ in the MF group (Fig. 4A). The results of the KEGG pathway analysis demonstrated that genes in the black module were mainly involved in ‘antigen processing and presentation’, ‘viral myocarditis’, ‘cell adhesion molecules cams’ and ‘allograft rejection’, among others (Fig. 5A).

Figure 4.

Gene Ontology analysis of the genes in the black and turquoise modules. (A) Gene Ontology analysis of the genes in the black module. (B) Gene Ontology analysis of the genes in the turquoise module. BP, biological process; CC, cellular component; MF, molecular function.

Figure 5.

KEGG analysis of the genes in the black and turquoise modules. (A) KEGG analysis of the genes in the black module. (B) KEGG analysis of the genes in the turquoise module. KEGG, Kyoto Encyclopedia of Genes and Genomes.

The genes in the turquoise module were mainly enriched in ‘leukocyte differentiation’, ‘T cell activation’ and ‘regulation of lymphocyte activation’ in the BP group, ‘cell leading edge’, ‘lamellipodium’ and ‘cytoplasmic side of plasma membrane’ in the CC group and ‘nucleoside-triphosphatase regulator activity’, ‘GTPase regulator activity’ and ‘phospholipid binding’ in the MF group (Fig. 4B). The results of the KEGG pathway analysis demonstrated that the genes in the turquoise module were mainly involved in ‘chemokine signaling pathway’, ‘B cell receptor signaling pathway’ and ‘T cell receptor signaling pathway’, among others (Fig. 5B).

Genes with |MM+GS|=5% in the black and turquoise modules were selected as candidate prognostic genes, and all other genes were removed. A PPI network of all genes in the black and turquoise modules was constructed using Cytoscape. The network comprised 222 nodes and 1,416 edges according to the STRING database (Fig. 6). Among those, C-X-C motif chemokine receptor 4 (CXCR4), interleukin 7 receptor (IL7R) and phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit γ (PIK3CG) were positively associated with overall survival (Fig. 7D-F). Based on TCGA data, the expression levels of CXCR4, IL7R and PI3KG were upregulated in primary tumors compared with metastatic tumors (Fig. 7A-C). In addition, the ROC curves indicated that CXCR4, IL7R and PI3KG exhibited excellent efficacy for diagnosing primary and metastatic tumor tissues (Fig. 7G-I).

Figure 6.

Protein-protein interaction network of genes in the black and turquoise modules. The size of the circle represents the degree of the node; lines indicate interactions between genes.

Figure 7.

Validation of the expression of hub genes in primary and metastatic melanoma using TCGA database. Expression levels of (A) CXCR4, (B) IL7R and (C) PIK3CG in primary and metastatic melanoma. (D-F) Survival analysis of the hub genes in TCGA dataset for (D) CXCR4, (E) IL7R and (F) PIK3CG. Red lines represent low expression of the hub genes; blue lines represent high expression. (G-I) Receiver operating characteristic curves and AUC statistics were calculated to evaluate the capacity of distinguishing primary and metastatic melanoma of (G) CXCR4, (H) IL7R and (I) PIK3CG. TGCA, The Cancer Genome Atlas; CXCR4, C-X-C motif chemokine receptor 4; IL7R, interleukin 7 receptor; PIK3CG, phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit γ; AUC, area under the curve.

Using genes which were identified as hub nodes in the PPI network (degree >30) associated with prognosis (P<0.05) and metastasis (AUC >0.7) as potential targets in the drug-gene interaction analysis (Fig. 8), the top 15 drugs ranked by the proximity of genes and drugs were screened as possible treatments for melanoma. The screened drugs could be divided into several major categories, including tyrosine kinase inhibitors (TKIs), vascular endothelial growth factor receptor (VEGFR) inhibitors, estrogen receptor modulators, proteasome inhibitors, Burton's tyrosine kinase (BTK) inhibitors and Raf kinase inhibitors. The top 15 drugs are: Ponatinib, nintedanib, tamoxifen, framycetin, regorafenib, dasatinib, sunitinib, bosutinib, benzylpenicilloyl polylysine, ibrutinib, pazopanib, methyl aminolevulinate, bortezomib, sorafenib, lenvatinib.

Figure 8.

Drug-gene interaction network. Network analysis of hub genes and their target drug networks in melanoma. Purple squares represent predicted drugs that may treat melanoma, pink dots are drug targets, green dots are genes associated with tumor development that are also drug targets, and gray dots are genes associated with tumor development that are not targets of drugs. A purple line indicates an interaction between a drug and a target, and a gray line indicates an interaction between genes.

Discussion

Among skin tumors, melanoma is the most malignant (22). High recurrence and metastasis rates affect the efficacy of melanoma treatment (23). The effects of conventional chemotherapy, immunotherapy and targeted therapy remain limited. Thus, identifying novel molecular targets and exploring therapeutic drugs for melanoma is important. In the present study, the GEO database was used to obtain genetic and clinical information from patients with melanoma, construct a co-expression network, select the most significant module and identify three hub genes: CXCR4, IL7R and PIK3CG. TCGA, which was used for further verification, revealed that three aforementioned specific molecules: CXCR4, IL7R and PIK3CG identified in melanoma tissues were associated with prognosis and metastasis. In addition, the top 15 drugs ranked by the proximity of genes and drugs were screened using a network screening method, and a drug-gene network was constructed.

CXCR4, which is a receptor of C-X-C motif chemokine 12 (CXCL12), is located on the surface of >23 human tumors, for example breast cancer, ovarian cancer, glioma, pancreatic cancer and prostate cancer (24). CXCL12 binds to CXCR4, which activates several extra- and intracellular signaling pathways, including the nuclear factor κB, Ca2+-dependent protein tyrosine kinase 2β, PI3K-Akt and mitogen-activated protein kinase signaling pathways (25). In various types of cancer, such as oral (26), esophageal (27), gastric, colon, liver, pancreatic, thyroid and ovarian cancer (28), and leukemia (29), CXCR4 expression is strongly associated with chemotaxis, invasion, angiogenesis and cell proliferation, all of which are involved in tumorigenesis and cancer. However, the results of the present study indicated that, compared with primary tumors, CXCR4 is downregulated in metastatic tumors, and is therefore associated with good prognosis in patients with melanoma. Mitchell et al (30) demonstrated that most of melanoma cases with mitosis, ulceration and regression were CXCR4-negative. Patients with American Joint Committee on Cancer (AJCC) stage (31) I and II melanoma exhibit higher expression of CXCR4 compared with those with AJCC stages III and IV, and a proportion of patients with AJCC stage III–IV melanoma are CXCR4-negative (30). Therefore, the role of CXCR4 as a biomarker warrants further investigation.

IL7R, which is expressed in immune cells, is crucial for the survival, development and homeostasis of the immune system (32). IL-7Rα activates Janus kinases 1 and 3, promoting the function of signal transducer and activator of transcription 5, which leads to the modulation of gene expression, as well as the activation of anti-apoptotic and pro-survival signaling pathways (33). Thus, IL7R is classified as an oncogene associated with several tumors, including esophageal and prostate cancer (34). However, a bioinformatics study has demonstrated that patients with colon cancer lacking IL7R (two cases of mortality out of three cases) had a median survival time of 34 months compared with patients with normal IL7R status, whose survival time was 45 months (35). Studies on the association between IL7R and melanoma, as well as the association between IL7R and metastasis, are lacking.

The PI3K signaling pathway modulates various biological processes, including cell proliferation, survival, motility, death and metabolism (36). Aberrations in these processes are pivotal for the pathogenesis of cancer. Based on structural differences, PI3K can be divided into several subunits, including PIK3CA, PIK3CB, PIK3CD and PIK3CG (37). A previous study has revealed that PIK3CG is expressed at undetectable levels in glioblastoma cells, and that blocking this specific subunit does not cause cytotoxicity (38). Another study has demonstrated that PIK3CG is downregulated in colorectal cancer, whereas 12 other genes in the PI3K-AKT signaling pathway are upregulated (39). However, a bioinformatics-based study reported that PIK3CG is significantly associated with melanoma metastasis to regional lymph nodes, which contradicted the results of the present study, suggesting that further investigation may be required to clarify the role of PIK3CG in the metastasis of melanoma (40).

In the present study, the GEO database, which comprised 214 melanoma samples, and TCGA database, which included 417 patients, were selected to verify the roles of the identified genes. Double validation and a large number of samples contributed to the reliability of the candidate genes. However, a limitation of the present study was a lack of clinical or experimental validation. Further study is required to verify the role of CXCR4, IL7R and PI3K3CG in melanoma.

The analysis of the association between genes and FDA-approved drugs demonstrated that the top 15 drugs were TKIs, VEGFR inhibitors, estrogen receptor modulators, proteasome inhibitors, Bcr-Abl kinase inhibitors, BTK inhibitors, Raf kinase inhibitors, framycetin, benzylpenicilloyl polylysine and methyl aminolevulinate. TKIs that function by blocking the Bcr-Abl tyrosine-kinase included dasatinib, ponatinib and bosutinib, which are used to treat chronic myelogenous leukemia and acute lymphocytic leukemia (41). Other drugs, including nintedanib, regorafen, sunitinib, pazopanib, sorafenib and lenvatinib inhibit several receptor tyrosine kinases, including platelet-derived growth factors, VEGFR, fibroblast growth factor receptors and Raf family kinases, which inhibit tumor angiogenesis and tumor cell proliferation (42). Ibrutinib, a BTK inhibitor, is used to treat chronic lymphocytic leukemia (43). Tamoxifen, a selective estrogen receptor modulator, is used for the treatment and prevention of estrogen receptor-positive breast cancer (44). Bortezomib was the first therapeutic proteasome inhibitor to be tested in humans; it serves a role in cell cycle arrest and apoptosis, and is approved in the United States for the treatment of relapsed multiple myeloma and mantle cell lymphoma (45). Framycetin, which is an antibiotic, is used to treat leg ulcers and other conditions associated with wound healing (46). Benzylpenicilloyl polylysine is used as a skin-testing reagent for individuals with a history of penicillin allergy (47). Methyl aminolevulinate, which is metabolized into phototoxic compounds, such as protopophyrin IX, may represent a candidate for photodynamic therapy, as it can induce oxidative damage to the cell (48).

Angiogenesis is a hallmark of several types of tumor, including melanoma. The process of angiogenesis is crucial for tumor development and metastasis (49). VEGF is one of the most important cytokines responsible for tumor-mediated angiogenesis (50). VEGF is strongly expressed in melanoma and serves a critical role in the progression of the disease (51). In a phase II study of sunitinib in patients with advanced melanoma, 4/31 (13%) patients exhibited a partial response and 8 (26%) had stable disease (52). Pazopanib, a VEGF and platelet-derived growth factor inhibitor, has been used in combination with paclitaxel in a phase II study of patients with metastatic melanoma; the 6-month progression-free survival rate was 68%, and the 1-year overall survival rate was 48% (53). The objective response rate was 37%, comprising one complete and 20 partial responses (54). A phase Ib study using lenvatinib (E7080) in combination with temozolomide for the treatment of advanced melanoma indicated an overall objective response rate of 18.8% (six patients), comprising all partial responses (55). SRC proto-oncogene, non-receptor tyrosine kinase (SRC) is a promising target in the treatment of solid types of cancer, including human melanoma; bosutinib, a SRC inhibitor, which induces cell death via lysosomal membrane permeabilization in melanoma cells, is a promising therapeutic agent for melanoma treatment (56). SRC inhibitor Dasatinib specifically inhibits p53 phosphorylation in melanoma; however, a comprehensive validation is required (57). Ibrutinib, a BTK inhibitor, has been used to treat chronic lymphocytic leukemia/small lymphocytic lymphoma and subsequent melanoma that occurs following leukemia (58). Sorafenib, a Raf inhibitor, is a first-line therapeutic agent used in advanced melanoma (phase I and open-label phase II) trials with an overall response rate of 12% with one complete response and nine partial responses (59). Bortezomib administration reduces the levels of proangiogenic cytokines in plasma (60). A clinical trial has indicated that tamoxifen therapy is not effective for treating melanoma, and that the mode of action of antiestrogens in melanoma is unclear (61). To the best of our knowledge, ponatinib, nintedanib, regorafen, framycetin, benzylpenicilloyl polylysine and methyl aminolevulinate have not been used to treat melanoma.

The present study used WGCNA to construct a gene co-expression network in order to determine the associations between genes and modules and to explore the association between the gene modules and clinical characteristics. Two significant modules (black and turquoise modules) shown in Fig. 3, were identified to be associated with the progression of melanoma. GO and KEGG pathway analyses demonstrated that this module was mostly involved in functions associated with antigen presentation. In addition, three hub genes, CXCR4, IL7R and PI3K3CG, were identified and demonstrated to be associated with the progression and prognosis of melanoma. Analysis of the interaction between genes and drug targets of the top 15 drugs for melanoma enabled the construction of a network of drug-gene interactions. Ponatinib, regorafen, nintedanib, framycetin, benzyl penicilloyl polylysine and methyl aminolevulinate, which were among the 15 drugs not currently used to treat melanoma, may be potential novel therapeutic drugs for this disease.