Jiaqi Zhang
Hui Liu
Wenhao Zhang
Yinfang Li3
Hua Jiang
Judong Luo- 1Department of Radiotherapy, The Affiliated Changzhou No. 2 People's Hospital of Nanjing Medical University, Changzhou, China
- 2Department of Dermatology, Graduate School of Dalian Medical University, Dalian, China
- 3Aliyun School of Big Data, Changzhou University, Changzhou, China
- 4Department of Oncology, Affiliated 3201 Hospital of Xi'an Jiaotong University, Hanzhong, China
- 5Department of Oncology, The Affiliated Changzhou No.2 People's Hospital of Nanjing Medical University, Changzhou, China
Skin cutaneous melanoma (SKCM) is an aggressive form of skin cancer that results in high mortality rate worldwide. It is vital to discover effective prognostic biomarkers and therapeutic targets for the treatment of melanoma. Long non-coding RNA (lncRNA) has been verified to play an essential role in the regulation of gene expression in diseases and tumors. Therefore, it is significant to explore the function of lncRNAs in the development and progression of SKCM. In this paper, a set of differentially expressed lncRNAs (DElncRNAs) and mRNAs (DEmRNAs) were first screened out using 471 cutaneous melanoma samples and 813 normal skin samples. Gene Ontology and KEGG pathway enrichment analysis were performed to obtain the significant function annotations and pathways of DEmRNAs. We also ran survival analysis on both DElncRNAs and DEmRNAs to identify prognostic-related lncRNAs and mRNAs. Next, a set of hub genes derived from protein-protein interaction (PPI) network analysis and lncRNA target genes screened from starbase-ENCORI database were integrated to construct a lncRNA-mRNA regulatory module, which includes 6 lncRNAs 4 target mRNAs. We further checked the capacity of these lncRNA and mRNA in the diagnosis of melanoma, and found that single lncRNA can effectively distinguish tumor and normal tissue. Moreover, we ran CMap analysis to select a list of small molecule drugs for SKCM, such as EGFR inhibitor AG-490, growth factor receptor inhibitor GW-441756 and apoptosis stimulant betulinic-acid, which have shown therapeutic effect in the treatment of melanoma.
1. Introduction
Malignant melanoma is one of the most aggressive malignancies that have a strong tendency to metastasize during the early stage of the disease (Kremenovic et al., 2020). Although surgical removal of the primary melanoma can significantly improve the survival rate of patients, there are approximate 10% cases diagnosed at advanced stage, by then the tumor has become metastatic and unresectable and thus yields to poor prognosis (Domingues et al., 2018; Kremenovic et al., 2020). Some therapeutic treatments, such as targeted therapy and immunotherapy, have been verified to somewhat improve the prognosis of patients with metastatic melanoma (Bai et al., 2019). However, most melanoma patients still suffer from poor prognosis with only 5-year relative survival rate. Therefore, it is urgent to uncover effective biomarkers for early diagnosis and novel therapeutic targets of melanoma.
The initiation and progression of melanoma arise due to the separate or combined action of genetic and epigenetic factors (Fischer et al., 2018). With the development of high-throughput sequencing technique, there has been significant advance in exploring the underlying molecular mechanisms, as well as the identification of crucial signal transductions and pathways in SKCM (Ko and Fisher, 2011). Many protein-coding genes have been shown to correlate with melanoma. Several growth factors, including basic fibroblast growth factor (bFGF), epidermal growth factor (EGF), platelet-derived growth factor (PDGF), transforming growth factors (TGF), and insulin-like growth factors 1 and 2 (IGF-1 and IGF-2), have been reported to be significantly up-regulated in melanoma cells (Polsky and Cordon-Cardo, 2003; Menezes et al., 2018). Mutations in the CDKN2A gene are the most common alteration in hereditary melanoma (Aoude et al., 2015). This gene encodes the p16 protein, which inhibits cyclin-dependent kinase (CDK) 4 and 6, and reduces cell cycle arrest and apoptosis (Soura et al., 2016). Moreover, it is reported that SOX10 is as equally specific as S100 for the detection of melanoma metastases in sentinel lymph node (Szumera-Ciekiewicz et al., 2020). Besides, some well-known melanoma markers HMB-45 and Melan-A are often used in the clinical diagnosis.
Meanwhile, more and more studies have shown that long non-coding RNAs (lncRNAs) play an important role in various biological processes (Chen et al., 2019). LncRNA is a type of RNA transcripts with non-coding potentials more than 200 nucleotides, which exhibits histone modification functions similar to protein-coding genes. LncRNAs comprise the majority of transcripts in the mammalian genomes, however, their functions still remain largely unknown (Ramilowski et al., 2020). Huang et al. reported that lncRNAs may participate in every step of gene expression via transcriptional, histone modification, post-transcriptional, and/or translational regulation (Huang et al., 2018). Some lncRNAs take effect in gene regulation and thus influence various facets of cellular homeostasis, including proliferation, survival, migration, and genomic stability (Ramilowski et al., 2020). Besides, they can modulate gene expression through transcriptional activation/repression and RNA editing/splicing/degradation (Moran et al., 2012).
Therefore, an increasing number of studies have poured attention to the correlation between lncRNA and cancer. For example, it has been reported that lncRNAs play crucial functional roles in tumorigenesis, including melanoma (Huarte, 2015). In an early study that identifies the lncRNA signature of melanoma metastasis, several metastasis-associated lncRNAs were tested for their possible function in lymphatic metastasis of melanoma. As a result, relative to matched primary melanomas, lncRNA HOTAIR was found to be highly expressed in lymph node metastases (Tang et al., 2013). Moreover, it has been shown that MALAT1 is intensively related to lymph node metastasis, and UCA1 is correlated with advanced melanoma (Tian et al., 2014), which indicates that lncRNAs may be effective biomarkers for identifying metastasis of SKCM. As lncRNA functions through regulating mRNA processing and post-transcriptional regulation (Geisler and Coller, 2013), lncRNA and mRNA may influence each other, thereby constitute complex regulation layers. The process of transcription and/or splicing of the lncRNAs confers a gene-regulation functionality. Transcription of lncRNA may influence local chromatin states and transcription factor (TF) binding on promoter and enhancer regions (Kopp and Mendell, 2018). For example, it is reported that impaired transcription or splicing of Blustr changed the chromatin state of the Sfmbt2 promoter, as well as decreasing RNA polymerase occupancy at the transcription start site and within the gene body of Sfmbt2. Transcription and splicing of the Blustr RNA is closely related to the expression of neighboring gene (Engreitz et al., 2016). However, the interactions between lncRNAs and mRNAs still remain unknown in melanoma to date. Investigation of the underlying molecular mechanisms of lncRNA participating in melanoma is helpful to preventing neoplasm metastasis and improve prognosis of melanoma patients.
In this paper, we set about to investigate the interactions between aberrantly expressed lncRNA and mRNA to uncover their regulation functions in SKCM. We collected genome-wide expression profiles from TCGA and GTEx database, and conducted bioinformatics analysis to explore the function of differentially expressed lncRNAs and mRNAs. A set of key lncRNAs and protein-coding genes were identified via survival analysis. GO functional annotations and KEGG pathways were subsequently identified via gene set enrichment analysis. By integrating lncRNA targets and hub genes of the PPI network, we established an essential lncRNA-mRNA regulatory module that underlie the occurrence and progression of SKCM, to provide new insight into the mechanism of lncRNAs and potential treatment of SKCM. Moreover, connectivity map (CMap) analysis yielded to several drugs that can reverse the expression profiles of the set of key genes and thus indicate potential treatment of melanoma.
2. Materials and Methods
2.1. Data Source
To evaluate the difference of mRNA and lncRNA expressions between the diseased and normal samples, we obtained the fragments per kilobase million (FPKM) values of mRNA and lncRNA expressions from UCSC Xena (Goldman et al., 2020), covering 471 SKCM samples and 1 normal sample from The Cancer Genome Atlas (TCGA) (Tomczak et al., 2015). Besides, we obtained 812 normal samples from Genotype Tissue Expression (GTEx) (Lonsdale et al., 2013).
As a result, a total of 471 tumor samples and 813 normal samples were collected. Clinical data with complete survival information of 463 SKCM patients were also retrieved from UCSC Xena. The expression profiles included 13,840 mRNAs and 2,249 lncRNAs, and the expression levels were calibrated and standardized.
2.2. Differential Expression Analysis
The Limma package (version 3.44.3) (Ritchie et al., 2015) was used to identify the differentially expressed lncRNA and mRNA. The log2 fold change >2 and p-value < 0.01 were used as the criteria to screen differentially express genes. The R packages ggplot2 and pheatmap were used to plot the volcano maps and heat maps of the two groups of DEGs.
2.3. GO and KEGG Pathway Enrichment Analysis
The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were conducted using the set of DEmRNAs for functional annotations and pathway enrichment analysis. The clusterProfiler (version 3.16.0) (Yu et al., 2012) was used to perform functional annotations and pathway analysis. Three type functional annotations, including biological process (BP), cellular composition (CC), and molecular function (MF), were covered. The threshold p-value < 0.05 was applied to select statistically significant functional terms and pathways.
2.4. Survival Analysis
To screen out key factors from the set of differentially expressed lncRNAs and mRNAs, survival analysis was conducted based on the melanoma cohort in TCGA. According to median levels, patients were divided into two groups with high and low expression levels. Overall survival (OS) rate analysis was then performed using the R package survival (Therneau and Lumley, 2014), based on which the OS differences between the two groups were estimated. Finally, the Kaplan Meier method (Efron, 1988) with a log-rank test was applied to evaluate the significance of difference, and the genes with p-value < 0.01 were considered statistically significant.
2.5. PPI Network and Hub Gene Identification
The protein-protein interactions (PPIs) were derived from STRING database (Szklarczyk et al., 2016). Cytoscape (Smoot et al., 2011) greatly facilitates the PPI analysis by readily visualizing the PPI network, in which each node represents a protein or a gene and an edge represent the interactions between two nodes. The degree of a node is the number of their connections with other nodes. Based on the PPI network, we used MCODE1 to screen out the important modules and then used cytoHubba2 to calculate the degree value. Only the nodes with degree score ≥20 were regarded as important ones.
2.6. Connectivity Map-Based Drug Screening
The connectivity map (CMap) is a large-scale online database that contains whole genomic expression profiles of cultured cancer cells treated with bioactive small molecules, including the drugs with clinical utility, known mechanism of action, or nomination from the NIH Molecular Libraries Program (Subramanian et al., 2017). By submitting the set of deregulated genes derived from PPI analysis to CMap web server, we obtained a list of drugs with connectivity scores scaled from –100 to 100, which indicate the functional connections among drugs, genes, and diseases. A positive score implies that an agent might facilitate the expression of the gene signature, while a negative score indicates that an agent might repress or reverse the expression of the gene signature (Zhang et al., 2020).
2.7. Construction of lncRNA-mRNA Network
Starbase-ENCORI (Li et al., 2014), an platform for studying the miRNA-ncRNA, miRNA-mRNA, and ncRNA-RNA interactions from CLIP-seq, degradome-seq and RNA-RNA interactions, was used to screen target genes of the DElncRNAs. The standard for selection of the genes was that the number of interactions is no less than 2 and at least 1 experiments had been performed. The overlapped mRNAs between the DEmRNAs and the lncRNA-targeted mRNAs were selected and used to construct the lncRNA-mRNA regulatory module. The Cytoscape was also employed to visualize the generation of the network.
3. Results
3.1. Differentially Expressed mRNAs Correlated to Melanoma
The differential expression analysis revealed 146 lncRNAs and 1,485 that are differentially expressed between tumor and normal samples. Among them, 100 lncRNAs and 514 mRNAs are up-regulated in SKCM, and 46 lncRNAs and 972 mRNAs are down-regulated. As is illustrated in Figure 1, the differentially expressed lncRNA and mRNAs are visually displayed in the form of volcano maps and heat maps. The top 20 significantly deregulated lncRNAs and mRNAs are shown in Tables 1, 2, respectively.
Figure 1. Differentially expressed genes between normal samples and melanoma patients (log2 fold change>2, p-value < 0.01). (A) Volcano plot of the p-value as a function of weighted fold change for DEmRNAs. (B) DElncRNAs, red dots represent significantly upregulated expressed genes and green dots represent significantly downregulated expressed genes (log2 fold change>2, p-value < 0.01). (C) Heat map for potential mRNAs (n = 1,485) showed significant expression changes, in which 971 were downregulated and 514 were upregulated. Red through green color indicates high to low expression level. (D) Heat map for potential lncRNAs (n = 146) showed significant expression changes, in which 46 were downregulated and 100 were upregulated.
Table 1. TOP 10 up-regulated/down-regulated lncRNAs in melanoma compared to normal samples.
Table 2. TOP 10 up-regulated/down-regulated mRNAs in melanoma compared to normal samples.
To uncover the role of the differentially expressed mRNAs (DEmRNAs) in the pathogenesis of melanoma, GO enrichment analysis was carried out. As shown in Figure 2, DEmRNAs are significantly enriched in skin epidermis-related biological processes, such as epidermis and skin development, keratinocyte differentiation, cornification, and epidermal cell differentiation. As for cellular component, these DEmRNAs are significantly enriched in the cornified envelope, collagen-containing extracellular matrix, melanosome, and pigment granule. For molecular function, they are significantly enriched in peptidase and endopeptidase-related activity and structural constituent of skin epidermis, as shown in Figures 2A,B. These enrichment analysis results indicate that DEmRNAs associated with epidermal cell structure and differentiation might contribute to melanoma tumorigenesis. For example, Kodet et al. report that melanoma cells are able to influence locally the differentiation pattern of keratinocytes in vivo as well as in vitro. They have shown that FGF-2, CXCL-1, IL-8, and VEGF-A participate in the activity of melanoma cells on keratinocytes (Kodet et al., 2015). Moreover, the extracellular matrix, which was known as major component of the local microenvironment, is commonly deregulated and becomes disorganized in cancer (Lu et al., 2012).
Figure 2. GO and KEGG enrichment analyses. (A,B) show the bar plot and cluster plot of significant GO funtional items. (C,D) show the bar plot and cluster plot of significant KEGG pathways.
According to the KEGG enrichment analysis, top 25 significantly enriched pathways are involved in the process of cancerogenesis. More specifically, DEmRNAs are mainly enriched in pathway related to immune and inflammatory response, including Staphylococcus aureus infection, antigen processing and presentation, viral protein interaction with cytokine and cytokine receptor. They are also involved in the Th17 cell differentiation and rheumatoid arthritis, as shown in Figures 2C,D. These results strongly indicate that the DEmRNAs may play important role in melanoma by involving in the inflammation and immune response.
3.2. Differentially Expressed lncRNAs and mRNAs Associated to Survival Rate
We evaluated the relationship between DElncRNAs and the clinical outcome of patients with melanoma. By performing Kaplan Meier survival curves analysis based on the survival data from TCGA, We found more than 10 lncRNAs were statistically significant in association with patient survival rate. The survival curves of the most significant 4 lncRNAs are illustrated in Figure 3. LIMD1-AS1 (p-value = 0.00137) signature has been shown to predict prognostic survival and direct clinical risk-specific treatments in melanoma (Liu et al., 2019). LINC00518 (p-value = 0.00803) has shown to act as a competing endogenous RNA to promote the metastasis of malignant melanoma via miR-204-5p/AP1S2 axis (Luan et al., 2019). Also, LINC00520 (p-value = 0.00931) promotes the proliferation and metastasis of malignant melanoma by inducing the miR-125b-5p/EIF5A2 axis (Luan et al., 2020). Another weighted gene co-expression network analysis show that the PI3K subunit PI3KCD exhibited excellent efficacy for diagnosing primary and metastatic tumor tissue (Wang et al., 2019).
Figure 3. Survival curves of four lncRNAs with the significant logrank p-value.
We also checked whether the DEmRNAs were statistically associated to patient survival rate. The top eight significant mRNAs are shown in Supplementary Figure 1. The expression profiles of GBP4 (p = 3.131e-08), FCGR2A (p = 5.825e-08), APOL1 (p = 1.408e-07), HLA-DRA (p = 4.035e-7), and HLA-DRB1 (p = 5.24e-07) are positively correlated with the survival rate, while NCCRP1 (p = 9.205e-08), RABIF (p = 3.505e-7), and FGD1 (p = 3.917e-07) are correlated with adverse prognosis. Among these genes, several have been reported to be associated with the prognosis of SKCM. Take GBP4 as an example, it has been verified that its high expression, together with GBP1, GBP2, GBP3, and GBP5, was correlated with favorable overall survival (OS) in the SKCM patients followed for a over 30-year study (Wang et al., 2018). Also, Brunner et al. identified GBP4 as one of the predictive genes whose expression are correlated with OS in cutaneous melanoma (Brunner et al., 2013). Also, the immune regulation pathway that FCGR2A participates in was reported to show significant association with cutaneous malignant melanoma risk (Yang et al., 2009). The functional mechanism study of FCGR2A shows that it binds more avidly to human IgG1 and IgG2 subtypes, thereby increases ADCC-mediated cell death and improve clinical outcomes (Vargas et al., 2018). The expression of class II MHC antigens, including HLA-DRA and HLA-DRB, has been chronically considered a crucial step in immune response toward colorectal carcinomas (Lee et al., 2012). Furthermore, the analysis showed that HLA-DR expression in melanoma cells may be a biomarker for tumors primed with activated T-cells and an appropriate IFNγ response to mediate sensitivity to PD-1/PD-L1 blockade, which is also associated with immune response (Johnson et al., 2016). For the genes related to adverse prognosis, NCCRP1 was also identified as an independent risk prognostic predictor for metastatic melanoma (Sun et al., 2019). FGD1 has been verified that it plays a direct role in the proliferation or invasion of melanoma cells (Hou et al., 2012), thereby the patients with high FGD1 expression frequently showed worse prognosis (Zeng et al., 2020). The complete results of the survival analysis for DElncRNAs and DEmRNAs are shown in Supplementary Table 1.
3.3. lncRNA mRNA Regulatory Module
PPI network analysis was conducted to reveal crucial genes that are systematically related to melanoma and interconnections among them. Based on the 332 DEmRNAs that have significance in survival and the PPIs derived from STRING, 252 genes were further selected out by using Cytoscape. As a result, a PPI network that contains 252 nodes and 1,566 edges was constructed, as shown in Figure 4A.
Figure 4. PPI network of differentially expressed mRNAs and module analysis. (A) PPI network contains 272 nodes and 1,577 edges. (B,C) display two significant modules identified by MCODE, respectively.
Next, we chose hub genes that have associations with many other genes and play key role, as they are considered as potential drivers of the development of diseases (Xu et al., 2020). This step yielded to 39 hub genes, including 22 up-regulated genes and 17 down-regulated genes. Meanwhile, two modules were found out. The first module consists of 22 nodes and 215 edges (score = 20.476), and the second consists of 17 nodes and 103 edges (score = 12.875), as shown in Figures 4B,C, respectively. Note that the darker a node color, the higher its degree.
To further explore the key regulatory factors in SKCM, target genes of the DElncRNAs were predicted by using Starbase-ENCORI database. We got 3,069 unique genes after removal of duplicates, and found four genes (CD52, CD53, DSC3, and LAPTM5) overlapped with the set of DEmRNAs. Subsequently, we constructed a lncRNA mRNA regulatory module that consists of 10 nodes (6 lncRNAs and 4 mRNAs) and 8 edges, as is shown in Figure 5A. The survival curves of the four key genes are shown in Figure 5B. The expression levels of CD52, CD53, and LAPTM5 are positively correlated with the survival rate, while DSC3 is correlated with adverse prognosis.
Figure 5. lncRNA-mRNA regulatory network and survival curves of four hub genes. (A) LncRNA-mRNA regulatory module, in which the red nodes represent up-regulated genes and green nodes represent down-regulated genes, triangles represent lncRNAs and circles represent mRNAs. (B) Survival analysis of the four overlapping mRNAs.
We have conducted extensive literature search to verify the lncRNA mRNA regulatory module. The lncRNA NEAT1 has been reported to facilitate melanoma cell proliferation, migration, and invasion via regulating miR-495-3p and E2F3 (Xia et al., 2019). In fact, NEAT1/miR-23a-3p/KLF3 has been found to form a novel regulatory axis in melanoma cancer progression (Ding et al., 2019). FOXD2-AS1 has also been found to be upregulated in cutaneous melanoma tissue specimens and cell lines compared in normal tissue and normal human epidermal melanocyte. Inhibition of FOXD2-AS1 can suppress cutaneous melanoma cell proliferation, migration and invasion through regulating phospho-Akt expression (Ren et al., 2019). Furthermore, MALAT1 downregulates miR-34a expression level in melanoma cells and tumor specimens (Li et al., 2019), thereby it is a potential biomarker for many human cancer diagnosis as well as prognosis (Amodio et al., 2018; Li et al., 2018).
In this regulatory module we found that the expression of most key lncRNAs and their target mRNAs are negatively regulated. First, lncRNA is often located in the upstream promoter region of coding gene and inhibit the expression of adjacent genes. Second, Trans-acting lncRNAs may also function by modulating the activity or abundance of proteins or RNAs to which they directly bind. Third, lncRNA exists on the antisense strand of the coding gene, which can interfere with the splicing of mRNA, thereby produce different forms of isoforms. lncRNA may also produce endogenous siRNA under the action of dicer enzyme to regulate gene expression level. In addition to RNA-binding proteins, lncRNAs also have the ability to regulate the abundance or activity of microRNA (miRNA), a category of transcripts termed competing endogenous RNAs (ceRNAs) (Kopp and Mendell, 2018).
Although lncRNAs mainly inhibit the expression level of their target mRNAs, our results demonstrated that not all the key lncRNAs and their target mRNAs are negatively correlated. For example, in our lncRNA-mRNA regulatory module, lncRNA MALAT1 can regulate the downstream gene LAPTM5 that is actually up-regulated, meanwhile MALAT1 may also down-regulated gene DSC3. In addition, it has been shown that the down-regulation function of MALAT1 restrain the development of uveal melanoma via inhibition of HOXC4. Therefore, we suppose that the regulatory relationship between lncRNA and mRNA is not absolutely positive or negative.
Also, previous studies have shown that some of the deregulated mRNAs are involved in melanoma. For example, DSC3 was upregulated in primary melanoma than in metastatic melanoma, its high expression level was correlated with adverse prognosis. and thus considered as new biomarkers in the therapeutics of metastatic melanoma (Sheng et al., 2020). Salerno et al. reported that DSC3 is associated with the lack of Th1 immune signatures in human melanoma metastases (Salerno et al., 2016). CD53 has significant homology to melanoma-specific antigen ME491 sequence (Bellacosa et al., 1991). Also, the expression of LAPTM5 closely related to tumorigenesis in many human cancers (Nuylan et al., 2016). For example, Somura et al. report that a 3-gene predictor which contains LAPTM5 may predict the early intrahepatic recurrence of hepatocellular carcinoma (Somura et al., 2008).
3.4. Key lncRNAs Function as Biomarker of Melanoma
To reveal the role of the key lncRNAs in the pathogenesis of melanoma, enrichment analysis was carried out by target genes of them. In terms of GO enrichment analysis results, the set of target genes are participated in biological processes such as SRP-dependent cotranslational protein targeting to membrane and viral gene expression. For cellular component, these target genes are significantly enriched in the cytosolic ribosome, large ribosomal subunit and MHC protein complex. As for molecular function, they are significantly enriched in structural constituent of ribosome and MHC class II protein complex binding, as shown in Figures 6A,B.
Figure 6. GO and KEGG enrichment analyses of key lncRNAs targeting mRNAs. (A,B) show the bar plot and cluster plot of significant GO funtional items. (C,D) show the bar plot and cluster plot of significant KEGG pathways.
As for KEGG enrichment analysis, top 25 significantly enriched pathways are shown in Figures 6C,D. LncRNA target genes are significantly involved in the process of antigen processing and presentation, ribosome, phagosome and allograft rejection pathway. The enrichment analysis results powerfully indicate that these genes may contribute to melanoma tumorigenesis. There is evidence for dysregulation of ribosome biogenesis in cancers. Ribosome expression is closely connected to cell growth and proliferation. The number of ribosomes per cell is proportional to the growth rate of that cell (Scull et al., 2019). MHC-II play an important part in antigen presentation to CD4+ T-lymphocytes and anti-tumor immunity, increasing evidence indicates that tumor-specific MHC-II related to favorable outcomes in patients with cancer, including those treated with immunotherapies, and with tumor rejection in murine models (Axelrod et al., 2018). These results strongly indicate that the DEmRNAs may play important role in cancerogenesis.
To demonstrate the effectiveness of six lncRNAs as biomarker for diagnosis of melanoma, we used the expression levels of each lncRNA to distinguish tumor and normal tissues. With the varying threshold, we can build the ROC curve with respect to each lncRNA and computed the AUC values, as shown in Figure 7.
Figure 7. ROC curves and the corresponding AUC values of the 6 lncRNAs and 4 mRNAs as biomarker of SKCM.
In fact, the average AUC of these regulatory factors reaches 0.97, among which FOXD2-AS1 (AUC: 0.973), MALAT1 (AUC: 0.989), NEAT1 (AUC: 0.988), AC245595.1 (AUC: 0.985), and SNHG5 (AUC: 0.964) are the most significant lncRNAs to predict diagnosis in SKCM. Previous studies have shown that these factors involved in melanoma. For example, Ren et al. have revealed that melanoma cell proliferation, migration and invasion can be suppressed by the inhibition of FOXD2-AS1 that regulates phospho-Akt expression in cutaneous melanoma. This loss-of-function study indicated that FOXD2-AS1 functions as oncogenic lncRNA in cutaneous melanoma cells (Ren et al., 2019). Also, the role of MALAT1 in metastasis of melanoma is a hot topic that raises much attention. Tian et al. have verified that knockdown of MALAT1 can attenuate the migrational ability of melanoma cells via in vitro studies as early as 2014, indicating the correlation between MALAT1 and melanoma metastasis (Tian et al., 2014). Thereafter, Luan et al. reported that MALAT1 promotes malignant melanoma growth and metastasis by sponging miR-22 (Luan et al., 2016). Sun et al. subsequently showed that MALAT1 can promote the proliferation, invasion and migration of melanoma cells by suppressing miR-140 expression (Sun et al., 2016). Another more recent study revealed that MALAT1 promotes melanoma development by downregulating miR-23a (Wang et al., 2020). NEAT1 was also proved to facilitate the melanoma cell proliferation, migration, and invasion via regulating miR-495-3p and E2F3 (Xia et al., 2019). SNHG5 plays an analogous role in melanoma by regulating the miR-26a-5p/TRPC3 pathway (Gao et al., 2019). In summary, all these findings strongly support these key lncRNAs are closely related to melanoma, and can be used as biomarker for diagnosis of melanoma.
3.5. CMap Analysis Reveals Potential Drugs
We went further to explore the small molecule drugs that can reverse the expression levels of 39 identified hub genes so that the development and progression of melanoma can be potentially inhibited. We submitted the set of genes to the Connectivity Map web server and choose skin melanoma cell line A375 to screen drugs. The returned drugs were prioritized according to their scores, and negative score indicates that drugs might repress or reverse the expression of the key genes. The top 15 small molecules with negative scores are listed in Table 3.
Table 3. Predicted small molecule drugs by CMAP analysis.
Expectedly, quite a few small molecules have been approved to inhibit the development of SKCM. For instance, EGFR inhibitor AG-490 can effectively suppress tumor cell proliferation by limiting the expression of cyclin D1 (Kamran et al., 2013). Betulinic acid has been demonstrated to induce programmed cell death with melanoma cells. Betulinic acid treatment of cultured UISO-Mel-1 (human melanoma cells) leads to the activation of p38 and stress activated protein kinase/c-Jun NH(2)-terminal kinase which widely accepted as proapoptotic mitogen-activated protein kinases (MAPKs) (Tan et al., 2003). In addition, HDAC and PARP inhibitors are also well-known as tumor-targeted drugs. More importantly, the potential effectiveness of these molecules implies that the set of hub mRNAs, based on which the CMAP analysis was conducted, are highly relates to the incidence and progression of SKCM, indicating the rationality of the bioinformatics pipeline presented in this paper.
4. Conclusion and Discussion
In this paper, we aim at exploring the underlying regulation network regarding to lncRNA and mRNA in SKCM using bioinformatics analysis. Firstly, a total of 146 lncRNAs and 1,485 mRNAs that are differentially expressed in SKCM tissues are identified. We use the set of DEmRNAs to perform Gene Ontology and KEGG enrichment analysis for exploring relevant biological functions at transcriptional level. The functional and pathway annotations demonstrate that the aberrantly expressed genes participate in melanoma-related biological processes, such as epidermis and skin development, keratinocyte differentiation, and cornification. Also, cellular compositions and molecular functions are related to cornified envelope, collagen-containing extracellular matrix, melanosome and peptidase inhibitor activity (Gajewski, 2007).
Next, survival analysis reveals a number of DEmRNAs that exhibit significant effects on melanoma patient survival rate, and PPI network analysis yields to 39 hub genes and two significant modules are identified. Finally, a central lncRNA mRNA module that may underlie the complex regulatory relationship between lncRNAs and mRNAs in melanoma is constructed. The lncRNA mRNA regulatory module including six key lncRNAs and four key genes are presented. We have found that the lncRNAs can be used as biomarker that effectively distinguish the tumor and normal tissue. Also, the connectivity map analysis screens a few small molecule drugs that can reverse the expression profiles of the 39 hub genes. Furthermore, our extensive literature mining also verifies that these lncRNAs and mRNAs essentially participate in the pathogenesis and prognosis of melanoma. These results demonstrate the effectiveness of our constructed lncRNA-mRNA regulatory network. We believe this regulatory module would provide novel insights of the biological mechanism in SKCM.
A limitation of our study mainly lies in that most of the data used in the bioinformatic analysis are obtained from public datasets, without verification of in vitro or in vivo biochemical experiments. Despite this limitation, we are the first to identify differentially expressed genes in skin melanoma by integrating the data of TCGA and GTEx databases. With the constructed lncRNA mRNA regulatory module, we can further comprehensively characterize the lncRNAs and mRNAs associated with melanoma, and provide new insight into the regulatory mechanism for melanoma.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author/s.
Author Contributions
JZ and HL conceived the main idea and analyzed the data. JZ and WZ drafted the manuscript. HL reviewed drafts of the paper and improved the manuscript. YL collected the data and performed the statistical analysis. ZF provided statistical advice. HJ and JL supervised the study and provided funding. All authors read and commented on the manuscript.
Funding
This work was supported by National Natural Science Foundation of China (82073339; 62072058; 81773224; 61672113), and the Scientific Projects of Jiangsu Province (BE2018643; BK20191157).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2020.615671/full#supplementary-material
Footnotes
References
Amodio, N., Raimondi, L., Juli, G., Stamato, M. A., Caracciolo, D., Tagliaferri, P., et al. (2018). Malat1: a druggable long non-coding RNA for targeted anti-cancer approaches. J. Hematol. Oncol. 11:63. doi: 10.1186/s13045-018-0606-4
Aoude, L. G., Wadt, K. A. W., Pritchard, A. L., and Hayward, N. K. (2015). Genetics of familial melanoma: 20years after cdkn2a. Pigment Cell Melanoma Res. 28, 148–160. doi: 10.1111/pcmr.12333
Axelrod, M. L., Cook, R. S., Johnson, D. B., and Balko, J. M. (2018). Biological consequences of major histocompatibility class-II expression by tumor cells in cancer. Clin. Cancer Res. 25, 2392–2402. doi: 10.1158/1078-0432.CCR-18-3200
Bai, X., Fisher, D. E., and Flaherty, K. T. (2019). Cell-state dynamics and therapeutic resistance in melanoma from the perspective of mitf and ifnγ pathways. Nat. Rev. Clin. Oncol. 16, 549–562. doi: 10.1038/s41571-019-0204-6
Bellacosa, A., Lazo, P. A., Bear, S. E., and Tsichlis, P. N. (1991). The rat leukocyte antigen mrc ox-44 is a member of a new family of cell surface proteins which appear to be involved in growth regulation. Mol. Cell. Biol. 11, 1864–1872. doi: 10.1128/MCB.11.5.2864
Brunner, G., Reitz, M., Heinecke, A., Lippold, A., Berking, C., Suter, L., et al. (2013). A nine-gene signature predicting clinical outcome in cutaneous melanoma. J. Cancer Res. Clin. Oncol. 139, 249–258. doi: 10.1007/s00432-012-1322-z
Chen, X., Sun, Y.-Z., Guan, N.-N., Qu, J., Huang, Z.-A., Zhu, Z.-X., et al. (2019). Computational models for lncRNA function prediction and functional similarity calculation. Brief. Funct. Genomics 18, 58–82. doi: 10.1093/bfgp/ely031
Ding, F., Lai, J., Gao, Y., Wang, G., Shang, J., Zhang, D., et al. (2019). Neat1/mir-23a-3p/klf3: a novel regulatory axis in melanoma cancer progression. Cancer Cell Int. 19:217. doi: 10.1186/s12935-019-0927-6
Domingues, B., Lopes, J. M., Soares, P., and Pópulo, H. (2018). Melanoma treatment in review. Immunotargets Ther. 7:35. doi: 10.2147/ITT.S134842
Efron, B. (1988). Logistic regression, survival analysis, and the kaplan-meier curve. J. Am. Stat. Assoc. 83, 414–425. doi: 10.1080/01621459.1988.10478612
Engreitz, J. M., Haines, J. E., Perez, E. M., Munson, G., Chen, J., Kane, M., et al. (2016). Local regulation of gene expression by lncRNA promoters, transcription and splicing. Nature 539, 452–455. doi: 10.1038/nature20149
Fischer, G. M., Vashisht Gopal, Y., McQuade, J. L., Peng, W., DeBerardinis, R. J., and Davies, M. A. (2018). Metabolic strategies of melanoma cells: mechanisms, interactions with the tumor microenvironment, and therapeutic implications. Pigment Cell Melanoma Res. 31, 11–30. doi: 10.1111/pcmr.12661
Gajewski, T. F. (2007). Failure at the effector phase: immune barriers at the level of the melanoma tumor microenvironment. Clin. Cancer Res. 13, 5256–5261. doi: 10.1158/1078-0432.CCR-07-0892
Gao, J., Zeng, K., Liu, Y., Gao, L., and Liu, L. (2019). LncRNA snhg5 promotes growth and invasion in melanoma by regulating the mir-26a-5p/trpc3 pathway. Oncotargets Ther. 12:169. doi: 10.2147/OTT.S184078
Geisler, S., and Coller, J. (2013). RNA in unexpected places: long non-coding RNA functions in diverse cellular contexts. Nat. Rev. Mol. Cell Biol. 14, 699–712. doi: 10.1038/nrm3679
Goldman, M. J., Craft, B., Hastie, M., Repečka, K., McDade, F., Kamath, A., et al. (2020). Visualizing and interpreting cancer genomics data via the xena platform. Nat. Biotechnol.38, 675–678. doi: 10.1038/s41587-020-0546-8
Hou, P., Liu, D., Dong, J., and Xing, M. (2012). The brafv600e causes widespread alterations in gene methylation in the genome of melanoma cells. Cell Cycle 11, 286–295. doi: 10.4161/cc.11.2.18707
Huang, Q., Yan, J., and Agami, R. (2018). Long non-coding RNAs in metastasis. Cancer Metast. Rev. 37, 75–81. doi: 10.1007/s10555-017-9713-x
Huarte, M. (2015). The emerging role of lncRNAs in cancer. Nat. Med. 21, 1253–1261. doi: 10.21037/atm.2016.06.27
Johnson, D. B., Estrada, M. V., Salgado, R., Sanchez, V., Doxie, D. B., Opalenik, S. R., et al. (2016). Melanoma-specific mhc-II expression represents a tumour-autonomous phenotype and predicts response to anti-pd-1/pd-l1 therapy. Nat. Commun. 7, 1–10. doi: 10.1038/ncomms10582
Kamran, M. Z., Patil, P., Shirsath, K., and Gude, R. P. (2013). Tyrosine kinase inhibitor ag490 inhibits the proliferation and migration and disrupts actin organization of cancer cells. J. Environ. Pathol. Toxicol. Oncol. 32, 361–371. doi: 10.1615/JEnvironPatholToxicolOncol.2013010051
Ko, J. M., and Fisher, D. E. (2011). A new era: melanoma genetics and therapeutics. J. Pathol. 223, 242–251. doi: 10.1002/path.2804
Kodet, O., Lacina, L., Krejčí, E., Dvořánková, B., Grim, M., Štork, J., et al. (2015). Melanoma cells influence the differentiation pattern of human epidermal keratinocytes. Mol. Cancer 14:1. doi: 10.1186/1476-4598-14-1
Kopp, F., and Mendell, J. T. (2018). Functional classification and experimental dissection of long noncoding RNAs. Cell 172, 393–407. doi: 10.1016/j.cell.2018.01.011
Kremenovic, M., Schenk, M., and Lee, D. (2020). Clinical and molecular insights into BCG immunotherapy for melanoma. J. Intern. Med. doi: 10.1111/joim.13037. [Epub ahead of print].
Lee, J., Li, L., Gretz, N., Gebert, J., and Dihlmann, S. (2012). Absent in melanoma 2 (aim2) is an important mediator of interferon-dependent and-independent hla-dra and hla-drb gene expression in colorectal cancers. Oncogene 31, 1242–1253. doi: 10.1038/onc.2011.320
Li, F., Li, X., Qiao, L., Liu, W., Xu, C., and Wang, X. (2019). Malat1 regulates mir-34a expression in melanoma cells. Cell Death Dis. 10, 1–11. doi: 10.1038/s41419-019-1620-3
Li, J.-H., Liu, S., Zhou, H., Qu, L.-H., and Yang, J.-H. (2014). starbase v2. 0: decoding miRNA-ceRNA, miRNA-ncRNA and protein–RNA interaction networks from large-scale clip-seq data. Nucleic Acids Res. 42, D92–D97. doi: 10.1093/nar/gkt1248
Li, Z.-X., Zhu, Q.-N., Zhang, H.-B., Hu, Y., Wang, G., and Zhu, Y.-S. (2018). Malat1: a potential biomarker in cancer. Cancer Manage. Res. 10:6757. doi: 10.2147/CMAR.S169406
Liu, N., Liu, Z., Liu, X., and Chen, H. (2019). Comprehensive analysis of a competing endogenous rna network identifies seven-lncRNA signature as a prognostic biomarker for melanoma. Front. Oncol. 9:935. doi: 10.3389/fonc.2019.00935
Lonsdale, J., Thomas, J., Salvatore, M., Phillips, R., Lo, E., Shad, S., et al. (2013). The genotype-tissue expression (gtex) project. Nat. Genet. 45, 580–585. doi: 10.1038/ng.2653
Lu, P., Weaver, V. M., and Werb, Z. (2012). The extracellular matrix: a dynamic niche in cancer progression. J. Cell Biol. 196, 395–406. doi: 10.1083/jcb.201102147
Luan, W., Ding, Y., Ma, S., Ruan, H., and Lu, F. (2019). Long noncoding RNA linc00518 acts as a competing endogenous RNA to promote the metastasis of malignant melanoma via mir-204-5p/ap1s2 axis. Cell Death Dis. 10:855. doi: 10.1038/s41419-019-2090-3
Luan, W., Ding, Y., Yuan, H., Ma, S., Ruan, H., Wang, J., et al. (2020). 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. 39:96. doi: 10.1186/s13046-020-01599-7
Luan, W., Li, L., Shi, Y., Bu, X., Xia, Y., Wang, J., et al. (2016). Long non-coding rna malat1 acts as a competing endogenous RNA to promote malignant melanoma growth and metastasis by sponging mir-22. Oncotarget 7:63901. doi: 10.18632/oncotarget.11564
Menezes, M. E., Talukdar, S., Wechman, S. L., Das, S. K., Emdad, L., Sarkar, D., et al. (2018). Prospects of gene therapy to treat melanoma. Adv. Cancer Res. 138, 213–237. doi: 10.1016/bs.acr.2018.02.007
Moran, V. A., Perera, R. J., and Khalil, A. M. (2012). Emerging functional and mechanistic paradigms of mammalian long non-coding RNAs. Nucleic Acids Res. 40, 6391–6400. doi: 10.1093/nar/gks296
Nuylan, M., Kawano, T., Inazawa, J., and Inoue, J. (2016). Down-regulation of laptm5 in human cancer cells. Oncotarget 7:28320. doi: 10.18632/oncotarget.8614
Polsky, D., and Cordon-Cardo, C. (2003). Oncogenes in melanoma. Oncogene 22, 3087–3091. doi: 10.1038/sj.onc.1206449
Ramilowski, J. A., Yip, C. W., Agrawal, S., Chang, J.-C., Ciani, Y., Kulakovskiy, I. V., et al. (2020). Functional annotation of human long noncoding RNAs via molecular phenotyping. Genome Res. 30, 1060–1072. doi: 10.1101/gr.254219.119
Ren, W., Zhu, Z., and Wu, L. (2019). Foxd2-as1 correlates with the malignant status and regulates cell proliferation, migration, and invasion in cutaneous melanoma. J. Cell. Biochem. 120, 5417–5423. doi: 10.1002/jcb.27820
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007
Salerno, E. P., Bedognetti, D., Mauldin, I. S., Deacon, D. H., Shea, S. M., Pinczewski, J., et al. (2016). Human melanomas and ovarian cancers overexpressing mechanical barrier molecule genes lack immune signatures and have increased patient mortality risk. Oncoimmunology 5:e1240857. doi: 10.1080/2162402X.2016.1240857
Scull, C. E., Zhang, Y., Tower, N., Rasmussen, L., and Schneider, D. A. (2019). Discovery of novel inhibitors of ribosome biogenesis by innovative high throughput screening strategies. Biochem. J. 476, 2209–2219. doi: 10.1042/BCJ20190207
Sheng, Z., Han, W., Huang, B., and Shen, G. (2020). Screening and identification of potential prognostic biomarkers in metastatic skin cutaneous melanoma by bioinformatics analysis. J. Cell. Mol. Med. 24, 11613–11618. doi: 10.1111/jcmm.15822
Smoot, M. E., Ono, K., Ruscheinski, J., Wang, P.-L., and Ideker, T. (2011). Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics 27, 431–432. doi: 10.1093/bioinformatics/btq675
Somura, H., Iizuka, N., Tamesa, T., Sakamoto, K., Hamaguchi, T., Tsunedomi, R., et al. (2008). A three-gene predictor for early intrahepatic recurrence of hepatocellular carcinoma after curative hepatectomy. Oncol. Rep. 19, 489–495. doi: 10.3892/or.19.2.489
Soura, E., Eliades, P. J., Shannon, K., Stratigos, A. J., and Tsao, H. (2016). Hereditary melanoma: Update on syndromes and management: genetics of familial atypical multiple mole melanoma syndrome. J. Am. Acad. Dermatol. 74, 395–407. doi: 10.1016/j.jaad.2015.08.038
Subramanian, A., Narayan, R., Corsello, S. M., Peck, D. D., Natoli, T. E., Lu, X., et al. (2017). A next generation connectivity map: L1000 platform and the first 1,000,000 profiles. Cell 171, 1437–1452. doi: 10.1016/j.cell.2017.10.049
Sun, L., Guan, Z., Wei, S., Tan, R., Li, P., and Yan, L. (2019). Identification of long non-coding and messenger RNAs differentially expressed between primary and metastatic melanoma. Front. Genet. 10:292. doi: 10.3389/fgene.2019.00292
Sun, L., Sun, P., Zhou, Q.-Y., Gao, X., and Han, Q. (2016). Long noncoding RNA malat1 promotes uveal melanoma cell growth and invasion by silencing of mir-140. Am. J. Transl. Res. 8:3939.
Szklarczyk, D., Morris, J. H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., et al. (2016). The string database in 2017: quality-controlled protein–protein association networks, made broadly accessible. Nucleic Acids Res. 45, D362–D368. doi: 10.1093/nar/gkw937
Szumera-Ciekiewicz, A., Bosisio, F., Teterycz, P., Antoranz, A., and Massi, D. (2020). Sox10 is as specific as s100 protein in detecting metastases of melanoma in lymph nodes and is recommended for sentinel lymph node assessment. Eur. J. Cancer 137, 175–182. doi: 10.1016/j.ejca.2020.06.037
Tan, Y., Yu, R., and Pezzuto, J. M. (2003). Betulinic acid-induced programmed cell death in human melanoma cells involves mitogen-activated protein kinase activation. Clin. Cancer Res. 9, 2866–2875.
Tang, L., Zhang, W., Su, B., and Yu, B. (2013). Long noncoding RNA hotair is associated with motility, invasion, and metastatic potential of metastatic melanoma. BioMed. Res. Int. 2013:251098. doi: 10.1155/2013/251098
Tian, Y., Zhang, X., Hao, Y., Fang, Z., and He, Y. (2014). Potential roles of abnormally expressed long noncoding RNA uca1 and malat-1 in metastasis of melanoma. Melanoma Res. 24, 335–341. doi: 10.1097/CMR.0000000000000080
Tomczak, K., Czerwińska, P., and Wiznerowicz, M. (2015). The cancer genome atlas (tcga): an immeasurable source of knowledge. Contemp. Oncol. 19:A68. doi: 10.5114/wo.2014.47136
Vargas, F. A., Furness, A. J., Litchfield, K., Joshi, K., Rosenthal, R., Ghorani, E., et al. (2018). Fc effector function contributes to the activity of human anti-ctla-4 antibodies. Cancer Cell 33, 649–663. doi: 10.1016/j.ccell.2018.02.010
Wang, L., Wei, C., Xu, Y., Deng, X., Wang, Q., Ying, J., et al. (2019). Prognostic genes of melanoma identified by weighted gene coexpression network analysis and drug repositioning using a networkbased method. Oncol. Lett. 18, 6066–6078. doi: 10.3892/ol.2019.10961
Wang, P., Hu, L., Fu, G., Lu, J., Zheng, Y., Li, Y., et al. (2020). LncRNA malat1 promotes the proliferation, migration, and invasion of melanoma cells by downregulating mir-23a. Cancer Manage. Res. 12:6553. doi: 10.2147/CMAR.S249348
Wang, Q., Wang, X., Liang, Q., Wang, S., Xiwen, L., Pan, F., et al. (2018). Distinct prognostic value of mrna expression of guanylate-binding protein genes in skin cutaneous melanoma. Oncol. Lett. 15, 7914–7922. doi: 10.3892/ol.2018.8306
Xia, Y., Zhou, Y., Han, H., Li, P., Wei, W., and Lin, N. (2019). lncRNA neat1 facilitates melanoma cell proliferation, migration, and invasion via regulating mir-495-3p and e2f3. J. Cell. Physiol. 234, 19592–19601. doi: 10.1002/jcp.28559
Xu, B.-F., Liu, R., Huang, C.-X., He, B.-S., Li, G.-Y., Sun, H.-S., et al. (2020). Identification of key genes in ruptured atherosclerotic plaques by weighted gene correlation network analysis. Sci. Rep. 10, 1–10. doi: 10.1038/s41598-020-67114-2
Yang, X. R., Pfeiffer, R. M., Wheeler, W., Yeager, M., Chanock, S., Tucker, M. A., et al. (2009). Identification of modifier genes for cutaneous malignant melanoma in melanoma-prone families with and without cdkn2a mutations. Int. J. Cancer 125, 2912–2917. doi: 10.1002/ijc.24622
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterprofiler: an R package for comparing biological themes among gene clusters. Omics 16, 284–287. doi: 10.1089/omi.2011.0118
Zeng, Y., Guo, Z., Hu, Z., Liu, M., Chen, Y., Chen, S., et al. (2020). Fgd1 exhibits oncogenic properties in hepatocellular carcinoma through regulating cell morphology, autophagy and mitochondrial function. Biomed. Pharmacother. 125:110029. doi: 10.1016/j.biopha.2020.110029
Keywords: skin cutaneous melanoma, long non-coding RNA, differential expression, gene regulation, messenger RNA (mRNA), survival analysis, prognostic biomarker
Citation: Zhang J, Liu H, Zhang W, Li Y, Fan Z, Jiang H and Luo J (2020) Identification of lncRNA-mRNA Regulatory Module to Explore the Pathogenesis and Prognosis of Melanoma. Front. Cell Dev. Biol. 8:615671. doi: 10.3389/fcell.2020.615671
Received: 09 October 2020; Accepted: 24 November 2020;
Published: 17 December 2020.
Edited by:
Liang Cheng, Harbin Medical University, ChinaReviewed by:
Jiajie Peng, Northwestern Polytechnical University, ChinaJingpu Zhang, Henan University of Urban Construction, China
Copyright © 2020 Zhang, Liu, Zhang, Li, Fan, Jiang and Luo. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Hua Jiang, lygcancer@126.com; Judong Luo, judongluo@163.com