Secondary data mining of GEO database for long non-coding RNA and Competing endogenous RNA network in keloid-prone individuals

Introduction

Both hypertrophic scars and keloids are pathological scars. Due to abnormal proliferation of dermal fibrous tissue, excessive deposition of collagen during the wound healing process often affects the appearance of patients and even affects local function [1]. Although many studies have addressed the mechanism underlying keloid occurrence and development and the treatment of keloids, the therapeutic effect of existing treatments on keloids is still not ideal [2]. Gene chips and sequencing are powerful tools currently used for studying gene expression profiles; however, due to limitations of the manufacturing process and algorithms for analysis of gene chips during the early stage of this technology, the required data could not be acquired completely [3]. As bioinformatic analysis methods have continually developed and the understanding of epigenetics has increased, many scholars have mined new transcriptome information from old microarray or sequencing data by improving the original algorithms and upgrading the databases to offer new ideas for clinical and scientific research [3–7].

In this study, sequencing and microarray data related to keloid formation were downloaded from the open-source National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database, and secondary data analysis was conducted by combining the R language with network database bioinformatic analysis technology. This analysis aimed to identify the coding RNAs (messenger RNAs, mRNAs) and non-coding RNAs (ncRNAs), including microRNAs (miRNAs) and long non-coding RNAs (lncRNAs), expressed during keloid formation and to predict their expression and relationships. Through this research, we intended to verify the possibility of secondary bioinformatic analysis of sequencing and microarray data to extract lncRNAs and construct a competing endogenous RNA (ceRNA) network and simultaneously mine these data for information related to regulation of the pathophysiological process of keloid formation. This study offers new research methods for more comprehensive exploration of keloid pathogenesis and treatment.

Results

Discussion

Keloid scarring is a clinical condition that is difficult to cure. Although keloids are benign tumours, they often extend past the extent of the damage, invading the surrounding normal skin, and seriously affecting the appearance of the skin after healing [4]. Keloids are most common in people between the ages of 10 and 30 years, are especially common in individuals of African, Hispanic, and Asian ethnicity, and frequently occur on the chest, earlobes, shoulders, and back [26]. Although different treatment methods, such as radiotherapy, hormone therapy, and surgical resection are available, keloids still have a relatively high recurrence rate [8]. Many scholars believe that fibroblasts are a main participant in the occurrence and development of keloids. After skin injury, a complicated regulatory signalling network is activated to regulate the proliferation, migration, and secretion of fibroblasts. Therefore, the biological behaviour of fibroblasts regarding the mechanism of keloid formation is a trending research topic [9]. Gene chips and sequencing are important methods for studying gene transcription profiles and are widely used in regenerative medicine [30], disease [27, 10–11], and cancer [7] research. However, the focus of transcriptome research has recently shifted gradually from protein-coding genes to epigenetics involving ncRNAs. In 2018, Onoufriadis et al. [26] compared patients with keloids with patients without keloids. They sampled wounds on the buttocks of all participants and sampled the untreated wounds 6 weeks later. Two specimens from each participant were obtained for transcriptome analysis, and the results revealed differences in mRNA and miRNA expression between the 2 groups and identified the miRNA-mRNA regulatory relationships. In the current study, the open-source GEO database to conduct secondary mining of ncRNAs and mRNAs through bioinformatics algorithms, aiming to identify suitable ncRNA molecules for further exploration of the epigenetic regulation mechanism of keloids.

In this study, the genes of keloid-prone patients were separated from those of control individuals and were analysed separately according to the wound healing process. Considering that the wound healing process is common to all populations, we excluded the intersecting set of genes between the keloid-prone patients and the control individuals, retaining only the genes differentially expressed in keloid-prone patients for further study. Because differential gene screening is generally based on counts, errors often occur. To reduce the false positive and false negative rates, the most recent version of DESeq2 software was used herein to obtain the required DEmRNAs and DElncRNAs in order to ensure the accuracy of the data analysis [12]. For typical prediction of miRNA target genes, more than 2 miRNA target gene prediction software programs are selected for prediction, and the intersection of the prediction results is then obtained to identify the most reliable miRNA target gene prediction results [26]. In this study, we used the StarBase database for analysis. The miRNA data in this database were verified by Ago2-CLIP-Seq [31] with a relatively high degree of reliability.

GO is an international standardized gene functional classification system that provides a set of dynamically updated controlled terms to comprehensively describe the attributes of genes and gene products in organisms. GO has 3 ontologies, which individually describe the MF, CC, and BP of genes. The basic unit of GO is the term, and each term corresponds to an attribute. Each gene is associated with 1 or more GO terms (aka GO functions) [13]. KEGG pathway analysis is a gene annotation database-based method for detecting pathways significantly enriched with differentially expressed genes. Therefore, a comprehensive database and complete pathway annotations are the keys to pathway analysis. Pathway analysis results reveal direct interactions between genes in the actual signalling pathways. Compared with GO analysis, pathway analysis allows researchers to study the biological effects of differentially expressed genes; moreover, it is complementary to GO analysis [30].

