2.1Screening of differentially expressed genes
Volcanic maps of gene expression were plotted for the mRNA datasets GES114445 and GSE7553, and the microRNA datasets GSE18512, GSE35579, and GSE24996, respectively (Figs. 1A–E). Data were analysed using GEO2R with screening conditions at values p < 0.05 and |logFC| > 1. In total, the GES114445 data set comprised 839 differentially expressed genes, including 528 up-regulated and 311 down-regulated genes. Moreover, the GSE7553 data set consisted of 1776 differentially expressed genes, including 551 upregulated and 1225 down-regulated genes. Furthermore, the GSE18512 data set comprised 22 differentially expressed genes, including 14 up-regulated and 8 down-regulated genes. Additionally, the GSE35579 data set consisted of 52 differentially expressed genes, including 18 up-regulated and 34 down-regulated genes. Eventually, the GSE24996 data set comprised 48 differentially expressed genes, including 20 up-regulated and 28 down-regulated genes.
An article entitled “Epigenetic silencing of CDR1as drives IGF2BP3-mediated melanoma invasion and metastasis” was retrieved from PubMed, and the circRNA differential expression data (the GSE138711 chip set from the GEO database) were downloaded. In total, 688 differentially expressed circRNAs were identified, including 572 differentially down-regulated and 116 up-regulated genes. The 10 most differentially expressed circRNAs were screened as the target molecules for subsequent research. Among these, the up-regulated circRNAs were LIFR, ZC3H13, ARID1A, XRN2, and ITCH, whereas the down-regulated circRNAs were SIPA1L1, WDR37, PAG1, ARHGEF12, and CCSER2 (Fig. 2, Table 1).
The intersection of the above-mentioned differential genes was obtained by drawing Venn diagrams of 283 and 26 common differentially expressed mRNAs and microRNAs, respectively (Figs. 3A, B).
Venn of differentially expressed gene intersection. A the intersection of GSE114445 and GSE7553. B the intersection of GSE18512、GSE35579 and GSE24996 .(the genes contained in at least two data sets were selected as the target genes)
2.2Acquisition of Hub genes
A network of 283 differentially expressed mRNAs was constructed in the STRING database to obtain a PPI network with 184 nodes and 490 edges; thereafter, the PPI network file was downloaded. Excel was used to collate the data in the network file, and two files ‘Node to Node Network’ and ‘degree’, were obtained. These files were imported into Cytoscape software and the network was enhanced, thereby obtaining the PPI network diagram as illustrated in Fig. 4. Moreover, 59 genes with degree value > mean degree value were screened as the Hub genes, and a sub-network was established (Fig. 5).
The subnetwork of Hub genes. Icons with different colors represent the degree value of the gene. The ligther the color, the smaller the degree value of the gene; Different icons represent the expression trend of the gene, the triangle the down-regulation of the gene expression in skcm, and the arrow symbol represents the up-regulation of the gene expression in skcm; The letters on the icon represest the name of the gene.
2.3GO and KEGG enrichment analysis of Hub genes
The GO and KEGG enrichment analyses were performed for 59 Hub genes in the DAVID database with screening conditions of p < 0.05 and count value > mean count value. Biological effects were mainly enriched in inflammatory response, immune response, innate immune response, positive regulation of ERK1 and ERK2, and drug response. Molecular functions were mainly enriched in protein binding. Moreover, cell localisation was mainly enriched in plasma membrane, extracellular exosomes, components of plasma membrane, cell surface, and desmosomes. Furthermore, the signalling pathways were mainly enriched in Staphylococcus aureus infection, osteoclast differentiation, cytotoxicity mediated by natural killer cells, influenza A, and in the chemokine signalling pathway. Figure 6A and 6B illustrate the results of GO analysis with more than three genes and the KEGG pathway enrichment analysis, respectively.
The results of GO and KEGG enrichment analysis. A the results of GO enrichment analysis(The GO_BP,GO_MF and GO_CC were selected with thresholds of P < 0.05 and count > 3);B the results of KEGG enrichment analysis༈The KEGG signaling pathways were selected with threshold of P < 0.05༉
2.4Survival analysis for Hub genes and identification of key genes
In order to screen the Hub genes and clarify their significance in melanoma, the 59 Hub genes screened above were imported into the GEPIA database for survival analysis (P < 0.05). In total, 31 key genes correlating with the overall survival rate of melanoma were screened, among which 25 genes were differentially up-regulated and 6 genes were differentially down-regulated (Table 1, Figs. 7).
Table 1
Correlation between mRNA and survival rate of SKCM
gene name
|
variation trend
|
P value
|
FBXW7
|
Down
|
3.30E-02
|
DSG1
|
Down
|
9.30E-03
|
EVPL
|
Down
|
2.60E-05
|
PKP1
|
Down
|
3.50E-04
|
DSC3
|
Down
|
6.50E-03
|
EPPK1
|
Down
|
6.10E-03
|
CXCR4
|
Up
|
2.10E-02
|
AURKA
|
Up
|
1.40E-02
|
CYBB
|
Up
|
3.40E-05
|
CCR1
|
Up
|
5.70E-05
|
CCL5
|
Up
|
4.50E-06
|
ICAM1
|
Up
|
2.60E-04
|
CD86
|
Up
|
1.40E-06
|
LCP2
|
Up
|
2.00E-06
|
TLR7
|
Up
|
4.00E-04
|
TYROBP
|
Up
|
1.30E-03
|
FCGR2A
|
Up
|
4.30E-06
|
C1QC
|
Up
|
9.30E-06
|
CD53
|
Up
|
5.00E-05
|
FCER1G
|
Up
|
9.90E-05
|
KPNA2
|
Up
|
1.80E-03
|
LAPTM5
|
Up
|
1.50E-06
|
DTL
|
Up
|
1.00E-03
|
FPR3
|
Up
|
2.50E-03
|
FYB
|
Up
|
6.20E-05
|
NCF1
|
Up
|
3.90E-06
|
SERPINA1
|
Up
|
2.90E-04
|
HAVCR2
|
Up
|
8.70E-06
|
SRGN
|
Up
|
8.20E-07
|
FCGR1B
|
Up
|
8.70E-06
|
GBP5
|
Up
|
3.20E-07
|
Correlation between mRNA and survival rate of SKCM. The mRNAs were selected with threshold of Logrank P < 0.05
2.5Construction of the ceRNA network
The miRNA–mRNA column in the starBase database was used to predict the corresponding relationships between the miRNAs and key genes, and the co-expression analysis of miRNAs and their corresponding mRNAs was performed in the starBase database. Based on the contradictory changing trends of miRNAs and mRNAs in the ceRNA network, the screening condition (coefficient-R < 0) was set for filtration, and the microRNA–mRNA docking relationships were obtained (Table 2).
In the miRNA–circRNA column of starBase database, the corresponding relationships of miRNAs with circRNA LIFR, ZC3H13, ARID1A, XRN2, ITCH, SIPA1L1, WDR37, PAG1, ARHGEF12, and CCSER2, were predicted. As circRNA often acts as a miRNA sponge in the ceRNA network to competitively inhibit the binding of miRNAs and mRNAs, the expression trends of circRNAs in the circRNA–miRNA–mRNA relationship are consistent with those of mRNAs, but contradictory to those of miRNAs. Based on the above-mentioned principles, we screened the following miRNA–circRNA corresponding relationships (Table 3). The miRNA–mRNA pairing relationships were integrated with the miRNA–circRNA pairing relationships, thereby obtaining the ceRNA network diagram (Fig. 8).
Table 2
Correlation between mRNA and miRNA
mRNA
|
Variation trend
|
microRNA
|
Coeffient-R
|
P value
|
FBXW7
|
Down
|
Hsa-miR-155
|
-0.071
|
1.31E-01
|
|
|
Hsa-miR-21
|
-0.156
|
8.90E-04
|
|
|
Hsa-miR-196a
|
-0.196
|
2.77E-05
|
|
|
Hsa-miR-425
|
-0.225
|
1.53E-06
|
|
|
Hsa-miR-142
|
-0.019
|
6.89E-01
|
DSG1
|
down
|
Hsa-miR-130b
|
-0.118
|
1.26E-02
|
|
|
Hsa-miR-142
|
-0.054
|
2.58E-01
|
PKP1
|
down
|
Hsa-miR-155
|
-0.017
|
7.23E-01
|
|
|
Hsa-miR-424
|
-0.078
|
9.92E-02
|
|
|
Hsa-miR-185
|
-0.043
|
3.60E-01
|
|
|
Hsa-miR-425
|
-0.024
|
6.14E-01
|
DSC3
|
down
|
Hsa-miR-106b
|
-0.107
|
2.35E-02
|
|
|
Hsa-miR-130b
|
-0.107
|
2.36E-02
|
|
|
Hsa-miR-155
|
-0.041
|
3.88E-01
|
|
|
Hsa-miR-20b
|
-0.065
|
1.67E-01
|
|
|
Hsa-miR-142
|
-0.011
|
8.20E-01
|
EPPK1
|
down
|
Hsa-miR-425
|
-0.132
|
5.11E-03
|
|
|
Hsa-miR-424
|
-0.105
|
2.62E-02
|
|
|
Hsa-miR-185
|
-0.095
|
4.38E-02
|
ICAM1
|
Up
|
Hsa-miR-214
|
-0.095
|
4.35E-02
|
|
|
Hsa-miR-125b
|
-0.057
|
2.30E-01
|
|
|
Hsa-miR-125a
|
-0.087
|
6.56E-02
|
LAPTM5
|
Up
|
Hsa-miR-185
|
-0.061
|
1.95E-01
|
DTL
|
Up
|
Hsa-miR-204
|
-0.104
|
2.70E-02
|
|
|
Hsa-miR-214
|
-0.021
|
6.50E-01
|
AURKA
|
Up
|
Hsa-miR-23b
|
-0.034
|
4.68E-01
|
|
|
Hsa-Let-7b
|
-0.115
|
1.44E-02
|
CXCR4
|
up
|
Hsa-miR-211
|
-0.238
|
3.39E-07
|
CYBB
|
up
|
Hsa-miR-196b
|
-0.110
|
2.02E-02
|
FAPTM5
|
up
|
Hsa-miR-27b
|
-0.042
|
3.79E-01
|
ceRNA Network Diagram. This figure clearly expressed the interaction among differentially expressed circRNA, microRNA and mRNA. The rectangular icon on the left represents differentially expressed mRNA; the oval icon on the right represesnts differentially expressed circRNA; the diamond icon in the middle represesnts the expressed microRNA. The color of the icon represents the expression trend of the gene in melanoma; Red is up-regulated and green is down-regulated.
Table 3
miRNA-circRNA docking relationship
microRNA
|
circRNA
|
Hsa-miR-211
|
ZC3H13
|
Hsa-miR-204
|
ZC3H13
|
Hsa-miR-214
|
ARID1A
|
Hsa-miR-185
|
ARID1A
|
Hsa-miR-125b
|
ARID1A
|
Hsa-miR-125a
|
ARID1A
|
Hsa-Let-7b
|
XRN2
|
Hsa-miR-27b
|
ITCH
|
Hsa-Let-7b
|
ITCH
|
Hsa-miR-23b
|
ITCH
|
Hsa-miR-214
|
ITCH
|
Hsa-miR-425
|
SIPA1L1
|
Hsa-miR-424
|
SIPA1L1
|
Hsa-miR-185
|
SIPA1L1
|
Hsa-miR-106b
|
SIPA1L1
|
Hsa-miR-20b
|
SIPA1L1
|
Hsa-miR-425
|
ARHGEF12
|
Hsa-miR-424
|
ARHGEF12
|
Hsa-miR-106b
|
ARHGEF12
|
Hsa-miR-20b
|
ARHGEF12
|
2.6ceRNA network verification
The results of mRNA and miRNA co-expression in the ceRNA network were retrieved from the starBase database. The conditions for ceRNA network verification were as follows: (1) P < 0.05 and (2) at least two common miRNAs between the mRNAs and circRNAs; satisfying one of the above conditions was necessary. The following target ceRNA corresponding relationships were screened based on the above conditions (Table 4).
The target ceRNAs and their corresponding GO/KEGG enrichment analysis results were plotted in Cytoscape (Fig. 9).
Differential analysis of mRNAs from the target ceRNA network was carried out in the GEPIA database, and the pan-cancer expression patterns obtained are illustrated in Fig. 10. Compared with other tumours, the expression of DSC3, PKP1, and EPPK1 in melanoma and normal skin tissue was specific. The correlation between target mRNA in the ceRNA network and the staging of melanoma was analysed in the GEPIA database (Fig. 11). Among FBXW7 (P = 0.0402), PKP1 (P = 0.000201), DSC3 (P = 0.000175), EPPK1 (P = 0.0243), ICAM (P = 0.423), DTL (P = 0.102), AURKA (P = 0.0534), and CXCR4 (P = 4.94e-06), only FBXW7, PKP1, DSC3, EPPK1, and CXCR4 were correlated with the tumour staging in melanoma (P < 0.05).
Survival analysis of miRNAs in the target ceRNA network was carried out in the UALCAN database (Fig. 12). Among has-miR-106b (P = 0.049), has-miR-20b (P = 0.027), has-miR-185 (P = 0.65), has-miR-214 (P = 0.84), has-miR-125b (P = 0.086), has-miR-204 (P = 0.056), has-miR-211 (P = 0.19), has-miR-425 (P = 0.60), has-miR-424 (P = 0.79), and has-let-7b (P = 1), only has-miR-106b and has-miR-20b were correlated with the overall survival rate of melanoma (P < 0.05).
Eventually, we screened the SIPA1L1-has-miR-106b/has-miR-20b-DSC3 regulatory axis and presumed that SIPA1L1 can regulate DSC3 expression by competitively binding to HAS-miR-106b and HAS-miR-20b, thereby affecting the development and progression of melanoma.
Table 4
mRNA-miRNA-circRNA docking relationship
序号
|
circRNA
|
miRNA
|
mRNA
|
1
|
SIPA1L1
|
hsa-miR-425
|
FBXW7
|
2
|
SIPA1L1
|
Hsa-miR-424/hsa-miR-425
|
PKP1
|
3
|
SIPA1L1
|
Hsa-miR-106b/hsa-miR-20b
|
DSC3
|
4
|
SIPA1L1
|
Hsa-miR-425/hsa-miR-424
|
EPPK1
|
5
|
ARID1A
|
Hsa-miR-214/hsa-miR-125b/hsa-miR-125a
|
ICAM1
|
6
|
ZC3H13
|
Hsa-miR-214
|
DTL
|
7
|
XRN2
|
Hsa-let-7b
|
AURKA
|
8
|
ZC3H13
|
Hsa-miR-211
|
CXCR4
|
Survival analysis of miRNA in skin cutaneous melanoma. It can be seen from the figure that has-miR-106b and has-miR-20b have related to the survival rate of skin cutaneous melanoma(The miRNAs were selected with threshold of P < 0.05).