{'fraction_non_alphanumeric': 0.031622935669224896, 'fraction_numerical': 0.021778438542210794, 'mean_word_length': 5.250125565042692, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Integrated analysis of circRNA-miRNA-mRNA regulatory network identifies potential diagnostic biomarkers in diabetic foot ulcer

Diabetic foot ulcer (DFU) is a common and serious complication of diabetes mellitus, which influences patients’ quality of life. Recently, circRNA regulated the mRNA levels by functioning as miRNA sponge in various disease, including diabetes mellitus. Nevertheless, the circRNA-miRNA-mRNA regulatory network involved in DFU remains obscure. The aim of this study is to construct a competing endogenous RNA (ceRNA) network and screen biological indicators as diagnostic factors in DFU. All the diﬀerentially expressed circRNAs, miRNAs and mRNAs were derived from Gene Expression Omnibus database. Furthermore, circRNAs identified by cytoHubba analysis and miRNAs obtained by human miRNA-disease database were used to construct DFU-specific ceRNA network with intersection of mRNAs. Functional enrichment analysis displayed the function and pathway of dysregulated mRNAs. Hub genes with high diagnostic value were screened by ClusterONE, GO semantic similarity and receiver operating characteristic (ROC) curve. Here, the ceRNA network consisted of 8 circRNAs, 11 miRNAs and 91 mRNAs. Functional enrichment analysis demonstrated diabetic complications-related pathway including TGF-beta, FoxO and Wnt signaling pathway. GO semantic similarity and ROC curve analysis showed 6 hub genes with high diagnostic value (the area under the ROC curve ≥ 0.8) in patients with DFU, including BCL2, CCND1, IRAK4, SMAD4, SP1 and SUFU, which were identified as potential target genes for DFU diagnosis. In conclusion, the present study looked at a circRNA-miRNA-mRNA regulatory network with DFU and screened the potential function of mRNA, then identified novel diagnostic biomarkers and therapeutic targets for patients with DFU.


Introduction
Diabetic foot ulcer (DFU) is one of the most common diseases of lower extremity complications in diabetes mellitus, resulting in a heavy socioeconomic burden in families and societies. The patients with DFU usually increase costs and risk because of required long-term medical treatments and poor availability of tools for effective debridement [1]. The average medical cost of patients with diabetic foot infection in China was about $2685 per person in 2012, which is significantly higher than that in 2004 [2]. In addition, the enhancement of wound healing represents a major clinical challenge during the diabetic foot disease [3]. It has been showed that a large number of physiologic factors involved in the wound healing, including circular RNAs (cir-cRNAs) [4,5].
CircRNA is an emerging class of noncoding RNA molecule of lacking 5′ caps and 3′ tails, which is consisted of a circular configuration through the 5′ to 3' phosphodiester bond [6]. Early studies indicated that the circRNA may be a byproduct resulted from the process of missplicing [7]. With the development of next generation sequencing technology, accumulating evidence has emphasized the importance of the circRNA to discuss the roles in the competition of splicing machinery between linear RNA and circRNA [8]. Recent reports demonstrated that the hsa_circ_0067301 that competes for binding to miR-141-5p strongly suppressed miRNA activity and results in increased levels of Notch-1 during human epithelial-mesenchymal transition [9]. An increasing number of studies showed that circRNAs functioned as miRNA sponge in the competing endogenous RNA (ceRNA) network and played an important role in the disease process [10]. For instance, circRNA_010567 in ceRNA network has been shown to promote progression of myocardium and cardiac fibroblasts (CFs) cells by sponging https://doi.org/10.1016/j.ncrna.2020.07.001 Received 18 June 2020; Received in revised form 28 July 2020; Accepted 28 July 2020 miR-141 in diabetic mice [11]. Recent researches investigated the role of circRNA-associated ceRNA network in diabetes complications including diabetic nephropathy [12], diabetic cardiomyopathy [13] and DFU [4]. However, only one study that investigated the functional effects of hsa_circ_0084443 on DFU has been reported up to now [4]. The circRNA-mediated ceRNA regulatory network and biological function correlated to DFU remains obscure.
To investigate the effects of circRNAs on DFU, we collected the expression files of circRNA, miRNA and mRNA from Gene Expression Omnibus (GEO) datasets. Differentially expressed circRNAs, miRNAs and mRNAs were identified by "limma" package in R software. Next, the regulatory network of circRNA-miRNA-mRNA was constructed using bioinformatics approaches in accordance with the ceRNA theory. The main functional pathways of DFU were evaluated using gene ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. The subnetwork of circRNA-miRNA-hub genes was established using the hub genes that were identified by protein-protein interaction (PPI) network and GO semantic similarity. The receiver operating characteristic (ROC) curve was used to identify the potential diagnostic biomarkers in the occurrence and development of DFU.

Data collection
The microarray data were acquired from the Gene Expression Omnibus (GEO) database. The circRNA expression data was obtained from GSE114248 (5 DFU tissues and 5 non-DFU tissues). The miRNA data was derived from GSE84971 (3 DFU tissues and 3 non-DFU tissues). The mRNA data was obtained from GSE80178 (6 DFU tissues and 3 non-DFU tissues). Diabetes mellitus-related miRNAs were collected from the HMDD v3.0 database [14].

Differential expression analysis of circRNAs, miRNAs, and mRNAs
The analysis of differentially expressed of cirRNAs, miRNAs and mRNAs were performed using "limma" package in R software. The differentially expressed circRNAs (DECs), differentially expressed miRNAs (DEMs) and differentially expressed genes (DEGs) were screened with P-value < 0.05 and |logFC| > 1 as the threshold.

Construction of the circRNA-miRNA-mRNA regulatory network
The Information about circRNAs were available on the CircInteractome (http://circinteractome.nia.nih.gov) and circBase database (http://www.circbase.org/). In order to identify the key circRNAs participating in the ceRNA network, the closeness centrality, degree and betweenness centrality of cicrRNAs were calculated using cytoHubba plugin based on cytoscape software [15]. The target miRNAs of circRNAs were obtained using the Encyclopedia of RNA Interactomes (ENCORI, http://starbase.sysu.edu.cn/index.php). Furthermore, the intersections between the DEMs and the miRNAs that were collected using HMDD database were considered as candidate miRNAs to construct the ceRNA network. Then, the TargetScan [16], miRDB [17] and miRWalk [18] were used to predict the relationships among DEMs and their target DEGs. The intersecting miRNAs negatively regulated by circRNAs and mRNAs were used to construct the ceRNA network according to the ceRNA theory. The visualization of regulatory network was performed by cytoscape 3.7.2.

Functional enrichment analysis and construction PPI network
The miRNA KEGG pathway analysis were performed using DIANA-mirPath 3.0 [19]. Moreover, KEGG functional analysis was conducted using DEGs in the ceRNA network with "clusterProfiler" package in R. Gene Ontology (GO) analysis was conducted using KOBAS 2.0. The false discovery rate (FDR) of less 0.05 was set as the cut-off criterion. GO semantic similarity analysis was performed to compute semantic similarity among GO terms by "GOSemSim" package in R software.
The protein-protein interaction (PPI) network was conducted using the Search Tool for the Tetrieval of Interacting Genes (STRING) database (http://string-db.org). The ClusterONE algorithm was used to screen modules of hub genes from PPI network. The visualization was performed by cytoscape software. The ROC curve were performed by calculating the sensitivity and specificity of target expression level to screen the potential diagnostic target genes in patients with DFU.

Convergence of gene expression signatures across different studies from GEO database
To investigate differentially expressed genes in each data sets, the "limma" package in R software were used for differentially expressed genes with P-value < 0.05 and |LogFC| > 1 as the threshold. The GSE114248 data set contained 945 DECs including 568 upregulated circRNAs and 645 downregulated circRNAs (Fig. 1A). The integrated analysis of GSE84971 dataset identified 30 DEMs, which included 29 upregulated miRNAs and 1 downregulated miRNAs (Fig. 1B). Furthermore, the DEGs were identified using the GSE80178 dataset. A total of 1715 DEGs were obtained, 430 of which were upregulated and 1285 of which were downregulated (Fig. 1C).

Construction of the ceRNA network
The role of circRNAs as a miRNA sponge is the main mechanism of circRNA function in disease. Therefore, we explore the role of circRNA through establishing a ceRNA network in DFU. The circRNA-miRNA-mRNA ceRNA network were constructed using a combination of circRNA-miRNA pairs and miRNA-mRNA pairs. In order to screen the novel circRNAs of these DECs in DFU, the betweenness centrality, degree and closeness centrality of cicrRNAs were performed by using the cytoHubba plugin (Table 1). We predicted the target miRNAs of cir-cRNAs using the ENCORI database, which has the advance that the results of the database could be filtered for experimentally-validated target miRNAs.
In addition, previous studies were performed using the open webresource HMDD v3.0 to investigate the claim that miRNA are truly specific for DFU. The HMDD database collected miRNAs that are associated with the development of a broad range of diseases [20]. Among miRNAs included in the database, we selected the miRNAs associated with the spectrum of diabetes according to the categories: "diabetes mellitus", "diabetes foot". Finally, 19 diabetes-associated DEMs were selected to use for further analysis with intersections of 30 DEMs. For miRNA-mRNA pairs, the TargetScan, miRDB and miRWalk programs were used to predict the target mRNAs with intersections of 1715 DEGs. Taken together, we established the circRNA--miRNA-mRNA ceRNA network, including 8 DECs, 11 DEMs, 91 DEGs (Supplementary Table 1). The ceRNA network was plotted using cytoscape 3.7.2 software (Fig. 2).

Functional enrichment analysis
In order to explore the altered biological processes of ceRNA network, the basic characteristics for the 8 DECs were performed using circInteractome and circBase database, as described in Table 2. Three of DECs (hsa_circ_0089761, hsa_circ_0089762 and hsa_circ_0001222) were not found in these database. The role of circRNAs as a miRNA sponge is the molecular mechanism of circRNA function in diseases. Therefore, the combinatorial effect of miRNAs in pathways was investigated using DIANA-miRPath v3.0. The results demonstrated two novel pathway including TGF-beta and FoxO signaling pathway (Fig. 3A), which was verified in KEGG pathway analysis of DEGs.
In addition, the functional enrichment analysis result of identified DEGs were performed by using KOBAS online database. The KEGG pathway analysis showed that the dysregulated mRNAs were enriched in diabetic complication-related pathways, including Wnt, TGF-beta, PI3K-Akt, p53, NF-kappa B and FoxO signaling pathway (Fig. 3B, Table 3). Furthermore, GO enrichment analyses were performed using KOBAS with P-value < 0.01 as the cutoff criteria. The top 20 terms were associated with multiple biological progress including cellular metabolic process, cell proliferation and protein metabolism (Fig. 3C, Table 4).

Construction of the PPI network and identification of the hub genes
A total of 80 protein-protein interactions among 91 DEGs in ceRNA network were screened by applying STRING database (Fig. 4A). With the help of ClusterONE, we identified 10 PPI modules with P-value < 0.05 as the cutoff thresholds (Fig. 4B). Subsequently, in order to investigate the function of modules in PPI network, GO enrichment analysis results of module 8, 10 and 11 indicated enrichment in biosynthetic process. The results of module 6 and 9 suggested that these dysregulated genes were enriched in canonical Wnt signaling pathway and ephrin receptor signaling pathway, respectively. The enriched most significant function for the differentially expressed genes associated with the 10 modules were listed in Table 5.
The nodes with high degree are considered as important genes in network. In our study, the top genes were selected in the PPI network with degree > 4 as the cutoff criteria. In addition, in order to investigate the function analysis of target genes. 16 dysregulated genes were enriched in diabetic complication-related signaling pathways. The intersection between the top genes in the network and genes that were associated with diabetic complication-related signaling pathways were AXIN2, BCL2, CCND1, CXCL12, EPHA4, FGF7, IRAK4, MDM4, MEIS1, PTEN, SMAD4, SMAD6, SP1, SSH2, SUFU and WIF1, which were considered as the hub genes of DFU (Fig. 4C). Furthermore, a circRNA-miRNA-hub gene subnetwork were construction according to the hub gene, as described in Fig. 4D.

ROC and GO semantic similarity analysis
To explore specific genes and potential biomarker that might be capable of discriminating the patients with DFU from controls, the ROC analysis were performed for 16 hub genes in ceRNA subnetwork using the GSE80178 dataset as training set and GSE132187 dataset as test set. The best performance demonstrated to be the discriminating biomarkers with area under ROC curve (AUC) > 0.8 as the threshold, including BCL2, CCND1, IRAK4, SMAD4, SP1 and SUFU (Fig. 5A). SMAD4, SUFU and IRAK4 displayed the highest discriminatory power as a single marker, exhibiting a sensitivity of 100%, and specificity of 100%. Furthermore, in order to measure the strength of the relationship between each protein and its partners by considering functions and location of proteins, we ranked these 16 genes in the subnetwork with GO semantic similarity. The result showed that SMAD4 was considered as the most significantly strength the relationship between each protein and its partners (Fig. 5B). The intersection between ROC analysis and GO semantic similarity were SMAD4, which was regarded as the most potential and robust biomarker candidates for DFU. The KEGG analysis showed that SMAD4 was enriched in Wnt, FoxO, TGF-beta and AGE-RAGE signaling pathway (Fig. 5C). Finally, these findings provided evidence that SMAD4 that was identified by DFU-associated ceRNA might play an important role as potential diagnostic biomarker in patients with DFU.

Discussion
CircRNAs are ignored in classical transcriptome researches for a long time because of the lack of 5′ caps and 3′polyadenylated tails [21]. With the development of the technologies of next-generation sequencing (NGS) and the methods of biochemical and computational biology, a large number of circRNAs were studied in various tissues and cells [10]. An increasing number of studies have revealed the important role of circRNA-associated ceRNA network in screening diagnostic and prognostic biomarkers in a myriad of human diseases, such as cancer and diabetes mellitus [10,22]. In hepatocellular carcinoma, circMTO1 as a prognosis predictor suppressed hepatocellular carcinoma progression by acting as the sponge of oncogenic miR-9 to promote p21 expression [23]. Furthermore, study have been demonstrated the role of circHIPK3 in retinal vascular dysfunction induced by diabetes mellitus. The circHIPK3 expression was significantly upregulated in diabetic retinas and retinal endothelial cells following stressors related to diabetes mellitus [24]. Although several circRNAs have been screened as participating in the pathogenesis of complications in diabetes mellitus, the circRNA-associated ceRNA network and biological function corrected to the pathogenesis of extremity complications in diabetes mellitus including diabetic foot remain obscure. It is still necessary to conduct a circRNA-miRNA-mRNA regulatory network and advance our understanding of the molecular mechanism of DFU. In the present study, we constructed a ceRNA regulatory network in DFU consisting of 8 differentially expressed circRNAs (DECs), 11 differentially expressed miRNAs (DEMs) and 91 differentially expressed genes (DEGs). In order to explore the interaction mechanism of ceRNA in patients with DFU, the miRNAs of circRNAs downstream were analyzed. Our results showed 11 DEMs were sponged by circRNAs involving in the ceRNA network. Among 11 DEMs, hsa-miR-146a-5p was identified as one of the highly up-regulated miRNA in the regulatory ceRNA network. In patients with type 1 diabetes mellitus, hsa-miR-146a-5p at 3  months correlated with residual beta cell function after diagnosis [25]. Furthermore, miR-146a as a potential molecular target contributed to inhibiting inflammation and apoptosis in the diabetic retina [26], suggesting miR-146a-5p might play an important role in pathogenesis and progression of diabetic complications. These findings were in consistence with our results that miR-146a-5p as a biological indicator was significantly increased in patients with DFU. Previous studies showed that circRNA was considered for its active role as miRNA sponge to indirectly regulate the function of mRNA in ceRNA network [27]. Studies have showed that dysregulated mRNAs involved in diabetic complications [28]. In diabetic foot ulcer, pro-inflammatory S100A8 and IL-8 proteins contributed to impaired wound healing and caused persistent inflammation [29]. In our analysis, 91 differentially expressed mRNAs were involved in the ceRNA network. Furthermore, we performed numerous bioinformatics analysis to identify potential novel genes based functional enrichment, PPI network, GO semantic similarity and ROC analysis. The results indicated that BCL2, CCND1, IRAK4, SUFU, SP1 and SMAD4 were found to be the key differentially expressed genes. These dysregulated mRNAs might act as the essential molecules that mediated the progression of DFU.
BCL2 (BCL2 apoptosis regulator) encodes an integral outer mitochondrial membrane protein that blocks the apoptotic death of cells [30]. Recent study showed that Bcl-2 modified adipose-derived stem cells effectively contributed to the wound healing process in diabetic mice [31]. Similarity, our results showed the effects of BCL2 on DFU. Studies have showed that CCND1 (Cyclin D1) were diagnosed as the top hub genes in diabetes mellitus to insulin resistance [32]. IRAK4 (Interleukin 1 Receptor Associated Kinase 4) are serine-threonine kinases involving in essential regulators of IL-1R and Toll-like receptor-mediated signaling pathways [33]. It has been reported that inhibition of IRAK-4 activity contributed to attenuating the progression of diabetic nephropathy under diabetic conditions [34]. SUFU (Suppressor of Fused) are a negative regulator of the hedgehog signaling pathway [35]. Studies have showed that SUFU could inhibit the activation of the insulin promoter mediated by Glis3, which has been linked to both type I and type II diabetes [36]. SP1 (Specificity protein-1) is a zinc finger gene that binds to GC-rich motifs of many promoters and acts in many cellular processes, including cell differentiation, cell growth, apoptosis KEGG pathway enrichment analysis were performed using 91 differential target genes that were negatively associated with 11 miRNAs in ceRNA network. C. Gene Ontology (GO) analysis of differential expressed mRNAs in ceRNA network. and immune responses [37]. Our result identified SP1 as a diagnostic gene in patients with DFU, which was in consistence with a recent report that Sp1 acted in diabetic wound healing by mediating matrix metalloproteinase-9 expression [38]. SMAD4 (SMAD family member 4) is a member of the Smad family of signal transduction proteins and is as an essential regulator of the transforming growth factor β (TGF-β) signaling pathway [39,40]. Previous researches showed that Smad4 was upregulated in human and mouse podocytes during diabetic nephropathy [41]. Depletion of Smad4 in podocytes protected mice from glomerulosclerosis in type 2 diabetes models [41]. Interestingly, GO semantic similarity analysis showed SMAD4 was the highest score among all of the hub genes in our results. ROC curve analysis showed that SMAD4 was a specific diagnostic biomarker in patients with DFU. These findings provided evidences that the candidate key genes played an important role in diagnosis and prognosis of diabetic complications. Our results identified these candidate key genes as diagnostic target indicators in patients with DFU. Although the molecular mechanism between some candidate key genes (CCND1, IRAK4, SUFU and SMAD4) and DFU have not been reported, further researches are needed to explore the molecular mechanism of novel target genes and implement combination therapy for DFU.
In order to investigate the potential causative roles of dysregulated genes in ceRNA network, the dysregulated biological pathways were identified by enrichment analysis. Our results showed that these dysregulated mRNA were enriched multiple signaling pathways, such as TGF-beta, FoxO and Wnt signaling pathway. Studies have found that the role of TGF-beta signaling pathway has been examined for the treatment of DFU [42]. TGF-β that was released from vicenin-2 film enhanced cell proliferation and wound contraction through the TGF-β signaling pathways [43]. In consistence with these findings, both functional enrichment analysis of dysregulated miRNAs and abnormal mRNAs in ceRNA network identified the novel TGF-β signaling pathways as an important regulator in patients with DFU in our results. Furthermore, Studies have shown that FOXO1 acted in the activation and regulation of key leukocytes needed to respond to bacterial challenge in oral mucosal tissues [44]. Although the effect of FoxO pathway on DFU has not been elucidated until now, the FoxO that are involved in diverse cellular function and clinical significance including immunity, wound healing and diabetic complications has been confirmed [45]. In addition, studies have showed that action of the Wnt/β-catenin pathway contributed to the inflammatory responses to diabetic ulcers and wound proliferation [46]. Another report demonstrated that inhibition of Wnt/β-catenin signaling was associated with impaired wound healing in patients with DFU [47]. Further study investigating strategies for molecular mechanism of wound healing in the DFU are needed. In our study, the comprehensive ceRNA regulatory network analysis of DFU was performed by using multiple bioinformatics analysis. However, there are some limitations to be acknowledged in our study. The biological experiments are still required to confirm the functions of the genes. Further studies warranted on the roles of hsa_-circ_0089761, hsa_circ_0089763, hsa-miR-146a-5p and SMAD4 as potential diagnostic biomarkers and therapeutic targets for DFU.
In summary, the present study showed a comprehensive network analysis and key communication of circRNA, miRNA and mRNA using bioinformatics approach. Our study provided novel insights into circRNA-associated ceRNA network and identified potential diagnostic biomarkers in patients with DFU.

Declaration of competing interest
The authors declared no conflicts of interest. Fig. 4. Construction of protein-protein interaction (PPI) network and ceRNA subnetworks. A. The PPI network was constructed for 80/91 differentially expressed mRNAs in ceRNA network. B. The modules were screened from the PPI network using "ClusterONE" algorithm containing 11 modules. Red dots represent upregulated mRNAs; blue dots represent downregulated mRNAs; grey dots represent mRNAs with no significant difference in expression. C. Interaction network of 16 hub genes. D. Sankey diagram for the circRNA-miRNA-hub gene subnetwork in diabetic foot ulcers, which was constructed using intersection between the top genes in the PPI network and potential genes that were associated with diabetic complication-related signaling pathways. Each rectangle indicates a gene. The size of the rectangle indicates the connection degree of each gene. Table 5 The enriched most significant function for the differentially expressed genes involved in 11 modules of PPI network.

Availability of data and materials
The sequencing data have been deposited in the Gene Expression Omnibus (GEO) database under the accession number GSE114248, GSE84971, GSE80178 and GSE132187.

CRediT authorship contribution statement
Shuping Liao designed experiments, collected data and helped analyze the data. Xiaolan Lin analyzed data and wrote the original draft. Weiqing Wu performed project administration and supervision. All authors read and approved the final manuscript.