Bioinformatics Prediction for Network-Based Integrative Multi-Omics Expression Data Analysis in Hirschsprung Disease

Hirschsprung’s disease (HSCR) is a rare developmental disorder in which enteric ganglia are missing along a portion of the intestine. HSCR has a complex inheritance, with RET as the major disease-causing gene. However, the pathogenesis of HSCR is still not completely understood. Therefore, we applied a computational approach based on multi-omics network characterization and clustering analysis for HSCR-related gene/miRNA identification and biomarker discovery. Protein–protein interaction (PPI) and miRNA–target interaction (MTI) networks were analyzed by DPClusO and BiClusO, respectively, and finally, the biomarker potential of miRNAs was computationally screened by miRNA-BD. In this study, a total of 55 significant gene–disease modules were identified, allowing us to propose 178 new HSCR candidate genes and two biological pathways. Moreover, we identified 12 key miRNAs with biomarker potential among 137 predicted HSCR-associated miRNAs. Functional analysis of new candidates showed that enrichment terms related to gene ontology (GO) and pathways were associated with HSCR. In conclusion, this approach has allowed us to decipher new clues of the etiopathogenesis of HSCR, although molecular experiments are further needed for clinical validations.


Introduction
Hirschsprung's disease (HSCR) is a neurocristopathy, caused by defective migration, proliferation, differentiation and/or survival of neural crest cells, leading to gut aganglionosis.It is a rare congenital disease, with an incidence of 1 in 3500-5000 live births [1].The disease is characterized by high heritability (>80%), significant sex bias (male: female ratio, 4:1) and high sibling and recurrence risk (3 to 17%) [2].The combination of all these distinctive characteristics are the hallmark of a typical multifactorial genetic disorder [3].Genetic studies on HSCR have identified more than 20 disease-associated genes, mainly RET (ret proto-oncogene) and EDNRB (endothelin receptor type B) [1].However, all these genes cumulatively explain <10% of cases, suggesting that additional genes may also be involved in the etiology of HSCR [3,4].
On the other hand, increasing evidence indicates that numerous microRNAs (miR-NAs) are altered in HSCR [4][5][6].MiRNAs are a class of small noncoding RNAs with function in posttranscriptional regulation of gene expression, thereby modulating diverse biological processes, such as neuronal cell differentiation, development, plasticity and survival [7,8].These regulators are involved in modulating the pathogenesis of HSCR by directly suppressing a variety of functional targets [4][5][6].MiRNAs have also been observed to be dysregulated in a variety of human pathologies, including cancer, neurodegenerative diseases and cardiovascular diseases [9][10][11].The ability of chosen miRNAs to target various diseaseassociated mRNAs makes them interesting candidates as targets of therapeutics (in the form of anti-miRNAs) or as therapeutics (in the form of miRNA mimics) [12].
A common challenge in rare disease research is to identify and understand the cellular and molecular mechanisms underlying the pathophysiology of the disease.The knowledge underlying normal gut development and motility by identification of dysregulated pathways in HSCR and their genetic causes may provide new possible treatment strategies.Over the past decades, many different approaches have been applied to human genetic studies-beginning with classical positional cloning, linkage analyses and genome-wide association studies (GWAS)-to the recent next-generation sequencing (NGS) studies, which have allowed us to identify novel genetic variants and genes linked to HSCR [13][14][15][16][17].
Studies implementing large-scale techniques (omics technologies) have increased our understanding of disease mechanisms and led to the discovery of new biological pathways, genetic loci underpinning disease progression, biomarkers, and targets for therapeutic development [18][19][20][21].Based on the hypothesis that molecular features within a system establish functional connections or are part of modules to carry out biological processes [21], advances in omics technologies allow for holistic studies into biological systems [21].These approaches further can help us to identify possible connections (e.g., genotype-phenotype correlations) and/or subnetworks (e.g., biological pathways) that are informative of an observed phenotype [22].Biological systems can be abstracted as a series of networks of interconnected molecular entities (genes, proteins, lncRNAs, miRNAs, etc.) [23], whereas a network is a representation of relationships (edges) between entities (nodes) [24].In general, network-based approaches obtain disease-related information from complex topological patterns [24].Networks commonly used in systems biology include gene regulatory networks, protein-protein interaction (PPI) networks, literature-curated networks and hybrid networks [25].
This study aims to identify new HSCR-associated modules, genes and miRNAs by applying a designed workflow based on a system biology approach, which integrates multiomics information and combines different biological networks to illustrate associations among multiple molecular features and to discover new biomarkers for this disease.

Materials and Methods
An overview of the different methods applied in this study is illustrated in Figure 1.In summary, a total of four major steps have been executed as follows: (1 • ) collecting and preprocessing expression data; (2 • ) identification of likely disease-associated modules and new disease candidate genes; (3 • ) prediction of disease-associated miRNAs; (4 • ) miRNA biomarker identification.Moreover, we included two additional working packages: (5 • ) functional enrichment analysis and (6 • ) disease relevance evaluation, which were conducted after steps 1, 2 and/or 3, as depicted in Figure 1.

Figure 1.
Flowchart containing the main steps followed on this study.The diagram is divided into four major sequential blocks, listed one through four.Additionally, the complementary working packages "5º Functional enrichment analysis" and "6º Disease relevance evaluation" are shown in yellow.Arrows show the direction of the workflow.Processes are represented by a rectangle, while the ovals correspond to the input/output data of each process.HSCR, Hirschsprung disease; Om-icsDI, Omics Discovery Index; GEO2R, interactive web tool of Gene Expression Omnibus (GEO); miRNA, microRNA; HDEGs, known HSCR DE genes; ODEGs, other DE genes; ORA, over-representation analysis; CDT, Comparative Toxicogenomics Database; BP, biological process; MF, molecular function; CC, cellular component; PPI, protein-protein interaction; HIPPIE, Human Integrated Protein-Protein Interaction reference database; SScore, significance score; ROC, receiver operating characteristic; MTI, MiRNA-target interaction; MRM, miRNA-regulatory module; TSR, total relevance score; miEAA, the miRNA Enrichment Analysis and Annotation Tool.

