Abstract

   Cutaneous melanoma is the most aggressive type of skin cancer which its
   incidence has significantly increased in recent years worldwide. Thus,
   more investigations are required to identify the underlying mechanisms
   of melanoma malignant transformation and metastasis. In this context,
   long noncoding RNAs (lncRNAs) are a new type of noncoding transcripts
   that their dysregulations are associated with almost all cancers
   including melanoma. However, the precise functional roles of most of
   the significantly altered lncRNAs in melanoma have not yet been fully
   inspected. In this study, a comprehensive list of lncRNAs was
   interrogated across cutaneous melanoma samples to identify the
   significantly altered/dysregulated lncRNAs. To this end, lncRNAs were
   filtered in several steps and the selected lncRNAs projected to a
   bioinformatic and systems biology analysis using several publicly
   available databases and tools such as GEPIA and cBioPortal. According
   to our results, 30 lncRNAs were notably altered/dysregulated in
   cutaneous melanoma most of which were co-expressed with each other.
   Also, co-expression/alteration and differential expression analyses led
   to the selection of 12 out of these 30 lncRNAs as cutaneous melanoma
   key lncRNAs. Furthermore, functional demonstrated that these 12 lncRNAs
   might be involved in melanoma-relevant biological processes and
   pathways. In addition, the end result of our analyses demonstrated that
   these lncRNAs are associated with the clinicopathological features of
   melanoma patients. These 12 lncRNAs need to be further investigated in
   future studies to characterize their exact roles in melanoma
   development and to identify their potential for being used as drug
   targets and/or biomarkers for cutaneous melanoma.

   Keywords: Cutaneous melanoma, long noncoding RNAs (lncRNAs),
   bioinformatics, systems biology

Introduction

   Melanoma is the most lethal type of skin cancer whose prevalence has
   significantly increased during the last decades.^[27]1 Among all
   invasive melanoma cases, almost 5% are associated with inheritance and
   familial setting, which indicates that low-penetrance/high-prevalence
   genes do not have much of a role in the incidence of melanoma rate in
   the world. However, environmental factors, especially sun exposure,
   play a significant role in melanoma initiation and progression.^[28]2
   Although sun exposure is the major risk factor for cutaneous melanoma,
   specific genetic alterations are also associated with melanoma
   risk.^[29]3 Cutaneous melanoma is hard to treat, and the available
   drugs and recent advancements in different treatment strategies
   including targeted therapy have not significantly increased the
   survival rate and period of melanoma patients.^[30]4,[31]5 However,
   recent advancements in melanoma immunotherapy have been considerably
   promising.^[32]6,[33]7

   Long noncoding RNAs (lncRNAs) are a new type of nonprotein coding
   transcripts that are present in almost all organisms and have unique
   features that differentiate them from messenger RNAs (mRNAs).^[34]8
   These long transcripts are usually more than 200 nucleotides, do not
   code for protein, and are involved in a variety of regulatory
   mechanisms including gene expression.^[35]9 Actually, lncRNAs play
   their roles, in relation to other molecules, as scaffolds, signals,
   guides, and decoys.^[36]10 Cumulative evidence are reporting the
   association of lncRNAs with developmental processes and different
   diseases, but the precise mechanism of action of most of them is not
   still elucidated.^[37]11 Also, advancements in high-throughput
   sequencing technologies have led to the discovery of an increasing
   number of novel transcripts including lncRNAs, most of which are not
   still fully annotated.^[38]12

   LncRNAs play key roles in a variety of diseases including cancer. In
   fact, they involve in the epigenetic regulation of other macromolecules
   through chromatin remodeling, as well as transcriptional and
   post-transcriptional regulation, and, consequently, their
   dysregulations lead to different diseases such as cancer.^[39]13 In
   other words, lncRNAs interact with other biological macromolecules
   including DNA, protein, and RNA and mediate several cancer
   manifestations.^[40]14 Frequent studies are reporting that a
   significant proportion of the genomic mutations in cancer occur within
   nonprotein coding regions, usually inside the lncRNA loci. Also, the
   association between aberrant expression of thousands of lncRNAs with
   different cancers has been demonstrated. However, only a small
   proportion of these lncRNAs have been functionally annotated. Notably,
   it has been indicated that these annotated lncRNAs involve in several
   cellular processes including proliferation, migration, or genomic
   stability and consequently their dysregulations drive malignant
   transformation.^[41]15,[42]16

   Since lncRNAs have a higher tissue specificity than mRNAs, they might
   adopt a dual functionality in cancer, executing either as an oncogene
   or tumor suppressor, according to their tissue of origin.^[43]17 As
   stated earlier, the association between thousands of lncRNAs and cancer
   has been reported. Among all, HOTAIR is the most known lncRNA in the
   context of cancer which plays an oncogenic role in many of the human
   malignancies including melanoma.^[44]9 Several other lncRNAs such as
   BANCR, MALAT1, UCA1, ANRIL, and SNHG5 are also connected with cutaneous
   melanoma and were reviewed by Hulstaert et al,^[45]5 some of which
   could be considered as biomarkers or therapeutic targets for malignant
   melanoma.^[46]5 Also, several studies have reported the association of
   individual lncRNAs such as LINC00673,^[47]18 GAS5,^[48]19 UCA1,^[49]20
   MEG3,^[50]21 and PVT1^[51]22 with melanoma. Even though, there are
   quite a number of lncRNAs that are significantly altered in cutaneous
   melanoma at the genomic and/or transcriptomic levels but are not still
   functionally annotated. In a similar study by Wang et al^[52]1 on the
   association of lncRNAs with the tumorigenesis and metastasis of
   cutaneous melanoma, 295 lncRNAs were identified as DE-lncRNAs between
   normal skin and metastatic melanoma tumor of which 4 lncRNAs were
   common with our results. However, we could not find common results
   between our study and another similar study on the association of
   lncRNAs and melanoma done by Zhang et al.^[53]23

   Bioinformatics and systems biology, in its simplest description, is the
   development and application of different computational methods and
   tools to convert biological “omic” data into the information that their
   interpretation could enhance the personalized treatment of
   patients.^[54]24,[55]25 In this study, we adopted an integrative
   bioinformatic and systems biology approach to examine if any of the
   known lncRNAs have significant associations with cutaneous melanoma. To
   this purpose, the processed data of TCGA Skin Cutaneous Melanoma study
   were interrogated. We also investigated the biological processes and
   pathways these lncRNAs might participate in the context of cutaneous
   melanoma. Overall, several novel candidate biomarkers and drug targets
   for cutaneous melanoma were presented in this study.

