{'fraction_non_alphanumeric': 0.03858750419885791, 'fraction_numerical': 0.02544507893852872, 'mean_word_length': 4.997733568370688, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 3, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of Hub Genes Associated with Hypertension and Their Interaction with miRNA Based on Weighted Gene Coexpression Network Analysis (WGCNA) Analysis

Background Hypertension is one of the most widespread health conditions in the world, and the molecular mechanism of it is still unclear. In this study, we identified the hub genes (hub miRNA genes) associated with hypertension and explored the relationship between hypertension miRNA-gene by constructing a mRNA co-expression network and a miRNA co-expression network, which can help to reveal the mechanism and predict the prognosis of hypertension progression. Material/Methods Based on gene expression profile data of hypertensive samples from the Gene Expression Omnibus database, WGCNA was used to detect hypertension-related biomarkers and key mRNA and miRNA modules. Then, DAVID was used to perform gene-annotation enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) and miRPath were used for pathway analysis of mRNA and miRNAs genes. Results We identified 3 key modules relating to hypertension, 2 mRNA modules named Msaddlebrown and Mgreenyellow and 1 miRNA module named Msalmon. In addition, 12 hub genes (RPL21, RPS28, LOC442727/PTGAP10, LOC100129599/RPS29P14, TBXAS1, FCER1G, CFP, FURIN, PECAM1, IGSF6, NCF1C, and LOC285296/UNC93B3) and 7 hub miRNAs (hsa-miR-1268a/b, hsa-miR-513c-3p, hsa-miR-4799-5p, hsa-miR-296-3p, hsa-miR-5195-5p, hsa-miR-219-2-3p, and hsa-miR-548d-5p) relating to hypertension were identified. HIF-1 signaling pathway and insulin signaling pathway were closely related to the 3 key modules. We also discovered 4 miRNAs (hsa-miR-548am-3p, hsa-miR-513c-3p, hsa-miR-182-5p, and hsa-miR-548d-5p) and 6 genes (IGF1R, GSK3B, FOXO1, PRKAR2B, HIF1A, and PIK3R1) were the core nodes in the hypertension-related miRNA-gene network, and hsa-miR-548am-3p was at the center of the network. Conclusions These findings will help improve the understanding of the pathogenesis of hypertension, and the discovered genes can serve as signatures for early diagnosis of hypertension.


Background
Hypertension, also known as high blood pressure, is one of the most widespread health conditions in the world. It is one of the most dangerous factors affecting cardiovascular death. Every year, more than 17 million people die from cardiovascular disease. Hypertension is responsible for more than 45% of deaths due to heart disease, and 51% of deaths due to stroke [1]. The risk factors for hypertension include age, ethnicity, weight, diet, alcohol and tobacco use, gender, as well as existing health condition such as diabetes, high cholesterol levels, and chronic kidney disease. The pathogenic mechanism is even more complex, involving a variety of molecules and pathways [2].
Currently, hypertensive patients are treated with drugs such as diuretics, vasodilators, rapid-acting intravenous antihypertensive agents, and other drugs. Unfortunately, hypertension medications have many side effects such as cough, diarrhea, dizziness, drowsiness, depression, ulcers, and more. Since hypertension is a multiple-factor disease, further research is needed to reveal the molecular mechanism of hypertension and new biomarkers need to be discovered. The increased availability of high-throughput sequencing and microarrays along with the development of bioinformatics tools and algorithms has allowed for the discovery of several biomarkers that regulate blood pressure. For example, the anti-inflammatory cytokine interleukin-10 (IL-10) has a strong antihypertensive effect [3]; ENOS gene can increase plasma nitric oxide levels to reduce blood pressure [4]; and renin-angiotensin system (RAS) improves insulin resistance and prevents the development of renal hypertension [5]. Many miRNAs related to hypertension have also been discovered [6], such as miR-155 [7], miRNA-126 [8], miR-124 [9], and miR-150 [10]. Most of these studies only focused on one single gene or miRNA, and the identification of these gene targets are limited to just differential expression. Very few studies focused on expression profiles of multiple genes, and there is also insufficient attention to the high degree of interconnection between genes.
The network of interactions between biomolecules provides an important basis for systematic research on disease. WGCNA is an R package, which is based on the similarity between genes to construct a weighted correlation network [11]. It has unique advantages in handling complex data with multiple samples, which is a powerful method to uncover basic mechanism of gene-disease relationships [12]. Using WGCNA, Zhang et al. discovered ten hub genes that could be used as biomarkers for oral squamous cell carcinoma tumors [13]. In another study, six hub genes were found to regulate the signaling pathway of clear cell renal cell carcinoma (ccRCC) [14]. Wu et al. applied WGCNA to the identification of potential therapeutic targets for angiotensin II (Ang II) induced hypertension [7].
Here, we utilized the WGCNA method to construct gene and miRNA modules, identified new biomarkers related to hypertension, and explored the relationship between hypertension and marker genes.

