Identification of dysfunctional modules and disease genes in congenital heart disease by a network-based approach

Identified dysfunctional modules

The work flow of our method of identifying dysfunctional modules is shown in Figure 1. To capture the information flow from causal genes to target genes and to identify dysfunctional modules from these causal paths, we first identified 85 target genes which are defined as those differentially expressed (DE) in sufficient proportion of patients (Additional File 1), and then connected each known causal genes of CHD with these target DE genes via shortest paths shown in Figure 1A and 1B. We model the protein interaction network as an electrical circuit where correlation coexpression of two end nodes of an interaction is used as the conductance of a resistor, and biological signals from disease genes propagate through PPI edges to responsive genes just like electrical current flows through resistors [11, 15]. Information from source genes will propagate its effect via protein-protein interactions, and DE genes which cover the majority of patients represent common dysregulated pathways in CHD. In the third step, we merged the paths from one causal gene to all target genes into a subnetwork, and computed the current flow for each gene to evaluate its importance in this local subnetwork as shown in Figure 1C. To measure the importance of one node in conducting electrical current, we computed the current flows through the node using the electronic laws [14], and defined the information flow score of the node as the sum of current through the node among all pair-wise combinations of the source node (the causal disease) and all target nodes (all responsive differentially expressed genes). Since the causal subnetwork can overlap, i.e. a gene can have several current flow scores, we assigned the gene to the subnetwork in which its current flow is maximum to derive mutually exclusive modules. The modules can also be thought of as an information-processing unit [5, 17], therefore we put the gene into the module where the largest amount of information passes through it. At last, reasoning that highly connected genes are more likely to coordinate and perform the same biological functions, we pruned the modules and iterated this process until each of them was connected as shown is Figure 1D (further details can be referred to Methods).

Totally, our method resulted in a connected CHD subnetwork consisting of 498 nodes and 2413 edges (Figure 2A). We identified 12 major disjoint modules by assigning genes to the group which maximizes the information flow scores while keeping each module inter-connected the CHD subnetwork (Figure 2B-E). The size of each module is shown in Table 1 and the list of all genes in each module is shown in Additional File 2.

Modules as disease-markers

To quantitatively investigate the relationships of the identified modules with disease phenotype, we computed mutual information (MI) of each module with sample phenotype. The results are shown in Table 1. We used the microarray data, which represents mRNA extracted from right ventricles (RV) of 16 children with Tetralogy of Fallot (TOF, the most common type of CHD) at the time of reconstructive surgery with 5 controls (NCBI GEO Accession ID GSE26125), to compute the activity vector of each module across all 21 samples. We used MI as the discriminative score to assess whether each module activity vector has significant relationship with phenotypes. Among 12 modules, 10 have significant mutual information (empirical p-value<0.05) compared with 1000 random gene sets of the same size (shown in Figure 3). Most modules are significantly correlated with disease phenotype, reflecting the synergistic differential expression within them.

To further evaluate the quality of our modules in terms of distinguishing disease from control samples, we computed the activity matrix of 12 modules (modules versus samples) and build logistic regression models [18] in a five-fold cross-validation approach. For comparison, we computed activity matrix of 12 KEGG pathways that are most enriched in the CHD subnetwork (Additional File 3). We also used 12 random gene sets of the same size as each module for control, and repeated this process 100 times. For experiments within the gene expression dataset of GSE26125, modules were first ranked by MI, and features were sequentially added to the classifier. Pathways were ranked by enrichment p-value, and were added sequentially. The result shows that modules are consistently better in the classification accuracy than pathways and random sets, while pathways are better than random sets only for the first five features (Figure 4A). To test the robustness of our module classifier, we also performed an independent cross-dataset experiment, in which classifiers were trained on GSE26125 and validated on GSE14970 (Figure 4B). In cross-dataset validation, both modules and pathways achieved similar high performance of classification accuracy, and clearly are much better than the random sets. Therefore, our modules are better classifiers of phenotype than pathways, and that their high classification performances are reproducible across independent microarray dataset. This provides evidence for the dysfunctional implication and importance of these identified modules.

Within module coexpression and cross module coordination