This study showed that during wound healing, DEmRNAs expressed only in keloid-prone individuals were mainly enriched in various ion channels in GO MF categories. KEGG pathway analysis further showed that these DEmRNAs were enriched mainly in the calcium signalling pathway. Liang et al. reported a similar finding in a microarray study focusing on DElncRNAs and DEmRNAs between keloid tissue and normal skin tissue and found that the lncRNA CACNA1G-AS, an antisense RNA to CACNA1G, which encodes a calcium channel subtype, was upregulated in keloid tissue in independent pairs of samples [14]. Calcium channels are associated with reactions to mechanical forces in different cell types, and elevated intracellular calcium levels were observed in fibroblasts and keratinocytes under application of hydraulic pressure [15]. Promotion of calcium channel expression could positively affect the migration of keloid fibroblast cells [16]. L-type calcium channel blockers, such as verapamil, have already been used in clinical treatment for keloids and have been proven to be safe and effective [17]. These results collectively suggest that dysregulation of a mechanical signal-related network, with mediation of calcium channel activation, could be pivotal in keloid development. KEGG pathway analysis also showed that the glutamatergic synapse pathway might also play an important role. Evidence has shown immunoreactivity of L-glutamate and its ionotropic receptor NMDAR in the human epidermis [18–20]. These effects may cooperatively regulate the calcium levels in related tissues and cells, as previously reported [21]. Moreover, our findings showed that neuroactive ligand-receptor interactions might also contribute. A 2018 bioinformatic analysis of the molecular mechanism of pathological scarring by Zhang et al. noted that genes expressed in scar tissues were enriched in this channel [22]. The neuroactive ligand-receptor interaction pathway includes hundreds of genes involving many types of neuroreceptors, suggesting a broader, yet-to-be-discovered, complex neuroactive ligand-receptor mediated network; therefore, further investigations are needed. Moreover, GO analysis showed that the genes were enriched not only in the BP terms muscle contraction, presynapse assembly, and regulation of membrane potential but also CC terms related to different structural parts of synapses. Last, by KEGG analysis, we also found DEmRNAs enriched in the nicotine addiction pathway, which has not been reported in previous studies and awaits experimental verification.

ncRNAs are RNAs that are not directly translated into proteins and were once thought to be non-functional components in gene expression and transcription [10]. As epigenetics and other research methods for genes and proteomic approaches have recently been developed, ncRNAs have been found to be associated with gene expression, not only regulating gene transcription and post-transcriptional modifications and translation but also forming ceRNA regulatory networks to participate in mutual reciprocal regulation, which in turn affects the biological functions of cells, tissues, and organisms [20]. According to their length, ncRNAs can be divided into small non-coding RNAs (sncRNAs), lncRNAs, and circular RNAs (circRNAs). Wang et al. studied the expression and effects of lncRNA-H19 in keloid fibroblasts and found that H19 regulated the viability and apoptosis of fibroblasts through the miR-29a/COL1A1 axis [19]. However, the functions of most lncRNAs are still unclear. Therefore, we used the existing gene data for keloid fibroblasts to identify the ncRNA molecules differentially expressed during keloid formation and to more comprehensively explore the epigenetic regulatory mechanism in keloids.

We found that a total of 6 lncRNAs were involved in the network: ADIRF-AS1, DLEU2, AP000317.2, AC006333.2, AL137127.1, and LINC01725. Among these DElncRNAs, DLEU2 and AP000317.2 were upregulated, and the others were downregulated. The ceRNA sub-network in which these lncRNAs are involved is shown in Figure 4A–4C. The experimental design for acquisition of the original GEO data involved invasive surgery on uninjured skin of patients with keloid scars. Therefore, no matching data can be found in public databases. To verify the expression of the DElncRNAs, we used pathological specimens of keloid scars as a substitute to compare with normal skin. The results showed that only the expression of DLEU2 and ADIRF-AS1 in keloid fibroblasts showed significant difference and that the expression of other genes was not changed significantly. We cannot exclude the possibility that our results were affected by differences in the sequencing models or qPCR analysis. Although these lncRNAs have not been reported in previous keloid-related literature, molecules such as DLEU2 and ADIRF-AS1 have been detected in other tissues and diseases [23, 24]. DLEU2, a gene located on chromosome 13q14, has frequently been observed to be deleted or epigenetically suppressed in leukaemia and is considered a tumour suppressor [24, 25]. Wu et al. found that the pattern of high lncRNA DLEU2 expression or low miR-30a-5p expression was an unfavourable prognostic factor for survival and tumour recurrence in patients with non-small cell lung cancer and showed that the lncRNA DLEU2 accelerates the development of cancer by sponging miR-30a-5p [24]. In this study, we predicted that DLEU2 acted on target genes, including development-related genes (such as LIMS3), apoptosis-related genes (such as SKIL), differentiation-related genes (such as PELO), and protein transporter-related genes (such as SLC35B4), through miR-30a-5p or miR-30b-5p, suggesting that DLEU2-associated ceRNAs may play an important regulatory role in keloid formation and wound healing. This possibility awaits verification by more basic and clinical experimental data.

Through data mining of the open-source GEO database, we acquired new regulatory information from previous sequencing and microarray data. However, our study has limitations. First, our study was based on secondary analysis of online sequencing data with limited amounts of original clinical samples. More clinical samples are needed to provide supportive evidence. The differential expression of lncRNAs and their targets in our ceRNA network also requires confirmation via qRT-PCR. Second, our hypotheses were inspired by similar findings in tumour studies, because keloids are thought—but have not been confirmed—to be tumour-like. Therefore, further well-designed experiments are needed to verify these hypotheses.

Conclusion

This study used bioinformatics methods to perform in-depth data mining of the ncRNA data in GEO microarray datasets. We found that the calcium signalling pathway, glutamatergic synapse pathway, neuroactive ligand-receptor interaction pathway and nicotine addiction pathway may have potential regulatory effects that influence keloid development. We also identified lncRNAs and a ceRNA regulatory network that exist only in keloid-prone patients during the wound healing process, possibly providing novel intervention targets and therapeutic options for patients with keloid scars.

Materials and Methods