Material and Methods
Data sources and preprocessing GSE75360, GSE75670, and GSE117261 datasets were downloaded from Gene Expression Omnibus (GEO) (http://www. ncbi.nlm.nih.gov/geo/). GSE75360 contained 10 hypertensive and 11 normal human gene expression data, (Illumina HumanHT-12 v.4.0 Expression BeadChip). GSE75670 contained 6 hypertensive and 6 normal human miRNA expression data (Exiqon mercury™ LNA™ microRNA array, 7 th generation [miRbase v18]). GSE117261 contained 58 pulmonary arterial hypertension and 25 normal human gene expression data ([HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array [transcript (gene) version]). Data preprocessing procedures such as correction of expression matrix, quality evaluation of expression data, and sample clustering were conducted in R version 3.5.2. Network analysis with WGCNA software package was also conducted in R version 3.5.2. The overall analysis workflow is showed in Figure 1.

Module construction base on WGCNA algorithm
WGCNA explores the complex relationships between genes and phenotypes by constructing scale-free co-expression networks. It transforms gene expression data into co-expression modules and then identifies hub genes in the modules. In this research, we constructed 2 weighted co-expression networks based on mRNA and miRNA expression data separately. The construction processes were the same except that some parameter values were different, so we only introduced the construction process of gene co-expression network. First, the correlation of all gene pairs was calculated to construct a similarity matrix. Second, the soft threshold b was a weighted parameter of the adjacent function, which the optimal value was obtained by the pickSoftThreshold function in the R package WGCNA [15]. Third, the TOM similarity function was used to convert the adjacency value into a TOM matrix. Then, using dissimilarity matrix dissTOM=1-TOM, we clustered the genes into the hierarchy to get the system clustering tree [13]. Fourth, mRNAs with similar expression profile were divided to the same module. According to the number of the genes and miRNAs, the minModuleSize of the mRNA was set to 50 and the minModuleSize of the miRNA was set to 30 [16]. Finally, we calculated the differences of the modules eigengenes, and set an appropriate cutline for the modules dendrogram to merge highly similar modules.

Identification of key modules in co-expressed networks
Two methods were used to identify key modules related to hypertension. The first method was to calculate the Pearson correlation coefficient and significance P-value of module eigengenes (MEs) and hypertension trait. Here, MEs represents the overall expression level of the gene module [17]. The second method calculates the gene significance (GS) and the module significance (MS). Here, GS is the correlation between a gene and the clinical features. MS is the average GS of all genes in a module [13]. Generally, the higher the absolute value of MS and GS, the more relevant the gene module is to hypertension.

Enrichment analysis of key modules
To further understand the function of a key module and its biological significance, we used the online functional annotation database DAVID (https://david. ncifcrf.gov/). For the mRNA modules, we used Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis to identify significantly enriched pathways (P-value <0.05). For miRNA modules, the online database miRPath V.3 (http://snf-515788.vm.okeanos.grnet.gr/) was used to predict miRNAs target genes for KEGG pathway enrichment [18].

Identify hub genes and hub miRNAs
This study used 2 methods to identify hub genes (miRNA) from key modules: 1) Importance threshold, and 2) MCC (maximal clique centrality function) algorithm which used the cyto-Hubba plugin. In the first method, hub genes are defined as genes with high GS value, high modular membership (MM) value, and low weighted P-value associated with genes and hypertension (P.weighted). MM was used to measure the importance of genes in modules [19]. The P-weighted value was calculated by the networkScreening function in WGCNA. A P.weighted value less than 0.05 was considered biologically significant [20]. In the second method, we used the MCC algorithm in cytoHubba plug-in from the Cytoscape software to identify hub genes [21].

