Qualitative reasoning of dynamic gene regulatory interactions from gene expression data
© Chen et al. 2010
Published: 2 December 2010
Skip to main content
© Chen et al. 2010
Published: 2 December 2010
A gene regulatory relation often changes over time rather than being constant. But many gene regulatory networks available in databases or literatures are static in the sense that they are either snapshots of gene regulatory relations at a time point or union of successive gene regulations over time. Such static networks cannot represent temporal aspects of gene regulatory interactions such as the order of gene regulations or the pace of gene regulations.
We developed a new qualitative method for representing dynamic gene regulatory relations and algorithms for identifying dynamic gene regulations from the time-series gene expression data using two types of scores. The identified gene regulatory interactions and their temporal properties are visualized as a gene regulatory network. All the algorithms have been implemented in a program called GeneNetFinder (http://wilab.inha.ac.kr/genenetfinder/) and tested on several gene expression data.
The dynamic nature of dynamic gene regulatory interactions can be inferred and represented qualitatively without deriving a set of differential equations describing the interactions. The approach and the program developed in our study would be useful for identifying dynamic gene regulatory interactions from the large amount of gene expression data available and for analyzing the interactions.
Many mechanisms of biological processes are controlled by complex regulatory interactions between genes rather than by a single gene . Therefore, identifying the gene regulatory interactions is essential to improving our understanding of biological processes. A gene regulatory relation often changes over time rather than being constant. However, many gene regulatory networks available in databases or literatures are static in the sense that they are either snapshots of gene regulatory relations at a time point or union of successive gene regulations over time. Static gene regulatory networks are simpler and easier to construct and understand than dynamic networks, but temporal aspects of gene regulations such as the order of the gene regulatory interactions and the pace of the interactions are ignored in static networks.
A gene involved in regulatory interactions with others has at least one activator or inhibitor. An activator initiates the transcription of the gene, so high level expression of the gene is not possible without an activator . Thus, identifying genes and their activators or inhibitors is the key to constructing gene regulatory networks. Silvescu et al.  characterize the gene regulatory network in a Boolean model with multiple-time delays. But the Boolean model is restricted to logical relationships between variables. Probabilistic Boolean networks  and dynamic Bayesian networks  can reconstruct longitudinal regulatory networks from a set of mathematical equations if the equations precisely specify the networks, but fail when the underlying model is not correct .
In general dynamic relations are best represented by a system of differential equations, but differential equations are not typically used to represent dynamic gene regulatory relations. This is because dynamic gene regulatory interactions are not understood fully enough to derive differential equations despite the large amount of gene expression data available today. Even if differential equations are derived, they are often hard to solve. As shown later in this paper, we have developed a qualitative method for inferring dynamic gene regulatory interactions and visualizing them without deriving or solving a set of differential equations.
This paper presents a computational approach to uncovering gene regulatory relations and their temporal properties from a time-series gene expression data using a modified Pearson correlation coefficient and a new score scheme. For the temporal properties of gene regulatory relations, we infer the order of the gene regulatory interactions and the pace of the interactions. The identified gene regulatory interactions and their temporal aspects are stored in the regulation list and visualized as a gene regulatory network. All the algorithms have been implemented as a program called GeneNetFinder (http://wilab.inha.ac.kr/genenetfinder/) and tested on several gene expression data. The rest of this paper presents the algorithms and their experimental results.
In Equation 1, N is the total number of time points contained in the time span, X k and Y k are the expression levels of genes X and Y at time k, and and are the average gene expression levels at all time points of the time span. The R1 score is in the range of [-1, 1]. Among the total i × p candidate regulations, the regulation with the maximum absolute value of R 1(X, Y, i, p) is selected as the regulatory relation between genes X and Y. If the expression level of gene X increases before that of Y increases, X is a candidate activator of gene Y; if the expression level of gene X increases before that of Y decreases, X is a candidate inhibitor of Y.
In Equation 2, and are the average gene expression levels at all time points in the time span. X max is the maximum gene expression value of gene X. and are the maximum and minimum gene expression value of gene Y, respectively.
Both R 1 and R 2 scores are intended to represent a relation between a regulator gene X and its potential receiver gene Y. The regulatory relation between X and Y is not symmetric, so R 1(X, Y) ≠ R 1(Y, X) and R 2(X,Y) ≠ R 2(Y,X) . An interesting observation from the actual data is that two genes with R 2 score < 3 tend to have an inductive relation, whereas those with R 2 scores > 6 tend to have an inhibitory relation. For example, in the dataset of 30 yeast genes, 89.2% of the activations have R 2 scores < 3, and 91.4% of the inhibitions have R 2 scores > 6. In an extended dataset of 6,177 yeast genes, 80.4% of the activations have R 2 scores < 3, and 92.1% of the inhibitions have R 2 scores > 6. Hence, we consider a gene regulation with R 2 score < 3 as an activation, and that with score > 6 as an inhibition.
+A(t 1) → +B(t 2): up-regulation of A at time t 1 is followed by up-regulation of B at time t 2 (t 2 > t 1).
-A(t 1) → +B(t 2): down-regulation of A at time t 1 is followed by up-regulation of B at time t 2 (t 2 > t 1).
+A(t 1) - > -B(t 2): up-regulation of A at time t 1 is followed by down-regulation of B at time t 2 (t 2 > t 1).
-A(t 1) → -B(t 2): down-regulation of A at time t 1 is followed by down-regulation of B at time t 2 (t 2 > t 1) .
The regulatory relation of gene A with gene B is determined by the sign of the R 1 score of the genes. A relation with a positive R 1 score implies that gene A activates gene B whereas a regulation with a negative R 1 score implies that A inhibits B. The R 1 score of each gene regulation is iteratively calculated using Algorithm 1. For genes A and B, the regulation with the largest absolute R1 score is chosen for the regulation between the genes and represented as R 1(A, B, t 1, p). Algorithm 1 provides the top-level description of the method for inferring gene regulations and constructing a list with the regulations.
Construct a regulation list
Compute R 1(A, B, t 1, p) between gene A at time point t 1 and gene B at time point t 1 + p for all pairs of genes.
Select the regulation with the largest absolute value of R 1(A, B, t 1, p) .
If 0 < p < 6, classify the regulation into one of the four types, +A(t 1) → + B(t 1 +p) , - A(t 1) → + B(t 1 +p) , +A(t 1) → -B(t 1 +p) , -A(t 1) → - B(t 1 +p) , and add it to the regulation list.
If p = 0, two genes are co-expressed or co-inhibited, and such gene regulation is not added to the regulation list.
If the new gene regulation is already in the regulation list, merge it with the previous regulation.
Go to step 2 to find the next gene regulation until no more regulation found.
A time point with the minimum expression level of the regulator gene becomes a splitting point of the time span.
Each sub-timespan starts with at least 3 consecutive time points that have a positive slope of a curve representing gene expression levels, and ends with at least 3 consecutive time points with a negative slope.
Each sub-timespan encompasses at least 6 time points, including the start and ending time points.
For example, the time span CDC 15 of Figure 3A is the longest one in the gene regulatory relation +CLB1(T)-> +SWI5(T+1), and split into 3 sub-timespans. The first sub-timespan of CDC15 has less than 6 time points, so the R 2 score is not computed for the first sub-timespan. The R 2 scores for the second and third sub-timespans are 0.24 and 0.38, respectively. Figure 3B is another example of a gene regulatory relation +CLB6(T) -> -CLB1(T+1). The time span CDC15 is split into 4 sub-timespans, and the third sub-timespan has only 3 time points. So, the R 2 score is computed for the remaining three sub-timespans, which are 27.02, 20.34, and 9.87, respectively.
All gene regulations identified in the previous step are visualized as a 2-dimensional gene regulatory network, in which a node represents a gene. Edge types and edge labels of the network represent gene regulatory relations. Arrows represent inductive interactions (relations+A(t1) -> +B(t2) and-A(t1) -> +B(t2)) and blocked arrows represent inhibitory interactions (relations +A(t1) -> -B(t2) and -A(t1) -> -B(t2)). The regulator gene, type of regulation (+ for induction and - for inhibition), and time delay of the regulation are annotated as edge labels. Each edge is labeled with R/s/T to indicate a regulator gene R, sign s of the log-ratio of the expression level of R, and the time delay T of the regulation.
The layered layout algorithm, described in Algorithm 3, puts all nodes to horizontal layers. The node with the maximum degree is assigned to the top layer, and the nodes connected to the node are put in the next layer. If a layer has two nodes connected to each other, it makes a new layer above the layer and moves the node with a smaller degree to the new layer. The layered layout usually takes more time than the grid layout.
Find the node with the highest degree, and place it at the center grid point.
If there is a tie, select a node with a higher out-degree. Position all nodes connected to the center node in the adjacent grids. Nodes with a higher degree are positioned earlier than those with a lower degree in the east, north, west, south, northeast, northwest, southwest, and south east grid of the center node. If more than 8 nodes are connected to the center node, put the 9th node to the east of the 1st node, the 10th node to the north of the 2nd node, and so on.
Repeat step 2 for nodes connected to nodes that are already positioned. If an adjacent grid is occupied, move to the next possible position until it is available.
If there are disconnected nodes, repeat step 1 for the nodes and put them to the right of the previous subgraph.
Put the node with the maximum degree at layer 1. If there is a tie, select a node with a higher out-degree.
Assign the nodes connected the nodes at layer i to layer i+1.
Repeat steps 1 and 2 for the remaining nodes.
If two nodes at the same layer are connected to each other, make a new layer between the layer and the upper layer and move the node with a smaller degree to the new layer. Nodes with 0 out-degrees are also moved to the new layer.
Order the nodes in each layer by the Barycenter method .
If there are disconnected nodes, repeat steps 1-5 for the nodes and put them to the right of the previous subgraph.
We implemented the algorithms in a program called GeneNetFinder using Microsoft Visual C#. GeneNetFinder is executable on any Windows systems, and the program and sample data of GeneNetFinder are available at http://wilab.inha.ac.kr/GeneNetFinder. Given a time-series data of gene expressions in log-ratios, it identifies gene regulatory interactions and visualizes them. This section shows the experimental results with the gene expression data of yeast cell cycling and human cell cycling.
The dataset of the yeast cell cycle, shown in Figure 1, includes 30 genes of yeast cell cycle from the Saccharomyces cerevisiae cell cultures . The 30 yeast genes are known to be involved in the cell-cycle regulations. For the cell cycle genes, we first selected well-known genes and their regulator or regulated genes. CLB genes, for example, are known to promote cell cycle progression into mitosis . CLN genes were selected because they have regulatory relations with CLB genes. The remaining genes were chosen randomly. There are 18, 24, 10 time points in the time spans of alpha-factor, cdc15 and cdc28, respectively.
Gene regulations identified in the time-series expression data of yeast cell cycle.
+CLB1(T) → +SWI5(T+1)
+CLB1(T+1) → +CDC20(T+3)
+CLB1(T+2) → +CDC20(T+4)
-CLB1(T) → -CDC20(T+1)
+ CLB2(T+1) → +SICl(T+3)
+CLB2(T+2) → +SWI4(T+5)
+CLB2(T+3) → +SWI4(T+6)
+CLB6(T+1) → -CLB1(T+2)
+CLB6(T) → -CLB2(T+1)
+CLB6(T+1) → -CLB2(T+2)
-CLN1(T) → -CLB4(T+1)
-CLN1(T+1) → -CLB4(T+2)
-CLN2(T) → -SWI6(T+1)
-CLN2(T+2) → +SIC1(T+5)
+CLN3(T) → +SIC1(T+1)
+CLN3(T) → +CLB6(T+3)
+CDC5(T) → +CDC14(T+1)
+CDC5(T+1) → +CDC14(T+2)
+CDC5(T+3) → +CDC20(T+5)
+CDC5(T) → +CDC20(T+1)
+CDC5(T+1) → +CDC20(T+2)
-CDC14(T) → +CLN1(T+1)
-CDC14(T+1) → +CLN1(T+2)
-CDC14(T+2) → +CLN1(T+3)
+CDC28(T) → +SWI4(T+1)
-CDC34(T+2) → +CDC34(T+5)
+ CDC53(T) → -CLN3(T+1)
+CDC55(T) → +USV1(T+1)
+CDC55(T+1) → +USV1(T+2)
+CDC55(T+2) → +USV1(T+3)
+CDC55(T+3) → +USV1(T+5)
+CDC55(T) → +SWI5(T+1)
+CDC55(T+1) → +SWI5(T+2)
+CDC55(T+2) → +SWI5(T+3)
+HCM1(T) → -CLB5(T+1)
+HCM1(T+2) → -CLB5(T+4)
+HCM1(T) → -CLN1(T+1)
+MCM(T+1) → -MBP1(T+2)
-MEC1(T) → +CBF1(T+1)
-MEC1(T+3) → +CBF1(T+5)
+MGA1(T) → +CDC5(T+1)
+PDR3(T) → +SWI5(T+1)
+PDR3(T+1) → +SWI5(T+2)
+PDR3(T+2) → +SWI5(T+3)
-SKP1(T) → -SWI4(T+1)
-SKP1(T+3) → -MBP1(T+4)
-SKM1(T) → +CLB6(T+1)
-SKM1(T+1) → +CLB6(T+2)
-SKM1(T+2) → +CLB6(T+3)
-SKM1(T+3) → +CLB6(T+4)
-SIC1(T+1) → +SWI5(T+3)
-SIC 1(T+2) → +SWI5(T+4)
-SIC 1(T+3) → +SWI5(T+5)
-SIC1(T+1) → -CLN2(T+3)
-SIC 1(T+2) → -CLN2(T+4)
-SIC 1(T+3) → -CLN2(T+5)
+SWI5(T) → -CLNl(T+1)
+SWI5(T+1) → -CLN1(T+2)
-SWI5(T) → +CLB6(T+1)
-SWI6(T+1) → -SKP1(T+2)
-SWI6(T+3) → -CDC20(T+6)
There are several resources that can be used to assess the gene regulations inferred by the program. We searched the KEGG, SGD and CYGD databases and literatures to see whether the databases contain a gene regulation that agrees with the identified gene regulation. KEGG has 29,471 pathways whereas SGD (http://www.yeastgenome.org/) and CYPD (http://mips.gsf.de/genre/proj/yeast/) provide genetics and functional networks of the budding yeast Saccharomyces cerevisiae, respectively. Gene ontology  can also be used to obtain gene regulatory relations; if a gene has a GO annotation of 'regulates' or 'regulated by', the gene can be considered as being involved in a gene regulation.
A gene regulation is certain when it has at least two supporting evidences.
A gene regulation is possible when its R 2 score is either < 3 or > 6, and no other evidences.
A gene regulation is uncertain when it has no supporting evidence.
In the dataset of the yeast cell cycle, a total of 73 gene regulations were inferred by GeneNetFinder (Table 1). There were 48 certain gene regulations. 4 out of the 73 gene regulations, +CDC28->+SWI4(T+1),+HCM1(T) -> -CLB5(T+1), +HCM1 -> -CLN1(T+1) and +MGA1(T) -> +CDC5(T+1), had supporting evidences both in the databases and GO annotations. 43 out of the 73 gene regulations showed exact agreement with the known regulations in the databases. Only one gene regulation +USV1(T)->-SKM1(T+1) had the R 2 score > 6 and the GO annotation about the regulation. These 48 gene regulations are marked as underlined entries in Table 1, and more details are available in Additional file 1.
12 out of the 73 gene regulations, written in italics in Table 1, were possible regulations. They had the R 2 score either < 3 or > 6 with no further supporting evidence. These regulations should be verified by experimental evidence. The remaining 13 gene regulations had the R 2 score between 3 and 6 and no supporting evidence, so they are uncertain regulations. Even if the 13 regulations are false positives, at least 82.2% (60 out of 73) of the regulations inferred by GeneNetFinder agreed with known gene regulations. Table 1 shows the regulations for the first four time points only, and there are two more regulations identified at time points T+4 and T+5: -CLB6(T+4) -> +CDC28(T+5)  and -CDC20(T+4) -> -CLB6(T+5) .
For further evaluation, we selected some genes. Gene CCNA2 encodes proteins of the highly conserved cyclin family which plays an important role in the cell cycle. CCNA2 binds and activates CDC2 kinases, and thus promotes both cell cycle G1/S and G2/M transitions. Then CDC2 encodes proteins which are members of the Ser/Thr protein kinase family. This protein is a subunit of the highly conserved protein kinase complex and essential for G1/S and G2/M phase transitions. In the KEGG pathway database, we found that CDC2 interacts with E2F1 and SKT15, and E2F1 has direct regulatory relationships with CDC2, CCNA2 and BRCA1. The protein encoded by the gene E2F1 is a member of the E2F family of transcription factors. The E2F family plays a crucial role in the control of cell cycle and action of tumor suppressor proteins. In summary, CCNA2 interacts with E2F1, CDC2, and CDKN1A; CDC2 interacts with E2F1, CCNA2, CCNB1, CDC25A, CDC25B, CDC25C and CDKN1A; and CCNB1 interacts with CDC2, CCNF, BRCA1 and CDKN1A. All these regulations agree with the regulations identified by GeneNetFinder.
Effect of changes in genes on the prediction performance.
#genes in a dataset
#certain regulations (a)
#possible regulations (b)
#total regulations (c)
Accuracy (%) ((a+b)/c)
There are a few programs that can infer gene regulatory interactions from time-series gene expression data [18, 19]. ASIAN , for example, finds correlation relationships between gene clusters and visualizes them as an undirected graph. While gene regulatory networks visualized by GeneNetFinder represent gene regulatory interactions between individual genes along with temporal aspects of the interactions, networks visualized by ASIAN represent correlations between gene clusters. Thus, it cannot show regulatory interactions between individual genes nor the order or pace of the interactions. BioTapestry  is an interactive tool for building and visualizing gene regulatory networks. For visualizing gene regulatory networks BioTapestry uses different layout methods from GeneNetFinder, and temporal aspects of gene regulatory interactions are not automatically shown as edge labels of the networks. BioTapestry allows flexible edge labels, which can be annotated with any properties specified by the user.
Gene regulatory interactions usually change over time rather than being constant, but many databases or literatures provide static gene regulatory networks only. They are either snapshots of gene regulatory relations at a time point or union of successive gene regulations over time. Static gene regulatory networks are easier to construct and understand than dynamic networks, but cannot provide information on temporal aspects of gene regulations.
This article has presented an algorithm for qualitatively reasoning dynamic gene regulatory relations from gene expression data using two types of scores, R 1 and R 2 . The algorithm has been implemented in a program called GeneNetFinder. From the time-series data of gene expression, GeneNetFinder infers not only gene regulatory interactions but also the temporal aspects of the regulatory interactions. As for the temporal aspect of gene regulatory relations, it identifies the order of the gene regulatory relations and the pace of the relations. The identified gene regulatory interactions and their temporal aspects are stored in the regulation list and visualized as a gene regulatory network. In the network visualized, gene regulations and their temporal aspects are represented by edge types and edge labels.
We tested GeneNetFinder on several datasets, including the yeast cell cycle data and the human cell cycle data. Experimental results indicate that the dynamic nature of dynamic gene regulatory interactions can be identified and represented qualitatively without deriving or solving a set of differential equations describing the interactions. GeneNetFinder is yet a prototype, but the approach of our work would be useful for identifying dynamic gene regulatory interactions from the large amount of gene expression data available at the present time. In particular, the gene regulatory networks constructed by GeneNetFinder can be used to find new gene regulatory relations or to refine known regulatory relations.
Yu Chen developed the algorithm and prepared the first draft of the manuscript. Byungkyu Park helped Yu Chen develop the method. Kyungsook Han supervised the work and rewrote the manuscript.
This work was supported by the National Research Foundation of Korea (NRF) grant (2009-0091509) and by the Mid-career Researcher Program through NRF grant (2008-0058358) funded by the MEST.
This article has been published as part of BMC Genomics Volume 11 Supplement 4, 2010: Ninth International Conference on Bioinformatics (InCoB2010): Computational Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2164/11?issue=S4.
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.