Besides comparing module expression profiles under disease and control status to evaluate its synergistic differential expression, another way to utilize microarray data is to exploit gene coexpression relationships within modules and across modules. Under the well-accepted hypothesis that genes exhibiting similar expression patterns across sample status are likely to have functional relevance [19], we reasoned whether genes within modules are significantly coexpressed, and whether modules have intense interactions in terms of interacting partners that are highly coexpressed. To analyze coexpression within modules, we computed the average Pearson's correlation coefficient (PCC) of all interacting proteins in each module, and compared it with that of 1000 random selected protein pairs of the same size (Figure 5A, Table 1). Result shows that 9 modules have significant coexpression (Mann-Whitney test, p-value<0.05), and all of them also have significant MI at the same time. If all 12 modules are regarded as a whole and compared with 12 random gene sets, the result is also significant (Mann-Whitney test, p-value = 2.96e-06).

To analyze intra-modular connectivity, we detected the interaction among disjoint modules based on weighted protein interactions (Figure 5B). Module interaction network is visualized in Figure 5C, where only the 22 significant interactions (corrected p-value<0.05) are shown. For each module, we defined its intra-modular connectivity as the sum of the absolute value of correlation for all its outgoing edges, and used this score to evaluate its network centrality (Table 1). Module interaction network not only reveals connectivity patterns between modules and their positions in CHD subnetwork, but also sheds light on the implications of various biological processes and pathways represented by each module. Therefore, besides topological meanings of this score, we would identify its biological underpinnings after examining the compositions of modules using annotated gene ontologies in the next subsection.

Functional enrichment analysis of modules reveals crucial biological processes of CHD

We have shown that modules with highest scores, being in the central place of the subnetwork, are enriched in core processes related to heart development, and that CHD genes participate in all the top ranked GO processes in a module. We then investigated the modules with lowest scores in the marginal place of subnetwork and identified their relations with CHD. Module 4 (Additional File 5) is ranked as the last, and is enriched in phosphorus metabolic process. Module 9 (Additional File 5) is enriched in cell communication, response to steroid hormone stimulus, lipid metabolic process and positive regulation of transcription. Module 5 (Additional File 5) is enriched in regulation of macromolecule biosynthetic process and nitrogen compound metabolic process. Although these GO processes seem to have less relationship with CHD, these auxiliary roles are important to facilitate key processes of heart development. Considering the fact that fetal heart developmental program involves intense transcription regulation, ligand-receptor interactions and signaling transduction, various macromolecules including hormones, cytokines and growth factors in the circulation or in the extracellular space of the heart, acting as ligands, can stimulate receptors in the cell membrane of cardiac cell [2], it is reasonable to assume that cell communication, transcription regulation, steroid hormone stimulus and macromolecules metabolism are more generally associated with CHD and can actually facilitate all related processes. Since they are not specific to CHD, they appear on the marginal place of the subnetwork. The fact that these modules can also contain causal disease genes demonstrates their relevance to CHD. For example, GATA4 and TBX5 are in Module 9, and both of them are important transcription factors in developmental processes. TBX5 protein associates with cardiac transcription factors including GATA4 and NKX2-5, and they activate many downstream cardiac effector targets such as sarcomeric proteins and vasoactive proteins. Various kinds of mutations in GATA4 and TBX5 can lead to various subtypes of CHD [25, 26].

Module-pathway crosstalk analysis reveals hub modules regulating key pathways of CHD

