Qualitative reasoning of dynamic gene regulatory interactions from gene expression data

Scoring scheme for gene regulatory relationships

The gene expression data of m genes with n samples is represented as an m × n matrix, where rows represent genes and columns represent various samples such as experimental conditions or time points in a biological process. Each element of the matrix represents the expression level of a particular gene in a particular sample. Two genes with similar expression patterns tend to be co-expressed at different time points. Figure 1 shows an example of the gene expression data for yeast genes during the yeast cell cycle, obtained from the Yeast Cell Cycle Analysis Project [6].

The gene expression matrix is analyzed for similarity between gene expressions at different time points. Three metrics are often used to measure the similarity of genes: Pearson correlation coefficient [7], Euclidean distance [8] and Spearman correlation [9]. To evaluate the regulatory relation between two genes, we modified the Pearson correlation coefficient. R 1(X, Y, i, p) in Equation 1 represents the correlation between gene X at time point i and gene Y at time point i + p. p is the time span of the gene regulation.

(1)

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.

The modified Pearson correlation coefficient R 1 is useful for finding gene regulatory relations with a significant change in expression levels. But, it cannot distinguish gene regulatory relations with the same correlation but different gene expression levels (see Figure 2 for an example). Therefore, we use an additional score R 2, which is the Euclidean distance of the expression levels of the two genes. The score R 2 for the gene regulatory relation between X and Y is computed by Equation 2.

(2)

where

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.

Inferring gene regulatory relationships

In the microarray data for gene expression, we use the log-ratio (in base 2) of the red and green intensities. Thus, genes with mRNA abundance have positive log-ratios whereas those with absence of mRNA have negative log-ratios. From the gene expression profiles, we identify the gene regulations and include them in the regulation list. In the regulation list, +A(t) indicates that gene A is up-regulated at time t, and -A(t) indicates that gene A is down-regulated at time t. The symbol '→' represents a directional relationship between genes. There are four possible gene regulatory relations between two genes:

1. +A(t₁) → +B(t₂): up-regulation of A at time t₁ is followed by up-regulation of B at time t₂ (t₂ >t₁).

2. -A(t₁) → +B(t₂): down-regulation of A at time t₁ is followed by up-regulation of B at time t₂ (t₂ >t₁).

3. +A(t₁) -> -B(t₂): up-regulation of A at time t₁ is followed by down-regulation of B at time t₂ (t₂ >t₁).

4. -A(t₁) → -B(t₂): down-regulation of A at time t₁ is followed by down-regulation of B at time t₂ (t₂ >t₁).

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₁, p). Algorithm 1 provides the top-level description of the method for inferring gene regulations and constructing a list with the regulations.

Algorithm 1

Construct a regulation list

1: Compute R 1(A, B, t₁, p) between gene A at time point t₁ and gene B at time point t₁ + p for all pairs of genes.

2: Select the regulation with the largest absolute value of R 1(A, B, t₁, p).

3: If 0 <p < 6, classify the regulation into one of the four types, +A(t₁) → + B(t₁+p), - A(t₁ ) → + B(t₁+p), +A(t₁) → -B(t₁+p), -A(t₁) → - B(t₁+p), and add it to the regulation list.

4: If p = 0, two genes are co-expressed or co-inhibited, and such gene regulation is not added to the regulation list.

5: If the new gene regulation is already in the regulation list, merge it with the previous regulation.

6: Go to step 2 to find the next gene regulation until no more regulation found.

After we construct a regulation list, we compute the R 2 score for the gene pairs in the regulation list. Some local maximum or minimum values are ignored when computing the R 2 score in a long time span. For example, all the three time spans shown in Figure 3A include at least 10 time points (18 time points in alpha-factor, 24 in cdc15, and 10 in cdc28). Expression levels of gene CLB1 has several wave peaks in the time span of CDC15, but only the maximum value in sub-timespan 2 of CDC15 will be selected when computing the R 2 score in the time span of CDC15. Time spans are divided into smaller sub-timespans as follows, and the R 2 score is not computed for sub-timespans with less than 6 time points.

1. A time point with the minimum expression level of the regulator gene becomes a splitting point of the time span.

2. 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.

3. 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.

Visualization of gene regulatory networks

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.

For visualization of gene regulatory networks, two layout algorithms have been developed: grid layout and layered layout. As described in Algorithm 2, the grid layout algorithm positions all nodes at grid points. The node with the highest degree will be placed at the center grid point (node S in Figure 4A). Then, we position all nodes connected to the center node at adjacent grid points. Nodes with a higher degree are positioned earlier than those with a lower degree in the east, north, west, south, northeast, southwest, and southeast grid point of the current node (nodes 1-8 of Figure 4). Other nodes connected to the positioned nodes are placed in the same manner.

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.

Microarray data of 30 genes in the yeast cell cycle

The dataset of the yeast cell cycle, shown in Figure 1, includes 30 genes of yeast cell cycle from the Saccharomyces cerevisiae cell cultures [6]. 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 [11]. 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.

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 [13] 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.

Gene regulations identified by the program are classified as certain, possible and uncertain depending on the number of supporting evidences (agreement with the databases, agreement with GO annotations, R 2 score < 3 for activation and R 2 score > 6 for inhibition):

1. A gene regulation is certain when it has at least two supporting evidences.

2. A gene regulation is possible when its R 2 score is either < 3 or > 6, and no other evidences.

3. 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) [14] and -CDC20(T+4) -> -CLB6(T+5) [15].

The gene regulations should be consistent with the phase characteristics of the cell cycle. For example, CLB6 promotes progression of cells into the S phase and expresses periodically throughout the cell cycle [16]. CLB1 and CLB2 both promote cell cycle progression into mitosis and their transcripts accumulate during G2 and M. These biological processes are characterized by two regulations, +CLB6 -> -CLB1 and +CLB6 -> -CLB2, and these are included in the regulation list found by GeneNetFinder (Table 1). Figure 5 shows the gene regulatory network of all the gene regulations of Table 1.

Microarray data of 20 genes in the human cell cycle

This data set includes 20 genes in the human cell cycle [17]. The gene expression during the human cell cycle is synchronized by double thymidine blocking (Thy-Thy1, Thy-Thy2, and Thy-Thy3), thymidine-nocodazole blocking (Thy-Noc) and Mitotic selection (M). Additional file 2 shows the data of 20 genes in the human cell cycle. All regulations identified by GeneNetFinder are given in Additional file 2. 71 of 113 potential gene regulations were found in at least one of KEGG, Entrez Gene, BIND and PUBMED, and 44 of the 71 regulations had been determined by experimental methods. Thus, at least 62.8% of the identified regulations are in agreement with known regulations. Figure 6 shows the regulatory network of 15 human genes in the first time span along with the user interface of GeneNetFinder.

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.