Preeclampsia: a bioinformatics approach through protein-protein interaction networks analysis

Background In this study we explored preeclampsia through a bioinformatics approach. We create a comprehensive genes/proteins dataset by the analysis of both public proteomic data and text mining of public scientific literature. From this dataset the associated protein-protein interaction network has been obtained. Several indexes of centrality have been explored for hubs detection as well as the enrichment statistical analysis of metabolic pathway and disease. Results We confirmed the well known relationship between preeclampsia and cardiovascular diseases but also identified statistically significant relationships with respect to cancer and aging. Moreover, significant metabolic pathways such as apoptosis, cancer and cytokine-cytokine receptor interaction have also been identified by enrichment analysis. We obtained FLT1, VEGFA, FN1, F2 and PGF genes with the highest scores by hubs analysis; however, we also found other genes as PDIA3, LYN, SH2B2 and NDRG1 with high scores. Conclusions The applied methodology not only led to the identification of well known genes related to preeclampsia but also to propose new candidates poorly explored or completely unknown in the pathogenesis of preeclampsia, which eventually need to be validated experimentally. Moreover, new possible connections were detected between preeclampsia and other diseases that could open new areas of research. More must be done in this area to resolve the identification of unknown interactions of proteins/genes and also for a better integration of metabolic pathways and diseases.


Background
Preeclampsia is a pregnancy related disease associated with hypertension, proteinuria and increased maternal and perinatal morbidity and mortality, without known underlying mechanism and preventive treatment [1,2]. On the other hand, the future health or possible risks of women with previous history of preeclampsia are important areas of investigation. In this direction, it is well known the increased risk of future cardiovascular disease and renal dysfunction, however, other risks are also being discussed [1,[3][4][5]. Owing to the morbidity and mortality of this pregnancy related disease and the possible multifactorial causes involved [1][2][3][4][5], several experimental procedures have been applied by researchers in the last two decades, evidently, generating an elevated number of unprocessed information.
Although some bioinformatic analysis has been performed in particular microarray assays [6,7], an extensive data evaluation and processing has not yet been performed. Furthermore, the capabilities of bioinformatics tools for gene prioritization, network analysis, gene ontology and gene-disease relationships [8,9], together with all available data on protein/gene expression during preeclampsia bring an interesting and valuable opportunity for an in-silico study of the disease. Therefore, the present study is focused on two main areas: I) collection and basic analysis of the genes/proteins-diseases dataset, including, protein-protein interaction network and pathway enrichment analysis and II) exploration of the related gene-diseases in order to evaluate other genetic diseases possibly related with preeclampsia.

Protein-protein interaction network analysis
Preeclampsia PPI network topology reveals (Figure 1) a similar behavior with respect to general topology of PPI following a power law behavior [10] and therefore scalefree properties. These types of networks have the particular feature that some nodes are highly connected compared with others on the same network. These highly connected nodes (hubs) in general, represent important proteins/genes in biological terms and therefore are treated with special attention.
The top 50 genes with high scores and also present in the initial set (347) are shown in Table 1, however, other genes were found with high scores value but there are not part of the initial gene group. As expected some of the selected genes like FN1, FLT1, F2, VEGFA, PGF, TNF, NOS and INHBA, are well known preeclampsia relates genes (see discussion) and several of them are related with signaling pathways.
A total of 27 communities (k = 3) covering 161 genes were identified by communality analysis. In Figure 2 (Left) we represent those communities that are superimposed to form a large connected graph. The genes involved in communities overlapping are also highly represented in the Table 1 (and also the genes members of the large community). The model based clustering analysis reveal an optimal number of 8 clusters (BIC = 152192.4) with an ellipsoidal distribution with equal volume, shape and variable orientation. The genes are grouped in the clusters (C1. . .8) as follow: C1 (67), C2 (56), C3 (1806), C4 (59), C5 (133), C6 (23), C7 (95) and C8 (161). The C8 and C4 correspond with the highest mean scoring value: 393.3 and 348.9 respectively, and contain all the 100 genes with highest score values (part or not of the initial gene set). Furthermore, 161 genes of C8 are also the same genes detected in the communality analysis.
Gene ontology (GO) enrichment analyses were performed in all obtained clusters. However, for simplicity only C4 and C8 are presented (Figure 2 Right). The GO analysis reveals that C8 comprise several processes related with angiogenesis, apoptosis and cell proliferation and also shared with C4 several processess involved in cell activation and biological adhesion. The relation between these processes as well as the fact that both groups are representative of the highest scored genes could indicate a particular relevance of the clusters in terms of genes-disease relationship. On the other hand, also these processes are well known involved in preeclampsia and are also consistent with the pathway enrichment analysis.