Data Collection
Available omics expression data of HSCR patients were searched in Omics Discovery Index (OmicsDI-www.omicsdi.org(accessed on 2 January 2023) [26,27].The dataset inclusion criteria for the analysis were as follows: (1) presence of data on mRNA, miRNA and protein expression in colon tissue samples of HSCR patients; and (2) possibility of performing comparison groups of aganglionic and ganglionic segments of colon tissue samples.These requirements led to the following sets being selected for analysis: GSE98502 [28], GSE96854 [29], GSE77296 [30] and PXD021292 [31].A summary of these four individual studies is shown in Table S1.

Data Collection
Available omics expression data of HSCR patients were searched in Omics Discovery Index (OmicsDI-www.omicsdi.org(accessed on 2 January 2023) [26,27].The dataset inclusion criteria for the analysis were as follows: (1) presence of data on mRNA, miRNA and protein expression in colon tissue samples of HSCR patients; and (2) possibility of performing comparison groups of aganglionic and ganglionic segments of colon tissue samples.These requirements led to the following sets being selected for analysis: GSE98502 [28], GSE96854 [29], GSE77296 [30] and PXD021292 [31].A summary of these four individual studies is shown in Table S1.

Data Preprocessing and Differential Expression Analysis
Selected microarray datasets (GSE98502 [28], GSE96854 [29] and GSE77296 [30]) of HSCR patients in the GEO database [42] were individually investigated for differential expression (DE) patterns between ganglionic and aganglionic colon samples by using the same analytical tool, GEO2R [43] (GEO's online tool for analyzing GEO data, available at http://www.nci.nlm.nih.gov/geo/geo2r/(accessed on 2 January 2023)), to maintain consistency during individual analysis.To perform differential expression (DE) analysis, we selected Force normalization in the Options tab, the false discovery rate (FDR) p-value adjustment for multiple testing by Benjamini-Hochberg method [44], the log data transformation method to normalize the results and the cutoff value to filter out the results of FDR ≤ 0.05 and |logFC| ≥ 0.5.
GEO2R [43] provided a list of probes and corresponding gene/miRNA symbols ranked according to their degree of DE.In the second step, the gene/miRNA symbol annotation of probe sets was updated.If a particular gene/miRNA was mapped to several probe sets, the highest expression value was selected for further analysis, while probes without gene/miRNA annotations were removed.
On the other hand, the list of significant DE proteins (with a |fold change| ≥ 1.2 and a p-value ≤ 0.05) between the ganglionic and aganglionic regions of HSCR patients, which were in common among the male L-HSCR, male S-HSCR, and female S-HSCR assays, was downloaded from the originally published article [31].The "Retrieve/ID Mapping" tool of UniProt [45] was used to retrieve UniProt Knowledgebase (UniProtKB) proteins from the list of DE protein identifiers and to convert Uniprot protein ID into respective gene symbols.
Finally, the multi-symbol checker tool [46] was used to check if the gene nomenclature was HGNC-approved symbols.Otherwise, the miRBase miRNA-version and miRNA-toprecursor converter options, from the miRNA Enrichment Analysis and Annotation Tool (miEAA) miEAA 2.0 [47], were used to support research settings where older releases of miRBase [48] are in use.

Identification of Novel Candidate Disease Genes and Disease-Associated Modules
We employed a previously described method [49] for predicting disease risk genes, based on currently known disease-associated genes collected from DisGeNET [33] and differentially expressed genes according to omics expression data and PPI network analysis.
In the first step, we created the HSCR disease-relevant PPI network by selecting data from the Human Integrated Protein-Protein Interaction reference (HIPPIE) database [38] and determined high-density clusters in the PPI network according to DPClusO algorithm by using nine different and independent density values from 0.1 to 0.9 [49][50][51][52][53][54].
In a second step, cluster enrichment of DE disease-associated genes was determined by Fisher's exact test, and a significance score (SScore) [49], a measure of confidence of prediction of each gene based on the p-values of the clusters they belong to, was assigned to each gene.
In a third step, receiver operating characteristic (ROC) analysis [55,56] was used to evaluate the relevance of SScore [49] to predict DE disease-associated genes collected from DisGeNET [33] as well as the optimal cluster density value to discriminate among them.
Therefore, clustering stratification outcome (total number of clusters) and cluster membership (cluster composition) were optimized by identifying the cluster density which allow obtain the best HSCR disease-associated gene prediction.In this light, we searched the optimal "outcome-driven" stratification, measured by the known good-quality HSCR gene distribution across clusters obtained (based on area under the ROC curve values for outcome classification).
Finally, the statistically significant clusters (disease-associated modules), and the potential new disease-associated genes (comprised in these clusters) were determined after multiple testing correction by Benjamini-Hochberg method [44].

Prediction of Disease-Associated miRNAs
We applied a previously reported method [57] for predicting disease-associated miRNAs based on currently known disease genes.In summary, we built a set of genes related to HSCR by combining the reported HSCR disease-associated genes in the four databases previously mentioned [32-36] (Figure S1) together with the set of potential candidate disease genes identified in this study.This list was used as a query to find their miRNA targets and to create the MTI datasets.This information was collected from the three different validated online databases of MTI [39][40][41]58] with the same selection criteria described in the original method [57].
The interactions between miRNAs and gene targets obtained from each of the three interaction sources were individually represented as a bipartite graph, which is called MTI network.The biclusters (miRNA-regulatory modules, MRM) of these networks were determined by BiClusO algorithm [54,59,60].Then, we created HSCR-related sub-MRMs from the MRMs obtained by identifying the presence of HSCR genes.Finally, we assigned a relevance score (RS) [57] for each miRNA present in these sub-MRMs.The relevance score of ith miRNA (RS miRNA(i) ) is given by where NoofHSCR miRNA(i) is the number of HSCR genes attached to ith miRNA in the HSCR MRM set and Noofcluster miRNA(i) is the number of HSCR MRMs attached to ith miRNA.After normalizing the score of individual miRNAs in each dataset, we calculated the total relevance score (TSR) [57] for each miRNA as follows: where TRS miRNAi is the total relevance score of ith miRNA based on all datasets, RSni is the relevance score of ith miRNA in nth dataset, Cni the number of clusters in nth dataset and Eni is the Boolean value, measuring whether ith miRNA is in the nth dataset.Every miRNA with a total relevance score (TSR) > 0 was proposed in this study as a candidate disease-associated miRNA.

Identification of Potential miRNA Biomarkers in HSCR
To reveal potential candidate miRNA biomarkers in HSCR, we used the bioinformatics tool miRNA-BD [61].In this tool, the biomarkers are identified based on two parameters: the values of single-line regulation (NSR) and transcription factor percentage (TFP) on the condition-specific MTI network.
According to the user guide of this tool, miRNAs from the input data with significantly high NSR and TFP values (cutoff p-value ≤ 0.05 in Wilcoxon signed-rank test), calculated based on the disease-specific MTI network, were selected as a candidate disease biomarker.Moreover, we provided the disease-specific MTI networks, TF gene datasets and different kinds of input datasets (list of DE and/or disease-associated miRNA/mRNAs).
Thus, we identified potential miRNA biomarkers from (a) the list of predicted HSCR-associated miRNAs obtained in this study and (b) the list of DE miRNAs obtained from GSE77296 [30].
The HSCR-specific MTI network was constructed in the following way: (a) extracting from miRNet 2.0 [62] those target genes of DE miRNAs in HSCR and HSCR-known miRNAs according to three different databases (DIANA [39], miRecords [40] and miRTarBase [41]); and (b) extracting from miRNet 2.0 [62] the target miRNAs of DE genes in HSCR, known HSCR-associated genes and proposed new candidate HSCR genes obtained in this study.
Finally, the union of all these miRNA-target interaction pairs was selected to create the MTI network by the web-based tool OmicsNet 2.0 [63].NAP [64] was employed for visual analysis of the biological network.
Moreover, TF gene dataset was obtained from TRRUST v.2 [65], the reference database of human transcriptional regulatory interactions.

Functional Enrichment Analysis
Malachite [66], a Gene Enrichment Meta-Analysis (GEM) Tool for ToppGene [67], was used twice, to analyze and compare differences at once between (1) the DE gene datasets and (2) the statistically significant disease-associated modules identified.The enrichment categories tested were drugs, diseases, pathways and gene ontology terms.Moreover, the same enrichment categories were analyzed for the new candidate disease genes using ToppFun [67].Finally, the functional enrichment analysis of predicted miRNA was performed with miEAA 2.0 [47].

Disease Relevance Evaluation
The implication of candidate HSCR-related genes and miRNAs in enteric nervous system (ENS) and/or HSCR onset was assessed by performing an automated bibliography search using a previously published script [15].Moreover, the importance of the proposed HSCR-related genes in biological systems was explored by searching several databases [65,[68][69][70][71][72].The total number of HSCR-related genes attached to new predicted HSCR-related miRNAs in the MTI networks created for this study was determined.Predicted miRNAs interacting with the clinical gene panel for HSCR in ClinGen [73] were highlighted.TAM 2.0 database [74] was employed to identify specific colon tissue miRNAs among the predicted HSCR-related miRNAs.
Predicted miRNAs with altered expression in HSCR stenotic colon tissue was further studied.For this purpose, the HSCR miRNA expression profiles from the publicly available microarray dataset GSE77296 [30] was re-analyzed by GEO2R [43] and compared with our results.Moreover, the list of predicted miRNAs was compared with reported DE miRNAs obtained from the GSE77296 dataset [30] by other citing articles [30,[75][76][77].Additionally, we performed a manual search of altered expression in HSCR colon tissue [4,5,[78][79][80][81][82][83][84] for the predicted miRNAs previously associated with HSCR.On the other hand, a manual review of miRNAs previously described as potential peripheral HSCR biomarkers [5,75,81,84,85] was conducted.The miRNA disease network was created with miRNet [62] using as a query the list of predicted HSCR-related miRNAs to check if HSCR is present in it.

Statistical Analysis
We analyzed data statistically using SPSS 20.0 (IBM, Chicago, IL, USA).False discovery rate (FDR) correction by Benjamini-Hochberg method [44] of the uncorrected p-values was performed using the SDM FDR online calculator [91].The VENN DIAGRAMS tool [92] was used for comparing gene/miRNA lists.

Identification of DE Genes in HSCR and Gene Enrichment Meta-Analysis
The analysis of the two mRNA microarray datasets by GEO2R [39] allowed for the identification of 89 DE genes in GSE98502 (8 upregulated and 81 downregulated genes) and, 476 DE genes in GSE96854 (220 upregulated and 256 downregulated genes).
Moreover, we included a list of 651 DE proteins (coming from a total of 666 proteincoding genes) which were common among the male L-HSCR, male S-HSCR, and female S-HSCR assays previously published [31].This list included 48 DE proteins upregulated, 58 downregulated, and 545 with different tend of regulation based on HSCR phenotype and/or patient sex.
The Venn diagram of the overlapping DE genes between these three datasets is shown in Figure 2a.We took the union set of the DE genes from these three comparisons and combined these DE genes into a single set consisting of 1180 genes.We found that 31 of the DE genes were reported in the DisGeNET database as HSCR-related genes.We considered these 31 genes as "known HSCR DE genes" (HDEGs) and the other 1149 DE genes as "other DE genes" (ODEGs).
Interestingly, the DE genes are slightly overlapped (3.98%) in these three datasets.Moreover, overlapped genes are differentially expressed in the same direction of regulation across all datasets except for APOC3, GC and APOA1, which were downregulated at the mRNA level in GSE96854 but upregulated at the protein level in PXD021292 [31] of aganglionic bowel compared with ganglionic segments of colon samples.Figure 2b shows the complete details of the list of overlapped DE genes found in the three datasets.
We applied Malachite to these three datasets to conduct enrichment analyses on these HSCR omics expression data.The significant GO terms for "Biological Process" and "Cellular component" enriched in all three datasets are listed in Table S2.
Although no pathway was statistically significant enriched in all three datasets, we found that most GO terms' enrichment was implicated in nervous system developmental processes, namely axon, synapse, cell adhesion and neuron projection development, HSCRassociated GO terms.Moreover, we found 4 diseases and 29 drugs or small molecules also enriched in all three datasets (Table S2).All shared diseases presented neurological symptoms, while several of the small molecule compounds enriched in common were shown to alter the behavior of neurons, such as acetylcholine, dopamine, glutamate, GABA and serotonin, and/or there were substances involved in muscle contraction, such as Lant-6 and Lupex (Table S2).
Therefore, despite the limited number of common DE genes obtained from these three omics experiments, the number of significant shared enrichment terms in all three datasets is higher and comprises a bigger number of DE genes.

Overlapped sets
Total 1 Overlapped DEG With different tendencies 3 (2); 1 Total number of differentially expressed genes (DEG) present in the overlapped set. 2 Profile shows the tendency of regulation of DEG for GSE96854, GSE98502 and PXD021292 assay.The number in brackets shows gene counts. 3Gene with different tendencies of regulation at protein level on base of Hirschsprung patient's sex and disease phenotype (short-segment vs. long-segment). 4HDEGs. 5Gene downregulated at protein level 6 Gene upregulated at mRNA level 7 Gene upregulated at protein level.depicting the detailed information of the overlap genes that were significantly regulated in more than one sequencing experiment.

Identification and Enrichment Analysis of Candidate Disease-Associated Modules and Potential New Disease Genes
From the 1180 DE genes previously identified, we constructed the relevant HSCR PPI network based on the HIPPIE database.A total of 14,343 interactions were collected involving 3806 different proteins from 3850 protein-coding genes (29 HDEGs, 988 ODEGs and 2833 interacting genes named "other genes, OG").
The main topological characteristics of this PPI network obtained with NAP are shown in Table 1.Of note, this network is typically scale-free (fit power law TRUE), which gives our PPI network some important features such as stability, invariability to changes of scale and vulnerability to targeted attack.Network's diameter is large, with more than six-step separation, which means that proteins are not very closely linked in spite of the average path length being 3.83.The transitivity or clustering coefficient of our PPI network is low, showing that the network contains few communities or groups of nodes that are densely connected internally, which is also reflected by the low value of average number of neighbors.   11.00 Average path length 5  3.83 Clustering coefficient 6  0.05 Modularity 7  0.42 Number of self loops 8  269.00Average eccentricity 9  8.16 Average eigenvector centrality 10  0.03 Average number of neighbors 11  0.97 Centralization betweenness 12  0.18 Centralization degree 13  0.14 Fit power law 14  TRUE 1 PPI, Protein-protein interaction. 2Number of edges: shows the number of connections in the PPI network. 3Number of nodes: shows the number of vertices in the PPI network. 4Diameter: shows the length of the longest geodesic. 5Average path length: the average number of steps needed to go from a node to any other. 6Clustering coefficient: a metric to show if the network has the tendency to form clusters. 7 Modularity: this function calculates how modular is a given division of a graph into subgraphs. 8Number of self-loops: how many nodes are connected to themselves. 9Average eccentricity: the eccentricity of a vertex is its shortest path distance from the farthest other node in the graph. 10Average eigenvector centrality: this shows the influence of a node in a network. 11Average number of neighbors: how many neighbors each node of the network has on average. 12Centralization betweenness: betweenness centrality quantifies the number of times a node acts as a bridge along the shortest path between two other nodes. 13Centralization degree: This is defined as the number of links incident upon a node. 14Fit power law: the degree distribution of proteins approximates a power law P(k On the other hand, Table 2 provides the results of each DPClusO clustering and Figure 3 shows the area under the curve (AUC) for the ROC curves performed for each cluster density.We observed that a smaller density value resulted in a larger size and fewer number of clusters (Table 2).
The area under the ROC curves (AUC) described a poor discrimination with values in the range of 0.5 < AUC < 0.7, which may be due to incomplete information on known good-quality HSCR genes (Figure S1).The highest AUC (0.62) was obtained in the case of the cluster set generated using density = 0.5, for which this cluster dataset was selected to identify the candidate disease-associated modules and predict the new potential HSCR genes.A total of 55 out of the 563 initial statistically significant clusters were finally proposed as candidate disease-associated modules after adjusting the p-value for multiple testing (Table S3).These 55 clusters comprised a total of 195 different genes (17 HDEGs,59 ODEGs and 119 OGs).Therefore, after excluding the 17 HDEGs, these 178 other genes (with adjusted p-value < 0.05) were proposed as new HSCR candidate genes (Table S4).proximates a power law (k) = k −γ .
On the other hand, Table 2 provides the results of each DPClusO clustering and Figure 3 shows the area under the curve (AUC) for the ROC curves performed for each cluster density.We observed that a smaller density value resulted in a larger size and fewer number of clusters (Table 2).To validate our results, we searched how many of the predicted genes were exactly matched with reported HSCR genes by four reliable sources (Figure S1).
After considering overlapping between databases [32-36], we found that 14 out of 178 purpose genes (7.87%) matched with reported HSCR genes.Bearing in mind that our predictions are based only on a limited set of reported HSCR genes (a total of 527; Figure S1), the 7.87% coincidence with good-quality data is significant (p-value = 0.00027, calculated based on a hypergeometric distribution assuming a total number of human genes of 20,000).Even if we consider as universal just the total number of protein-coding genes in our PPI network (3850) and the set of included reported HSCR genes (196), the number of successes that HSCR genes found among our list of candidate HSCR genes is statistically significant (0.031).Moreover, after manually reviewing the results obtained from the automated bibliography search of the implication of these genes in enteric nervous system (ENS) and/or HSCR onset, we found that 39 of these 178 genes were associated with the ENS and/or HSCR development (Table S4).As expected, the 14 predicted HSCR-related genes that matched with reported HSCR genes in any of the four gene-disease databases (Figure S1) have been previously reported as implicated in the enteric nervous system (ENS) and/or HSCR onset.
Furthermore, we found that 25 of our candidate HSCR-related genes are described as essential, 17 as oncogenes, 25 as tumor suppressor genes, 44 as housekeeping genes and 8 as transcription factors (Table S4).
Gene enrichment analysis for all 178 genes was performed with ToppFun.The top 10 significant enrichment GO terms and pathways are shown in Table 3.According to the values observed both in the pathways and in the enriched GO terms, approximately 30% of the 178 candidate genes are related to cell adhesion processes (Table 3).Other cellular processes regulated by cell adhesion such as cell morphogenesis and tissue development in multicellular organisms were also found in the top 10 enriched GO: biological process terms (Table 3).Moreover, HSCR-related enriched terms were found both in the molecular function category (signaling receptor binding) and in the cellular component category (membrane raft and plasma membrane region) (Table 3).We also found several HSCR-related pathways among the Top10 most significant pathways (Table 3), for example, signaling events regulated by Ret tyrosine kinase obtained from Malacards [35,36], and developmental biology and signaling by interleukins obtained by filtering out from CTD using "Hirschsprung" as a keyword.
Cluster enrichment analysis for the 55 potential disease modules was performed with Malachite in a similar way to that which we mentioned above.REACTOME pathway analysis showed that several pathways related to signal transduction, such as the MAPK signaling pathway and Ras signaling pathway, were detected across more than 50% of potential disease modules (Figure 4).In this line, KEGG pathway analysis showed that the VEGF signaling pathway was also detected across more than 50% of potential disease modules.
Other HSCR-associated pathways (Figure 4) obtained by filtering out from CTD and Malacards were central carbon metabolism in cancer, thyroid cancer; B cell receptor signaling, signaling by ERBB4, GRB2 events in ERBB2 signaling, SHC1 events in ERBB2 signaling, SHC1 events in ERBB4 signaling and the Fc epsilon RI signaling pathway.
Finally, we proposed two new HSCR-associated pathways, since they were observed in more than 50% of potential disease modules (Figure 4).The proposed new HSCR-associated pathways were the prolactin signaling pathway and aldosterone-regulated sodium reabsorption (Figure 4b).
Malachite in a similar way to that which we mentioned above.REACTOME pathway analysis showed that several pathways related to signal transduction, such as the MAPK signaling pathway and Ras signaling pathway, were detected across more than 50% of potential disease modules (Figure 4).In this line, KEGG pathway analysis showed that the VEGF signaling pathway was also detected across more than 50% of potential disease modules.The only enrichment GO term for "Biological Process" found in more than 50% of significant clusters was peptidyl-tyrosine autophosphorylation (GO:0038083).The wide functional diversity of disease modules is reflected by the high number of significant singleton enrichment GO terms and pathways, which was also remarkable.In this line, we found more than 2371 different and unique GO (BP) terms.As expected, among the enrichment diseases, we found HSCR (C0019569) present in 35 clusters.

Prediction of HSCR-Associated miRNAs
According to our workflow (Figure 1), the candidate HSCR genes obtained in this study together with the previously reported HSCR genes obtained from four different sources (Figure S1) were combined in a unique dataset of 691 genes, which was employed to construct the HSCR MTI network from three different experimental validate databases (DIANA, miRecords and miRTarBae), whose individual characteristics are shown in Table S5.
The number of MRMs, the HSCR-related sub-MRMs and the number of comprised miRNAs and genes present in the sub-MRMs for each MTI network created, after bicluster analysis by BiclusO, are shown in Table 4. 1 MTI, miRNA-target interaction. 2MRM, miRNA-regulatory module. 3sub-MRM, Hirschsprung related bicluster. 4miRNA, microRNA comprised in sub-MRMs. 5Gene, number of genes comprised in sub-MRMs.
We calculated the total relevance score (TSR) for each miRNA, predicting 137 HSCRrelated miRNAs (Table S6) (with a TSR > 0) out of the initial 430 miRNAs, as shown in Figure 5.As expected, the sum of number of miRNAs encompassed in sub-MRMs obtained from each MTI network (Table 4) was larger than the total number of predicted HSCRrelated miRNAs after biclustering and TSR calculation, since a miRNA might be attached to more than one MTI network and because the relevance score (SR) of miRNAs encompassed in a sub-MRM without any direct interaction with a HSCR-related gene is equal to zero.

Assessment of the Relevance of the Identified Candidate miRNAs in HSCR
To validate our results, we searched how many of the predicted miRNAs are exactly matched with well-curated known HSCR miRNAs in five public databases (Figure S2).
Of 137 predicted miRNAs, 17 (12.41%)matched with good-quality known HSCR miRNAs, which, based on a hypergeometric distribution, is significant (p-value= 9.01× 10 −8 ) if assuming a total number of human mature miRNAs of 2300 [93].Even if we consider as universal just the total 430 miRNAs present in different MTI datasets before biclustering and TSR calculation, which included 30 HSCR miRNAs, the probability of pulling 17 good-quality known HSCR miRNAs among the list of 137 predicted HSCR-associated miRNAs is also statistically significant (p-value = 0.0022).According to literature reports, at least 35 of our predicted miRNAs have been previously described as associated with HSCR onset (Table S6).Remarkably, among this group of 35 HSCR-related miRNAs, we found 15 of the 17 well-curated known HSCR miRNAs in five public databases (Figure S2).Furthermore, we studied the relevance of the 137 predicted miRNAs in HSCR disease by different aspects.
First, we identified the total number of HSCR genes attached to them in the MTI networks created for this study (Table S6).Interestingly, it appears that HSCR-related genes are enriched in the top 10 miRNAs.The total number of HSCR genes attached to the top 10 miRNAs was 196, whereas the total number of HSCR genes attached to all 137 miRNAs was 325.Thus, an approximate ratio of 8.27:1 was achieved in terms of attachment to the HSCR genes for the top 10 miRNAs (Figure 6).

Assessment of the Relevance of the Identified Candidate miRNAs in HSCR
To validate our results, we searched how many of the predicted miRNAs are exactly matched with well-curated known HSCR miRNAs in five public databases (Figure S2).
Of 137 predicted miRNAs, 17 (12.41%)matched with good-quality known HSCR miRNAs, which, based on a hypergeometric distribution, is significant (p-value= 9.01× 10 −8 ) if assuming a total number of human mature miRNAs of 2300 [93].Even if we consider as universal just the total 430 miRNAs present in different MTI datasets before biclustering and TSR calculation, which included 30 HSCR miRNAs, the probability of pulling 17 good-quality known HSCR miRNAs among the list of 137 predicted HSCR-associated miRNAs is also statistically significant (p-value = 0.0022).According to literature reports, at least 35 of our predicted miRNAs have been previously described as associated with HSCR onset (Table S6).Remarkably, among this group of 35 HSCR-related miRNAs, we found 15 of the 17 well-curated known HSCR miRNAs in five public databases (Figure S2).Furthermore, we studied the relevance of the 137 predicted miRNAs in HSCR disease by different aspects.
First, we identified the total number of HSCR genes attached to them in the MTI networks created for this study (Table S6).Interestingly, it appears that HSCR-related genes are enriched in the top 10 miRNAs.The total number of HSCR genes attached to the top 10 miRNAs was 196, whereas the total number of HSCR genes attached to all 137 miRNAs was 325.Thus, an approximate ratio of 8.27:1 was achieved in terms of attachment to the HSCR genes for the top 10 miRNAs (Figure 6).
Fifth, after having created the miRNA disease network with miRNet, we found that 68 out of the 137 predicted miRNAs were associated with 268 diseases.As expected, we identified HSCR interacting with two miRNAs (hsa-miR-107 and hsa-miR-429) (Figure S3).From the network (Figure S3), the diseases with highest node connectivity degree (≥10) were related both to different kinds of cancer (e.g., hepatocellular carcinoma, lung cancer, prostate cancer and colorectal cancer), and other diseases, such as cardiac hypertrophy and Duchenne muscular dystrophy (both related to muscular dystrophies), Alzheimer disease (a neurodegenerative disease associated with colonic dysmotility/inflammation in its earliest stages) and ulcerative colitis (a chronic inflammatory bowel disease).
Fifth, after having created the miRNA disease network with miRNet, we found that 68 out of the 137 predicted miRNAs were associated with 268 diseases.As expected, we identified HSCR interacting with two miRNAs (hsa-miR-107 and hsa-miR-429) (Figure S3).From the network (Figure S3), the diseases with highest node connectivity degree (≥10) were related both to different kinds of cancer (e.g., hepatocellular carcinoma, lung cancer, prostate cancer and colorectal cancer), and other diseases, such as cardiac hypertrophy and Duchenne muscular dystrophy (both related to muscular dystrophies), Alzheimer disease (a neurodegenerative disease associated with colonic dysmotility/inflammation in its earliest stages) and ulcerative colitis (a chronic inflammatory bowel disease).
Seventh, we found 27 HSCR-related lncRNA-miRNA interactions involving 27 of the miRNAs predicted as HSCR-associated in this study after constructing a miRNA-HSCRrelated lncRNA network in miRNet (Table S8).Interestingly, MEG3 and other dysregulated circRNAs identified in HSCR were predicted to act as sponges of hsa-miR-326.
On the other hand, 65 HSCR-associated pathways from CTD were enriched terms in our predicted miRNAs (Table 5).Interestingly, 11 out of these HSCR-associated pathways (Table 5) were in the Top 20 of the most significant enriched terms in the categories KEGG (miRPathDB) and Reactome (miRPathDB), to know cellular responses to stress, cellular senescence, cytokine signaling in immune system, disease, diseases of signal transduction, immune system, microRNAs in cancer, pathways in cancer, PIP3 activates AKT signaling, signal transduction and signaling by interleukins.Moreover, we looked for the prolactin signaling pathway (with 37 observed miRNAs and an adjusted p-value = 4.68 × 10 −13 ), and aldosterone-regulated sodium reabsorption (with nine observed miRNAs and an adjusted p-value = 0.002) among the statistically significant pathways for the category KEGG (miRPathDB), confirming that these proposed HSCR-related pathways were also enriched in the predicted miRNAs.Additionally, Table 6 shows the enriched terms significantly over-represented obtained for the category Annotation (gene ontology).Remarkably, 14 of these significant enriched GO terms were associated with HSCR.  2  1.99 × 10 −3 8 Extracellular exosome GO:0070062 2  3.66 × 10 −3 15 Negative regulation of cardiac muscle cell apoptotic process GO:0010667 2  4.25 × 10 −3 10 Table 6.Cont.

Potential miRNA Biomarker Identification in HSCR Based on miRNA-Target Regulatory Network Analyses
We generated a reference HSCR-specific MTI network based on DE miRNAs/genes in HSCR and HSCR-known associated miRNAs/genes from five different databases [32][33][34][35][36][37].The reference HSCR-specific MTI network contained 129,139 microRNA-target regulatory pairs.After calculating NSR and TFP parameters, miRNAs with significantly high values were highlighted (Table 7).Overall, 12 predicted miRNAs in this study were screened as candidate biomarkers.Except hsa-miR-30a-5p, all of them have been previously described with a dysregulated expression in HSCR colon tissue.Among them, three were present as DE miRNAs in GSE77296 (hsa-miR-200b-3p, hsa-miR-148a-3p and hsa-miR-429).However, when we used as input the list of DE miRNAs in GSE77296 to screen their potential as biomarkers, we only found two candidates (hsa-miR-146a-5p and hsa-miR-200b-3p).Therefore, the only overlapped miRNA was hsa-miR-200b-3p.

Discussion
Network-based computational approaches are being used to identify disease-associated subnetworks and for efficient disease-gene association prediction [94].Biological network studies are based on two underlying assumptions: Firstly, human diseases are the consequences of disruption in molecular networks [95], and secondly, genes associated with the same or similar diseases tend to reside in the same neighborhood of a PPI network [96,97] or a pathway [98].
Network-based methods have been successfully applied in the study of rare diseases [99][100][101][102], and specific methodology has been developed for searching and ranking disease-causing genes in this kind of disorders [103][104][105].According to Orphanet information, 71.9% of rare diseases are genetic and 69.9% are exclusively pediatric onset [106], such as HSCR.However, HSCR and other rare conditions are characterized by a wide diversity of symptoms and signs that vary from patient to patient, and their pathogenesis is not fully understood because of our partial knowledge of the molecular basis and causative genes [4].For this reason, in this study, we aimed to test if reported methods [49,57] previously validated in the context of inflammatory bowel disease (IBD, a non-rare disease) can be generalized to find disease modules and disease-associated genes for HSCR.
For the employed of the first method [49], a condition-specific PPI network was constructed by mapping the gene expression data obtained in transcriptome and proteome experiments onto the global PPI network, and by locating well-known disease genes (deposited in DisGeNET [33], which collects disease-gene association obtained from experimental assays such as GWAS, NGS, etc.) in the network.This strategy allowed us to propose up to 178 HSCR-related genes and validate this approach after identifying 14 well-curated HSCR genes among our candidates.Therefore, these results could involve an increase in the number of genes potentially involved in HSCR of up to 33.78% considering the 527 HSCR genes according to the four used databases, allowing us to guide and focus subsequent studies.This value is considerably higher than that obtained by other strategies also aiming to identify new HSCR-associated genes [14 -17], which shows the efficacy of this cost-free study, even without being exempt from certain limitations and weaknesses, such as (i) the lack of consideration of expression differences between the postnatal stenotic tissue (our input) and the developing fetus in which the neurogenesis occurs [28]; (ii) the limited number of published high-throughput omics expression studies in HSCR and the absence of these types of data in public repositories [107][108][109][110]; (iii) the absence of demographic information and clinical disease parameters from all participants of the included studies [28][29][30][31]; and (iv) the employment of only some of the available integrative database resources for disease-gene associations [32][33][34][35][36]111] both for gene-disease module identification and result validation.
This enrichment analysis revealed that a high percentage of our candidate genes were involved in previously HSCR-associated pathways and GO terms [16,17].In this line, it is noteworthy that approximately 30% of the 178 candidate genes were related to cell adhesion processes (Table 3).This biological functionality is related to HSCR across different patient populations [16].In a biological context, cell adhesion molecules (CAMs) and fibroblast growth factors (FGFs) are known to stimulate neurite growth through the activation of various FGF receptors (FGFRs) on neurons [112].The marked decrease in CAM expression in the stenotic intestine of HSCR patients suggests that this alteration would lead to defects in the migration of enteric neuroblasts [112].Furthermore, approximately 20% of these candidate genes were described as tumor suppressor and/or oncogenes in any context (Table S4), which is interesting considering that all clinical gene panels for HSCR in ClinGen have been suggested as important players in the development of human carcinomas [113][114][115][116].In this line, HSCR has been reported in the co-occurrence of several kinds of cancers, such as familial medullary thyroid carcinoma, multiple endocrine neoplasia type 2A and type 2B [117], neuroblastoma [118,119] hepatoblastoma [120,121], central nervous system tumors [122] and retinoblastoma [123,124].
Moreover, this approach also allowed us to identify 55 significant disease gene modules enriched in known HSCR DE genes in patient colon tissue samples.As expected, comparative enrichment analysis showed that several of the significantly enriched pathways present in more than 50% of these clusters (Figure 4) were related to previously HSCR-associated pathways, such as VEGF signaling, which results in the activation of multiple downstream pathways, including the Ras/MAPK pathways [125,126] (pathways systematically associated with HSCR [16,17]), regulating cell proliferation and gene expression; and/or the PLCγ pathway, controlling vascular permeability [127].In this line, despite there being no reports of increased vascular permeability in the gastrointestinal tract (GIT) of HSCR patients, a recent study has suggested that vascular permeability may be etiologic for postoperative Hirschsprung-associated enterocolitis (HAEC) [128].In addition, vascular endothelial growth factor A (VEGF-A) mediates angiogenesis, and altered vascular density has also been reported in GIT of HSCR patients and in some mouse models of HSCR [129].However, contradictory results were observed when comparing mouse and human data, which might be due to several reasons, to know that (a) the human data were not standardizable; (b) for the analysis of the vascular density in colon samples of HSCR patients, only one area per sample (which exhibited the highest numbers of microvessels) was investigated; and (c) the statistical analysis was performed using a Mann-Whitney test on two dependent samples because analyzed ganglionic and aganglionic colon tissue samples belong for the same group of HSCR patients.Therefore, in spite of the relevance of these preliminary results, their specific clinical importance in HSCR is as yet unclear [128].
Additionally, other two KEGG pathways-the prolactin signaling pathway and aldosterone-regulated sodium reabsorption pathway-were proposed as candidate HSCRrelated pathways.Prolactin-like material (PLM) has been previously identified in the gut mucosa, but is absent for aganglionic segment of HSCR patients [130].In spite of the lack of experimental validation assays, it has been suggested that, as pituitary prolactin, PLM might be involved in metabolism and osmoregulation [130], like the aldosterone-regulated sodium reabsorption pathway (Figure 4b).This pathway is an electrogenic sodium transport through the amiloride-sensitive epithelial sodium channel (ENaC) [131].Along the intestine, electrogenic absorption through EnaC is limited to surface epithelial cells of the distal colon and rectum [132].The distal aganglionic colon segment in HSCR patients leads to abdominal distention with intraluminal loss of chloride, potassium, and sodium, resulting in metabolic alkalosis with secondary hyperreninemia and hyperaldosteronism [133].However, after proctocolectomy, EnaC starts to be expressed in the distal part of the small intestine.Therefore, possible symptoms of the pseudo-Bartter syndrome in neonates with HSCR completely disappeared after a pull-through operation with the removal of the causal factor [133].The exact importance of this pathway in HSCR has not been clarified yet but it is known that Endothelin receptor type B (EDNRB), the second most mutated gene in Hirschsprung (5% cases) [1], is an antagonist of EnaC activity [134].Moreover, the downregulation of EnaC with a reduction in sodium reabsorption in the colon may contribute to diarrhea associated with IBD [131], a disease with relatively low co-occurrence in HSCR patients with a similar clinical presentation to HAEC [135,136].
Recently, miRNAs have received widespread attention due to their role in the regulation of cellular processes [137] and the increase in experimental evidence postulating that their aberration can be among the triggering causes of pathological phenomena, including HSCR [85, [138][139][140].Additionally, miRNAs are endogenous, stable, easily accessible in the extracellular space and detectable in patients' tissue and blood specimens [85].Therefore, several authors have suggested that they may have utility as diagnostic and prognostic biomarkers for HSCR [5].For these reasons, our second aim was to decipher key miRNAs related to HSCR by using computational tools to unravel miRNA-disease associations and to identify miRNA biomarkers.
The application of biclustering algorithms can solve the traditional problems associated with the discovery of regulatory modules that control gene transcription in biological model [141].The biclustering algorithm of BiClusO [54,59,60], as well as the one employed in ComiRNet [142], allowed for the efficient discovery of overlapping and highly cohesive biclusters.Then, an overall relevance score (TSR) [57] based on the exploration of only those biclusters containing genes of interest (known disease-associated genes) was applied to predict Hirschsprung-related miRNAs in this study.The use of this second previously reported method, validated for IBD, is suitable for isolating miRNAs for similar types of diseases [57].Although HSCR was not identified in its miRNA-disease network, as we mentioned above, the co-occurrence of IBD, especially Crohn's disease, has been recently reported in children with HSCR [133,134], and common pathogenesis and hub gene biomarkers have been also identified [143].
In fact, the application of this method allowed us to predict 137 HSCR-related miRNAs and validated the method, which is especially relevant for a rare disease with a limited knowledge of its genetic landscape.Remarkably, we found five RET-regulated miRNAs (hsa-miR-15a-5p, hsa-miR-27b-3p, hsa-miR-195-5p, hsa-miR-218-5p and hsa-miR-128-3p) and two L1CAM-regulated miRNAs (hsa-miR-1-3p and hsa-miR-34a-5p) among our predicted HSCR-related miRNAs.Moreover, we identified three miRNAs with specific colon tissue expression among our predicted miRNAs, with has-miR-194-5p being the only one without a detailed study of its implication in HSCR to our knowledge [75].
However, when we tried to evaluate the relevance of the predicted miRNAs in the disease context, we noticed some methodological considerations that should be mentioned: Firstly, a majority of microRNA-mRNA interactions remain unidentified [144].However, to avoid a high rate of false positives among our predicted miRNAs, we refused the use of predictive algorithms, and only considered experimentally validated miRNA-target interaction [58].Secondly, other molecules such as lncRNAs and circRNAs can affect the functions of miRNAs in gene expression and play an important role in regulating cell physiological functions and diseases [145].Therefore, an integrative analysis of the lncRNA/circRNA-miRNA-mRNA ceRNA network would be required.However, the current availability of lncRNA and circRNA expression data in HSCR patients is not enough to carry it out, although our results also pointed to the interconnection of these molecules in regulatory modules.Concretely, we found that six of our predicted miRNAs were predicted targets of dysregulated circRNA in HSCR.Nevertheless, only hsa-miR-142-3p and hsa-miR-338-3p have been described as significantly upregulated in stenotic segment tissues of HSCR patients [30].Of note, several experimental validated target mRNAs of these two miRNAs were also dysregulated in intestine samples of HSCR patients and were related to the pathophysiology of this disease, such as ATP2A2, CSE1L, LRP8 and THBS4 (target of hsa-miR-142-3p) and CDH2 and PKM (target of hsa-miR-338-3p).Finally, the software tool employed for miRNA biomarker identification is based on just two significant parameters, but it lacks a disease-specific signal in this approach [61].Nonetheless, all putative biomarkers obtained (Table 7) from our list of predicted HSCR-associated miRNAs were involved in disease pathogenesis according to well-curated known HSCR miRNAs in five public databases and literature reports [30] (Table S6).
On the other hand, functional analysis of predicted miRNAs revealed the meaningful results obtained in this study.These results will serve as a valuable resource for further exploration by experiments.For example, it is potentially important that the putative biomarkers originally identified were able to be validated in blood samples, since blood could be easily obtained for clinical translation [5,84].

Conclusions
In summary, large-scale high-dimensional omics expression data analyses have been successfully applied to the discovery of functional connections among diseases, genetic perturbation, and drug action.However differential expression analysis alone does not provide any functional insights.For this reason, we designed the bioinformatics framework proposed in this study, which has allowed us to decipher putative modules, genes and miRNAs in HSCR, which could open new lines of research on this disease.In this line, this study has identified new candidate genes and miRNAs not previously associated with HSCR, but which are involved in important biological processes related to enteric nervous system development or HSCR onset.The employed approach in this study also was useful for amplification and facilitating research into HSCR specific-functional connections such as the proposed HSCR pathways: the prolactin signaling pathway and aldosterone-regulated sodium reabsorption pathway.Therefore, we consider that this methodology can be also successfully applied to other rare complex diseases for which conducting large, randomized trials is difficult.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biom14020164/s1, Figure S1: Venn diagram representing the overlaps among Hirschsprung-related genes collected from four different databases; Figure S2: Venn diagram showing overlapping between Hirschsprung-related miRNAs collected from five different sources; Table S1: Characteristics of the selected datasets for this study; Table S2: Terms enriched in all three datasets comparing stenotic colon to ganglionic colon samples; Table S3: List of potential diseaseassociated modules; Table S4: List of predicted Hirschsprung-related genes; Table S5: Hirschsprung miRNA-target interaction information collected from three different databases; Table S6: List of predicted Hirschsprung-related miRNAs; Table S7: Dysregulated CircRNAs identified in Hirschsprung's disease with predicted target miRNAs proposed in this study; Table S8: Hirschsprung-related lncR-NAs with predicted target miRNAs proposed in this study; Table S9: Results obtained for Gender and Age category after ORA analysis performed on predicted HSCR-related miRNAs.Funding: This work was supported by the Instituto de Salud Carlos III (ISCIII), Spanish Ministry of Economy and Competitiveness, Spain and co-funded by the European Union (ERDF, "A way to make Europe") (PI19-01550, PI22-01428), The strategic plan for the Precision Medicine Infrastructure associated with Science and Technology-IMPaCT (IMP-0009), IMPaCT-Data (IMP-00019) and Regional Ministry of Health and Families of the Autonomous Government of Andalusia (PEER-0470-2019) and co-funded by European Union.H.L.-P. is supported by fellowship CD20/00171 (Sara Borrell) from ISCIII and co-funded by European Union (ERDF/ESF, "A way to make Europe"/"Investing in your future").

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Figure 1 .
Figure 1.Flowchart containing the main steps followed on this study.The diagram is divided into four major sequential blocks, listed one through four.Additionally, the complementary working packages "5 • Functional enrichment analysis" and "6 • Disease relevance evaluation" are shown in yellow.Arrows show the direction of the workflow.Processes are represented by a rectangle, while the ovals correspond to the input/output data of each process.HSCR, Hirschsprung disease; OmicsDI, Omics Discovery Index; GEO2R, interactive web tool of Gene Expression Omnibus (GEO); miRNA, microRNA; HDEGs, known HSCR DE genes; ODEGs, other DE genes; ORA, over-representation analysis; CDT, Comparative Toxicogenomics Database; BP, biological process; MF, molecular function; CC, cellular component; PPI, protein-protein interaction; HIPPIE, Human Integrated Protein-Protein Interaction reference database; SScore, significance score; ROC, receiver operating characteristic; MTI, MiRNA-target interaction; MRM, miRNA-regulatory module; TSR, total relevance score; miEAA, the miRNA Enrichment Analysis and Annotation Tool.

Figure 2 . 2 .
Figure 2. Gene overlapping bioinformatic analysis.(a) Venn diagram showing the number of common differentially expressed genes among the ganglionic and stenotic colon segment groups from Figure 2. Gene overlapping bioinformatic analysis.(a) Venn diagram showing the number of common differentially expressed genes among the ganglionic and stenotic colon segment groups from the transcriptomic (GSE96854 and GSE98502) and proteomic (PXD021292) experiments; (b) Tabledepictingthe detailed information of the overlap genes that were significantly regulated in more than one sequencing experiment.

Figure 3 .
Figure 3. Plot of AUCs corresponding to nine density sets of clusters after ROC analysis.Figure 3. Plot of AUCs corresponding to nine density sets of clusters after ROC analysis.

Figure 3 .
Figure 3. Plot of AUCs corresponding to nine density sets of clusters after ROC analysis.Figure 3. Plot of AUCs corresponding to nine density sets of clusters after ROC analysis.

Figure 6 .
Figure 6.Total score (TRS) of top 10 predicted HSCR-related miRNAs with the number of attachments to different databases.

Figure 6 .
Figure 6.Total score (TRS) of top 10 predicted HSCR-related miRNAs with the number of attachments to different databases.

Table 1 .
Topological parameters of the protein-protein interaction (PPI) network constructed in the present study obtained with NAP.

Table 2 .
Results of each DPClusO clustering based on the HSCR relevant PPI network.

Table 2 .
Results of each DPClusO clustering based on the HSCR relevant PPI network.

Table 3 .
Top 10 significant pathways and GO terms obtained after enrichment analysis of predicted Hirschsprung-related genes with ToppFun.

Table 4 .
Summary of BiclusO clustering results obtained in this study for each Hirschsprung MTI network.
1Total number of observed miRNAs.

Table 6 .
Enriched terms significantly over-represented obtained for the category Annotation (gene ontology) of predicted HSCR-associated miRNAs with miEAA.

Table 7 .
Prediction values for candidate miRNA biomarkers identified for Hirschsprung disease (HSCR) based on the model with the reconstructed HSCR miRNA-mRNA reference network.Values of single-line regulation (NSR) and transcription factor percentage (TFP) parameters were calculated to characterize the regulatory pattern of the miRNAs in the constructed reference network and to identify miRNAs that could have a crucial role in the progression of HSCR for each (a) HSCR miRNA predicted in this study; (b) differentially expressed (DE) miRNA in GSE77296.