Gene selection and hub gene verification
We used limma R package to analyze differently expressed genes (DEGs) between normal samples and hypertension samples in the dataset GSE75360 and set the cutoff value to log-2FC< |0.182| and P-value <0.05. The volcanic map and hierarchical clustering analysis were performed by "ggplot2" and "pheatmap" package of R, respectively. Then, we used the jvenn (http://jvenn.toulouse.inra.fr) to draw Venn diagrams to overlap the genes in DEGs and hub genes [22]. We verified our hub genes using a different DEG dataset (GSE117261). We queried the role of hub genes through The Human Protein Atlas database (https://www.Proteinatlas.org/) and used NCBI (https://www.ncbi.nlm.nih.gov/gene/) to verify whether these hub genes could be used as biomarkers of hypertension.
Construction of miRNA-gene interactive network DAVID and miRPath were used to identify enriched mRNA and miRNA pathways, respectively. Then, the overlapping pathways were chosen for network construction. Cytoscape was used to construct miRNA-gene interaction networks.

Data preprocessing
After preprocessing, an expression matrix containing 29 595 genes was obtained. We calculated standard deviation (SD) of all genes, then, ranked the SD from large to small, and selected the top 6000 genes as the input data for the construction of the gene expression network. We performed the same preprocessing procedures on the miRNA expression matrix. A total of 1916 miRNAs were selected for co-expression network construction of miRNAs.

Construction of weighted co-expression network
To obtain a network that meets the scale-free topology criterion, we calculated network structures by using different softthresholding power range from 1 to 20. The scale-free topological fitting index of mRNA co-expression network reaches 0.8, when the soft-thresholding power was 14 (Supplementary Figure 1A, 1B), which met the scale-free network criterion. Then, the genes were clustered into modules by hierarchical clustering according to expression value, and the most similar modules were merged by setting the MEDissThres cutting line to 0.2 (Supplementary Figure 1C). Finally, 19 mRNA gene modules were identified ( Figure 2). For miRNA co-expression networks, the soft-thresholding power was set at 9 and the scale-free topological fitting index at 0.88 (Supplementary Figure 2A, 2B). The MEDissThres cutting line was set to 0.1 (Supplementary Figure 2C), and 14 miRNA gene modules were identified ( Figure 3). The relationship for each mRNA module was analyzed by drawing a network heat map (Supplementary Figure 3). The network heat map was produced with the modules in the miRNA co-expression network ( Figure 4). The different colors in the vertical and horizontal axes stand for different miRNA modules. The yellow color in the middle area indicates a degree of connection for each miRNA module. The figure shows no significantly differences in the interaction between the modules, indicating that the miRNA modules had relatively high independence.

e923514-4
Identifying of key modules We used 2 methods to identify key modules. The first method calculate the Pearson correlation and significance P-value between MEs and hypertension, and the other method calculates GS. We found 2 mRNA modules with significant correlations to hypertension, named M saddlebrown and M greenyellow . Compared to other mRNA modules, the correlation coefficient of theese two modules were the largest ( Figure 5A). In the miRNA coexpression network, the M salmon module was highly related to hypertension ( Figure 5B). In order to ensure the identified modules were significantly associated with hypertension, we calculated the GS in the modules and verify the key modules again ( Figure 6A, 6B). For the mRNA co-expression network, the M saddlebrown and M greenyellow modules had the largest GS scores at 0.451 and 0.410, respectively. For miRNA co-expression networks, the absolute GS value of M salmon was 0.398. These GS values all indicate that these modules were significantly associated with hypertension.

Enrichment analysis of key modules
We performed GO analysis and KEGG analysis on the genes in the M saddlebrown and M greenyellow modules to explore their biological significance (Supplementary Tables 1, 2). The M saddlebrown module was related to RNA transcriptional translation, such as  The cluster dendrogram of mRNA in mRNA expression data, each branch represents a gene, and each color below represents a co-expression module. The first ribbon represents the module detected by dynamic tree cutting, and the second ribbon represents the module after merging the similar module.

DATABASE ANALYSIS
SRP-dependent co-translational protein targeting to membrane, translational initiation, and RNA binding. The M greenyellow module was involved in a variety of immune and metabolic processes, such as positive regulation of IL-6 production, positive regulation of transcription from RNA polymerase II promoter, glycoside catabolic process, and so on. Furthermore, KEGG pathway analysis showed enrichment in ribosome, HIF-1 signaling pathway, and osteoclast differentiation pathways. The significant pathways in M salmon module were Hippo signaling pathway, adherens junction, and proteoglycans in cancer (Supplementary Table 3). Among them, the HIF-1 signaling pathway and the insulin signaling pathway were shared by both mRNA and miRNA modules.

Identification of hub genes and miRNAs
According to the definition of module connectivity, we calculated the MM and GS of the genes (miRNAs) in each of the key modules to select the hub genes (miRNAs). Then, we used the networkScreening function to obtain the P.weighted of each gene (miRNA). |MM| >0.8, |GS| >0.2 and P.weighted <0.05 were used as the identification criteria. Finally, 26 genes were obtained in the M saddlebrown module (Supplementary Table 4), 53 genes were obtained in the M greenyellow module (Supplementary Table 5), and 22 miRNAs were obtained in the M salmon module (Supplementary Table 6). We also import the files of these 3 key modules into Cytoscape software and used the plugin cy-toHubba to identify and visualize pivotal genes and miRNAs (Supplementary Figure 4A-4C). We finally determine that there were 4 hub genes in M saddlebrown module, 8 hub genes in M greenyellow module, and 7 hub miRNAs in M salmon module (Table 1 and Supplementary Figure 5A-5C).

The verification and functional analysis of hub genes
Since M saddlebrown and M greenyellow modules were negatively correlated with hypertension status, we wondered if the hub genes in these 2 modules were also negatively correlated with hypertension. We verified this by correlating hub genes expression values with hypertension status. There was a significant difference in hub genes (miRNAs) expression levels between normal and hypertensive group Supplementary Figure 6A-6C. We detected the expression of 12 hub genes by DEGs on original dataset GSE75360. Figure 7 shows the DEGs in a volcano map. The hierarchical clustering heat map of DEGs is shown in Supplementary Figure 7. We used the jvenn tool to overlap the hub genes in the M saddlebrown and M greenyellow module and DEGs, respectively, and found 11 hub genes were both in the DEGs gene list and the hub genes lists (Supplementary Figure 8A, 8B). We used the same method to verify the hub gene in another dataset GSE117261 from the GEO database, and found TBXAS1, FCER1G, and IGSF6 were in the DEGs gene list and hub genes lists (Supplementary Figure 9). We compared the expression status of

DATABASE ANALYSIS
these three hub genes in normal and hypertensive patients from the GSE117261 dataset, and the results were consistent with those from the GSE75360 dataset (Supplementary Figure 10).

Analysis of miRNA-gene interaction networks
We discovered that the HIF-1 signaling pathway and the insulin signaling pathway were also found in the miRNA networks. In order to better understand the regulatory relationship between genes and miRNAs, the miRNA-gene interaction network was constructed, which was based on genes and miRNAs involving the same pathway. There was a total of 46 nodes (21 genes and 25 miRNAs) and 112 pairs of interactions in the miRNAgene interaction network ( Figure 8). We used the MCC algorithm in cytoHubba plugin to screen the top 10 miRNA-genes in the network. Four miRNAs (hsa-miR-548am-3p, hsa-miR-513c-3p, hsa-miR-182-5p, and hsa-miR-548d-5p) and 6 genes (IGF1R, GSK3B, FOXO1, PRKAR2B, HIF1A, and PIK3R1) were

DATABASE ANALYSIS
the core nodes of the network (Table 2, Figure 9), and hsa-miR-548am-3p was considered as the core regulator because it targeted these 6 genes.

Discussion
In this study, mRNA and miRNA co-expression networks of hypertensive patients samples were constructed by using the WGCNA method. Out of the 19 identified mRNA modules, M saddlebrown and M greenyellow had the most significant correlations to hypertension. For miRNAs, we identified 14 modules, of which M salmon was the most significantly associated with hypertension. Finally, we identified 4 hub genes in the M saddlebrown module, 8 hub genes in the M greenyellow module, and 7 hub miRNAs in the M salmon module that were correlate with hypertension. Four miRNAs and 6 genes were also associated with the genetic susceptibility to hypertension. Our findings help us better understand the pathogenesis of hypertension, which in turn will provide us candidate biomarkers for clinical decision-making, potential therapeutic targets for accurate diagnosis, and treatment targets of hypertension.
In the M saddlebrown module, the GO analysis indicates that M saddlebrown mainly refered to membrane formation, translation, and ribosome and rRNA processing. For the M greenyellow module,    the insulin signaling pathway and the HIF-1 signaling pathway were enriched. The HIF-1 signaling pathway is reported to be important for development of pulmonary hypertension in chronic hypoxia [23]. The HIF-1 pathway is related to proliferation of pulmonary arterial smooth muscle cells (PASMCs), which is also a central pathological component for a kind of hypertension [24]. We found 12 hub genes associated with hypertension (RPL21, RPS28, LOC442727/PTGAP10, LOC100129599/ RPS29P14, TBXAS1, FCER1G, CFP, FURIN, PECAM1, IGSF6, NCF1C, and LOC285296/UNC93B3). Other studies have shown that some of these hub genes are related to hypertension, such as TBXAS1, FCER1G, FURIN, and PECAM1. The enzyme encoded by the TBXAS1 gene can catalyze many reactions involving cholesterol, steroids, drug metabolism and other lipid synthesis; and a study has reported that TBXAS1 is a potent vasoconstrictor [25]. FCER1G is the Fc fragment of FcepsilonRI (IgE) receptor Ig and IgE expression on the cell surface of human platelets. It has been reported that FCER1G has a regulatory effect on diabetic kidneys, and as hypertension is closely related to this disease, it can be inferred that it has a certain effect on hypertension [26]. FURIN encodes proteases, and Li et al. reported that FURIN may be a candidate gene for human hypertension [27]. PECAM1 encoded immunoglobulin are involved in processes of leukocyte migration and angiogenesis, and studies have shown that PECAM-1 expression is upregulated in pregnancy-induced hypertensive patients [28].
Some hypertension-related pathways were found in the M salmon module, which indicates that miRNA modules may also play a role in hypertension. The pathways include the adipocytokine signaling pathway, the thyroid hormone signaling pathway, the insulin signaling pathway, and the HIF-1 signaling pathway. The 3 key modules (M saddlebrown , M greenyellow and M salmon ) have significant relations to hypertension, so we speculate that there may be regulatory relationships of the miRNAs and genes in the same pathway that affect blood pressure. In order to find out their relationship, we constructed a miRNA-gene regulatory network (Figure 8). We found 10 core nodes of the network (hsa-miR-548am-3p, hsa-miR-513c-3p, hsa-miR-182-5p, hsa-miR-548d-5p, IGF1R, GSK3B, FOXO1, PRKAR2B, HIF1A, and PIK3R1), and hsa-miR-548am-3p was considered the core regulator. Although these 4 miRNAs have not been reported to regulate hypertension, their target genes have blood pressure regulation effect. It is found that IGF1R, the target gene of hsa-miR-548am-3p, hsa-miR-513c-3p, hsa-miR-182-5p, and hsa-miR-548d-5p, may be associated to the genetic susceptibility to hypertension [29]. Due to its function in lung development, IGF-1R may play a role in the pathogenesis of pulmonary hypertension [30]. GSK3B, a target gene of hsa-miR-548am-3p, hsa-miR-513c-3p and hsa-miR-182-5p, encodes for the enzyme glycogen synthase kinase 3 beta (GSK3B). Platelet-derived factor (PDGF) can affect the abnormal growth of pulmonary arterial smooth muscle cells (PASMC) in pulmonary hypertension
This study had several limitations. First, our sample size was small and may not fully represent hypertension patients. Second, there is no biological experimental verification of the hub genes (miRNAs) and relationship between genes and miRNA. In a follow-up study, the molecular verification experiment will be conducted to uncover the molecular level mechanisms of miRNA-gene interactions and their relations to hypertension. Verified molecular markers can be used as new diagnostic indexes of hypertension in the future.

Conclusions
This study established a WGCNA-based gene expression data process workflow, identified 2 mRNA modules and 1 miRNA module related to hypertension, and provided potential candidate biomarkers for hypertension treatment. Our analysis revealed novel miRNA-gene interactions as well as central miRNAs and genes that play critical roles in hypertension.

Conflict of interests
None. dataset. The red dots on the right corresponds to 2-fold up changes with P-value less than 0.05, and the blue dots on the left means 2-fold down changes with P-value less than 0.05.

Supplementary Data
Scale free topology model t, signed R 2

Scale independence
Soft threshold (power)     There is no significant difference in the interaction between the different modules, indicating a high degree of independence between these modules.