Diseases and metabolic pathway enrichment analysis
Several types of diseases were found statistically significant in the enrichment analysis; partial results are presented in Table 2. Obviously, preeclampsia and even hypertension have to be present in the analysis. In the GAD database there are several disease classes and beside the presented in Table 2, others like hematologic (pvalue = 6.1E-06) and renal diseases (p-value = 2.5E-04) were also significant. Even when we present only the results obtained with GAD database, analysis was also performed with OMIN database confirming the ovarian cancer (p-value = 0.019) and also indicating colorectal cancer as statistically significant (p-value = 0.011), however, better results were obtained with GAD and also more consistent with literature information and pathway enrichment analysis.
It is important to consider that several genes in the PPI network do not present a known relation with specific diseases, at least reported in the GAD or OMIN databases. Only around 30% of the 2400 genes were found in the databases. This difficulty means that we have to be cautious with the preeclampsia genesdiseases relationships and with the reliability of the statistical p-value, even when some important and significant inferences, can be made.
A similar situation occurred with the pathway enrichment analysis (Table 3). Even when the KEGG database is the most representative of our gene space and high coherence was noticed with the physiopathology of PE, the results only cover around 50% of the initial 2400 genes. A similar procedure was also performed with the Reactome and BioCarta databases with a less covering (37% and 27% respectively) and showing a high coherence with Table 3 results. These databases reveal other significant pathways like NGF, PDGF, BMP, EPO and EGFR signaling as well as apoptosis and hemostasis pathways (data not shown). Some cancerous pathways (i.e. prostatic, pancreatic and lung) were also found statistically significant in KEGG but were excluded from Table 3, in order to simplify, because many of them have similar reactions with the general cancer pathway, already presented in the Table 3.
In order to simplify and enhance the understanding of the involved pathways and their relationship with the selected hubs, a fusion between both was made ( Figure 3). However, it is important to exalt that from the 50 hubs previously selected; only 22 present some significant pathway association with Table 2. The genes: NDRG1, LGALS3BP, BANF1, SGTA, TRIM29, RGS20, PLEC, GRN, ST13, AKAP5, FSTL3, DST, PKIA, QKI, MLF2 and KRT19, for example, were not found in the KEGG database.