To analyze pathways represented by modules, and how they are coordinated by modules to account for the observed phenotypes, we went beyond identifying lists of pathways significantly enriched in each module, or simply counting the number of overlapping genes shared by the module and the pathway. As previously shown, there is a distinction between a pathway and a module. A pathway is a specific information-flow conduit, consisting of a series of molecular interactions while a module is an information-processing unit with self-contained cellular functions [17]. Therefore, modules can contain multiple pathways while pathways can be coordinated by various modules to allow inter-module connections. With this in mind, we implemented an analysis of module-pathway crosstalk similar to previous module interaction procedure. For a given module-pathway pair, we considered both common genes and weighted protein-protein interactions between them (See Methods for details of constructing module-pathway crosstalk network). Network view of module-pathway crosstalk is shown in Figure 6A, where circles represent modules, rectangles represent pathways, and edge width corresponds to strength of interaction between module and pathway. Heatmap of pathway-module crosstalk is shown in Figure 6B, where rows represent 85 pathways that are significantly (empirical p-value<0.01) influenced by at least one module, and columns represent modules and their influence to the 85 pathways. Module 10 is the most influential and many of its influenced pathways are closely related to heart progression, such as cardiac muscle contraction, Dorso-ventral axis formation, gap junction and regulation of actin cytoskeleton. Specifically, the gap junction is very important in cardiac muscle, through which the signal to contract is efficiently passed, allowing the heart muscle cells to contract in tandem. For example, GJA1 encodes a gap junction protein, and gene conversion in GJA1 has been found in patients with CHD [27]. In brief, these influenced pathways are consistent with major roles of Module 10 in biological processes like anatomical structure morphogenesis and cytoskeleton organization. Module 12 is the only one that can influence Notch signaling pathway, which is also consistent with its GO enrichment. Module 12 and module 2 are the only two modules that regulate ECM (extracellular matrix)-receptor interaction, a very important pathway in tissue and organ morphogenesis and in the maintenance of cell and tissue structure and function. The extracellular matrix protein, FLNA, for example, can cross-links actin filaments and participates in anchoring actin cytoskeleton to membrane proteins. Loss-of-functions mutations in FLNA are lethal to male and cause various CHD-related phenotypes in female [28, 29]. Although Module 12 and Module 2 do not affect as many pathways as Module 10 does, their dominant roles in regulating two essential pathways demonstrate their substantial importance during heart development. From Figure 6B, Module 7 also seems important because it interacts with many pathways, while we find that it has much less influence to pathways highly related to CHD. Compared with Module 10, it shows significant drop for many essential pathways related to CHD such as cardiac muscle contraction, Dorso-ventral axis formation, gap junction and regulation of actin cytoskeleton, while some pathways without direct relations to CHD such as Type I diabetes mellitus and Olfactory transduction, exhibit significant increase. Therefore, we reasoned that the importance of a module should be evaluated not by the number of pathway that it affected, but by its influence to important pathways implicated in CHD.

Prioritization of candidate disease genes by module analysis

HAND2 has significantly higher score than other candidate genes and ranks 1^st in our list. HAND2, and causal disease gene HAND1, whose somatic mutations have been reported to contribute to CHD [30], both belong to basic helix-loop-helix (bHLH) transcription factors, and are expressed in cardiac mesoderm during embryogenesis. HAND1 expression is limited in future left ventricle while HAND2 expression is limited in future right ventricle [31]. Targeted gene deletion of HAND2 in mouse embryos resulted in embryonic lethality from heart failure [32]. In a recent research, Shen et al. [33] screened for mutations in the HAND2 genes in 131 patients with various forms of CHD. Seven mutations in HAND2 were identified in 12 out of these patients. It is very likely that HAND2 is indeed a CHD disease gene. FOS is ranked 2^nd, and can form transcription factor complex AP-1 with proteins from the JUN family. FOS proteins are implicated as regulators of cell proliferation, differentiation and transformation. Literature text-mining using GeneCards Version 3 [34] reveals that FOS is involved in cardiac hypertrophy. NOTCH2 ranked 3^rd, and together with NOTCH1, whose autosomal dominant mutations can cause CHD [23], belong to single-pass transmembrane receptors that regulate many developmental pathways. Mutations in gene encoding NOTCH2 leads to Alagille syndrome [35], with cardiac phenotypes of peripheral pulmonary artery stenosis and septal defects. Besides genes with high score and well supported by literature search, our candidate gene list also contains genes with low scores. Currently, these genes have fewer, if any, favorable GO evidence similar to disease genes in the module, but some of them would be potential novel disease genes and worth further experimental research.

Datasets

GSE26125 consists of RNA extracted from right ventricle of 16 idiopathic TOF (Tetralogy of Fallot) patients and 5 controls obtained at the time of reconstructive surgery. The microarrays were CodeLink Human Whole Genome Bioarray, which contain > 54,000 probes. Raw data was preprocessed as described in the original paper [8, 9]. The microarrays of GSE14970 were Affymetrix Human Genome U133 Plus 2.0 Array. Raw data was preprocessed using RMA algorithm in R bioconductor [40]. Probe sets were mapped to NCBI entrez genes, and for genes with multiple probe sets, the average expression value of all corresponding probe sets was used. Together 8547 genes were expressed in the right ventricle tissues on GSE26125, which will be used for our study. 6818 genes were expressed in both GSE26125 and GSE14790, whose expression values will be used for cross-data validation. GSE14970 consists of 12 diseased samples (5 pre-operation acyanotic TOF and 7 pre-operation cyanotic TOF in the RV tissue) without controls but is the only alternative CHD-related dataset at this time, we added 5 control samples from GSE26125 and z-transformed the combined dataset.