Materials and Methods

Data preparation

   First, 3823 lncRNAs were retrieved from HGNC BioMart
   ([56]https://biomart.genenames.org/)^[57]26 (Supplementary file 1)
   according to the following criteria: HGNC Mart Genes Datasets; Filter
   by genes with RefSeqs accession; Status: Approved; Locus group:
   noncoding RNA; Locus type: RNA, long noncoding. Then, using cBioPortal
   database,^[58]27,[59]28 all lncRNAs were queried at both genomic and
   transcriptomic levels against the Skin Cutaneous Melanoma (TCGA,
   Provisional) study, comprising 478 samples, and lncRNAs with
   alterations (copy-number alteration [CNA] and expression dysregulation)
   in ⩾10% of samples were selected for downstream analyses (30 genes;
   Supplementary file 2). cBioPortal is a database that facilitate in
   visualizing, analyzing, and downloading of large-scale cancer genomics
   data sets. Finally, the Expression (RNA Seq V2 RSEM), Expression
   z-Scores (RNA Seq V2 RSEM), and relative linear copy-number values of
   all selected lncRNAs were imported from the cBioPortal database into R
   statistical software by means of CDGS-R package v1.2.8
   ([60]https://CRAN.R-project.org/package=cgdsr). These data have been
   preprocessed and normalized by the cBioPortal database. The CGDS-R is a
   package for R statistical software that provides a basic set of
   functions for querying the Cancer Genomic Data Server (CGDS).

Genotranscriptomic analysis

   The mean expression level of each lncRNA across all samples was
   calculated to compare the expression of these lncRNAs in cutaneous
   melanoma. Also, all 30 lncRNAs were interrogated across single-cell
   RNA-Seq data sets using scRNASeqDB
   ([61]https://bioinfo.uth.edu/scrnaseqdb/index.php) and the expression
   of lncRNAs with significant rank in melanoma data sets were retrieved.
   It is also worth mentioning that the expression of our
   under-investigation lncRNAs were tested in scRNASeqDB regardless of the
   library preparation method used for the generation of different data
   sets. These analyses help identify the lncRNAs with noticeable
   expression in cutaneous melanoma at both bulk and single-cell levels.
   Furthermore, the correlation between CNA and z-Score expression data of
   each lncRNA was analyzed using the Spearman method to identify the
   lncRNAs which their dysregulation could have been due to their genomic
   alterations.

Analysis of the co-alteration/expression of lncRNAs with each other

   The significant mutual exclusivity and co-occurrence of alterations
   among 30 lncRNAs were obtained from cBioPortal to interrogate the
   interactions of these lncRNAs (Supplementary file 3). Then, gene pairs
   with significant co-occurrent alterations were recruited to reconstruct
   an association network using Cytoscape v3.5.1
   ([62]http://www.cytoscape.org).^[63]29 Also, to identify the hub (most
   important) lncRNAs within the association network, the topology of the
   network was analyzed using the Analyzer tool of the Cytoscape software.
   Similarly, both expression and z-Score expression data were used to
   analyze the coexpression and z-Score coexpression (i.e. the correlation
   between z-Score expression data) of 30 lncRNAs, respectively. To this
   purpose, the correlation of both expression and z-Score expression data
   was calculated independently using the Pearson method, visualized as an
   integrated network and analyzed using Cytoscape v3.5.1.

Differential expression analysis

   GEPIA web server ([64]http://gepia.cancer-pku.cn)^[65]30 was used to
   investigate the differential expression of 30 lncRNAs in cutaneous
   melanoma according to the TCGA RNA-Seq data and using the LIMMA method.
   GEPIA is a web server that interactively analyzes the RNA-Seq
   expression data of tumor and normal samples from the TCGA and the GTEx
   projects. Furthermore, a comprehensive literature review was conducted
   to find if the same dysregulations of these lncRNAs have been
   previously reported in melanoma or other cancers.

Selection of key lncRNAs

   To focus our investigations on the most important lncRNAs, we combined
   the list of cutaneous melanoma differentially expressed lncRNAs
   (DE-lncRNAs) with the hub lncRNAs derived from the topological analysis
   of co-alteration/expression networks and named them as the “cutaneous
   melanoma key lncRNAs” (CMK-lncRNAs). Accordingly, the CMK-lncRNAs
   include 2 groups of lncRNAs: (1) those lncRNAs that are both
   significantly differentially expressed in cutaneous melanoma samples
   compared with normal ones and have genomic alterations in over 10% of
   cutaneous melanoma samples; (2) those lncRNAs that are both hubs of the
   cutaneous melanoma co-alteration/expression networks and have
   mutations/CNAs (genomic alterations) in over 10% of cutaneous melanoma
   samples. Also, all downstream analyses were implemented based on the
   CMK-lncRNAs.

Coexpression analysis of CMK-lncRNAs with other genes

   The coexpression of CMK-lncRNAs with each other and all other genes
   were interrogated to predict the functional similarity of these lncRNAs
   according to their associated genes (guilt by association principle).
   To obtain the coexpressed genes with the highest level of confidence,
   the preprocessed data of both cBioPortal and GEPIA databases were
   integrated. To this end, the coexpression of each CMK-lncRNA with other
   genes in the skin cutaneous melanoma (TCGA, Provisional) study was
   analyzed independently using cBioPortal, and genes with Pearson and
   Spearman correlation coefficient over 0.3 were outputted (Supplementary
   file 4). Besides, the coexpression of CMK-lncRNAs with other genes in
   the cutaneous melanoma TCGA study was evaluated via GEPIA web server,
   and gene pairs with Pearson correlation coefficient (PCC) over 0.3 were
   selected (Supplementary file 5). Then, the intersection of coexpression
   data derived from cBioPortal and GEPIA databases was identified.
   Moreover, to increase the sensitivity of our downstream analyses, the
   expression status (either dysregulation or normal expression) of
   coexpressed genes was imposed to that of CMK-lncRNAs and the
   inconsistencies were filtered out. To this purpose, the list of all
   differentially expressed genes in cutaneous melanoma TCGA study was
   retrieved from the GEPIA web server (Supplementary file 6) and
   subjected to the expression data; if the CMK-lncRNA was dysregulated,
   the coexpressed genes that were dysregulated were retained, and if the
   CMK-lncRNA was normally expressed, the coexpressed genes with normal
   expression were selected.

Functional prediction

   The coexpression data obtained in the previous step were used to
   perform a computational functional analysis. First, the coexpressed
   genes of CMK-lncRNAs were employed to perform a gene set enrichment
   analysis of gene ontology-biological process (GO-BP) terms using
   Enrichr web server
   ([66]http://amp.pharm.mssm.edu/Enrichr).^[67]31,[68]32 Next, all of the
   CMK-lncRNAs were interrogated via FuncPred database
   ([69]http://www.funcpred.com)^[70]33 to investigate the association of
   them with GO-BP terms in normal skin tissue based on the evolutionary
   conserved and tissue-specific expression data. Finally, first-ranked
   GO-BP terms of FuncPred (according to the lowest FDR) and Enrichr
   (according to the highest combined score) results as well as their
   intersection were selected as the most significant GO-BP terms.
   Finally, the GO-BP terms and their associated lncRNAs were used to
   reconstruct an association network and identify the most important
   CMK-lncRNAs within that using the Cytoscape. Besides, considering the
   modulated lncRNAs in the lncRNA-GO-BP network, the ConsensusPathDB
   database ([71]http://cpdb.molgen.mpg.de)^[72]34,[73]35 was used to
   figure out if coexpressed genes of CMK-lncRNAs in each module are
   overrepresented in any of the Reactome pathways. Also, a literature
   review was performed to investigate the association between identified
   biological processes and first-ranked pathways with cutaneous melanoma
   development.

Clinicopathological analysis

   At last, the association between the CMK-lncRNAs and the
   clinicopathological features of cutaneous melanoma patients was
   assessed. To this end, the GEPIA web server was used to analyze the
   prognostic value of all CMK-lncRNAs using the Kaplan–Meier estimates.
   Likewise, the differential expression of CMK-lncRNAs between different
   stages of cutaneous melanoma was investigated using the GEPIA web
   server. In both of these analyses, R statistical software was also used
   to perform a multiple testing correction and calculate the adjusted P
   values (p.adj) using the “Hochberg” method. Moreover, the CMK-lncRNAs
   that were differentially expressed ([74]Table 1) were focused to
   further investigate the clinicopathological significance of their
   dysregulations in cutaneous melanoma. To this purpose, the samples
   containing significant alterations in the expression (|expression
   z-Score| > 2) of any of the 9 focused CMK-lncRNAs were extracted from
   the skin cutaneous melanoma (TCGA, Provisional) data set. Then, the
   clinicopathological data of extracted samples were imported from the
   cBioPortal database into R statistical software by means of CDGS-R
   package v1.2.8. Finally, the overall survival (OS), disease-free
   survival (DFS), and pathological tumor stage of all extracted samples
   were studied against the selected genes.

Table 1.

   DE-lncRNAs in cutaneous melanoma.
   Gene symbol Gene symbol Dysregulation Log2FC Adj P^[75]* Previous
   confirmation
   LINC00518 long intergenic nonprotein coding RNA 518 Over-expressed 4.03
   2.22E−257 Melanoma^[76]36
   LINC00174 long intergenic nonprotein coding RNA 174 Under-expressed
   −1.404 9.00E−185 Glioma^[77]37
   LINC00467 long intergenic nonprotein coding RNA 467 Over-expressed 1.77
   1.78E−224 Lung adenocarcinoma^[78]38
   CASC15 cancer susceptibility 15 Over-expressed 1.674 5.88E−151
   Melanoma^[79]39
   LINC00964 long intergenic nonprotein coding RNA 964 Under-expressed
   −1.521 5.27E−215 Esophageal cancer^[80]21 and HNSCC[81]^**^[82]40
   LINC00265 long intergenic nonprotein coding RNA 265 Under-expressed
   −1.08 1.65E−121 Colorectal Cancer^[83]41
   MIR205HG MIR205 host gene Under-expressed −7.482 0.00E+00 Cervical
   cancer^[84]42
   HCP5 HLA complex P5 Under-expressed −1.219 4.76E−47 Melanoma^[85]43
   LINC00662 long intergenic nonprotein coding RNA 662 Over-expressed
   1.018 6.55E−134 Lung cancer^[86]44
   [87]Open in a new tab
   ^*

   Adjusted P value.
   ^**

   Head and neck squamous cell carcinoma.

Statistical and topological analysis

   All the statistical analyses were done by R statistical software, R
   Development Core Team (2014), freely available at
   [88]http://www.r-project.org. The absolute z-Score > 2 was considered
   as the significant expression dysregulation throughout the study. Also,
   the absolute Pearson and Spearman correlation coefficients >0.3 were
   considered as statistically significant. Moreover, the P-value and P
   adj < .05 was used as the significance threshold in all of the
   analyses. Regarding the topological analysis of networks, 2 different
   network metrics including degree (the number of edges incident to each
   node) and betweenness centrality (a measure of node centrality based on
   shortest paths) were simultaneously employed to identify the hub nodes.

Results

Remarkable expression of lncRNAs in cutaneous melanoma

   The genotranscriptomic analysis of all lncRNAs using cBioPortal
   demonstrated that 30 lncRNAs had significant alterations (CNA and
   expression dysregulation) in ⩾10% of samples ([89]Figure 1). Also, it
   is notable that 29 out of 30 lncRNAs had up-regulative (i.e.
   amplification and over-expression) alterations, while CDKN2A-AS1 lncRNA
   was mostly deleted in cutaneous melanoma. The expression data of 30
   lncRNAs, retrieved from cBioPortal, demonstrated that all of them
   except LINC00051 and HPYR1 have notable expression in cutaneous
   melanoma ([90]Figure 2A). Similarly, the single-cell RNA-seq analysis
   of lncRNAs demonstrated that 21 out of 30 lncRNAs were significantly
   expressed in a single-cell RNA-seq melanoma data set (GEO series
   [91]GSE72056). Although, the analyzed data of only 19 lncRNAs were
   accessible ([92]Figure 2B). Interestingly, HPYR1 lncRNA had
   considerable expression at the single-cell level. Furthermore, analysis
   of the correlation between CNA and z-Score expression data of lncRNAs
   demonstrated that the z-Score expression of 19 out of 30 lncRNAs is
   significantly correlated (rho > 0.3, P-value < 1.91E–09) with their
   copy-number values ([93]Figure 2C).

Figure 1.

   [94]Figure 1.
   [95]Open in a new tab

   30 lncRNAs with genotranscriptomic alterations in ⩾10% of samples of
   TCGA skin cutaneous melanoma.

   The percentages on the Y axis are indicative of the proportion of
   samples with genomic and transcriptomic alterations in the
   corresponding lncRNAs in cutaneous melanoma. Blue and red colors
   demonstrate the type of genomic or transcriptomic alterations, and the
   light gray color is indicative of samples with no alteration in the
   corresponding lncRNAs in TCGA skin cutaneous melanoma. This figure was
   achieved by means of cBioPortal database.

Figure 2.

   [96]Figure 2.
   [97]Open in a new tab

   The expression of selected lncRNAs in cutaneous melanoma: (A) Mean
   expression levels of selected lncRNAs in cutaneous melanoma according
   to the data retrieved from cBioPortal. Each color represents a
   different lncRNA and the mean expression levels are log2 transformed.
   (B) maximum expression levels of selected lncRNAs in melanoma cell
   groups according to the data obtained from scRNASeqDB database. Each
   color represents a different lncRNA and the maximum expression values
   of lncRNAs in melanoma cell groups are log^[98]2 transformed. (C) The
   correlation between z-score expression and copy-number alteration data
   of selected lncRNAs in cutaneous melanoma.

   All of these figures were drawn by R statistical software.

Co-alteration/expression of 30 lncRNAs

   The network of gene pairs with co-occurrent alterations illustrated
   that the selected lncRNAs have notable co-occurrent alterations with
   each other ([99]Figure 3A). According to the topological analysis of
   this network, the lncRNAs LINC00240, SNHG11, and HCG18 were identified
   as hub nodes of the network. Also, the coexpression analysis of lncRNAs
   indicated that 20 out of 30 lncRNAs are significantly coexpressed with
   each other (PCC > 0.3) according to both coexpression and z-Score
   coexpression data ([100]Figure 3B). Likewise, the topological analysis
   of coexpression network revealed the lncRNAs LINC00265, SNHG11, and
   HCG18 as hub nodes of the network.

Figure 3.

   [101]Figure 3.
   [102]Open in a new tab

   Co-expression/alteration of selected lncRNAs with each other: (A)
   co-occurrent alteration network of lncRNAs with each other in cutaneous
   melanoma. In this network, both genomic and transcriptomic alterations
   are considered. The intensity of purple node color is proportional to
   the betweenness centrality score; stronger purple color indicates
   higher betweenness centrality score. The node size is proportional to
   node degree; larger node indicates higher node degree. (B)
   co-expression network of lncRNAs with each other in cutaneous melanoma.

   Red and blue edges are indicative of coexpression and z-Score
   coexpression of lncRNAs with each other, respectively. Also, the edge
   width demonstrates the Pearson correlation coefficient; wider edge is
   indicative of higher Pearson correlation coefficient. The intensity of
   green node color is proportional to the betweenness centrality score;
   stronger green color indicates higher betweenness centrality score. The
   node size is proportional to node degree; larger node indicates higher
   node degree.

Differential expression of lncRNAs in cutaneous melanoma compared with normal
samples

   The differential expression analysis of 30 lncRNAs according to the
   TCGA and GTEx RNA-Seq data demonstrated that 9 out of 30 lncRNAs were
   differentially expressed (|Log2FC|> 1, q value < 0.05) in cutaneous
   melanoma 5 of which were under-expressed while others were
   over-expressed ([103]Table 1). Interestingly, the literature review
   demonstrated same dysregulations of 5 out of these 9 lncRNAs in
   melanoma or other cancers. Among all, LINC00518 lncRNA had the most
   remarkable differential expression (Log2FC = 4.03) in TCGA cutaneous
   melanoma samples compared with normal ones.

Functional roles of CMK-lncRNAs

   Twelve out of 30 lncRNAs were finally selected as the CMK-lncRNAs. The
   coexpression analysis of CMK-lncRNAs with other genes using both
   cBioPortal and GEPIA web server demonstrated that they are
   significantly coexpressed (PCC > 0.3) with other genes (Supplementary
   file 7). The coexpressed genes with any of the CMK-lncRNAs were used to
   perform a gene set enrichment analysis, which demonstrated that most of
   the CMK-lncRNAs might be associated with cancer-related biological
   processes ([104]Figure 4; Supplementary file 8). For instance, the
   potential of melatonin, as a potent antineoplastic factor, in the
   suppression of cell growth in some cutaneous and uveal melanoma cell
   lines have been previously demonstrated, which might be associated with
   the involvement of melatonin in the regulation of keratinocyte
   proliferation (GO:0010837; associated with CASC15).^[105]45 Also, it
   has been demonstrated that targeting of NGLY1, a protein
   deglycosylation (GO:0006517; associated with LINC00265) enzyme, could
   promote the efficiency of melanoma targeted therapy and chemotherapy by
   inducing stress signaling-associated apoptosis, pleiotropic responses,
   and cytokine surges.^[106]46 Similarly, the genes involved in the
   Regulation of Macrophage Chemotaxis (GO:0010758; associated with
   LINC00174) are notably downregulated in FoxO1-deficient macrophages
   which results in the transition of macrophages to an M2-like/TAM-like
   phenotype and development of melanoma.^[107]47 Moreover, it has been
   reported that genes involved in antigen receptor-mediated signaling
   pathway (GO:0050851; associated with LINC00265) could be used for
   predicting prognosis and immunotherapy response in melanoma.^[108]48
   Furthermore, Zhang et al^[109]49 showed that downregulated genes in
   uveal melanoma were involved in several pathways including regulation
   of Ras protein signal transduction (GO:0046578; associated with
   LINC00265 and LINC00174). Interestingly, several of these
   melanoma-associated biological processes were connected to lncRNAs of
   the module B of the lncRNA-GO-BP network ([110]Figure 4). In addition,
   according to the pathway enrichment analysis via ConsensusPathDB
   database, coexpressed genes of CMK-lncRNAs in Module B of the
   lncRNA-GO-BP network were significantly enriched (P-value and q value <
   .05) in different Reactome pathways (Supplementary file 9) of which the
   “Cytokine Signaling in Immune system” pathway was identified as the
   first-ranked result (P-value = 2.37E−11, q value = 3.41E−09). Notably,
   27 out of 217, 12.44%, of coexpressed genes of lncRNAs in Module B of
   the lncRNA-GO-BP network were enriched in this pathway. It was
   demonstrated that the “Cytokine signaling in Immune system” was among
   the 4 pathways most impacted by Rosiglitazone (RGZ) drug, the
   application of which promotes the growth and development of melanoma
   via activating the paracrine signaling.^[111]50 Also, Chen et
   al^[112]51 demonstrated that immune check-point blockade therapies
   promote melanoma treatment via influencing several immune-associated
   pathways including the “Cytokine signaling in Immune system.”

Figure 4.

   [113]Figure 4.
   [114]Open in a new tab

   The lncRNA-GO-BP association network.

   This network demonstrates the lncRNAs with different biological
   processes. Pink nodes represent CMK-lncRNAs. The intensity of blue node
   color indicates the betweenness centrality score; stronger blue color
   indicates higher betweenness centrality score. The edge color is
   representative of the source database of predicted GO term and the
   reference tissue type. The network was reconstructed using the
   Cytoscape software.

Clinicopathological significance of CMK-lncRNAs

   The survival analysis of CMK-lncRNAs across all TCGA cutaneous melanoma
   samples indicated that the expression level of 4 out of 12 CMK-lncRNAs
   was significantly associated (Logrank test P adj < .05) with the OS of
   cutaneous melanoma patients ([115]Figure 5A). More precisely, while the
   expression level of LINC00518 and MIR205HG was negatively associated
   with the OS of cutaneous melanoma patients that of HCP5 and LINC00964
   had positive correlation with the OS of patients. Similarly, the
   cutaneous melanoma pathological stage analysis of CMK-lncRNAs showed
   that the expression of 4 out of 12 CMK-lncRNAs is significantly altered
   (P adj < .05) in different stages of the cutaneous melanoma
   ([116]Figure 5B). Moreover, the clinicopathological investigation of
   the cutaneous melanoma samples with dysregulation in any of the
   differentially expressed CMK-lncRNAs demonstrated that the OS and DFS
   of more than 56% of the selected samples were lower than 5 years, and
   more than 55% of the melanoma samples were in stage III or IV
   (Supplementary file 10).

Figure 5.

   [117]Figure 5.
   [118]Open in a new tab

   The prognostic value of CMK-lncRNAs for melanoma: (A) the expression
   level of 4 out of 12 CMK-lncRNAs had significant associations with
   survival rate of melanoma patients based on Logrank test analyses. (B)
   The significant differential expression of 4 out of 12 CMK-lncRNAs
   between different stages of cutaneous melanoma.

   CMK-lncRNAs stands for cutaneous melanoma key lncRNAs. These plots were
   obtained using the GEPIA web server.

Discussion

   LncRNAs are a new type of noncoding transcripts that play pivotal
   regulatory roles in the cell. An increasing number of studies are
   reporting the association between lncRNA dysregulations and a variety
   of cancers. Cutaneous melanoma has not been an exception, and the
   involvement of a number of lncRNAs in the development of this
   malignancy has been previously demonstrated.^[119]5 However, since
   melanoma is the most lethal type of skin cancer in the world,
   investigation of the lncRNAs involved in the initiation and progression
   of this cancer should be highly regarded. In this study, we passed 3823
   lncRNAs through several filtration steps and conducted a bioinformatic
   and systems biology analysis to identify key lncRNAs associated with
   cutaneous melanoma and to characterize their functional roles. A
   similar study on the association of lncRNAs with survival of cutaneous
   melanoma patients^[120]52 reported 6 significant lncRNAs but only one
   of them, namely HCP5, was common with our results. In addition, while
   GEPIA and cBioPortal are well-known specialized cancer databases,
   significantly dysregulated lncRNAs obtained from these 2 databases did
   not include the previously reported melanoma-associated lncRNAs such as
   SPRY4-IT1 (SPRY4 intronic transcript 1).^[121]53,[122]54 This could
   possibly be referred to the differences between TCGA Research Network
   and individual/nonglobal research projects in the context of the cohort
   of melanoma patients/samples used and the adopted analysis protocol.

   According to our results, 30 lncRNAs were significantly altered in
   cutaneous melanoma and, thus, could be essential factors in cutaneous
   melanoma initiation/progression. However, the differential expression
   analysis of these lncRNAs using the GEPIA web server demonstrated that
   9 out of 30 lncRNAs were differentially expressed in cutaneous melanoma
   compared with normal samples. It is also worth mentioning that the
   dysregulation directionality (either up- or down-regulation) of only 4
   out of 9 differentially expressed lncRNAs, including LINC00518,
   LINC00467, CASC15, and LINC00662, was consistent with the data
   retrieved from the cBioPortal. This is possibly due to the integration
   of TCGA SKCM normal (#1) and GTEx normal (#557) skin samples for
   differential analysis by the GEPIA web server, while cBioPortal
   database relies only on the TCGA data. Also, the positive correlation
   between expression dysregulations and CNA of 19 out of 30 lncRNAs
   demonstrated that amplifications/deletions of these lncRNAs play a key
   role in the alteration of their expression level. Moreover, the
   expression analysis also demonstrated that most of these lncRNAs had
   notable expression level according to both single-cell and bulk RNA-seq
   data, which is confirming that these lncRNAs might play essential roles
   in cutaneous melanoma cells and possibly normal melanocytes.

   The network analysis of co-occurrent alteration and coexpression of
   lncRNAs demonstrated that most of the selected lncRNAs are
   significantly coexpressed/co-altered in cutaneous melanoma. This
   indicates that these lncRNAs are possibly involved in common biological
   processes and pathways. In addition, topological analysis of these
   networks identified lncRNAs SNHG11, HCG18, LINC00265, and LINC00240 as
   the hub lncRNAs of which LINC00265 and LINC00240 have confirmed
   associations with lung adenocarcinoma^[123]55 and esophageal squamous
   cell carcinoma,^[124]56 respectively. However, there is still no
   evidence of the association of SNHG11 and HCG18 with neither melanoma
   nor other malignancies.

   According to our applied methodology (see Methods), 12 lncRNAs,
   including, LINC00518, LINC00174, LINC00467, CASC15, LINC00964,
   LINC00265, MIR205HG, HCP5, LINC00662, HCG18, SNHG11, and LINC00240,
   were presented as CMK-lncRNAs. Although 3 out of 12 CMK-lncRNAs,
   including HCG18, SNHG11, and LINC00240, were not differentially
   expressed according to the GEPIA web server, these lncRNAs might still
   be pivotal factors in cutaneous melanoma progression. Actually, these
   lncRNAs might be considered as “Class II cancer genes.”^[125]57 Class
   II cancer genes are those cancer-related genes that, though, are
   neither frequently mutated nor differentially expressed in cancer, they
   could play a role as “signal linkers” in coordinating cancerous signals
   between differentially expressed and mutated genes.^[126]57

   The gene set enrichment analysis of coexpressed genes of CMK-lncRNAs
   indicated that most of the CMK-lncRNAs might be involved in
   cancer-associated biological processes. Among all, 3 biological
   processes including regulation of keratinocyte proliferation
   (GO:0010837), keratinocyte differentiation (GO:0030216), and epidermis
   development (GO:0008544) as well as immune-related biological pathways
   are the most remarkable ones in the context of cutaneous melanoma.
   Also, the pathway enrichment analysis of lncRNAs in Module B of the
   lncRNA-GO-BP network demonstrated that the coexpressed genes of these
   lncRNAs are most significantly enriched in the “Cytokine Signaling in
   Immune system” pathway. This result indicates that CMK-lncRNAs in
   Module B might collaboratively/synergistically function in the same
   pathway. Fleming et al^[127]58 previously demonstrated that miR-150
   activity in the Cytokine Signaling in Immune system is associated with
   recurrence of melanoma. Also, the relation between “Cytokine Signaling
   in Immune system” pathway and 2 other cancer types including
   hepatocellular carcinoma^[128]59 and prostate cancer^[129]60 has been
   reported. Moreover, the “Cytokine Signaling in Immune system” pathway
   is associated with macrophage polarization and plasticity,^[130]61
   which is an important factor in the development of various malignancies
   including cutaneous melanoma.^[131]62

   CMK-lncRNAs are also associated with clinicopathological features of
   cutaneous melanoma patients. According to the survival analyses, the
   expression level of 4 out of 12 CMK-lncRNAs, including HCP5, LINC00518,
   LINC00964, and MIR205HG, is associated with the OS of cutaneous
   melanoma patients. Thus, these lncRNAs could be collectively used for
   early prognosis and possibly diagnosis of cutaneous melanoma. In
   addition, LINC00518 and upstream negative regulators of 3 other lncRNAs
   could be considered as candidate drug targets for cutaneous melanoma
   patients. Similarly, the expression level of 4 out of 12 CMK-lncRNAs,
   including HCP5, CASC15, SNHG11, and MIR205HG, is associated with
   different stages of the cutaneous melanoma. Accordingly, the expression
   level of these lncRNAs could be considered for characterization of
   different stages of melanoma tumors. Furthermore, based on our results,
   most of the cutaneous melanoma samples with dysregulation in any of the
   differentially expressed CMK-lncRNAs belong to the stage III or IV
   melanoma patients with less than 5 years OS/DFS.

Conclusion

   To sum up, this is a thorough study with the purpose of identifying
   most essential lncRNAs in the development of cutaneous melanoma.
   According to our results, among 3823 lncRNAs, 30 lncRNAs were proposed
   as the ones with most significant associations with cutaneous melanoma.
   Also, 12 out of these 30 lncRNAs were labeled as CMK-lncRNAs, and our
   downstream analyses demonstrated that most of these CMK-lncRNAs are
   associated with cancer-related biological processes. We also
   interrogated the functional roles of CMK-lncRNAs and demonstrated their
   associations with clinicopathological features of cutaneous melanoma
   patients. In short, the CMK-lncRNAs could be considered as candidate
   diagnostic/prognostic biomarkers and/or drug targets for cutaneous
   melanoma patients. Considering all of our results simultaneously, the
   lncRNA LINC00518 is the most remarkable one among all 12 CMK-lncRNAs
   and might be an essential oncogenic factor in the cutaneous melanoma
   development. However, further experimental assays are required to
   evaluate the potential of these lncRNAs in practice and to further
   characterize their precise roles in cutaneous melanoma tumorigenesis
   and metastasis.

Supplemental Material

   sj-pdf-1-bbi-10.1177_1177932220988508 – Supplemental material for
   Functional Prediction of Long Noncoding RNAs in Cutaneous Melanoma
   Using a Systems Biology Approach
   [132]Click here for additional data file.^ (246.6KB, pdf)

   Supplemental material, sj-pdf-1-bbi-10.1177_1177932220988508 for
   Functional Prediction of Long Noncoding RNAs in Cutaneous Melanoma
   Using a Systems Biology Approach by Mozhdeh Shahmoradi and Zahra
   Rezvani in Bioinformatics and Biology Insights

Acknowledgments