PPI analysis
The top 50 genes selected were manually analyzed with the scientific literature in order to validate its connection with PE or even changes during pregnancy and we corroborate that several of them like FLT1, FLT4, VEGFA, PGF, TNF, FN1, F2 and NOS3 can be related to the main lines of research in preeclampsia pathogenesis hypothesis and specially angiogenesis [11][12][13][14][15][16]. Moreover, PGF, INHBA and related inhibin as well as activin proteins have been considered in several predictive model of PE [17,18]. The presence of those genes in our selection could validate the method applied and increase our confidence with respect to those genes that are poorly explored or unknown in their association with preeclampsia. In the latter group, we have for example the genes LYN, PDIA3, NDRG1 and TBK1. The PDIA3 genes encode an endoplasmic protein highly related with folding processes and in its relation with pregnancy only one article was found [19] revealing that PDIA3 gene intervene in trophoblast invasion via interleukin (IL) 11. Similarly, TBK1 could also intervene during the first moments of pregnancy to secure the implantation in relation with the nuclear factor kappa B [20]. However, there are not articles of PDIA3, TBK1 or LYN related with preeclampsia. Moreover, only one article was found describing an increased expression of NDRG1 during preeclampsia [21]. A similar problem is found with other genes such as IQGAP1, DNM1, SAT1, MEN1 and SH2B2 that also have been little explored during pregnancy.
Cluster analysis indicates that mainly C8 but also C4 probably embrace the most significant genes, at least related with centrality. All genes highest scored are part of these clusters as well as the large community graph and as we can notice the genes that lead to a community superposition are also highly scored in Table 1 ( Figure 2). On the other hand there are other genes like ENG, VEGFB and INHA that are well known related to preeclampsia and are also part of the C8 cluster [11][12][13][14][15][16][17][18]. It is important to notice that there are other genes not The table shown the first 50 genes selected according to the score values and arranged in descendent order. We also provide the gene symbol the associated Entrez Gene identifier. Notes: 1) Entrez gene ID.
shown in Table 1 because were not present in the initial gene set. In this group, EGFR and GRB2, are both with highest scores and there are both well related with preeclampsia [22,23]. Thus, even when the work presented focused on the analysis of the initial group, it is possible that relevant genes were identified by PPI topology but not included in the initial subset. Moreover, the GO enrichment analysis clearly reveals that C8 and C4 clusters are related with angiogenesis, apoptosis and biological adhesion, which are also biological processes obtained with the pathways enrichment analysis. Angiogenesis and related processes (i.e. vascular growing, cell differentiation and hypoxia) are considered as central processes related to preeclampsia [24] and therefore it could support the reliability of the procedure and also the necessity to increase the study of the C8 and C4 genes. Future combination of centrality indexes and specific clusters selection together with machine learning procedures or genetic algorithm optimization based on groups differentiability or external prediction subset, could reduce even more the final gene space and we are currently working in this direction.

