Identification of highly synchronized subnetworks from gene expression data

Simulation study

Simulation utilized the sample expression data gal80R given in Cytoscape (http://cytoscape.org/). There are 331 genes and 361 interactions in this network. Within it, we randomly selected subnetworks at three different sizes n (n = 40, 60, 80), as condition-responsive. In each responsive subnetwork m% (80%, 90%, 100%) of genes are defined to be active. The significance values of active genes were assigned randomly with top $n\times m\%$ significance values in gal80R, and that of the other genes were randomly sampled from the rest of the significance values. The phase locking index λ (see 2.3) of the interactions in the predefined responsive subnetwork were sampled from $N\left(0.8,0.5\right)$, i.e. a normal distribution with μ = 0.8, σ = 0.5; while λ for the remaining edges were sampled from $N\left(0.4,0.3\right)$. The choice of these values was based on the distribution of the λ values of gene pairs in protein complexes and of randomly selected gene pairs. For protein complexes we used the MIPS annotation (http://mips.helmholtz-muenchen.de/genre/proj/yeast) edited by Gerstein Lab (http://www.gersteinlab.org/proj/bottleneck/mips.txt).

We used the average sensitivity, specificity and F score to measure the performance of TopoPL. The performance is also evaluated with Receiver Operating Characteristic (ROC) curve, a plot of the true positive rate against the false positive rate [11].

Gene expression and protein-protein interaction data

Gene expression data was downloaded from EMBL's Huber group (http://www.ebi.ac.uk/huber-srv/scercycle/). It is a time course study of yeast cell cycle, where cells were arrested using alpha factor or cdc28. The alpha factor dataset contains 41 time points and the cdc28 dataset contains 44 time points, both at 5-minute resolution. These datasets provide strand-specific profiles of temporal expression during the mitotic cell cycle of S. cerevisiae, monitored for more than three complete cell divisions [14]. Yeast PPI data were downloaded from BioGRID (thebiogrid.org, version 3.1.69).

Phase locking analysis

In a perfect locking $\lambda =\left|\mathsf{\text{exp}}\left(i{\mathrm{\Psi }}_0\right)\right|=1$, and $\lambda \to 0$ when $\mathrm{\Psi }\left(t\right)$ is randomly distributed. λ offers a new measure to infer potential interaction between gene pairs [11].

TopoPL

$z_A^{TopoPL}$ captures the dynamic topological property of the subnetwork, and hub genes (genes with high network degree) contribute more to this metric. $|z_i{|}^{0.5}*a_{ij}*|z_j{|}^{0.5}*sgn\left(z_i+z_j\right),i\ne j$ can be regarded as the "activity measurement" of the interaction. Gene pairs with significant and synchronized expression changes, and whose gene products interact, contribute more to the activity of the subnetwork.

This metric is an improved version over the PCI that we previously proposed to identify active pathways from gene expression data [13]: $PCI={\sum }_{i,j}|x_{is}{|}^{0.5}*a_{ij}*|x_{js}{|}^{0.5}*sgn\left(x_{is}+x_{js}\right)$, where $x_{is}$ is normalized log expression measurement of gene i in sample s, and $(a_{ij})$ is the adjacency matrix of the PPI network of genes in the pathway. The merit of PCI has been demonstrated in previous works [13]. To reduce the potential impact on the network measure from residual inter-sample and inter-array biases after normalization, here we adopted the non-parametric measure $z_i$ in place of $x_{is}$. A similar metric to Eq. (5) was developed recently by us to predict candidate disease genes for type 1 diabetes, where $z_i$ is the z-score of disease relevance of gene i. There again we demonstrated the advantage of incorporating network structural information [16].

We implemented the searching procedure based on simulated annealing. The pseudocode of the algorithm is described below:

Input: the entire network $G_0=\left(V,E\right)$; a set of parameters for running simulated annealing: start temperature $T_{start}$ (= 1 in this study), end temperature $T_{end}$ (= 1e-8 in this study), number of iterations $N$.

Output: the subnetwork with the highest score.

Steps: initialize each node with its expression significance score $z_i$ and each edge with its phase locking index; select the largest connected component (subnetwork) $G_{out}$ from top 10% significant nodes of $G_0$; calculate score of $G_{out}$ and obtain its score $z_{out}^{TopoPL}$; then run the following:

For i = 1 to N, Do

Calculate the current temperature $T_i=T_i*0.8^{1/N}$; $G_{try}\leftarrow G_{ou{t}^{\prime }}$

Exit loop if $T_i<T_{end}$

Randomly pick a node $n\in V$

IF ($n\in G_{try}$), remove n from $G_{try}$;

ELSE add n to $G_{try}$;

Calculate score $z_{try}^{TopoPL}$for the largest connected component of $G_{try}$;

Calculate $\mathrm{\Delta }=z_{try}^{TopoPL}-z_{out}^{TopoPL}$;

IF Δ> 0, then$G_{out}\leftarrow G_{try}$;

ELSE, accept $G_{out}\leftarrow G_{try}$ with the probability$p=e^{\mathrm{\Delta }/T_i}$;

END

These steps can be iterated to identify subnetworks with the next highest scores and so on.

Simulation study

Using the simulated yeast gene expression data, we compared TopoPL with two other methods: (1) Additive scoring method (see definition Eq. (7) in Methods); and (2) TAPPA (see definition Eq. (8) in Methods) [13]. Additive does not use any structural information of the network, TAPPA uses only predefined static network structure ignoring the dynamic, condition-specific changes in interaction patterns. Figure 1 summarizes the average sensitivity, precision and F score from all simulated data: 10 replicates each of three network sizes (n = 40, 60, 80), at three states of activity (m = 80%, 90%, 100%). Though the three methods have similar sensitivity, the precision of TopoPL is higher. F scores showed that TopoPL performs better than TAPPA and Additive. The ROC curves also indicate that TopoPL performs better than the other two approaches, with the highest Area Under Curve (AUC), as shown in Figure 2.

Yeast cell cycle data

Presently, there is no "gold standard" to evaluate the biological relevance of network modeling algorithms. Here we investigated the functional enrichment of the proteins in the identified subnetworks [9], and compared to that obtained using Additive and TAPPA. The p values (Bonferroni corrected) of the top 2 terms are 3.33E-13 and 6.5E-12 with TAPPA, and 3.05E-8 and 3.13E-8, with Additive, respectively. TAPPA's are slightly larger than TopoPL, but Additive gave much larger p values. This indicates that including interaction structure, especially its dynamics, improves the sensitivity at identifying biologically relevant gene subnetworks.

The top 30 high-degree and high-betweenness nodes from the identified subnetwork and their interactions are presented in Figure 3. We hypothesize that they constitute a relevance core to yeast cell cycle, and provide a holistic picture of the primary molecular basis of cell cycle. In the core there are 18 genes annotated with GO:0007049 cell cycle (round rectangles), this rate (18 out of 39) is higher than that of the whole identified subnetwork (128 out of 524, a 1.9 fold enhancement, p = 0.11), and that of all genes in yeast (612 out of 5286, p = 0.00013). These results suggest that degree and betweenness can be utilized to further improve the performance of functional gene module identification.

We investigated the distribution of the phase locking index within the identified subnetwork. Clearly on average there is a higher degree of phase locking in it than in the whole PPI network (Figure 4). Interestingly the synchronization in the core is even higher, indicating that these core genes may work more closely in a coordinated fashion than others in the identified subnetwork.

Highly synchronized protein complex

We further examined the highly synchronized regions in the network core. Figure 5 shows the top 20 most synchronized interactions (corresponding ~1% of interactions in the identified subnetworks), MAK21 (NOC2) is at the center of this region. MAK21 is involved in preribosome export from the nucleus to the cytoplasm. Though it is not annotated with cell cycle GO term, but its homologue, SWA2 likely plays a role in ribosome biogenesis that is essential for the coordinated mitotic progression [23].

In protein complexes, the core components, which consist of two or more proteins that are present in most complex isoforms, are often regarded as functional units as they show surprisingly high degree of functional, essentiality, and localization homogeneity [24, 25]. We therefore also surveyed protein complexes and core components in the identified subnetwork. We found that all core components in complex 56 are in our core subnetwork, and they are shown in Figure 6. Interestingly all six genes show extremely high synchronization (0.976±0.006, see Figure 4). Their expression profiles are given in Figure 7. We also included their expression profiles in the cdc28 dataset; again high synchronization in expression is evident. This means that they are coordinated to work closely during cell cycle. This is not surprising as a large percentage of protein pairs within the core subnetwork were coexpressed at the same time during cell cycle [24]. Our algorithm is naturally good at finding highly synchronized genes pairs, therefore tends to include more core components from the same complexes.

Interestingly all six genes are annotated with GO:0042254 (ribosomal chaperone activity), it is defined as "A cellular process that results in the biosynthesis of constituent macromolecules, assembly, and arrangement of constituent parts of ribosome subunits; includes transport to the sites of protein synthesis".

Transcription factor binding motif analysis

FKH1 and MCM1 are well studied cell cycle related transcriptional factors [27]. TOD6 (Pbf1) and DOT6 (Pbf2) as PAC-binding factors, important in the regulation of ribosome biogenesis. Existing ChIP-chip studies suggest that genes have the highest occupancy by TOD6 and DOT6 are highly enriched for the GO Biological Process ''ribosome biogenesis" [28].

Agreement between the datasets

A good algorithm should be efficient at uncovering the true biology underlying different datasets, which should be consistent. In this study, we identified 484 genes with the cdc28 dataset, and 524 genes with the alpha factor dataset. There are 156 (~31%) overlapping genes in them (p < 0.00001, Fisher Test). In contrast, there are only 87 (~17%) overlapping genes with the Additive method (alpha: 501 genes; cdc28: 509 genes), and 145 (~29%) with TAPPA (alpha: 499 genes; cdc28: 503 genes). This indicates that incorporating network structural and dynamic information can generate robust results.