|Year : 2018 | Volume
| Issue : 7 | Page : 1644-1649
Personalized discovery of disrupted pathways and significant genes in preeclampsia based on accumulated normal tissue data
Ying Luo, Xiao-Chen Ma, Qing Gao, Lu-Quan Cao
Department of Prenatal Diagnosis, Jinan Maternity and Child Care Hospital, Jinan 250001, PR China
|Date of Web Publication||19-Dec-2018|
Department of Prenatal Diagnosis, Jinan Maternity and Child Care Hospital, No. 2 Jianguo Xiaojing San Road, Jinan 250001
Source of Support: None, Conflict of Interest: None
Purpose: This study was designed to identify disrupted pathways in an individual with preeclampsia (PE) using accumulated normal sample data based on individualized pathway aberrance score (iPAS) method.
Materials and Methods: Pathway data were obtained from the Reactome database. Next, the average Z algorithm was utilized to compute the iPAS. The disrupted pathways in a PE sample were identified by means of t test according to the pathway statistics values of normal and PE samples. In addition, we screened the differential expressed genes (DEGs) using SAMR package and constructed the differential co-expression network comprising DEGs. Subsequently, topological analysis for the co-expression network was conducted to identify hub genes.
Results: Under the threshold of false discovery rate <0.05, 69 disrupted pathways were selected. Among them, formation of tubulin-folding intermediates by containing t-complex polypeptide 1 (CCT)/TCP1 ring complex (TriC) was the most remarkable pathway. Degree analysis for co-expression network of DEGs suggested that there were several hub-disrupted pathway-related genes, for instance, TCP1 and TUBA1A. More importantly, these two hub genes were enriched in the most significant pathway of formation of tubulin-folding intermediates by CCT/TriC.
Conclusion: The iPAS method is suitable for identifying disrupted pathways in PE. Pathway of formation of tubulin folding intermediates by CCT/TriC might play important roles in PE.
Keywords: Differentially expressed genes, disrupted pathways, hub genes, individualized pathway aberrance score, preeclampsia
|How to cite this article:|
Luo Y, Ma XC, Gao Q, Cao LQ. Personalized discovery of disrupted pathways and significant genes in preeclampsia based on accumulated normal tissue data. J Can Res Ther 2018;14:1644-9
|How to cite this URL:|
Luo Y, Ma XC, Gao Q, Cao LQ. Personalized discovery of disrupted pathways and significant genes in preeclampsia based on accumulated normal tissue data. J Can Res Ther [serial online] 2018 [cited 2019 Dec 6];14:1644-9. Available from: http://www.cancerjournal.net/text.asp?2018/14/7/1644/203603
| > Introduction|| |
Preeclampsia (PE) is a disorder of pregnancy, characterized by high blood pressure and a large amount of protein in the urine. It is a major contributor to maternal and neonatal morbidity and mortality  and occurs in 2%–4.5% of all pregnancies. Moreover, PE can result in severe complications including seizures (eclampsia), multi-organ failure, stroke, and death. Because the pathogenesis of PE is not completely understood, effective strategies for preventing PE are lacking. Thus, a deeper understanding of the molecular mechanism is helpful in guiding prevention and treatment of PE.
As reported, in addition to obesity, prior hypertension and diabetes mellitus,, several genes and pathway alterations have been identified to predict response to therapeutic regimens of PE. For instance, it has been verified that imprinted gene PHLDA2 and MTHFR  are related to PE. Moreover, Delta-like ligand 1 and Notch 1 receptor have also been demonstrated to play important roles in PE. Zhu et al. have indicated that ErbB signaling pathway may be relevant to the epigenetic pathogenesis of PE. However, the potential mechanisms for PE development and progression are still not fully understood.
Currently, pathway analysis has become the first option to identify abnormal pathways to explain the potential biology of genes, due to the advantage of reducing complexity and increasing explanatory power. Deficiently, traditional pathway analyses are mainly focused on extracting disrupted pathways between normal and disease sample, but these traditional methods are not fit for identifying the pathway aberrance which might occur in an individual sample. Moreover, the traditional method of defining significant biological pathways is to analyze the expression of thousands of genes. Therefore, it is crucial to extract the pathway aberrance of individual sample based on the accumulated normal tissue data through analyzing the pathway level (individualized pathway aberrance score, iPAS), instead of analyzing the expression of genes in a brute-force manner.
Herein, in this article, a novel method was used to quantify the aberrance of pathways in an individual PE sample based on the comparison of one PE sample with many accumulated normal samples (reference meant the accumulated normal samples). Briefly, microarray profile of E-GEOD-25906 about PE and pathway data from Reactome database were obtained. Next, gene-level statistics, pathway-level statistics, and a significant test were successively implemented to screen out the altered pathways. Subsequently, co-expression network was constructed using Empirical Bayes (EB) method based on differentially expressed genes which were identified by means of significance analysis of microarrays (SAM) package. Afterward, topological analysis was conducted for co-expression network to identify hub genes and significant pathways.
| > Materials and Methods|| |
In this article, we downloaded the microarray profile of E-GEOD-25906 about PE from ArrayExpress database (http://www.ebi.ac.uk/arrayexpress/) to perform our analysis. E-GEOD-25906 included 23 PE and 37 control individuals based on the platform of A-MEXP-930 - Illumina Human-6 v2 Expression BeadChip. The gene probe IDs were converted into gene symbol level. Duplicated symbols in matrix were removed.
Herein, all human pathways were downloaded from Reactome pathway database (http://www.reactome.org/) which is an online-curated resource for human pathway data and provides infrastructure for computation across the biologic reaction network. Pathways having small number of genes are more easily understood by researchers. Thus, we removed the pathway with gene size >100. In addition, by intersecting with gene expression data, the pathways with the intersection of 0 between genes in pathway and original genes obtained in our analysis were also excluded. Finally, we obtained 1012 pathways (4648 genes).
Individualized analysis for pathways
Data pretreatment and gene-level statistics
In the present analysis, accumulated normal samples were utilized as a reference and the expression values of genes were calculated through comparing one PE sample with accumulated normal samples according to RMA package.
Specifically, the genes in the collected normal samples were normalized to the reference one by one, and average value as well as standard deviation (SD) of the expression level was counted. For individual PE samples, normalization was implemented using quantile package after combining the single PE microarray with all reference samples. Assuming that genes with multiple probes, gene expression level was calculated through averaging the expression level of probes. Gene-level statistics of each gene in a single PE sample was standardized according to the average and SD of reference.
The formula was defined as follows:
Where, mean (Nj) and SD (Nj) represented the average expression value and SD of expression value of j-th gene of the normal, respectively; Tij symbolized the expression value of j-th gene of a single PE sample; and Zij stood for the gene level statistics of each gene in a single PE sample.
Pathway-level statistics based on average Z method
Average Z method was applied to evaluate iPAS to examine an individual PE sample's pathway aberrance using the accumulated normal samples. In this experiment, we extracted the gene-level statistics value of all genes in one pathway and defined the average of the gene-level statistics value as the pathway statistics value. The formula was defined as follows:
Where, iPAS represented the expression status of a pathway, Zi was the gene level statistics of i-th gene, n denoted the number of genes in each pathway. Ultimately, we obtained the pathway statistics value of each pathway in an individual PE sample.
Significant analysis of the disrupted pathways
To analyze the alterations of the pathways in individual PE samples, Wilcoxon-test was applied to build the cluster heat map of pathways. Then, we employed t-test to measure the pathway statistics of each pathway in normal and PE samples. All the collected normal samples were compared with the reference one by one to generate the null distribution of pathway statistics. Afterward, P value was produced according to the comparison between the null distribution and a statistic of an individual PE sample. Subsequently, P values were adjusted for multiple hypothesis testing by means of false discovery rate (FDR) control by Benjamini and Hochberg. Disrupted pathways were identified when FDR was set as 0.05.
Validation of disrupted pathways
In an attempt to identify more significant genes and pathways and to further validate the disrupted pathways in PE samples, we performed differential expression analysis and co-expression network analysis.
Detection of differential expressed genes
At present, the common methods of screening of differential expressed genes (DEGs) mainly include limma, SAM, and siggenes. In this article, we utilized SAM package, where it correlated a large number of features (for example, genes) with an outcome variable, such as a group indicator, quantitative variable, or survival time  to screen the DEGs between PE patients and control placentas. Specifically, each gene was assigned to a score on the basis of the change in gene expression, relative to the SD of repeated measurements for this gene. Genes with scores greater than a liminal value were potentially significant. The percentage of falsely significant genes relative to the significant genes was termed as FDR. In an attempt to increase the stringency for significant difference in gene expression, delta value was computed using the function of SAMR. Compute.delta.table is a function, which is used to defined the delta value (a SAM parameter). DEGs in PE were extracted based on the delta cutoff value of 1.44. Taking the intersection of all genes in 1012 pathways with DEGs, we got the distribution of DEGs in pathways.
Construction of co-expression networks using Empirical Bayes method
Differential co-expression genes refer to gene clusters that have the same properties and work together in gene expression profiles under different conditions and could reflect the interaction relationship among genes under the different conditions. EB approach  is proven to be a useful complement to existing differential expression methods by simulations and case studies. In this section, we used EB method to achieve the co-expression relationship pairs between the DEGs. Furthermore, we constructed the co-expression network for these relationship pairs based on Cytoscape. Moreover, the genes in differential pathways obtained through iPAS approach were mapped to the co-expression network, and those mapped genes were extracted for subsequent analysis.
To explore the biological significance of nodes in the differential co-expression network, the topological centrality associated with the local scale (degree) was analyzed to describe the importance of nodes. Degree is the most evident index and determined as the number of interactions of a given node through totalizing the number of its adjacent genes. In our work, the nodes with degrees more than 18 were extracted as the hub genes. Significantly, the interactions among genes in altered pathways and the degree distribution of these genes were obtained. In our work, the hub genes were determined as the nodes with degrees not <18 in the co-expression network and, incidentally, were the top 20% genes of the highest connectivity degree. Then, we performed the comprehensive analysis on the hub genes and significant pathways.
| > Results|| |
Identification of disrupted pathways
In this article, we conducted the quantile normalization for individual PE samples to evaluate their gene-level statistics, making use of accumulated normal data. After the gene-level statistics of genes was transformed into the pathway-level statistics value of each pathway, we calculated the pathway-level statistics about each of 1012 identified pathways based on the mean value of gene-level statistic values of all genes enriched in this pathway. Based on FDR value <0.05, a total of 69 disrupted pathways were extracted. Then, we performed the cluster analysis to discover the alterations of these 69 pathways, and the heat map is shown in [Figure 1]. More significantly, these pathways were sorted in ascending order according to FDR value, and the top ten significant pathways were displayed in [Table 1]. The most significant disrupted pathway was the formation of tubulin-folding intermediates by containing t-complex polypeptide 1 (CCT)/TCP1 ring complex (TriC) (FDR = 3.2E-03).
|Figure 1: Cluster individual pathway aberrance score of preeclampsia dataset. Pathways (n = 69) and samples were clustered based on individualized pathway aberrance score. The colors represented the pathway statistics values. Horizontal axis was samples; vertical coordinate was disrupted pathways|
Click here to view
|Table 1: The top ten disrupted pathways based on the false discovery rate value <0.05|
Click here to view
Gene compositions of disrupted pathways
A pathway was made up of several genes which worked together to regulate biological functions or participate in cellular processes. With the goal of further revealing the properties and functions of disrupted pathways, the gene compositions of these disrupted pathways were studied on gene expressed levels. Taking the pathway of formation of tubulin-folding intermediates by CCT/TriC, for example, this pathway consisted of 17 genes (CCT6A, CCT8, TUBA4A, CCT5, TUBB6, TUBB2B, TCP1, TUBB3, TUBA3C, CCT2, CCT3, TUBB1, CCT7, TUBA1C, TUBA1A, TUBA1B, and CCT4). We calculated the gene-level statistics of these 17 genes in this pathway across PE and normal samples, as shown in [Figure 2]. From this figure, we discovered that the gene expression level in PE was disturbed comparing with that in normal sample. Moreover, in normal samples, gene levels for the 17 genes were similar. Of note, the change of the expression level of TUBA1C was remarkably obvious. Thus, we speculate that the different gene levels result in the generation of disrupted pathways in PE compared with the accumulated normal data.
|Figure 2: Expression pattern of genes in the pathway of formation of tubulin-folding intermediates by containing t-complex polypeptide/TCP1 ring complex. Each line stands for sample (red: preeclampsia; blue: normal)|
Click here to view
Differential co-expression network construction and topological analysis
With the delta = 1.44, a total of 414 DEGs were identified. Furthermore, taking the intersection of all genes in 1012 pathways with 414 DEGs, a total of 167 DEGs were distributed in the 424 different pathways. Based on the percentage of DEGs enriched in pathways >0.2, there were 9 pathways. Among these 9 pathways, the pathway of formation of tubulin-folding intermediates by CCT/TriC had the largest number of DEGs.
Moreover, all 414 DEGs were used to construct differential co-expression network comprising 405 nodes and 2493 co-expression interactions. After mapping the genes of the 69 disrupted pathways to the co-expression network, we found that a total of 30 genes in disrupted pathways were aligned to the co-expression network. Moreover, the subnetwork of the 30 genes is shown in [Figure 3]. Furthermore, we analyzed the degree distribution of the nodes in co-expression network. Based on the degree ≥18, 81 hub genes were extracted and 10 hub genes were disrupted pathway-enriched genes, as listed in [Table 2]. Afterward, we performed the comprehensive analysis on the ten genes and significant pathways. The pathway results of the top ten DEGs were shown in [Table 3]. We found that three pathways including formation of tubulin-folding intermediates by CCT/TriC, G-protein beta: gamma signaling, and postchaperonin tubulin-folding pathway were simultaneously identified based on iPAS method and comprehensive analysis. According to this result, we infer that iPAS is available to predict marker pathways for PE.
|Figure 3: The sub co-expression network that comprised the intersection of disrupted pathway-related genes and DEGs. Red nodes represented upregulated genes and green nodes were downregulated genes|
Click here to view
|Table 2: The hub genes enriched in disrupted pathways in the co-expression network|
Click here to view
Obviously, hub gene TCP1 was enriched in the pathway of formation of tubulin-folding intermediates by CCT/TriC and folding of actin by CCT/TriC. Moreover, another hub gene TUBA1A was enriched in the most significant pathway of formation of tubulin-folding intermediates by CCT/TriC. Significantly, these two hub genes TCP1 and TUBA1A were downregulated.
| > Discussion|| |
To obtain a deeper understanding of molecular mechanism of PE, a straightforward and novel pathway analysis approach was employed to screen out altered pathways in PE, which was in accordance with the comparison of a PE sample with the accumulated normal samples as “reference.” In the present paper, a total of 69 disrupted pathways were identified based on FDR <0.05. Among them, formation of tubulin-folding intermediates by CCT/TriC was the most remarkable pathway. Degree analysis for co-expression network of DEGs suggested that there were several hub-disrupted pathway-related genes, for instance, TCP1 and TUBA1A. More importantly, these two hub genes were all enriched in the most significant pathway of the formation of tubulin-folding intermediates by CCT/TriC.
The cytosolic chaperon in CCT is also known as the TRiC. Significantly, CCT/TRiC has been reported to mediate the folding of tubulins. Currently, few data are available about the clinical role of tubulin status in PE. Excitedly, a former study has indicated that the promoter of β-tubulin gene includes transcription factors which participate in the hypoxic response., Moreover, hypoxia-induced factors (for instance, HIF1α, and HIF2α) are needed for placental development and trophoblast differentiation. Interesting enough, hypoxia has been indicated to be a major cause of PE. In the current study, the pathway of formation of tubulin-folding intermediates by CCT/TriC was the most significant one. Altogether, these observations strongly suggest that the formation of tubulin-folding intermediates by CCT/TriC pathway might play keys role in PE development, partially through regulating the hypoxic response.
Furthermore, we obtained several hub genes from the disrupted pathways that were likely crucial to maintain functions and coherence of metabolic mechanisms. For example, hub genes TCP1 and TUBA1A all existed in the disrupted pathway of formation of tubulin-folding intermediates by CCT/TriC, which were downregulated. T-complex protein 1 subunit alpha is a protein that is encoded by the TCP1 gene in humans. Although TCP1 and TUBA1A have not been reported to be related with PE, a previous study has indicated that TCP-1 complex is a component for tubulin which is involved in cell cycle progression. Moreover, TUBA1A mutation has been reported to result in the defects of tubulin folding. Unek et al. have also indicated that placental alterations of PE might be connected to the cell cycle arrest. Moreover, the microinjection of an anti-CCTε antibody (that was shown to change the rate of and tubulin interactions with TCP upon translation in rabbit reticulocyte lysate) into cultured cells led to a delay in G1/S phase transition. Based on these, we infer that hub genes TCP1 and TUBA1A might be important in the development and progression of PE through the regulation of cell cycle. Our contributions provide the first clues for the molecular mechanisms of TCP1 and TUBA1A in PE.
In the current analysis, personalized-based pathway method was used to identify the significant pathways through creating the concept of the comparison of a single PE sample with accumulated normal samples. The clinical significance of our method is that it can diagnose a PE case in an individual patient. In light of our findings, we indicated that pathway of formation of tubulin-folding intermediates by CCT/TriC and hub genes (TCP1 and TUBA1A) might play significant roles in PE and they will be potentially preventive and prognostic markers in PE. However, we must take into consideration the limitations of our study. This analysis was performed based on the bioinformatics methods, but the findings have not been proved by animal experiments or clinical researches yet. Hence, further subtle experiments are warranted to explore the changes of the altered pathways in the understanding of PE pathology using animal experiments or PE patients' tissues.
Financial support and sponsorship
Conflicts of interest
There are no conflicts of interest.
| > References|| |
Eiland E, Nzerue C, Faulkner M. Preeclampsia 2012. J Pregnancy 2012;2012:586578.
Bijlenga D, Koopmans CM, Birnie E, Mol BW, van der Post JA, Bloemenkamp KW, et al.
Health-related quality of life after induction of labor versus expectant monitoring in gestational hypertension or preeclampsia at term. Hypertens Pregnancy 2011;30:260-74.
Palmsten K, Setoguchi S, Margulis AV, Patrick AR, Hernández-Díaz S. Elevated risk of preeclampsia in pregnant women with depression: Depression or antidepressants? Am J Epidemiol 2012;175:988-97.
Sibai B, Dekker G, Kupferminc M. Pre-eclampsia. Lancet 2005;365:785-99.
Al-Jameil N, Aziz Khan F, Fareed Khan M, Tabassum H. A brief overview of preeclampsia. J Clin Med Res 2014;6:1-7.
World Health Organization. WHO Recommendations for Prevention and Treatment of Pre-eclampsia and Eclampsia: Evidence Base. Geneva: World Health Organization; 2011.
Jin F, Qiao C, Luan N, Shang T. The expression of the imprinted gene pleckstrin homology-like domain family a member 2 in placental tissues of preeclampsia and its effects on the proliferation, migration and invasion of trophoblast cells JEG-3. Clin Exp Pharmacol Physiol 2015;42:1142-51.
Ge J, Wang J, Zhang F, Diao B, Song ZF, Shan LL, et al.
Correlation between MTHFR gene methylation and pre-eclampsia, and its clinical significance. Genet Mol Res 2015;14:8021-8.
Shimanuki Y, Mitomi H, Fukumura Y, Makino S, Itakura A, Yao T, et al.
Alteration of delta-like ligand 1 and notch 1 receptor in various placental disorders with special reference to early onset preeclampsia. Hum Pathol 2015;46:1129-37.
Zhu L, Lv R, Kong L, Cheng H, Lan F, Li X. Genome-wide mapping of 5mC and 5hmC identified differentially modified genomic regions in late-onset severe preeclampsia: Apilot study. PLoS One 2015;10:e0134119.
Glazko GV, Emmert-Streib F. Unite and conquer: Univariate and multivariate approaches for finding differentially expressed gene sets. Bioinformatics 2009;25:2348-54.
Ahn T, Lee E, Huh N, Park T. Personalized identification of altered pathways in cancer using accumulated normal tissue data. Bioinformatics 2014;30:i422-9.
Tsai S, Hardison NE, James AH, Motsinger-Reif AA, Bischoff SR, Thames BH, et al.
Transcriptional profiling of human placentas from pregnancies complicated by preeclampsia reveals disregulation of sialic acid acetylesterase and immune signalling pathways. Placenta 2011;32:175-82.
Croft D, O'Kelly G, Wu G, Haw R, Gillespie M, Matthews L, et al.
Reactome: A database of reactions, pathways and biological processes. Nucleic Acids Res 2011;39:D691-7.
Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, et al.
Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics 2003;4:249-64.
Gehan EA. A generalized wilcoxon test for comparing arbitrarily singly-censored samples. Biometrika 1965;52:203-23.
Green GH, Diggle PJ. On the operational characteristics of the Benjamini and Hochberg false discovery rate procedure. Stat Appl Genet Mol Biol 2007;6:1-23.
Dawson JA, Kendziorski C. An empirical Bayesian approach for identifying differential coexpression in high-throughput experiments. Biometrics 2012;68:455-65.
Smoot ME, Ono K, Ruscheinski J, Wang PL, Ideker T. Cytoscape 2.8: New features for data integration and network visualization. Bioinformatics 2011;27:431-2.
Haythornthwaite C. Social network analysis: An approach and technique for the study of information exchange. Libr Inf Sci Res 1996;18:323-42.
Gao Y, Thomas JO, Chow RL, Lee GH, Cowan NJ. A cytoplasmic chaperonin that catalyzes beta-actin folding. Cell 1992;69:1043-50.
Yaffe MB, Farr GW, Miklos D, Horwich AL, Sternlicht ML, Sternlicht H. TCP1 complex is a molecular chaperone in tubulin biogenesis. Nature 1992;358:245-8.
Kaluz S, Kaluzová M, Stanbridge EJ. Expression of the hypoxia marker carbonic anhydrase IX is critically dependent on SP1 activity. Identification of a novel type of hypoxia-responsive enhancer. Cancer Res 2003;63:917-22.
Wagner KD, Wagner N, Wellmann S, Schley G, Bondke A, Theres H, et al.
Oxygen-regulated expression of the Wilms' tumor suppressor Wt1 involves hypoxia-inducible factor-1 (HIF-1). FASEB J2003;17:1364-6.
Cowden Dahl KD, Fryer BH, Mack FA, Compernolle V, Maltepe E, Adelman DM, et al.
Hypoxia-inducible factors 1α and 2α regulate trophoblast differentiation. Mol Cell Biol 2005;25:10479-91.
Redman CW, Sargent IL. Latest advances in understanding preeclampsia. Science 2005;308:1592-4.
Willison K, Kelly A, Dudley K, Goodfellow P, Spurr N, Groves V, et al.
The human homologue of the mouse t-complex gene, TCP1, is located on chromosome 6 but is not near the HLA region. EMBO J 1987;6:1967-74.
Tian G, Jaglin XH, Keays DA, Francis F, Chelly J, Cowan NJ. Disease-associated mutations in TUBA1A result in a spectrum of defects in the tubulin folding and heterodimer assembly pathway. Hum Mol Genet 2010;19:3599-613.
Unek G, Ozmen A, Mendilcioglu I, Simsek M, Korgun ET. The expression of cell cycle related proteins PCNA, Ki67, p27 and p57 in normal and preeclamptic human placentas. Tissue Cell 2014;46:198-205.
Grantham J, Brackley KI, Willison KR. Substantial CCT activity is required for cell cycle progression and cytoskeletal organization in mammalian cells. Exp Cell Res 2006;312:2309-24.
[Figure 1], [Figure 2], [Figure 3]
[Table 1], [Table 2], [Table 3]