Diseases and pathway analysis
We manually evaluate through scientific literature the connection of PE with the identified metabolic pathways and diseases. The relationship between preeclampsia and inflammatory, immunologic, angiogenic and hemodynamic responses are very well documented [22][23][24][25][26][27][28] and therefore are expected not only in the gene space but also in the metabolic pathways analysis. The metabolic pathway analysis (Table 3) indicates a strong significance of the cancer pathway that is consistent with the disease analysis (Table 2) and also with the VEGF signaling, that is present in several related processes like cytokinecytokine receptor, angiogenesis and also cancer pathway.  The table provides the p-value obtained from disease enrichment analysis according to the Disease Class and also the sub-classes diseases following the GAD database structure. The results showed are partials as discussed in the respective section. Figure 3 Genes-metabolic pathways interaction. The genes and pathways were selected after hubs detection (see Table 1) and enrichment analysis respectively (see Table 3). Progesterone-mediated oocyte maturation 2.50E-07 Renin-angiotensin system 1.80E-02 The table provides the p-value obtained from pathway enrichment analysis in the KEGG database. The pathways names correspond with the KEGG nomenclature.
The results showed are partials as discussed in the respective section.
With the procedure carried out in the present study a simplification of the metabolic pathways was possible, however, more needs to be done in this area in order to better integrate not only the hubs genes but also the comprehensive data created by the PPI interaction.
Considering the associated metabolic pathways we can notice that several signaling pathways are statistically significant and some of them are connected with PIK3CB, LYN and linked in several cases with TNF (Figure 3) that is a central gene affecting several processes and also widely studied in its relation with preeclampsia [28]. Similarly, SH2B2 that encodes an adapter protein of several tyrosine kinase receptors is also connected with metabolic pathways indirectly related with apoptosis. Contrasting PDIA3 and NDRG1 are not present in the Figure 3 connected with any metabolic pathway, but the recent relation found between PDIA3 and IL11 could open a relationship with cytokine-cytokine receptors [19], specifically through the hematopoietins receptor family. We can also notice ( Figure 3 and Table 3) that in the pathways associated with central proteins (the hubs), those highly connected are the cytokine-cytokine receptors, focal adhesion and apoptosis pathways and they contain almost the complete space of genes mainly studied in preeclampsia.
The elevated number of signaling pathways that we found statistically significant, together with the hubs distribution detected in Table 1, seems to point out the idea of a signaling related disease similar to cancer, however, the apoptosis and angiogenic mechanism had been related previously with PE and are also highly represented in our study.
On the other hand, the relationship between PE and cardiovascular diseases is well known. Women with previous history of PE or even pregnancy hypertension present an increased risk of future development of cardiovascular disease or hypertension [29][30][31]. This is clearly expressed in our results (Table 2). Moreover, our results also suggest a significant relationship between cancer and PE (at this point it is not possible to say if this correlation is positive or negative) indicated by diseases enrichment analysis and also by the metabolic pathways. However, the experimental and epidemiological evidence that could support the influence of PE in future cancer development is for now inconclusive, even when could be reasonably expected by the wide kind of genes that are actually shared between both diseases (i.e. angiogenesis related genes). Several articles report a reduced risk of future cancer development in women with previous history of PE, but others could not find any statistical relationship [32][33][34][35].
There are neither clinical signs in common between PE and Alzheimer and nor epidemiological studies. However, the connection between both diseases has been questioned before (at least in late-onset Alzheimer's disease) [36] and was significant in the present study. Therefore, although it is not possible to validate our findings experimentally with the current information, it actually opens new possibilities for future works. A similar problem concerns to other ageing related factors like atherosclerosis, arthritis or longevity (that for obvious reasons will be difficult to explore in PE during long-term studies).
The present study points out several advantages of the bioinformatics analysis but some limitations were also found. The detection of genes that are very well known to be involved in PE thought the applied methodology as well as the consistence of the results across metabolic information could support the novel candidates found, however, it is necessary to reduce even more the genes space applying other methodologies as well as to design new experimental experiences. On the other hand, the limitation of the human protein interaction information suggests that also orthologous genes should be needed in order to increase the PPI covering of the initial dataset and to increase the capabilities of the metabolic pathways and disease enrichment analysis.

Conclusions
In the present study we applied several bioinformatics tools in order to create a comprehensive initial database of genes statistically related to preeclampsia and a further expansion through the construction of related protein-protein interaction network. Using several centrality indexes for hubs detection, some well know preeclampsia related genes like FLT1, TNF, VEGFA and PGF were detected as well as other genes with high scores, like PDIA3, NDRG1, TBK1, LYN, IQGAP1, DNM1, SAT1, MEN1 and SH2B2 that have been poorly explored or unknown in the current state of the art of preeclampsia physiopathology.
Through disease enrichment analysis we corroborated the well know connection between PE and cardiovascular disease, but we also found a possible link of PE with aging and cancer. Moreover we also found the cancer pathway, focal adhesion, ECM-receptor interaction, apoptosis and cytokine-cytokine receptors interactions metabolic pathways highly represented by the PPI network that are in agreement with some of the preeclampsia-diseases relationship found, as well as with the central topics of the current preeclampsia investigations.

Dataset construction
Experimental proteomic/genomic data comparing normal (N) and preeclamptic (PE) pregnancies was obtained analysing the Gene Expression Omnibus (GEO) database [37]. The datasets considered are represented in Table 4.
Each experiment was analyzed independently in order to reduce the number of genes. In our case we considered an adjusted p-value ≤0.05 and a fold expression ≥2 as discussed elsewhere [6,7,[25][26][27]38,39]. Initially the p-value was obtained by a bootstrapping procedure with 1000 or 10000 iterations (depending on size of the sample) obtaining 645 statistically significant modulated genes, however, applying the false discovery rate (FDR) correction by the Benajmini-Hochberg method [40], this sample was reduced to 330 genes.
In addition, several text mining exploration tools were used to complement the GEO results. There are several tools to perform a text data mining analysis but several of them require extra information (i.e. chromosome region) instead phenotype or diseases notation (i.e. diseases name or related keywords). In our case we choose those methods that do not require previous genetic knowledge of the disease [8]. Moreover, the text mining procedures usually could provide several false positive associations and therefore those tools which also combine text-mining with other data sources in the analysis are preferred [8,41]. Considering these aspects, we used the following tools: PolySearch [42], Candid [43] and Phenopred [44]. Candid and PhenoPred use several heterogeneous data sources to overcome bias while in Poly-Search analyse was restricted to PubMed publications. Obviously many other algorithms could be used in alternative. In order to reduce the risk of including biased relationships, the top 10-20 genes/proteins with the highest scores were selected and individually analyzed considering the preeclampsia related scientific publications. Some of the top genes were also present in the previous dataset (GEO), therefore, the final dataset contained 347 genes.

Protein-protein interaction network (PPI)
The proteins associated with the previous 347 genes were identified and cross-referenced with the IRefIndex (v1.16) [45] and a signaling curated databases [46] that were used to create the protein-protein interaction (PPI) network. The IRefIndex database provide an index of protein interactions available in several databases like: BIND, BioGRID, DIP, HPRD that simplify the task consuming process of inter-database mapping and lead to a comprehensive covering of the available known protein interactions space. On the other hand, this PPI database is easily integrated in Cytoscape. Additionally, many diseases are related with signaling pathways modifications and therefore the inclusion of this interaction database considerable improve the PPI space. The interaction search was restricted to Homo Sapiens and includes all kind of experimental procedure as well as some predictive interactions (mainly from the OPHID database). The curation of the final database was performed both, manually and using home-made software to remove duplicate interactions and unify isoforms notation with unique genes. We obtained our final PPI network with 3279 interactions and 2400 nodes.
Some of the proteins present in our initial dataset had not any known experimental interaction (at least in humans) and therefore the 2400 nodes cover only 234 (67.45%) genes of the initial set (347). The network visualization and network topology indexes, calculated in the hubs detection process, were carried out using Cytoscape 2.8.2 and CytoHubba [47,48].
Several methodologies are available for hubs and essential genes identification, and all of them with the respective advantage and limitations [47,[49][50][51][52][53][54]. Some strategies are the use of genetic algorithm or machine learning procedures [49,50], however, the centrality approaches are by far the most applied procedures even by simplicity and because several studies had being pointed out its applicability [47,51,52]. Therefore, several centrality indexes were evaluated: Betweenness, bottleneck, density of maximum neighborhood component (DMNC), node degree, edge percolation component (EPC), eccentricity, maximal clique centrality (MCC), maximum neighborhood component (MNC), radiality and stress [47]. On the other hand, to obtain a scoring index we created the measurement (Score I) as follow: Where Ic i is the values of centrality indexes and i = 1. . .Nc, and is the number of calculated centrality indexes (Nc = 10). As we can note, score I is the sum of all the indexes percent value after individual normalization and therefore is restricted to a maximal value of 100×Nc which simplify even better the top genes selection. With the normalized centrality indexes we also performed a model-based clustering analysis using R-package [53] in The Table shown the GEO data sources, references and complementary information used in the study. Notes: 1) The 10 samples were collected for placenta biopsy during delivery. Five of the preeclamptic samples were collected before 31 weeks of gestation and the rest after this time.
2) The 23 samples were collected for placenta biopsy during delivery, however, because in GEO appear two platform types, both were analysed.
order to study hubs distribution with respect to centrality ranks. We also performed a communality (or cliques) network analysis by clique percolation method using CFinder [54]. The communality analysis provides a better topology description of the network including the location of highly connected sub-graphs (cliques) and/or overlapping modules that usually correspond with relevant biological information.

Pathway and diseases enrichment analysis
The pathways and diseases enrichment analysis were performed through the DAVID bioinformatics resource 6.7 [55], exploring the well know databases: KEGG, BioCarta and Reactome (pathways related) as well as OMIN and Genetic Association Database (GAD) (diseases related analysis). This online resource (DAVID) integrate, in a faster computational analysis, a wide range of enrichment analysis thought different databases providing also a substantial statistical description. The analysis was carried out considering the complete gene space of the PPI network. We also used DAVID in order to perform a gene ontology enrichment analysis in the obtained clusters.