J Cancer 2021; 12(10):2921-2932. doi:10.7150/jca.51851 This issue Cite

Research Paper

Bioinformatics analysis of lncRNA‑associated ceRNA network in melanoma

Yi Ding1, Min Li1, Tuersong Tayier2, MeiLin Zhang3, Long Chen1, ShuMei Feng1 Corresponding address

1. Department of Histology and Embryology, School of Basic Medical Sciences, Xinjiang Medical University, Urumqi, Xinjiang, China.
2. Department of Pharmacology, Pharmacy College, Xinjiang Medical University, Urumqi, China.
3. Xinjiang Urumqi City Center Blood Station, Urumqi, China.

Citation:
Ding Y, Li M, Tayier T, Zhang M, Chen L, Feng S. Bioinformatics analysis of lncRNA‑associated ceRNA network in melanoma. J Cancer 2021; 12(10):2921-2932. doi:10.7150/jca.51851. https://www.jcancer.org/v12p2921.htm
Other styles

File import instruction

Abstract

Graphic abstract

Melanoma is an extremely malignant tumor with early metastasis and high mortality. Little is known about the process of by which melanoma occurs, as its mechanism is very complex and only limited data are available on its long non-coding RNA (lncRNA)-associated competing endogenous RNAs (ceRNAs). The purpose of this study was to screen out potential prognostic molecules and identify a ceRNA network related to the occurrence of melanoma. We screened 169 differentially expressed mRNAs (DEmRNAs) from E-MTAB-1862 and GSE3189; gene ontology (GO) enrichment analysis showed that these genes were closely related to the development of skin. In the protein-protein interaction network, we screened out a total of 19 hub genes. Furthermore, we predicted the microRNAs (miRNAs) that regulate hub genes using the miRWalk database and then intersected these with GSE35579, resulting in nine DEmiRNAs. We also predicted the lncRNAs that regulate the miRNAs using the LncBasev.2 database. According to the ceRNA hypothesis, and based on the intersection of the DElncRNAs with merged GTEx and TCGA data, we obtained 20 DElncRNAs. A total of four DEmRNAs, nine DEmiRNAs, and 20 DElncRNAs were included in the ceRNA network. Based on Cox stepwise regression and survival analysis, we identified five biomarkers, ZSCAN16-AS1, LINC00520, XIST, DTL, and let-7a-5p, and obtained risk scores. The results showed that most of the differentially expressed genes were related to epithelial-mesenchymal transition (EMT) in melanoma. Finally, we obtained a LINC00520/let-7a-5p/DTL molecular regulatory network. These results suggest that ceRNA networks have an important role in evaluating the prognosis of patients with melanoma and provide a new experimental basis for exploring the EMT process in the development of melanoma.

Keywords: Bioinformatics analysis, ceRNA, differentially expressed gene, melanoma, survival analysis.

Introduction

Melanoma is the most frequently occurring malignancy in the elderly and has high mortality rates [1]. In recent decades, the number of patients with melanoma has increased with the aging of the population [2]. The development of melanoma is caused by the progression of a benign nevus, sometimes accompanied by sun exposure, physical stimulation, genetic variation, and other factors [3]. Nevi are benign products of normal melanocytes in the skin. In general, nevi will not develop into melanoma [4]. However, constant rubbing or physical stimulation of the skin in the affected area increases the probability of malignant transformation. At present, treatments for melanoma include surgery, radiotherapy, chemotherapy, and immunotherapy, but their efficacy is not optimal [5-7]. The transformation of a normal nevocytic nevus into melanoma is regulated by a complex molecular network. Therefore, it is of great importance to clarify the molecular mechanism of melanoma formation and find new molecular therapeutic targets. This study aimed to unravel the pathogenesis, prognostic factors, and potential molecular therapeutic targets of melanoma. First, we used multiple public databases to screen out differentially expressed genes (DEGs), mRNAs (DEmRNAs), microRNAs (DEmiRNAs), and long non-coding RNAs (DElncRNAs). Then, we constructed a competing endogenous RNA (ceRNA) network based on gene enrichment, target gene prediction, survival analysis, and multivariate Cox regression analysis. Finally, five predictive genes, ZSCAN16-AS1, LINC00520, XIST, DTL, and let-7a-5p were identified; these genes are closely related to the occurrence of melanoma. This study provides new molecular therapeutic targets and lays the foundation for a better understanding of the mechanism of melanoma formation.

Materials and Methods

Data preparation

