{'fraction_non_alphanumeric': 0.039480460033186476, 'fraction_numerical': 0.027178577559077644, 'mean_word_length': 5.15639309616062, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 2, 'https://': 5, 'lorem ipsum': 0, 'www.': 4, '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 and validation of key non-coding RNAs and mRNAs using co-expression network analysis in pre-eclampsia

Abstract Pre-eclampsia (PE) is a common complication of pregnancy, associated with maternal and fetal morbidity and mortality. In this study, we aimed to explore important long non-coding RNAs (lncRNAs) and their possible mechanisms in PE. GSE60438 expression profile including 25 PE samples and 23 normal samples were obtained from gene expression omnibus (GEO) database. After normalization with betaqn package in R, differentially expressed lncRNAs (DElncRNAs) and mRNAs (DEmRNAs) were identified using the limma package. Gene Ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway were analyzed using DAVID 6.7 and GSEA 3.0. LncRNAs-mRNAs coexpression was implemented using weighted gene co-expression network analysis (WGCNA). MicroRNAs linked with these DElncRNAs and DEmRNAs were predicted and a competitive endogenous RNA (ceRNA) network was built. A total of 53 DElncRNAs and 301 DEmRNAs were identified between control and PE samples. These DEmRNAs were enriched into pathways such as protein digestion and absorption, osteoclast differentiation. WGCNA constructed a lncRNA-mRNA coexpression network, among which SUMO1P3, NACAP1, NCF1C, ANXA2P1, GTF2IP1, NAPSB, OR7E37P were hub genes. ceRNA network was constructed together with microRNAs (miRNAs), and functional analysis indicated cellular membrane and sugar binding were involved in PE progression. Five lncRNAsANXA2P1, GTF2IP1, NACAP1, NCF1C and OR7E37P were successfully validated in our clinical specimens. The DElncRNAs, including ANXA2P1, GTF2IP1, NACAP1, NCF1C and OR7E37P might play important roles in PE. However, the exact mechanism of these lncRNAs in prediction and diagnosis of PE should be further explored.


Introduction
Preeclampsia (PE) is characterized by new onset hypertension with other maternal organ dysfunction and fetal growth restriction. [1] As a highly variable and heterogeneous syndrome, PE is often associated with new-onset pretension, thrombocytopenia, impaired liver function, renal insufficiency, pulmonary edema, and cerebral disturbances. [2] It is estimated that 2% to 8% of pregnancies are affected by PE, which is evaluated to cause over 5000 maternal deaths worldwide per year and is considered as the leading cause of maternal and morbidity and mortality. [3] Delivery of the fetus and placenta is the only known cure for PE. [4] The pathogenesis of PE is involved abnormal placentation and the development syndrome, but the exact underlying pathogenesis is largely unknown. [5] In recent years, long non-coding RNAs (lncRNAs) have been found play crucial roles in many diseases. Based on high-throughput techniques, several lncRNAs that might play important roles in PE have been identified in recent years. For example, He et al revealed that lncRNAs including LOC284100, LOC391533, and CEACAMP8 were aberrantly expressed in PE and might contribute to PE pathogenesis. [6] Long et al found lncRNA RP11-465L10.10 associated with the MMP9 gene was downregulated in PE patients. [7] Tong et al identified lncRNAs HK2P1, BNIP3P1, and PGK1P1 were core regulatory genes in PE. [8] In addition, microRNAs (miRNAs) have been reported as biomarkers in many diseases including cancer, cardiovascular disease, and diabetes mellitus. [9] Up-regulation of circulating miR-516-5p, miR-517-5p, miR-520a-5p, miR-525, and miR-526a is observed in PE. [10] Increasing evidence suggested that lncRNA might functions as miRNA spongers to regulate expression of target genes, which was termed as "competitive endogenous (ceRNA)". Although aberrantly expressed mRNAs, lncRNAs, and miRNAs have been identified extensively in PE, the lncRNA-miRNA-mRNA network-ceRNA network has been poorly reported.
Weighted Gene Co-expression Network Analysis (WGCNA) draws gene expression network, in which pairs of genes were identified based on the correlated expression across samples. [11] WGCNA analysis has been widely applied in recent years to identify core genes involved in disease occurrence and development. [12,13] In this study, we identified the differentially expressed genes (DEGs), including differentially expressed mRNAs (DEmRNAs) and differentially expressed lncRNAs (DElncR-NAs) between PE samples and normal samples. WGCNA was implemented on DEGs to construct coexpression network. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway analyses were performed for the DEmRNAs. Additionally, the hub lncRNAs were verified in clinical samples. Our research provided newly candidate genes that are important in PE and might have potential to serve as biomarkers in PE diagnosis.

Data source
The raw data of expression profile of GSE60438 based on platforms of GPL6884 and GPL10558 were obtained from the gene expression omnibus (GEO) database (https://www.ncbi. nlm.nih.gov/geo/). This dataset included genome-wide transcriptome sequencing data on 65 normotensive and 60 PE patients and was contributed by Yong et al. [14] The sequencing data based on the platform of GPL6884, which included 48 samples from 25 PE samples and 23 normal samples, were extracted for further analysis.

Identification of DEGs
The sequencing data were downloaded and genes were annotated into mRNAs, and lncRNAs based on GENCODE database (release 26). Data were normalized with betaqn package in R (version 1.16.0; https://www.rdocumentation.org/packages/wa teRmelon/versions/1.16.0/topics/betaqn-exprmethy450-meth ods). Further, differential gene expression analysis was implemented using the limma package (version 3.10.3) in R. The resulting P value was adjusted using Benjamini & Hochberg method. The mRNAs with adjusted P value <.05 and jlog fold change (FC)j > 0.5 and the lncRNAs with adjusted P value <.05 were regarded as differentially expressed.

Functional enrichment analysis
The GO was performed using DAVID version 6.7 (https://davidd.ncifcrf.gov/) and visualized by GO plot in R. [15] KEGG pathway enrichment was implemented using Gene set enrichment analysis version 3.0 (GSEA) with the adjusted P value <.05 (https://www.genome.jp/kegg/). The pathway with normalized enrichment score (NES) >0 represented the activated pathway while pathway with NES <0 was suppressed.

Weighted correlation network analysis
WGCNA (version 1.61) in R was used to construct coexpression network. The input genes were first filtered by removing those with median absolute deviation (MAD) <0.01 and among the later 25% MAD. The genes with missing value were also removed. Then, the Pearson correlation coefficient (PCC) between each pair of genes was calculated, then the adjacency function was defined, and module partitioning was performed based on the threshold of minimal module size of 30 genes. In addition, correlation between gene modules and clinical phenotypes including age, gestational age, and infant weight was calculated and clinical phenotyperelated modules were identified.

Construction and functional annotation of ceRNA network
The PCCs between DElncRNAs and differentially expressed mRNAs were calculated using psych package (version v1.8.12; https://www.rdocumentation.org/packages/psych/versions/ 1.8.12) in R based oncorr. test method. Coexpressed pairs were screened based on criteria of adjusted P value <.05 and jrj ≥ 0.75 and coexpression network was visualized using Cytoscape (version 3.7.0).
MiRNAs that could target the coexpressed lncRNAs and mRNAs were predicted from starbase database (http://starbase. sysu.edu.cn/starbase2/) and mirwalk database (http://mirwalk. umm.uni-heidelberg.de/). LncRNA-mRNA pairs regulated by the same miRNAs were integrated into lncRNA-miRNA-mRNA network and visualized by Cytoscape. The function of the mRNAs in these RNAs network was analyzed for GO and KEGG pathway enrichment using DAVID with P value <.05. [16] 2.6. Validation of hub lncRNAs in the ceRNA network by qRT-PCR The decidua basalis was obtained from 5 PE patients and 5 normotensive maternal women who were delivered in the department of Obstetrics and Gynecology, Shanxi Bethune Table 1 The sequences of primers for qRT-PCR.

Gene
Direction Hospital. This study was approved by the ethic committee of Shanxi Dayi Hospital. Total RNA of decidua basalis was extracted using TRIzol reagent (Thermo Fisher). Qualified RNA was reverse-transcribed into cDNA using the Prime Script RT master mix (TaKaRa, Shiga, Japan). The amplification was performed on ViiA7 Realtime PCR System (Thermo Fisher) using Power SYBR Green PCR master (Thermo Fisher). The relative expression level of lncRNA was calculated using the 2 ÀDDCt method using GAPDH as internal control. The primers used are listed in Table 1.

Statistical analysis
Each reaction of qRT-PCR was repeated 3 times and data were represented as mean ± standard deviation. Comparisons between groups were conducted by student's t test in Graphpad Prism 5.0. P < .05 was regarded as significant.

DEGs identification in PE samples
A total of 14337 mRNAs and 11091 lncRNAs were included in the dataset of GSE60438. Among them, 96 upregulated mRNAs and 205 downregulated mRNAs as well as 26 upregulated ncRNAs and 27 downregulated ncRNAs were filtered out ( Fig. 1A and B). Hierarchical clustering analysis was performed and the top 10 up-regulated and down-regulated genes were illustrated in heatmap (Fig. 1C).

WGCNA analysis of DEGs
By integrating the clinical symptom and the gene expression data, the expression profile of 305 DEGs in 48 samples were obtained.
WGCNA analysis further filtered out 228 DEGs for further analysis. Most of these genes were significantly classified into 3 modules, including 131 genes in ME turquoise module, 60 genes in ME blue module, and 34 genes in ME brown module (Fig. 4A). Genes in the same module showed significant correlation, while revealed poor correlation other modules ( Fig. 4A and 4B). Among the 3 identified modules, genes enriched in ME turquosis module showed higher degree, compared with those in the ME blue and ME brown modules. These data indicated that genes in ME turquosis demonstrated more function than genes in the other modules (Fig. 4C). In this network, high degree genes including TMEM71, SELL, Rgr, RGS18, P2RY13, NALP12, MCEMP1, LST1, LRG1, and LILRA2 were filtered out. Clinical symptoms such as disease status, gestational age and infant weight showed great correlation with these 3 modules. However, maternal age did not reveal significant difference ( Table 2).

Construction and functional annotation of ceRNA network
To further analyze the genes in the 3 modules, paired genes with adjusted P value <.05 and jrj ≥ 0.75 were calculated and a lncRNA-mRNA coexpression network was constructed (Fig. 5). There were 74 nodes and 102 edges including 58 mRNAs and 16 ncRNAs. High connect nodes were identified such as SUMO1P3, NACAP1, NCF1C, ANXA2P1, GTF2IP1, NAPSB, and OR7E37P.  We screened ncRNA-mRNA pairs regulated by a common miRNA, which was predicted from starbase and miRwalk databases. These ncRNA-mRNA pairs together with miRNAs were integrated into a ceRNA network including 25 mRNAs, 7 ncRNAs and 43 miRNAs (Fig. 6). Functional analysis of the network revealed that genes in this network were involved in molecular function of carbohydrate binding and sugar binding and were mostly located on cellular membrane (Table 3).

Discussion
PE is defined as hypertension associated with one or more newonset conditions including proteinuria, renal insufficiency, liver involvement, neurological complications, haematological complications, or uteroplacental dysfunction. [5] Various pathways have been implicated in the pathogenesis. In this study, we found DEGs were enriched in oxygen transport, immune response, inflammatory response, cell adhesion, angiogenesis, and platelet degranulation function. These genes were involved in immune related pathways and many other signaling pathways. Next, WGCNA analysis was performed to explore the hub genes combined with clinical syndrome. Further analysis revealed that disease, gestational age and infant weight were closely associated with identified modules. We identified 7 hub lncRNAs: SUMO1P3, NACAP1, NCF1C, ANXA2P1, GTF2IP1, NAPSB and OR7E37P. Five of these hub lncRNAs, including ANXA2P1, GTF2IP1, NACAP1, NCF1C, and OR7E37P were successfully validated in our clinical specimens.
The placental and maternal dysfunctions contribute to the pathogenesis of PE. Several genetic, angiogenic and other pathways have been showed in PE. Trophoblast invasion impairment would lead to imbalance between pro-and antiangiogenic factors. [17,18] Endothelial cells count a lot in Figure 7. Validation of hub lncRNAs in the ceRNA network by qRT-PCR. Comparisons between groups were conducted by student's t test in Graphpad Prism 5.0. * * indicated P < .01, * * * indicated P < .001.  inflammatory response, which would up-regulated the cytokines and adhesion factors secretion in PE. [19,20] Additionally, placental oxygenation, redox and immune tolerance have also been involved in PE. [21,22] Our research is consistent with previous studies that DEGs enriched in oxygen transport, immune response, inflammatory response, cell adhesion and angiogenesis, platelet degranulation function. All these pathways have been implicated in PE. Notably, the "osteoclast differentiation" was found to be significantly suppressed in PE. During pregnancy, bone turnover increases significantly in order to meet the demands from the fetus. [23] Evidences suggested that bone metabolism is altered in hypertensive pregnancies compared with normotensive pregnancies and bone resorption is increased and bone formation is decreased in PE patients. [24,25] Vitoratos et al demonstrated that PE patients exhibit higher OPG (osteoprotegerin) levels, a key factor in inhibiting bone resorption, which might be compatible with lower bone turnover. [26] Consistent with these studies, our study supports that bone turnover in PE was suppressed in PE compared with control pregnancies. Various investigations have been implemented for possible biomarkers in PE. Angiogenic factors emerged as important biomarkers in PE. Imbalance of these factors located in the central position of pathogenesis of PE. The expression of PIGF, sFLT1, and sENG differs significantly between women with PE and normotensive pregnancies. The ratio of sFLT1 to PIGF and PIGF to sENG is characterized with good performance in prediction. [27][28][29][30] In addition, ICAM-1 and VCAM-1 in endothelial dysfunction, cell free fetal DNA marker HYP2, NK cells in immune reaction, ROS in oxidative stress and other biomarkers including ADAM-12, PAPP-A, and PP-13 have been implicated their clinical function in PE. [31][32][33][34] In our research, 7 hub genes: SUMO1P3, NACAP1, NCF1C, ANXA2P1, GTF2IP1, NAPSB, and OR7E37P were identified as the diagnostic biomarker in prediction of PE. All these genes including protein coding gene and non-coding gene have not been reported before.
Small ubiquitin-like modifier 1 pseudogene 3 (SUMO1P3) is a pseudogene-expressed lncRNA which served as the oncogenic lncRNA in many kinds of human malignancy. Differently expressed SUMO1P3 has been reported in gastric cancer, nonsmall cell lung cancer, bladder cancer and hepatocellular carcinoma, which revealed its diagnostic function. [35][36][37][38][39] Annex in A2 pseudogenes (ANXA2P1) is differently expressed pseudogenes between glioblastomas and normal control tissues. [40] Deletion breakpoints mapping to GTF2IP1 is associated with Williams-Beuren syndrome. [41] Over-expressed NAPSB has been identified in carcinoma of the uterine cervix and this result indicated its correlation with CACX. [42] To our surprise, the function of NACAP1, NCF1C and OR7E37P have yet to be reported. All these genes were firstly reported in pathology of PE, thus providing the underlying clinically predictive significance. Other points we should pay attention to are all these identified genes were pseudogenes and these remind us that pseudogenes could be new biomarkers in PE progress.
There are some limitations in this study. First, though 5 of these hub lncRNAs, including ANXA2P1, GTF2IP1, NACAP1, NCF1C, and OR7E37P were successfully validated in our clinical specimens, the sample size is relatively small. Further studies involved in a large sample size are still warranted. Second, though the key lncRNAs involved in PE were identified, the exact mechanism of these lncRNAs involved in PE should be further studied.
In conclusion, we identified 355 DEGs between PE and normal samples. Functional analysis showed that immune pathway, inflammation pathway, and cell adhesion pathway were involved in development of PE. Further WGCNA analysis constructed ncRNA-mRNA network. Three modules were constructed and hub genes were filtered out. The hub ncRNAs combined with physical status might not only mediate the relationship between biological pathways, but also offer novel insights into the diagnosis and pathogenesis of PE. Combined with miRNAs, ceRNA network was constructed. Hub lncRNAs including SUMO1P3, NACAP1, NCF1C, ANXA2P1, GTF2IP1, NAPSB, and OR7E37P were identified and 5 lncRNAs ANXA2P1, GTF2IP1, NACAP1, NCF1C, and OR7E37P were successfully validated in our clinical specimens. However, further validation in a large sample size set and molecular mechanism of these key lncRNAs should be investigated in future.