Conditional-specific interaction network

We utilized protein-protein interaction data from several curated interaction databases in human, i.e., HPRD, BIND, BioGrid, IntAct, and MINT [41–45]. Those interactions contained in at least three of the five databases were selected. Only genes both in the PPI network and in the microarray were used in the following study, and we extracted the maximum connected component of the network, which consists of 4761 nodes and 18084 edges. To construct a context-specific protein interaction network, the weight of each interaction between protein u and v, as well as the conductance of each edge e(u,v), was defined as their correlation value, i.e. w(u,v) = |corr(u,v)|

Module identification

where I _kNei is the current between node k and its neighbor, ∑ _Nei is the sum of all neighbors connected to node k.

where P is a symmetric (N-1)×(N-1) conductance matrix, V is a (N-1)×1 vector of node voltages and I is a (N-1)×1 vector of currents passing through the nodes, the row and column of ground node is removed since its voltage is zero.

Such a linear system considered all interactions in the global network and required that every interaction can have a regulatory role for the expression of the target gene. We therefore implemented a straightforward heuristic approach: we identified the nodes that are on the maximum weighted shortest paths from a source gene to all target genes and performed the calculations of information flow score on this subnetwork. The modules were identified using the following algorithms.

The source code for module identification and related data in this study can be downloaded at http://doc.aporc.org/wiki/CHD.

Mutual information

Here, a' is obtained by discretizing a into [log2(# of samples)+1] = 6 equally spaced bins [13], and the midpoint of each bin is used as the value for a'. × and y enumerate values of a and c, p(x,y) is the joint probability density function (pdf) of a' and c, and p(x) and p(y) are the marginal pdf's of a and c.

where N is set as N = 1000.

Phenotype classification comparison

To evaluate the ability of our modules for discriminating phenotypes, we performed both within-dataset and cross-dataset classification validation. To allow for comparison with pathways and random controls, we aggregated the expression values of DE genes in each module and in each pathway from the CHD subnetwork. We also aggregated the expression of random gene sets with the same size as the module. The activity vector aggregation method is the same as in the former Mutual Information section. Modules were ordered in decreasing significance of MI, pathways were ordered in decreasing significance of enrichment and random gene sets were ordered the same sequence as modules. Then, logistic regression models were trained on the module activity matrix (modules versus samples), pathway activity matrix (pathways versus samples) and random gene sets (gene sets versus samples). For within-dataset experiments, each time 4/5 of the samples were implemented as the training set to build the classifier and the remaining 1/5 samples were used as testing set to evaluate the performance (five-fold cross validation). Assuming that there are K modules, then for each k ≤ K, select the first k_th modules to train the classifier. Pathways and random gene sets were also used in similar manner. The final classification performance was reported as the Area Under ROC Curve (AUC [47]) on the testing set using the classifier optimized from the validation set. For cross-dataset experiment, all 12 modules/pathways/random gene sets were trained on GSE26125 and validated on 12 disease samples from GSE14970 combined with 5 controls from GSE26125. All preprocessed microarray expressions were z-transformed before activity vector is aggregated. For both within- and cross-dataset experiments, random controls were repeated 100 times.

Correlation and crosstalk

where ∑ _{x ∈ M} ,_y∈P|weight(e(x,y))| corresponds to the sum of the weight of all edges, and {e(x,y):x∈M,y∈P} are those edges that connect a pathway P and a module M. $\sum _{x\in M\cap P}^{}1$ corresponds to the total number of common nodes shared by M and P, where each node is given the weight equals to 1. $\sqrt{\left|M\right|}\sqrt{\left|P\right|}$ corresponds to the square root product of the number of nodes in module M and pathway P. The summation is such divided in order to normalize the size of module and pathway.

where N₂ is set as N₂ = 1000.