We downloaded two expression profile datasets, E-MTAB-1862 and GSE3189, from the European Bioinformatics Institute ArrayExpress database (http://www.ebi.ac.uk/arrayexpress) and the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database, respectively. The E-MTAB-1862 dataset consisted of 11 nevus samples and 21 melanoma samples [8], based on the A-AFFY-44-Affymetrix GeneChip Human Genome U133 Plus 2.0 array. The GSE3189 dataset, submitted by Talantov et al., included 18 nevi samples and 45 melanoma samples [9] and was based on the GPL96 [HG-U133A] Affymetrix Human Genome U133A array. In addition, miRNA sequencing data for melanoma was downloaded as GSE35579; this dataset incorporates 11 benign nevi and 20 primary melanoma samples, sequenced on the GPL15183 CRUK/Melton lab-Human melanoma-71-v2-microRNA expression profiling platform. Finally, we obtained lncRNA sequencing data from the merged GTEx and The Cancer Genome Atlas (TCGA) databases (https://www.cancergenome.nih.gov/). The GTEx melanoma data was downloaded from UCSC Xena (http://xena.ucsc.edu/) [10]. As of June 20, 2020, a total of 813 normal samples and 471 melanoma samples were merged.

Identification of DEGs in melanoma

After downloading the mRNA expression profiles of E-MTAB-1862, we used the EdgeR package in R version 4.0.0 to identify DEmRNAs[11]. Then, we used the ggplot2 R package to draw a heat map and volcano plot of the DEmRNAs. Next, we filtered the DEmRNAs from GSE3189 using the GEO2R software [12]. To find overlapping DEmRNAs, we used the Venn software online to obtain the intersection of the two datasets. Moreover, we used GEO2R to screen DEmiRNAs from GSE35579 between nevi tissue and melanoma tissue. Finally, after merging the GTEx and TCGA data using R, we filtered out DElncRNAs with the EdgeR package, and generated the corresponding heat map and volcano map using ggplot2. The DEmRNAs and DEmiRNAs were then filtered using adjusted p<0.05 and |log fold change (FC)|>1, and the DElncRNAs with thresholds of |log FC|>1, p <0.05, and false discovery rate (FDR) < 0.05.

Functional enrichment analysis

GO annotation analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using the STRING (Search Tool from the Retrieval of Interacting Genes; https://string-db.org/) database (version 11.0) [13]. STRING can be used to predict protein-protein interactions (PPI) directly or indirectly. (p < 0.05 and FDR<0.05).

PPI analysis

We used the STRING database to visually evaluate DEmRNAs with a confidence score >0.4. Furthermore, we used the Cytoscape software (version 3.7.2) to visualize the PPI [14]. We used the Molecular Complex Detection (MCODE) application in Cytoscape to screen hub genes of the PPI network with degree >8 (degree of cutoff = 2, node score cutoff = 0.2, k-core = 2, maximal depth = 100).

Construction of ceRNA regulation network in melanoma

We used the miRWalk and miRDB databases to predict the miRNAs that regulate DEmRNAs with mini p-value >0.9 [15]. To find overlapping DEmiRNAs, we used Venn software online to identify the intersection between the predicted DEmiRNAs and those from the GSE35579 dataset. Furthermore, we predicted lncRNAs that regulated the DEmiRNAs using the LncBasev.2 database with threshold >0.9 [16]. Then, we identified the intersection of the DElncRNAs obtained after merging the GTEx and TCGA databases with the predicted lncRNAs. Finally, based on the ceRNA hypothesis, we constructed a ceRNA network and visualized the results using Cytoscape.

Survival analysis of DEGs

We used the GEPIA database (http://gepia.cancer-pku.cn/index.html) to analyze the expression of particular lncRNAs and mRNAs in melanoma [17]. Then, we carried out a survival analysis of specific lncRNAs, miRNAs, and mRNAs using the oncolnchttp website (http://www.oncolnc.org/) [18] with log-rank p <0.05.

Multivariate Cox regression model analysis

Only DEGs identified as meaningful in the survival analysis were considered for further analysis. Based on the clinical information of patients with survival status and duration in the TCGA database. First of all, we rule out those patients who have no survival status. Furthermore, we used the survival R package to perform multivariate Cox regression analysis using the six meaningful DEGs (LINC00520, ZNF561-AS1, XIST, MMP9, DTL, and let-7a-5p). We also performed a stepwise regression analysis of these six variables and obtained risk scores. To evaluate the prognostic value of the final gene combination in patients with melanoma, patients were divided into high- and low-risk groups based on the median risk score. Moreover, survival analysis and receiver operating characteristic (ROC) curve analysis of the final gene combination were carried out using the survival and survivalROC R packages, respectively. Finally, we calculated the area under the curve (AUC). The heat map of the six DEGs was constructed using ggplot2 (p<0.05). To better predict the prognosis of patients, we used melanoma patients with gender, stage, and TNM stages in the TCGA database as inclusion criteria. At the same time, we used a survival package combined with clinical information of patients and the risk model constructed in this experiment for univariate and multivariate COX analyses and obtained the related forest plots. Furthermore, the ROC curve related to clinical information was drawn to evaluate the prognostic value of the model. P-value < 0.05.

Results

Identification of DEGs in melanoma

Through analysis of gene expression profiles in E-MTAB-1862, we screened a total of 367 DEGs, of which 210 were upregulated and 156 were downregulated in melanoma. Then, we used a volcano plot and heatmap to map the top 100 DEmRNAs (Figure 1). We also screened 3090 DEmRNAs from the GSE3189 dataset using GEO2R, of which 1215 were upregulated and 1875 were downregulated in melanoma tissues. Finally, we took the intersection of the two datasets using the Venn software online and found that a total of 169 DEmRNAs appeared in both datasets, of which 74 were upregulated and 95 were downregulated in melanoma samples (Figure 2a). In addition, we analyzed the miRNA sequencing data from the GSE35579 dataset; 91 DEmiRNAs were screened, of which 31 were upregulated and 60 were downregulated. We identified nine miRNAs in the intersection of GSE35579 with the DEmiRNAs predicted from the miRWalk database (Figure 2b). We also screened DElncRNAs by combining the melanoma GTEx and TCGA data, and obtained 684 upregulated and 1160 downregulated DElncRNAs in melanoma. The ggplot2 R package was used to plot the heat map of DElncRNAs (Figure 3).

 Figure 1 

Distribution of DEmRNAs in melanoma. (A) Heat map of the first 100 DEmRNAs. (B) Volcano plot of DEmRNAs. Red indicates upregulation, green indicates downregulation, and black indicates normal expression in the volcano plot. |log FC|>1, adjusted p<0.05.

J Cancer Image
 Figure 2 

Venn diagram. (A) 169 DEmRNAs; (B) nine DEmiRNAs.

J Cancer Image
 Figure 3 

Heat map of DElncRNAs. |log FC|>1, p<0.05, and FDR< 0.05.

J Cancer Image

GO term and KEGG pathway analysis

Analysis of the top 10 GO enrichment results (Table 1, Figure 4a) showed that the DEmRNAs were mainly involved in tissue development, epidermis development, skin development, and cell differentiation with respect to biological processes. Regarding cellular components, the DEmRNAs were mainly enriched in the extracellular matrix, intercellular bridge, cornified envelope, and cell surface. Regarding molecular functions, the DEmRNAs were significantly involved in proteoglycan binding, structural constituent of the cytoskeleton, carboxylic acid binding, and DNA-binding transcription activator activity. Furthermore, KEGG pathway analysis showed that the DEmRNAs were involved in bladder cancer, ECM-receptor interaction, PPAR signaling pathway, pathways in cancer, IL-17 signaling pathway, and microRNAs in cancer (Table 1, Figure 4b).

PPI analysis

The PPI network contained 169 nodes and 407 edges, including 74 upregulated genes and 95 downregulated genes (Figure 5a). We used Cytoscape to visualize the PPI network. Furthermore, we used MCODE to identify two modules comprising 37 genes. Then, taking degree >8 as the screening index (Figure 5b), we selected 19 genes, of which the CD44 gene had low expression and the other genes had high expression in melanoma.

Constructing the ceRNA regulation network in melanoma

We used the miRWalk database to predict the miRNAs regulating 19 mRNAs (Table 2); these miRNAs were also found in the miRDB database. Then, we took the intersection of the miRNAs in the predicted miRNA database and GSE35579 dataset, resulting in nine miRNAs (Table 3). Furthermore, we used the LncBasev.2 database to predict the lncRNAs regulating the miRNAs (threshold >0.8), and used the Venn software online to intersect the different lncRNAs from the merged GTEx and TCGA data with the predicted lncRNAs. According to the ceRNA hypothesis, there is a positive regulatory relationship between lncRNAs and mRNAs, and a negative regulatory relationship between miRNAs and mRNAs. Finally, we obtained 20 lncRNAs with significant expression differences in melanoma (Table 4). We constructed the ceRNA network using Cytoscape (Figure 6).

Survival analysis of DEGs

We used GEPIA to analyze the expression of specific DEGs (Figure 7a). The results showed that compared with normal nevus tissues, LINC00520, ZSCAN16-AS1, MMP9, and DTL had higher expression in melanoma, whereas XIST had lower expression. Survival analysis using the oncolnchttp website showed that LINC00520, ZSCAN16-AS1, XIST, MMP9, DTL, and let-7a-5p were significantly correlated with overall survival (OS) (Figure 7b). Kaplan-Meier survival analysis showed that the survival time of patients with high expression of LINC00520, ZSCAN16-AS1, MMP9, and DTL was significantly shortened, whereas that of patients with high expression of XIST and let-7a-5p was significantly prolonged (log-rank p<0.05).

 Table 1 

GO and KEGG pathway enrichment analysis of DEmRNAs.

Term IDTerm descriptionObserved countsFDR
Biological processes
GO:0009888Tissue development461.95E-09
GO:0007275Multicellular organism development762.91E-06
GO:0008544Epidermis development192.91E-06
GO:0048856Anatomical structure development802.91E-06
GO:0048513Animal organ development554.84E-06
GO:0030855Epithelial cell differentiation236.93E-06
GO:0043588Skin development171.66E-05
GO:0060429Epithelium development291.66E-05
GO:0030154Cell differentiation592.35E-05
GO:0032501Multicellular organismal process902.35E-05
Cell component
GO:0005576Extracellular region503.20E-06
GO:0044464Cell part1611.60E-04
GO:0031012Extracellular matrix113.50E-03
GO:0071944Cell periphery693.70E-03
GO:0005886Plasma membrane676.40E-03
GO:0045171Intercellular bridge58.50E-03
GO:0001533Cornified envelope51.05E-02
GO:0099081Supramolecular polymer191.05E-02
GO:0009986Cell surface161.15E-02
GO:0005622Intracellular1412.49E-02
Molecular function
GO:0043394Proteoglycan binding59.20E-03
GO:0005200Structural constituent of cytoskeleton71.51E-02
GO:0005504Fatty acid binding41.51E-02
GO:0005515Protein binding821.51E-02
GO:0031406Carboxylic acid binding91.51E-02
GO:0033293Monocarboxylic acid binding51.64E-02
GO:0050786RAGE receptor binding31.64E-02
GO:0001228DNA-binding transcription activator activity121.87E-02
GO:0036041Long-chain fatty acid binding31.87E-02
GO:1990837Sequence-specific double-stranded DNA binding171.87E-02
KEGG
hsa05219Bladder cancer64.80E-04
hsa04512ECM-receptor interaction78.70E-04
hsa03320PPAR signaling pathway62.90E-03
hsa05200Pathways in cancer152.90E-03
hsa04657IL-17 signaling pathway66.50E-03
hsa05206Micrornas in cancer71.09E-02
hsa04610Complement and coagulation cascades51.72E-02
hsa04514Cell adhesion molecules (cams)62.77E-02
hsa04974Protein digestion and absorption52.77E-02
hsa05146Amoebiasis52.77E-02
hsa05222Small cell lung cancer52.77E-02
hsa00590Arachidonic acid metabolism43.20E-02
 Figure 4 

GO and KEGG pathway enrichment analysis of DEmRNAs. (A) GO enrichment analysis. BP, biological process; CC, cellular component; MF, molecular function. (B) KEGG pathway enrichment analysis.

J Cancer Image
 Table 2 

Expression of 19 mRNAs in melanoma.

EXPRESSIONGENESlogFCadj.P.Val
mRNAE-MTAB-1862GSE3189E-MTAB-1862GSE3189
UpregulatedDTL1.301.371.37E-036.89E-08
SHCBP11.221.302.73E-048.49E-07
KPNA21.511.136.81E-055.37E-12
DLGAP51.231.604.92E-031.06E-08
KIAA01011.121.521.29E-043.14E-09
CEP551.282.021.57E-031.39E-11
HMMR1.051.234.56E-034.80E-06
NDC801.402.155.32E-048.88E-07
RRM21.571.582.75E-031.02E-11
EZH21.371.019.25E-051.33E-04
KIF231.273.853.46E-038.07E-14
KIF20A1.161.981.40E-022.15E-10
PBK1.761.201.00E-032.09E-08
PTTG11.021.673.36E-034.74E-15
TPX21.071.525.92E-032.17E-12
NCAPG1.152.455.24E-031.87E-10
CKS21.931.888.21E-046.15E-12
MMP91.491.301.36E-021.10E-03
DownregulatedCD44-1.13-0.832.85E-051.23E-03
 Table 3 

Information about nine miRNAs in the ceRNA network.

miRNAlogFCadj.P.Val
hsa-miR-149-5p-2.344.04E-06
hsa-miR-429-1.654.38E-02
hsa-miR-1290-1.735.03E-03
hsa-miR-141-3p-2.073.18E-03
hsa-miR-200c-3p-1.472.39E-05
hsa-miR-200b-3p-3.432.87E-06
hsa-miR-200a-3p-2.691.21E-05
hsa-miR-20b-5p2.344.62E-04
hsa-let-7a-5p-1.001.57E-02

Multivariate Cox regression analysis

As different marker combinations appeared very promising as prognostic indicators, we further evaluated the prognostic value of six gene combinations in patients with melanoma. Using stepwise regression analysis of six variables, we obtained five variables by Cox regression. Risk scores were calculated as follows, risk score=0.68192*ZSCAN16-AS1+0.14792*LINC00520+0.74931*DTL-3.07155*let-7a-5p-0.12864*XIST. To further evaluate the prognostic value of the five gene combinations in patients with melanoma, patients were divided into high- and low-risk groups according to the median of the Cox regression model. Kaplan-Meier survival analysis showed that the survival time of patients in the high-risk group was significantly shorter than that in the low-risk group (p<0.05). ROC curve analysis showed that a combination of four genes had high accuracy and specificity in evaluating the prognosis of patients with melanoma, with AUC=0.716. The expression of five genes in the high and low-risk groups was visualized using a heat map (Figure 8). Also, univariate COX analysis of the risk model and clinical information showed that both stage and TNM stages could predict the prognosis of patients with melanoma (Supplemental Figure S1A). The forest plot results of multivariate COX analysis showed that T and N could be used as independent prognostic factors (Supplemental Figure S1B). However, the ROC results showed that the AUC of the risk model, T and N was 0.705, 0.679, and 0.667, respectively (Supplemental Figure S1C). Therefore, compared with TNM staging, the risk model we constructed is more accurate in predicting patient survival.

 Table 4 

The list of 20 DElncRNAs in the ceRNA network.

lncRNAlogFCpValueFDR
LINC008391.745.23E-871.16E-86
TRPM2-AS1.051.11E-722.14E-72
BOLA3-AS11.692.62E-1591.43E-158
LINC009202.793.17E-1843.2E-183
CASC151.295.5E-1281.85E-127
ZSCAN16-AS11.614.18E-1592.26E-158
HOTAIR-1.243.45E-1241.1E-123
AGAP11-1.246.98E-1816.38E-180
FAM182B-1.081.6E-1014E-101
XIST-1.734.98E-568.23E-56
NEAT1-4.941.57E-1943.9E-193
HAND2-AS1-1.443.14E-1853.45E-184
MEG3-5.042.2E-1903.33E-189
GABPB1-AS1-1.197.81E-1393.06E-138
HOTAIRM1-1.021.28E-872.86E-87
TTN-AS1-2.508.34E-1963.82E-194
SNHG14-1.521.57E-1711.12E-170
COX10-AS1-1.643.33E-1926.33E-191
LINC005202.135.8E-1131.63E-112
ST3GAL5-AS11.254.24E-1723.09E-171
 Figure 5 

PPI network, most significant DEmRNA module, and interaction network of hub genes. (A) PPI network of 169 DEmRNAs constructed in STRING. (B) PPI network of DEmRNAs constructed with Cytoscape. (C-D) Two modules obtained from the PPI network using MCODE Cytoscape plug-in. Red nodes represent upregulated genes, blue nodes represent downregulated genes.

J Cancer Image
 Figure 6 

ceRNA network in melanoma. Red nodes indicate increased expression, blue nodes indicate decreased expression; rectangles represent mRNAs, rhomboids represent miRNAs, and triangles represent lncRNAs.

J Cancer Image

Discussion

Epithelial-mesenchymal transformation (EMT) is the process by which epithelial cells transform into mesenchymal cells. Under normal physiological conditions, EMT is related to human embryonic development, wound healing, and tissue regeneration [19]. In the pathological state, EMT is associated with cancer progression and fibrosis, and is significantly related to the occurrence and development of cancer. During the development of cancer, because of gene mutations, adhesion molecules and connections between cells are reduced, the extracellular matrix is also reduced, and the cytoskeleton is reconstructed [20]. Therefore, cancer cells can migrate and invade, and eventually spread to other organs. Moreover, cancer cells further promote metastasis by changing the tumor microenvironment and secreting enzymes that degrade the extracellular matrix. EMT-treated cells also gain the ability to resist apoptosis. In this study, we screened 169 DEmRNAs, 91 DEmiRNAs, and 1844 DElncRNAs. The results of the GO analysis showed that the function of DEmRNAs was mainly related to epidermis development, proteoglycan binding, and structural constituents of the cytoskeleton; these factors are closely related to the occurrence of melanoma. Notably, the top 10 biological process terms were all related to skin development. With respect to molecular functions, most of the DEmRNAs were related to protein binding and structural constituents of the cytoskeleton. Therefore, we speculate that the formation of melanoma may be closely related to changes in the cytoskeleton structure of normal nevus cells.

In addition, the KEGG pathway analysis showed that the DEmRNAs were mainly enriched in bladder cancer, ECM-receptor interactions, the PPAR signaling pathway, and the IL-17 signaling pathway. Surprisingly, we found that most of the DEmRNAs were related to ECM-receptor interaction, which is closely related to cell adhesion and migration. Generally, EMT is related to the metastatic and invasive ability of cancer cells. During EMT [21], the integrity of normal epithelial cells is lost, the adhesion and connections between cells are reduced, and the cytoskeleton is reconstructed, which leads to greater movement, migration, and invasion of cancer cells. Furthermore, we found that most of the genes were enriched in the bladder cancer pathway, suggesting a significant relationship between bladder cancer and melanoma.

 Figure 7 

Expression levels and survival analysis of specific genes in melanoma compared with normal melanocytes. (A) Expression levels of specific genes in melanoma according to GEPIA (*p<0.05). (B) Survival analysis for specific genes by oncolnchttp (log-rank p<0.05).

J Cancer Image

After using the PPI network to cluster DEmRNAs, we found that among the 19 hub genes, only CD44 showed low expression in melanoma. CD44 can interact with moesin protein, which has been shown to promote the migration and infiltration of breast cancer cells by remodeling the cytoskeleton [22], thereby further promoting the metastatic ability of cancer cells [23, 24]. CD44 is a cell surface glycoprotein that is highly expressed in some cancer stem cells and is closely related to cell motility and cancer metastasis. In normal nevocytic nevus, CD44 shows strong staining. In melanoma, the expression of CD44 decreases gradually with increasing invasiveness [25]. DTL regulates the cell cycle and maintains the stability of the genome. DTL is highly expressed in invasive hepatocellular carcinoma and is closely related to tumor grade and prognosis [26]. Chen et al. found that knocking down the expression of DTL in liver cancer inhibited the growth and invasion of liver cancer cells, accelerated the apoptosis of cancer cells, and inhibited tumorigenesis [27]. The same effects were found in gastric cancer, breast cancer, and cervical cancer [28]. Through bioinformatics analysis, Zhou et al. identified DTL as a hub gene in hepatocellular carcinoma cells and found that it was highly expressed [26]. Through survival analysis, a close relationship was identified between DTL and the prognosis of patients with liver cancer. In addition, previous studies have reported that DTL is highly expressed in melanoma compared with benign melanocytic nevi and significantly correlated with OS [29]. However, the effects of DTL on the migration and invasion of melanoma have not been reported.

 Figure 8 

ROC curve and heat map of five gene combinations. (A) OS curves for five gene combinations in melanoma. (B) ROC curve and AUC. (C) Heat map of the five gene combinations.

J Cancer Image

MMP9 is a matrix metalloproteinase that promotes the migration and invasion of cancer cells by degrading the extracellular matrix [30]. MMP9 has been reported to play an important part in the EMT process of cancer cells. Its overexpression can increase the invasive ability of hepatocellular carcinoma cells and promote metastasis [31], leading to higher TNM staging and poor prognosis. High expression of MMP9 is often found in breast cancer [32]; in malignant breast cancer tissues, MMP9 shows strong positive expression. Overexpression of MMP9 can enhance the colony formation and migration ability of breast cancer cells, whereas inhibition of MMP9 leads to a significant decrease in invasive ability. Furthermore, survival analysis showed that MMP9 has a significant relationship with prognosis in breast cancer patients and could thus represent a prognostic biomarker for breast cancer [33]. In addition, MMP9 can promote angiogenesis; Li et al. showed that the angiogenic ability of melanoma cells could be inhibited by knocking down the expression of MMP9 [34]. There is a well known positive correlation between angiogenesis and invasion of cancer cells. Bakos et al. showed that the expression of MMP9 was low in normal melanocytes [35]. Therefore, we infer that MMP9 plays an important part in the transformation of normal melanocytes into melanoma.

In this study, we also constructed a lncRNA-related ceRNA network to identify related miRNAs in melanoma. The functions of lncRNAs have been implicated in a variety of cancers. LINC00520 has been widely reported in breast cancer, colorectal cancer, skin squamous cell carcinoma, melanoma, and other cancer cells, and shows increased expression in breast cancer tissue compared with normal breast tissue. Henry et al. showed that the migration ability of breast cancer cells could be blocked by knocking down the expression of LINC00520 [36]. LINC00520 can also promote the proliferation of cutaneous squamous cell carcinoma, as well as enhancing its migration and invasiveness via upregulating the expression of MMP9 [37]. Wen et al. found that LINC00520 had a similar function in melanoma. Compared with normal melanocytes, LINC00520 was highly expressed in melanoma [38], and overexpression of LINC00520 could promote the proliferation, migration, and invasion of melanoma cells. In addition, it was significantly associated with a shorter OS in patients with melanoma. The role of XIST in cancer has also been widely reported. Tian et al. showed that XIST was highly expressed in melanoma cells [39]. Overexpression of XIST can promote the proliferation, migration, and invasion of melanoma cells. Moreover, XIST is associated with TNM stage and lymphatic metastasis in patients with melanoma. We showed here for the first time that lncRNA ZSCAN16-AS1 was highly expressed in patients with melanoma and significantly correlated with poor OS. To the best of our knowledge, no role of ZSCAN16-AS1 has previously been reported in any disease. Let-7a-5p is a gene that has been studied extensively and is been found in a variety of cancer cells. In the current study, it was found to regulate downstream target genes. Functional experiments with let-7a-5p showed that its overexpression could inhibit the proliferation, migration, and invasion of lung cancer cells [40], as well as promoting autophagy and apoptosis. Li et al. found that let-7a-5p inhibited the proliferation of lung cancer cells by inhibiting the G1/S phase process [41]. Moreover, let-7a-5p was associated with shorter OS in patients with lung cancer. The miRNA sequencing analysis by Babapoor et al. showed that the expression of let-7a-5p in invasive melanoma was lower than that in normal melanocytes [42].

Different from previous studies, our study focused more on the relationship between the ceRNA network and melanoma epithelial-mesenchymal transformation (EMT). Also, the data sources of this study were more abundant. For example, mRNA and miRNA data come from GEO and ArrayExpress databases, and lncRNA data comes from TCGA and GTEx databases. Moreover, we not only constructed a melanoma ceRNA network but also carried out multivariate COX analysis of the key molecules in the ceRNA network. Finally, a risk model was constructed. To evaluate the clinical value of this risk model, through multivariate COX analysis, combined with the clinical data of melanoma patients in the TCGA database, we found that this risk model can predict the survival of melanoma patients more accurately than traditional TNM staging. To independently verify our risk model, we query the melanoma data from online databases such as the GEO database. Unfortunately, most data sets only sequence mRNA or miRNA and lack relevant lncRNA data, and some melanoma data lack patient survival data. Therefore, this study has not been able to independently verify this risk model. However, we will still strive to collect tissue samples and clinical information of melanoma patients in the future, and further, verify our risk model through gene chip sequencing data.

In the current study, we found that DTL, MMP9, LINC00520, and ZSCAN16-AS1 were highly expressed in melanoma compared with normal melanocytes, whereas XIST and let-7a-5p showed the opposite trend. The results of the survival analysis showed that these six genes were closely related to the prognosis of patients with melanoma. Furthermore, the results of our multivariate Cox analysis suggested that the combination of DTL, LINC00520, ZSCAN16-AS1, XIST, and let-7a-5p could identify high-risk melanoma patients. According to the ceRNA network, both LINC00520 and XIST can bind to let-7a-5p, and let-7a-5p can interact with its target gene DTL. However, previous studies showed that XIST was highly expressed in melanoma, contrary to our results. The role of XIST in melanoma has not been fully elucidated; we will investigate this in future research. Finally, we constructed a ceRNA network for LIINC00520/let-7a-5p/DTL.

In summary, we constructed a ceRNA network using miRNA, lncRNA, and mRNA expression profiles in melanoma and normal nevus. This will provide a basis for subsequent studies of the regulation mechanism of melanoma. In addition, we identified a combination of five genes that could represent a new signature for the diagnosis and prognostic analysis of patients with melanoma. In the future, we will further expand the collection of clinical data based on these data, and conduct a more detailed analysis of the pathogenesis of melanoma to clarify the role of these genes in the ceRNA network.

Supplementary Material

Supplementary figures and tables.

Attachment

Acknowledgements

This study was supported by the National Natural Science Foundation of China [grant nos. 81660324 and 81660521]; the Educational Research and Reform Project of Xinjiang Medical University, China [grant no. YG2019053]; the Key discipline Construction of the 13th Five-Year Plan in Xinjiang, China-Plateau Discipline Project to the article.

Competing Interests

The authors have declared that no competing interest exists.

References

1. Wei CY, Zhu MX, Lu NH. et al. Bioinformatics-based analysis reveals elevated MFSD12 as a key promoter of cell proliferation and a potential therapeutic target in melanoma. Oncogene. 2019;38:1876-91

2. Cheng PF. Medical bioinformatics in melanoma. Curr Opin Oncol. 2018;30:113-7

3. Sun C, Wang L, Huang S. et al. Reversible and adaptive resistance to BRAF(V600E) inhibition in melanoma. Nature. 2014;508:118-22

4. Damsky WE, Bosenberg M. Melanocytic nevi and melanoma: unraveling a complex relationship. Oncogene. 2017;36:5771-92

5. Lim SY, Lee JH, Diefenbach RJ. et al. Liquid biomarkers in melanoma: detection and discovery. Mol Cancer. 2018;17:8

6. D'Arcangelo D, Giampietri C, Muscio M. et al. WIPI1, BAG1, and PEX3 Autophagy-Related Genes Are Relevant Melanoma Markers. Oxid Med Cell Longev. 2018;2018:1471682

7. Roszik J, Markovits E, Dobosz P. et al. TNFSF4 (OX40L) expression and survival in locally advanced and metastatic melanoma. Cancer Immunol Immunother. 2019;68:1493-500

8. Eriksson J, Le Joncour V, Nummela P. et al. Gene expression analyses of primary melanomas reveal CTHRC1 as an important player in melanoma progression. Oncotarget. 2016;7:15065-92

9. Talantov D, Mazumder A, Yu JX. et al. Novel genes associated with malignant melanoma but not benign melanocytic lesions. Clin Cancer Res. 2005;11:7234-42

10. Goldman MJ, Craft B, Hastie M. et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol. 2020;38:675-8

11. Debroas G, Hoeffel G, Reynders A. et al. [Neuroimmune interactions in the skin: a link between pain and immunity]. Med Sci (Paris). 2018;34:432-8

12. Barrett T, Wilhite SE, Ledoux P. et al. NCBI GEO: archive for functional genomics data sets-update. Nucleic Acids Res. 2013;41:D991-5

13. Szklarczyk D, Gable AL, Lyon D. et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47:D607-D13

14. Shannon P, Markiel A, Ozier O. et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498-504

15. Sticht C, De La Torre C, Parveen A. et al. miRWalk: An online resource for prediction of microRNA binding sites. PLoS One. 2018;13:e0206239

16. Paraskevopoulou MD, Vlachos IS, Karagkouni D. et al. DIANA-LncBase v2: indexing microRNA targets on non-coding transcripts. Nucleic Acids Research. 2016;44:D231-D8

17. Tang Z, Li C, Kang B. et al. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Research. 2017;45:W98-W102

18. Anaya J. OncoLnc: linking TCGA survival data to mRNAs, miRNAs, and lncRNAs. PeerJ Computer Science. 2016 2

19. Nieto MA, Huang RY, Jackson RA. et al. Emt: 2016. Cell. 2016;166:21-45

20. Singh M, Yelle N, Venugopal C. et al. EMT: Mechanisms and therapeutic implications. Pharmacol Ther. 2018;182:80-94

21. Lamouille S, Xu J, Derynck R. Molecular mechanisms of epithelial-mesenchymal transition. Nat Rev Mol Cell Biol. 2014;15:178-96

22. Wang CC, Liau JY, Lu YS. et al. Differential expression of moesin in breast cancers and its implication in epithelial-mesenchymal transition. Histopathology. 2012;61:78-87

23. Ichikawa T, Masumoto J, Kaneko M. et al. Moesin and CD44 expression in cutaneous melanocytic tumours. Br J Dermatol. 1998;138:763-8

24. Martin TA, Harrison G, Mansel RE. et al. The role of the CD44/ezrin complex in cancer metastasis. Critical Reviews in Oncology/Hematology. 2003;46:165-86

25. Harwood CA, Green MA, Cook MG. CD44 expression in melanocytic lesions: a marker of malignant progression? Br J Dermatol. 1996;135:876-82

26. Zhou Z, Li Y, Hao H. et al. Screening Hub Genes as Prognostic Biomarkers of Hepatocellular Carcinoma by Bioinformatics Analysis. Cell Transplant. 2019;28:76S-86S

27. Chen YC, Chen IS, Huang GJ. et al. Targeting DTL induces cell cycle arrest and senescence and suppresses cell growth and colony formation through TPX2 inhibition in human hepatocellular carcinoma cells. Onco Targets Ther. 2018;11:1601-16

28. van Dam PA, Rolfo C, Ruiz R. et al. Potential new biomarkers for squamous carcinoma of the uterine cervix. ESMO Open. 2018;3:e000352

29. Yang L, Dai J, Ma M. et al. Identification of a functional polymorphism within the 3'-untranslated region of denticleless E3 ubiquitin protein ligase homolog associated with survival in acral melanoma. Eur J Cancer. 2019;118:70-81

30. Yang F, Yu N, Wang H. et al. Downregulated expression of hepatoma-derived growth factor inhibits migration and invasion of prostate cancer cells by suppressing epithelial-mesenchymal transition and MMP2, MMP9. PLoS One. 2018;13:e0190725

31. Church MK, Kolkhir P, Metz M. et al. The role and relevance of mast cells in urticaria. Immunol Rev. 2018;282:232-47

32. Dong H, Diao H, Zhao Y. et al. Overexpression of matrix metalloproteinase-9 in breast cancer cell lines remarkably increases the cell malignancy largely via activation of transforming growth factor beta/SMAD signalling. Cell Prolif. 2019;52:e12633

33. Ren F, Tang R, Zhang X. et al. Overexpression of MMP Family Members Functions as Prognostic Biomarker for Breast Cancer Patients: A Systematic Review and Meta-Analysis. PLoS One. 2015;10:e0135544

34. Li L, Fan P, Chou H. et al. Herbacetin suppressed MMP9 mediated angiogenesis of malignant melanoma through blocking EGFR-ERK/AKT signaling pathway. Biochimie. 2019;162:198-207

35. Bakos RM, Bakos L, Edelweiss MI. et al. Immunohistochemical expression of matrix metalloproteinase-2 and -9 in melanocytic nevi is altered by ultraviolet B. Photodermatol Photoimmunol Photomed. 2007;23:250-4

36. Henry WS, Hendrickson DG, Beca F. et al. LINC00520 is induced by Src, STAT3, and PI3K and plays a functional role in breast cancer. Oncotarget. 2016;7:81981-94

37. Mei XL, Zhong S. Long noncoding RNA LINC00520 prevents the progression of cutaneous squamous cell carcinoma through the inactivation of the PI3K/Akt signaling pathway by downregulating EGFR. Chin Med J (Engl). 2019;132:454-65

38. Luan W, Ding Y, Yuan H. et al. Long non-coding RNA LINC00520 promotes the proliferation and metastasis of malignant melanoma by inducing the miR-125b-5p/EIF5A2 axis. J Exp Clin Cancer Res. 2020;39:96

39. Tian K, Sun D, Chen M. et al. Long Noncoding RNA X-Inactive Specific Transcript Facilitates Cellular Functions in Melanoma via miR-139-5p/ROCK1 Pathway. Onco Targets Ther. 2020;13:1277-87

40. Duan S, Yu S, Yuan T. et al. Exogenous Let-7a-5p Induces A549 Lung Cancer Cell Death Through BCL2L1-Mediated PI3Kgamma Signaling Pathway. Front Oncol. 2019;9:808

41. Li JP, Liao XH, Xiang Y. et al. Hyperoside and let-7a-5p synergistically inhibits lung cancer cell proliferation via inducing G1/S phase arrest. Gene. 2018;679:232-40

42. Babapoor S, Wu R, Kozubek J. et al. Identification of microRNAs associated with invasive and aggressive phenotype in cutaneous melanoma by next-generation sequencing. Lab Invest. 2017;97:636-48

Author contact

Corresponding address Corresponding author: Prof. ShuMei Feng, Department of Histology and Embryology, Basic Medical College, Xinjiang Medical University. Email: 87391167com.


Received 2020-8-11
Accepted 2021-2-15
Published 2021-3-15


Citation styles

APA
Ding, Y., Li, M., Tayier, T., Zhang, M., Chen, L., Feng, S. (2021). Bioinformatics analysis of lncRNA‑associated ceRNA network in melanoma. Journal of Cancer, 12(10), 2921-2932. https://doi.org/10.7150/jca.51851.

ACS
Ding, Y.; Li, M.; Tayier, T.; Zhang, M.; Chen, L.; Feng, S. Bioinformatics analysis of lncRNA‑associated ceRNA network in melanoma. J. Cancer 2021, 12 (10), 2921-2932. DOI: 10.7150/jca.51851.

NLM
Ding Y, Li M, Tayier T, Zhang M, Chen L, Feng S. Bioinformatics analysis of lncRNA‑associated ceRNA network in melanoma. J Cancer 2021; 12(10):2921-2932. doi:10.7150/jca.51851. https://www.jcancer.org/v12p2921.htm

CSE
Ding Y, Li M, Tayier T, Zhang M, Chen L, Feng S. 2021. Bioinformatics analysis of lncRNA‑associated ceRNA network in melanoma. J Cancer. 12(10):2921-2932.

This is an open access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/). See http://ivyspring.com/terms for full terms and conditions.
Popup Image