PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of narLink to Publisher's site
 
Nucleic Acids Res. 2010 January; 38(1): e1.
Published online 2009 October 23. doi:  10.1093/nar/gkp822
PMCID: PMC2800212

Maximization of negative correlations in time-course gene expression data for enhancing understanding of molecular pathways

Abstract

Positive correlation can be diversely instantiated as shifting, scaling or geometric pattern, and it has been extensively explored for time-course gene expression data and pathway analysis. Recently, biological studies emerge a trend focusing on the notion of negative correlations such as opposite expression patterns, complementary patterns and self-negative regulation of transcription factors (TFs). These biological ideas and primitive observations motivate us to formulate and investigate the problem of maximizing negative correlations. The objective is to discover all maximal negative correlations of statistical and biological significance from time-course gene expression data for enhancing our understanding of molecular pathways. Given a gene expression matrix, a maximal negative correlation is defined as an activation–inhibition two-way expression pattern (AIE pattern). We propose a parameter-free algorithm to enumerate the complete set of AIE patterns from a data set. This algorithm can identify significant negative correlations that cannot be identified by the traditional clustering/biclustering methods. To demonstrate the biological usefulness of AIE patterns in the analysis of molecular pathways, we conducted deep case studies for AIE patterns identified from Yeast cell cycle data sets. In particular, in the analysis of the Lysine biosynthesis pathway, new regulation modules and pathway components were inferred according to a significant negative correlation which is likely caused by a co-regulation of the TFs at the higher layer of the biological network. We conjecture that maximal negative correlations between genes are actually a common characteristic in molecular pathways, which can provide insights into the cell stress response study, drug response evaluation, etc.

INTRODUCTION

A molecular pathway is referred to as a series of actions among molecules in a cell leading to a certain end point of cell function. Pathway identification is usually aimed to uncover all biological molecules participating in the same functional pipeline, which may include DNA/gene, miRNA, protein or metal ion, etc. As DNA and protein play the major roles in a pathway, gene and protein's; indirect relations are of paramount importance for detecting and analyzing molecular pathways.

Gene expression data, especially time-course gene expression data, have been widely used to explore various relationships of the genes in the pathways, with the particular focus on the positive correlations. For example, Segal et al. (1) proposed to identify new pathways by assuming that most genes in the same pathway can exhibit a similar gene expression profile, and their proteins often interact. Multiplicative patterns and scaling patterns have been also used to describe the expression profiles of the genes in the same pathway (2–4). Co-regulation patterns, additive expression patterns or shifting patterns, have been conceptualized to detect regulatory modules from gene expression data (5,6). Further, geometric patterns based on trigonometric functions are believed to be related to circular regulation processes (7). Here, concepts such as profile similarity, shifting pattern, scaling pattern or geometric pattern are all concentrated on positive correlations among genes, implying that genes with expression homogeneity are possible to have the same biological function.

In this work, we are interested in negative correlation. It is also an important relationship among genes, and it has been previously observed in many biological processes. Schmid et al. had a study on development expression patterns for large gene families of Arabidopsis thaliana (8); they highlighted two groups of genes showing an opposing expression trends from an early seed development stage to a late stage (Figure 1a). In a study of expression patterns in the chondrogenic differentiation, James et al. (9) had a careful analysis on a 15-day temporal gene expression data of a Mouse micromass culture system and they reported an interesting example on two groups of genes displaying an opposite expression pattern. In that example, transcription factors (TFs) Sox9 showed high expression levels before day 6, then their expression decreased by about half; while, Ibsp transcripts showed low expression levels until day 12 and then had 1200-fold up-regulation (Figure 1b). Chuang et al. (10) introduced the notion of complementary gene expression patterns for inferring time-lagged genetic interactions, which is intended to capture contrasting expression patterns like the one that: when one gene's; expression increases, the other gene's; decreases, or vice versa. As another example, Stekel et al. (11) proposed a method for modeling the so-called self-negative regulation of TFs, in which the product of one gene is assumed to regulate this gene's; TF in a feedback way.

Figure 1.
Notion of negative correlations in biological studies. (a) Reprinted by permission from Macmillan Publishers Ltd: Nature Genetics, 37(5):501–506, ©2005. (b) Reprinted by permission from American Society For Cell Biology: Molecular Biology ...

Based on above biological ideas and primitive observations, we propose to formulate and systematically investigate the problem of maximizing negative correlations. The objective is to discover all maximal negative correlations from time-course gene expression data. A maximal negative correlation is defined as a pair of correlated gene sets where genes between the two sets must be negatively correlated in their expression over a time segment, and the number of genes and the number of time points in the segment are required to be maximized. Therefore, maximization of negative correlations, like detecting positive correlations, can expand and deepen down the analysis on molecular pathways. This is in a good agreement with a biological fact that: when two genes have a negatively correlated interaction, then the two group genes locating at the downstream of their participating pathways will also form a negative relationship (12). To capture such a pipeline of negative correlations spanning multiple time points from time-course gene expression data, maximization of negative correlations is a novel and effective attempt.

Maximization of negative correlations is extendable and potentially applicable to a wide range of pathway studies where biological molecular relationships can be mapped to gene expression correlations. Segal et al. (13) were interested in motif extraction from sequences of promoters, and observed that the expression patterns of regulatory motifs with nucleotide variation were bi-polarly different. See one example of such expression patterns in Figure 1c. Shieh et al. (14) proposed to study transcriptional compensation interactions or synthetic lethal pairs with the idea that: following some gene's; mutation, its compensatory gene will be upregulated (or downregulated). Millar et al. (15) were interested in the whole-genome pattern of histone lysine acetylation and methylation in Yeast to confirm a hypothesis that different combinations of histone modification sites are likely associated with specific and contrasting transcription behaviors (Figure 1d). They also pointed out that these patterns can exist in other organisms such as Schizosaccharomyces pombe genome, portions of the Drosophila melanogaster even Human genomes. Recently, research results all show that such histone modification patterns are correlated with Human diseases (16–18). Therefore, through maximizing negative correlations, these biological applications can be certainly deepened down.

Given a gene expression data matrix, a maximal negative correlation can be viewed as an activation–inhibition two-way expression pattern (AIE pattern), where the two groups of genes exhibit such a behavior that when one group of genes is upregulated, the other group is low-expressed, or vice versa, consistently at a continuous range of time points. Identifying a complete set of significant AIE patterns from gene expression data is computationally expensive. We design a new graph-based method for an exact and complete enumeration of AIE patterns with high efficiency. Our algorithm combines two mining strategies: a suffix-tree structure and a bi-clique approach for efficient search of the AIE patterns. To our best knowledge, there is no algorithm that can be specialized to identify AIE patterns. Clustering methods may be easy to find out activation–inhibition relation (7,19), but there will be a lot of false positives, and local negative correlations under different specific time points cannot be identified. Biclustering can find gene expression patterns related to specific conditions, but it is difficult to mine large number of genes within negative relations (20,21). The so-called anti-correlated patterns (22,23) are closely related to our AIE patterns, however, their mining algorithm cannot produce exactly AIE patterns (see a detailed comparison later). In our in-silicon evaluation, our method has been successfully applied to Yeast time-course gene expression data to reveal negative correlations in the molecular pathways of Saccharomyces cerevisiae for increasing the understanding of its biological mechanisms.

METHODS

Let M be a time-course gene expression data set denoted as a triplet M = (G,C,d ), where G = {g1, g2, …, gn} is a set of genes (rows), C = {c1, c2, …, cm} is an ordered set of continuous time points (columns) and d : G × C [implies] R is the level function by which d(gi, cj) represents the expression level of gene gi at time point cj.

A continuous subset of C = {c1, c2, … , cm} is an ordered subset of C with continuous time points. In other words, if T = {t1, t2, … , tk} is a continuous subset of C, then An external file that holds a picture, illustration, etc.
Object name is gkp822i1.jpg and ti = ci + j, i = 1, … , k, for some j[set membership]{0,1, … , mk}.

Definition 1 (AIE pattern). Let X = (I,J,d) be a submatrix of a time-course gene expression data matrix M = (G,C,d), where I is a subset of G and J is a continuous subset of C, X is defined as an AIE pattern if I can be divided into I1 and I2 such that at every time point j [set membership] J, the expression levels of genes in I1 and those in I2 satisfy either

equation image

or

equation image

Suppose X = (I, J, d ) is an AIE pattern in a time-course gene expression data matrix. Usually, we require |I| and |J| that are both maximal. This means that there are no genes or no time points that can be added into X = (I, J, d ) to maintain the conditions of the above definition.

Our method to enumerate a complete set of AIE patterns from a gene expression data set consists of three computational steps. The first step is to construct a dichotomy matrix based on the original data set, which captures and discretizes the expression difference between every pair of genes at every time point. The second step is to transform the dichotomy matrix into another representation as suffix tree, and to extract the time-series biclusters from this tree in order to locate the time period when the genes possibly show a negative correlation. The third step is to construct a bi-partite graph according to the row ids (i.e. the gene pairs) from those biclusters, and to distinguish two groups of genes forming a bi-clique in such bi-partite graph.

Definition 2 (Dichotomy Matrix). Given a matrix M = (G, C, d), its dichotomy matrix DM|G|×(|G|–1)/2,|C|(M) is defined as:

equation image

where i1 < i2. For each row k in DM(M), if

equation image

then, the row ID of this row is assigned as <i1,i2>.

See an example of dichotomy matrix at the left panel of Figure 2b which is derived from the gene expression matrix shown in the left panel of Figure 2a. In fact, each row of the dichotomy matrix can be considered as a 0–1 sequence, so a suffix tree of all 0–1 sequences/rows (the right panel of Figure 2b) can be constructed in linear time (22). The depth of the nodes corresponds to the number of time points; the leaf nodes and the splitting nodes are marked with the row/gene pair IDs of the dichotomy matrix. Thus, every splitting node or leaf node is a maximal submatrix with every row identical in the dichotomy matrix. A proof about the relation between the nodes in a generalized suffix tree and the maximal biclusters with continuous columns can be found in (22).

Figure 2.
A suffix-tree and a bi-clique search method are combined for AIE pattern mining.

For every splitting or leaf node, we construct a bi-partite graph using the gene pair IDs stored at the node. Suppose a node contains k number of gene pair IDs: An external file that holds a picture, illustration, etc.
Object name is gkp822i2.jpg, then we denote An external file that holds a picture, illustration, etc.
Object name is gkp822i3.jpg as the nodes at one side of the bi-partite graph, and denote An external file that holds a picture, illustration, etc.
Object name is gkp822i4.jpg as the nodes at the other side of the bi-partite graph. Meanwhile, assign an edge between An external file that holds a picture, illustration, etc.
Object name is gkp822i5.jpg and An external file that holds a picture, illustration, etc.
Object name is gkp822i6.jpg for j = 1, 2, … , k. See an example of bi-partite graph at the left panel of Figure 2c. Then, we enumerate all maximal bi-cliques from this bi-partite graph. Assume I1 and I2 are the two vertex sets of such a maximal bi-clique, then I1 and I2 are exactly the two nonoverlapping gene groups for an AIE pattern whose time points are decided by the edge labels of the path leading to the splitting node or the leaf node from the root node in the suffix tree.

The pseudo-code of our algorithm is shown in Algorithm 1. The core subroutine of this algorithm for the bi-clique mining is taken from (24). The whole algorithm can be divided into two parts and their computational complexities are analyzed as follows. The first part is a determinant routine to solve P-problem that discovers maximal row-identical submatrices. The initiation of this algorithm (line 4) needs O(G²C) time and space; the suffix-tree construction (line 6) needs O(G²C) time and space (20); and the identification of maximal row-identical submatrices from the suffix-tree (line 8) needs O(G²C²) time. So the time and space cost of this determinant subroutine are O(G²C²) and O(G²C). While, the second part is an exhaustive pattern mining subroutine to solve an nondeterministic polynomial time (NP)-problem that discovers all bi-cliques. The size of bi-partite graph input into bi-clique mining (line 9–17) is no more than

Algorithm 1 Suffix-tree and bi-clique two-stage method for detecting AIE patterns in time-course data

Require: a matrix M of size G × C
Ensure: AIE patterns
 1: Int GPN = G(G–1)/2
 2: Int GP[GPN][C]
 3: /*Produce dichotomy matrix of M at first stage*/
 4: GP = Dichotomy(M)
 5: /*Construct Suffix-tree at second stage*/
 6: SuffixTree ST = SuffixTree(GP)
 7: /*Extract row-identical submatrix Bic(I*, J ) in suffix-tree*/
 8: Bic <ROW,COL> = ExtractBic(ST)
 9: for each <rows,cols> in <ROW, COL> do
10: /*Construct bi-partite graph B(I)*/
11: Gragh GG = GenesGraph(rows)
12: /*Bi-clique mining at third stage*/
13: Biclique <ROWG,COLG> = Biclique(GG)
14: for each <rowsG,colsG> in <ROWG,COLG>do
15:  Output <rowsG,colsG,cols> is an AIE pattern
16: end for
17: end for

G × C, so under the worst condition, its time complexity is O(G²N) and space complexity is O(G²) (24). Here, N is the number of all maximal bi-cliques, or the number of all AIE patterns. It should be noted that, in the worst situation, the distribution of all potential AIE patterns will be dense in the original data. That means the N is so large that the time cost of AIE pattern mining will increase tremendously. Therefore, in actual applications, we will use the size, up-down index, and differential gap (see their definitions later) to limit the number of potential AIE patterns. The proposed algorithm will perform an exhaustive search for all AIE patterns under such specific constraints.

By definition, the negative correlation in an AIE pattern may go like the way that: at all time points in J, the expression of the genes in I1 is always higher, or always lower, than those of genes in I2. This is an extreme case of negative correlation. For the other cases, I1's expression may be higher at the first time point, while turn to be lower at the second time point, and then come higher again than I2's expression. This up-down trend at multiple time points can be measured by an ‘up-down index’ value. This index can be used to categorize the expression trends exhibited by different AIE patterns. We also believe that this index value of an AIE pattern is related to the strength of a negative correlation.

Definition 3 (Up-down index). Let X = (I, J, d ) is an AIE pattern, its up-down vector U is a |J|×1 vector determined by:

equation image

An up-down index value UI of X is defined as:

equation image

Figure 3 shows three examples of AIE pattern where the two genes marked with solid line are in one group, while the other two labeled as dotted line are in the second group. The AIE pattern in Figure 3a has an up-down index value of 3/5. While, the AIE pattern in Figure 3b has an UI value as 1/5. In fact, the up-down index value of an AIE pattern is the up-down change frequency of one group's; expression against the other group's; over a continuous time points. For the case of negative correlation shown in Figure 3c, its up-down index value is 0. However, it still looks interesting. Therefore, we introduce a generalized Pearson's; correlation coefficient to measure the negative correlation between two sets of variables. This is a finer correlation measurement compared with the up-down index.

Figure 3.
Examples of AIE patterns and their up-down trends.

Let I1 and I2 be two sets of gene variables with no overlapping. The generalized Pearson's; correlation coefficient (noted as R-value) between I1 and I2 is denoted by R(I1,I2), and it is calculated by

equation image

where r(gi,gj) is the Pearson's; correlation coefficient between the two variables gi and gj. Note that the value of R(I1,I2) is between −1 and 1, exactly with the same range as the conventional Pearson's; correlation. The correlation is the most negative when the value −1 is reached, while on the other hand, it is the most positive when the value of 1 is approached. The generalized Pearson's; correlation coefficient can also apply to a set I of gene variables to measure its inherent negative correlation if I can be properly divided into two sub-groups I1 and I2. In this work, such I1 and I2 are obtained by a hierarchical clustering approach.

Given an AIE pattern, the number of genes in one subgroup (I1 or I2) can be sometimes much smaller than the other group. So, we also introduce a ‘group ratio’ index GR(I1,I2) to measure the size balance between I1 and I2 for an AIE pattern, as calculated by

equation image

To gain more insight into the negative correlation of an AIE pattern, we also examine how wide the expression between the two groups is. When the difference becomes wider, the expression behavior of the two groups are more distinct, thus the negative correlation may be more significant. To highlight this expression divergence between the two groups of genes I1 and I2, we define a ‘differential gap’ to average the minimum difference of expressions between the two groups of genes, namely An external file that holds a picture, illustration, etc.
Object name is gkp822i7.jpg over multiple time points. This differential gap is calculated by

equation image

RESULTS

The data used in our evaluation is a time-course gene expression data repository related to Yeast. The raw data were published by the Yeast cell cycle analysis project (25). This project had acquired the expression measurements of 6179 genes involved in the Yeast cell cycle under three different conditions: Alpha factor-based synchronization, Cdc15-based synchronization and Cdc28-based synchronization. These measurements are actually the relative expressions against the control/background data; the zero expression represents the control level. We denote these three data sets simply as Alpha factor, Cdc15 and Cdc28 data set in this article. The number of time points was set by the project as 18, 24 or 17, respectively, for the Alpha factor, Cdc15, or Cdc28 condition. The time point number 18 for the Alpha factor condition (similarly, 24 and 17 for the other two conditions) means that the RNA samples for the 6179 genes were collected at 18 time points starting at 0 to 7, 14, 21, …, till the 119th min, which covers two cell cycles. In our data preprocessing, genes that do not occur in any known pathways were removed according to the Saccharomyces Genome Database (SGD) (26), which contains 142 known Yeast pathways covering 515 genes. There were 502 genes finally left in the three expression data sets, which were subsequently used as input in our experiments. We note that each of these 502 genes participates into at least one pathway, and some of them participate into up to seven different pathways. The number of genes involved in these 142 pathways varies from 1 to 23 with 5 on average.

We compare the P-values (biological significance) and R-values (negative significance) of AIE patterns with those of the gene clusters found by the widely used clustering/biclustering methods. The purpose of this comparison is to confirm that many important negative correlations were unable to be identified by those conventional algorithms. We also present representative AIE patterns, and other statistics information of AIE patterns, including the up-down index values and the differential gap information. More importantly, we take case studies to illustrate how the AIE patterns are biologically interpreted for enhancing the analysis of molecular pathways.

Comparison by using P-values and R-values

Conventional clustering methods under our comparison include a hierarchical clustering method (HCL) (19,27), a K-means clustering method (K-means) (27), the Cheng & Church biclustering method (CC) (3,27), the order preserving sub-matrix algorithm (OPSM) (3,27) and e-CCC-Biclustering (e-CCC) (23). All of their implementation are available at BicAT (27) or at BiGGEsTS (28).

It's; specially noted that the notion of CCC-Biclusters (22) and its extension e-CCC-Biclusters are closely related to our concept of AIE patterns, in particular when the sign-change rule (U↔D) is combined to form CCC-Biclusters or e-CCC-Biclusters. By definition, a CCC-Bicluster has to satisfy the condition that every possible gene pair in this bicluster shares a positive (coherent) expression change behavior over the time. e-CCC-Bicluster extends CCC-Bicluster by allowing a certain degree of noise (measured by the parameter e) in a CCC-Bicluster, such that the coherent expression change behavior in an e-CCC-Bicluster may not be always the same on some time points. The sign-change rule introduced in the e-CCC-Bicluster mining algorithm (23) enriches the diversity of these biclusters, and it can be used to detect the so-called anti-correlated patterns.

As introduced, the definition of AIE patterns is simple and different. Two non-overlapping gene groups I1 and I2 can form an AIE pattern if and only if I1 and I2 have a negative expression change behavior over a time segment (Definition 1). That is, the genes in the same group (I1 or I2) are not necessarily required to have exact coherent expression change over the time segment. Even if the sign-change rule is applied to AIE patterns, the gene pairs within I1[union or logical sum]I2 may still not have much coherence. Though an anti-correlated pattern may sometimes become an AIE pattern, on the other hand, an AIE pattern usually does not satisfy the conditions required by an anti-correlated pattern. Therefore, they two are not equivalent. Another difference lies in the algorithms of mining e-CCC-Bicluster and AIE patterns. Although both of them take an exhaustive enumeration approach, e-CCC-Biclustering needs the parameter e to control the noise level of e-CCC-Biclusters and also needs sign-change to allow negative correlation, while AIE pattern mining is a parameter-free algorithm (many proposed indexes such as up-down index and differential gap are just used in the post-analysis, although these indexes can also be used as predefined parameters to reduce the running time of AIE pattern mining). We would also like to point out that both CCC-Biclusters and AIE patterns use the maximality rule to get rid of some redundancy in the patterns.

In our comparison, the input parameters of these conventional methods are all set as the default values as suggested in (23,27). In order to avoid possible comparison bias, some of the gene clusters obtained by the biclustering methods are filtered if the overlap degree between any two clusters is >25% as previously done by (22,23,27).

We first compare P-values (29) which can indicate a biological significance of a gene cluster indirectly. In this article, P-values are calculated through gene set enrichment analysis with Fisher's; exact test (29). A P-value <10−3 is widely accepted as a gold standard in most biological significance analysis. Figure 4 shows the bar charts representing proportions of gene clusters whose P-values are less than some thresholds. For example, from the Cdc28 data set, there are 57 592 AIE patterns identified. Of them, 151 are left after the filtering, in which 99% have a P-value less than <10−3, and 22% with a P-value <10−6. In many cases, AIE patterns have better proportions of biologically significant gene clusters than those of the gene clusters found by the conventional methods. It is worth noting that the OPSM method seems to perform better than the other biclustering methods. A possible reason is that it only outputs no more than 30 top clusters and no more than 5 clusters after filtering (shown in the second column of Table 1). CC and e-CCC-Biclustering are slightly better than our AIE method in terms of these P-values. We also note that when the number of clusters is set to be big by HCL or Kmeans, the genes contained in each cluster tend to be small. As the size of a cluster affects the calculation of P-values, the performance of HCL or Kmeans is influenced greatly by the pre-given number of clusters. See Table 1 for detailed information of the clusters under different settings.

Figure 4.
Biological significance comparison between our AIE patterns and gene clusters by the conventional methods. Here, HCL-x or Kmeans-x stands for x number of clusters being pre-set.
Table 1.
Different methods on negative correlation mining

We have seen that AIE patterns have close competitive P-values with those of the gene clusters generated by the conventional clustering methods. Next, we report the R-values and group ratio information of AIE patterns. Figure 5a, Figure 5c and Figure 5e shows a box-plot view of the R-values of the gene clusters from the Alpha factor, Cdc15 and Cdc28 data set, respectively, by different methods. We can see that the R-values of the AIE patterns are almost all less than 0. This result is in full agreement with our original notion of negative correlation. In particular, nearly half of the AIE patterns have a R-value lower than −0.5. However, among the gene clusters found by the traditional clustering methods and OPSM, usually only a small proportion of them have negative R-values. So combining the R-values and P-values, we can note that only CC, e-CCC-Biclustering and our AIE approach are good to detect biologically significant and negative correlations. We also consider the group ratio information GR(I1,I2) of the clusters in the comparison. The CC and e-CCC-Biclustering method both tend to find clusters with small group ratio values around 0.2, while our AIE patterns usually have a group ratio value higher than 0.4. This indicates that CC and e-CCC-Biclustering prefer negative correlations between unbalanced gene group pairs. However, our AIE approach can indeed find negative correlation spanning two size-balanced gene groups which can have strong biological significance. See Figure 5b, d and f for a detailed comparison of the group ratio information.

Figure 5.
The generalized pearson’s correlation coefficient and group ratio comparison between AIE pattern mining and the conventional methods.

Representative AIE patterns from the three expression data sets

Table 2 presents statistics information of all AIE patterns discovered from the three data sets when the gene number threshold for |I1| and |I2| is set as 5, and the minimal number of time points in |J| is set as 10. In this table, the column ‘No. of AIE Patterns’ refers to the total number of AIE patterns from one data set, together with their average number of genes and their average time points in one pattern; the column ‘Differential gap’ indicates the smallest, biggest and average expression differences between the two subgroups for the AIE patterns in each data set; the column ‘Up-down index’ indicates the smallest, biggest and average up-down index values for the AIE patterns in each data set; and, the column P-value shows the minimal, maximal and average P-values of the AIE patterns in each data set. From this table, we can see that the shape and property of AIE patterns can vary very much. There are many other choices to set the size threshold for |I1|, |I2| and |J|; readers are referred to use our web site http://sunim1.birc.ntu.edu.sg/~aie/ to get the statistics information when a different threshold is set.

Table 2.
Some statistics information on the complete sets of AIE patterns

Figure 6 displays three typical examples of negative correlation under the three experimental conditions. These negative correlations can be categorized into two kinds of behavior according to their up-down index values. Under the experiment condition Alpha factor or Cdc28, the expressions of the two gene groups do not up-down change frequently. Figure 6a shows such an expression trend where the expression up-down change happened only once which was at the time point 9.

Figure 6.
Representative examples of AIE pattern from the three data sets. Here, two groups of genes colored in red and blue, show negative correlation during the time points in orange area.

However, most AIE patterns under the condition Cdc15 can have up-down index values close to 1.0. See the fourth column of Table 2. This means that the gene groups change their expression up-down very frequently under this cell environment. Figure 6b shows a perfect example for this expression trend that continuously crosses 14 time points. Therefore, we can conjecture that negative correlation can behave in tremendously different ways in the Yeast cell cycle when different conditions are applied, meanwhile their P-values are all very significant (as shown in the fifth column of Table 2 as well as in Figure 4).

Biological interpretation of AIE patterns: three case studies

We take special case studies to illustrate how AIE patterns and negative correlations can be used to infer new modules of gene regulation, TFs and regulatory networks (RNs). Our first example is the representative AIE pattern identified from the Alpha factor data set, whose gene expression profile is shown in Figure 6a. Total 11 genes are involved in the two gene subgroups of this pattern that have a negative correlation spanning 14 continuous time points. One subgroup consists of six genes: we specially name it ‘red’ group and denote by Ir = {YIR038C (GeneID:854856), YNR001C (GeneID:855732), YKL127W (GeneID:853732), YOL126C (GeneID:853994), YJL068C (GeneID:853377), YLR328W (GeneID:851039)}. The other group consists of five genes: We specially name it ‘blue’ group and denote by Ib = {YDR234W (GeneID:851820), YDL182W (GeneID:851346), YDL131W (GeneID:851425), YBR265W (GeneID:852568), YNR050C (GeneID:855786)}. We also denote this pattern simply as AIE-Alpha-87. It is one of the most negative AIE patterns in Alpha factor according to their R-values; and its P-value is 2.9 × 10–10. The main covering pathway of AIE-Alpha-87 is lysine biosynthesis, which has seven genes known currently and four of them are contained in the blue group Ib.

A new regulatory module and its putative TF

Let us start the analysis on the five genes in the blue group Ib of AIE-Alpha-87. As mentioned, four genes in Ib are directly involved in the lysine biosynthesis pathway. In fact, the four genes are co-regulated by a known TF YEL009C. As all of the five genes in Ib are inherently co-expressed, we can infer that they altogether form a regulatory module with YEL009C as a co-TF. To confirm this hypothesis, we examined the whole upstream of each gene to identify their binding motifs by using YEASTRACT (30). Three possible YEL009C binding sites were identified (Table 3). We can infer that the first four genes in this table share a binding motif ‘TGACTGA’, while the last two genes YNR050C and YBR265W share a binding motif ‘TGACTMT’. Thus, these five genes are all likely to be co-regulated by YEL009C through its binding upon regulation segments in the upstream of the five genes. This is a new insight into the gene regulatory behavior of this module and its TF. This new understanding is mainly attributed to the positive relationship of the genes within the blue group.

Table 3.
Possible YEL009C binding sites at the upstream of the five genes in the blue group of AIE-Alpha-87

Building a tree-structure RN

On top of the idea of Boolean RNs (31), we incorporate our negative correlations and introduce a tree-structure RN for genes involved in an AIE pattern. We take two steps to complete the induction of these trees. The first step is to use YEASTRACT (30) to construct an initial RN by taking as input all the genes in an AIE pattern as well as all of their known TFs in Yeast; the second step is to trim the RN to eventually become a tree structure. The trimming constraints are set as follows: (i) all the leaf nodes of the tree are required to represent the genes in the AIE pattern, (ii) the inner nodes all represent the known TFs, and (iii) the tree can be decomposed into two sub-trees by removing one edge in such a way that each sub-tree exactly covers all genes in one group of the AIE pattern. This edge under the removal is critically important. The regulation represented by this edge should be either ‘enabled’ or ‘disabled’, and there should be at least one negative regulator to produce this negative correlation between the expressions of the two gene groups. Figure 7 shows a tree-structure RN built from AIE-Alpha-87.

Figure 7.
A subregulatory tree for AIE-Alpha-87.

This regulatory subtree possesses strong biological meanings. First of all, the regulatory relations indicated as edges in this tree all have direct/indirect evidence supported by known literatures (30). In particular, YLR451W is a negative regulator as reported and studied in (32), while other TFs are positive regulators. Next, we explain why genes in the blue group can show the negative correlation with the genes in the red group. This is likely due to the participation of the negative regulator YLR451W (32), which is the next TF after YEL009C on the regulatory paths from the root YHR084W to the genes in the blue group. As the TFs on the paths from the root to the red genes are completely different from the TFs on the paths from the root to the blue genes, it is most likely that the negative co-regulation between the genes in the two groups is caused by YHR084W with an auxiliary help from YLR451W. In fact, a recent study found that YHR084W is a specific Yeast cell cycle TF (33). Thus, we can see that an AIE pattern can not only identify two groups of genes with a negative correlation on gene expression profiles, but also can imply that one group genes possibly come from a same pathway, having similar expression behavior as reference to the genes in the other group. Therefore, such negative correlations can uncover the connections among TFs which are at a higher layer in the RN.

New pathway components

Through the above analysis, we have already understood that there exists a possible regulation pathway of YEL009C on the five genes in the blue group, and that the six genes in the red group have a negative correlation with lysine biosynthesis which is likely caused by the co-regulation of YHR084W and the negative regulation of YLR451W. We next study whether the 11 genes of AIE-Alpha-87 have any biological function relations with lysine biosynthesis.

A known lysine biosynthesis pathway diagram can be obtained from SGD (26). As mentioned, four (YDL131W, YDL182W, YDR234W and YNR050C) of the eleven genes of the AIE-Alpha-87 pattern participate directly in this pathway. To see whether the other seven genes have any function related to a biological molecular or a biochemistry reaction in the current lysine biosynthesis, we attempted a function annotation. A more complete diagram for the lysine biosynthesis pathway is shown in Figure 8, where the dashed line means a new component assessed by the annotation. In detail, genes YBR265W and YOL126C can affect on molecules NAD and NADPH through oxidation. Genes YNR001C, YKL127W, YIR038C and YJL068C are also indirectly related to this pathway through their biochemical functions on the direct partners of lysine biosynthesis. The most amazing gene is YLR328W. This gene is a key component in the pathway of NAD biosynthesis which produces NAD, an important molecule for the lysine biosynthesis. Therefore, all the 11 genes have a direct or indirect relationship with the lysine biosynthesis pathway. If taking only positive correlation for the study, many genes indirectly related to the pathway like those in the red group of AIE-Alpha-87 would be ignored. Therefore, through function analysis on AIE patterns, we can expand existing pathways by adding indirect biochemistry reactions or by linking to other biosynthesis pathways. This is the main reason why we say AIE pattern can enhance our understanding on molecular pathways.

Figure 8.
An expanded diagram for the pathway lysine biosynthesis after our functional annotation is used to derive new components (shown as dashed lines) based on the negative correlation of AIE-Alpha-87.

The second example of our case study is to show how a negative correlation of the genes in the same pathway can be used to identify gene targets in biological experiments for testing the negative mechanism.

Gene target identification for testing negative mechanism

We found that many of our AIE patterns consist of two groups of genes that are from the same pathways. One example is the representative AIE pattern in Cdc15 condition named AIE-Cdc15-273421, whose gene expression profile is shown in Figure 6b. Its P-value is 4.54E–12, R-value is about −0.55, up-down index is 1.0 and differential gap is about 0.068. A reason to choose this pattern for illustration is because of its extremely high up-down change frequency (i.e. 1.0). Of the total 14 genes in this pattern, there are six genes contained in the pathway ergosterol biosynthesis. Three of them [YML008C (GeneID:855003), YGL001C (GeneID:852883) and YLR100W (GeneID:850790)] are located at upstream, while the rest three [YMR202W (GeneID:855242), YGL012W (GeneID:852872) and YMR015C (GeneID:855029)] are at the downstream. According to a known structure of the pathway ergosterol biosynthesis partially shown in Figure 9, these six genes function sequentially one by one to achieve the final function of this pathway. Associating these genes' expressions with their roles in this pathway, it can be observed that when upstream genes are up- or down-regulated, the genes located at the downstream have simultaneous opposite expressions. Intuitively, we can infer that there exists some negative control mechanism in the path from fecosterol to episterol, which is the functional boundary between the upstream and downstream gene groups (Figure 9). Interestingly, it has been previously reported that ergosterol exerts a negative feedback on its own biosynthesis in S. cerevisiae, particularly at the C-24 methylation step involving the gene YML008C (34,35). Therefore, we can see that the negative correlation of the genes in the same pathway is potentially useful to identify gene targets in the biological experiments for testing negative mechanisms.

Figure 9.
A partial diagram of the pathway ergosterol biosynthesis.

Invariable negative correlations under different conditions

We present our third case study and examine hierarchical clusterings of the expression profiles of the genes in one AIE pattern under the three different cell cycle environments. Taking again AIE-Alpha-87 as an example, hierarchical clusterings of the expression profiles of the 11 genes in Alpha factor, Cdc15 and Cdc28 are shown in Figure 10. [The drawing was done by the software PermutMatrix (36).] In the case of Alpha factor, the 11 genes are nicely divided into two groups: one matches with the blue group, the other maps with the red group. However, the negative correlation under Alpha factor, disappeared in the other two environments of cell cycle. We can also see that the four known genes in the pathway lysine biosynthesis always have co-behaviors under different conditions of cell cycle, but the other genes' negative correlation cannot be always maintained.

Figure 10.
The expression profiles of the genes in AIE-Alpha-87 under three different conditions (Alpha factor, Cdc15 and Cdc28). The signs following each gene name have specific meanings. -B stands for this gene in the blue group; -BP stands for this gene in the ...

Thus, another perspective to understand AIE patterns is to see whether one pattern as a whole can be conserved during different biological environments. Of course, the above example is not the case. In fact, there is very little overlapping between any two AIE pattern pairs from different environments of cell cycle. This is because most environmental perturbations can cause big change in expression, resulting in alteration in the complex regulatory system. Ronen et al. (37) had designed a glucose plus experiment to detect affection of the small constraint perturbation, which is, on the other hand, helpful to understand stable negative correlations. We briefly explain two AIE patterns identified from their data set. One AIE pattern before glucose plus (denoted as AIE-StressBefore-4411) contains 31 genes whose expression profiles are shown at Figure 11a. An AIE pattern after glucose plus (denoted as AIE-StressAfter-3827) contains 30 genes whose expression profiles are shown at Figure 11c. Interestingly, these two patterns share six genes that maintain the same negative correlation. The six genes are: YDR050C (GeneID:851620), YJL200C (GeneID:853230), YCL018W (GeneID:850342) and YIL083C (GeneID:854726), YOR321W (GeneID:854499) and YOR136W (GeneID:854303). The expression profiles of this stable negative correlation are displayed in the Figure 11b and d. It looks that these six genes and their negative correlation are a stable genetic indicator during glucose plus. From the viewpoint of gene function, the other 25 genes in AIE-StressBefore-4411 have a significant function hit on transferase activity (10 out of 25 genes), while the other 24 genes in AIE-StressAfter-3827 have a significant function hit on oxidoreductase activity (10 out of 24 genes). It is also known that oxidoreductase is closely related to high-glucose ambience (38). Therefore, it is suggestive that this stable AIE pattern shared by these six genes is likely involved in both normal and stress activity of Yeast.

Figure 11.
An invariable negative correlation shared by six common genes of two AIE patterns (AIE-BeforeStress-4411 and AIE-AfterStress-3827).

CONCLUSION

The main contribution of this work is the formalization of the widely observed negative correlations in genes' functions within molecular pathways. Through our mining algorithm which uses a suffix-tree data structure and a bi-clique search idea, all possible AIE patterns in a time-course gene expression data set can be enumerated. As some of them are perhaps of less interest, we have suggested to use the size threshold, up-down index and R-value index to control the quantity and quality of AIE patterns in the post-analysis. Although pairs of gene clusters computed by the traditional clustering methods can find some negative correlations, they are unable to detect negative correlations shown in time segments as our AIE patterns can do. The biclustering methods can iteratively conduct clustering from both genes and time points, it is still hard to detect all negative correlation candidates in large data set. However, our mining algorithm can overcome this difficulty.

Our experimental results on three Yeast cell cycle expression data sets have demonstrated that maximal negative correlation can occur between pairs of large groups of genes, one group or both covering many genes from the same pathway. Based on existing knowledge about molecular pathways, negative correlations can be used to infer new gene regulation modules, TFs and RNs, as shown in our case studies. With these new elements, we are able to get a fuller and deeper picture about the direct and/or indirect relationships of all components in a molecular pathway. Besides, significant invariable negative correlations are found in both normal activity and stress activity of Yeast. All these ideas and results highlight that maximal negative correlation is an important characteristic in the gene expression profiles within pathways, which is expected to be useful in the cell stress response study, drug response evaluation, cancer-related pathways' detection, etc.

FUNDING

Tier-1 grant (RG66/07) awarded by Nanyang Technological University, Singapore. Funding for open access charge: Tier-1 grant (RG66/07) awarded by Nanyang Technological University, Singapore.

Conflict of interest statement. None declared.

REFERENCES

1. Segal E, Wang H, Koller D. Discovering molecular pathways from protein interaction and gene expression data. Bioinformatics. 2003;19(Suppl. 1):i264–i271. [PubMed]
2. Cho RJ, Campbell MJ, Winzeler EA, Steinmetz L, Conway A, Wodicka L, Wolfsberg TG, Gabrielian AE, Landsman D, Lockhart DJ, et al. A genome-wide transcriptional analysis of the mitotic cell cycle. Mol. Cell. 1998;2:65–73. [PubMed]
3. Madeira SC, Oliveira AL. Biclustering algorithms for biological data analysis: a survey. IEEE/ACM Trans. Computat. Biol. Bioinform. 2004;1:24–45. [PubMed]
4. Aguilar-Ruiz JS. Shifting and scaling patterns from gene expression data. Bioinformatics. 2005;21:3840–3845. [PubMed]
5. Segal E, Shapira M, Regev A, Pe’er D, Botstein D, Koller D, Friedman N. Module networks: identifying regulatory modules and their condition-specific regulators from gene expression data. Nat. Genet. 2003;34:166–176. [PubMed]
6. Cheng KO, Law NF, Siu WC, Liew AW. Identification of coherent patterns in gene expression data using an efficient biclustering algorithm and parallel coordinate visualization. BMC Bioinformatics. 2008;9:210–237. [PMC free article] [PubMed]
7. Kim J, Kim H. Clustering of change patterns using Fourier coefficients. Bioinformatics. 2008;24:184–191. [PubMed]
8. Schmid M, Davison TS, Henz SR, Pape UJ, Demar M, Vingron M, Scholkopf B, Weigel D, Lohmann JU. A gene expression map of arabidopsis thaliana development. Nat. Genet. 2005;37:501–506. [PubMed]
9. James CG, Appleton CT, Ulici V, Underhill TM, Beier F. Microarray analyses of gene expression during chondrocyte differentiation identifies novel regulators of hypertrophy. Mol. Biol. Cell. 2005;16:5316–5333. [PMC free article] [PubMed]
10. Chuang CL, Jen CH, Chen CM, Shieh GS. A pattern recognition approach to infer time-lagged genetic interactions. Bioinformatics. 2008;24:1183–1190. [PubMed]
11. Stekel DJ, Jenkins DJ. Strong negative self regulation of prokaryotic transcription factors increases the intrinsic noise of protein expression. BMC Syst. Biol. 2008;2:6–19. [PMC free article] [PubMed]
12. Missero C, Pirro MT, Di Lauro R. Multiple RAS downstream pathways mediate functional repression of the homeobox gene product ttf-1. Mol. Cell Biol. 2000;20:2783–2793. [PMC free article] [PubMed]
13. Segal L, Lapidot M, Solan Z, Ruppin E, Pilpel Y, Horn D. Nucleotide variation of regulatory motifs may lead to distinct expression patterns. Bioinformatics. 2007;23:i440–i449. [PubMed]
14. Shieh GS, Chen CM, Yu CY, Huang J, Wang WF, Lo YC. Inferring transcriptional compensation interactions in yeast via stepwise structure equation modeling. BMC Bioinformatics. 2008;9:134–143. [PMC free article] [PubMed]
15. Millar CB, Grunstein M. Genome-wide patterns of histone modifications in yeast. Nat. Rev. Mol. Cell Biol. 2006;7:657–666. [PubMed]
16. Esteller M. Cancer epigenomics: DNA methylomes and histone- modification maps. Nat. Rev. Genet. 2007;8:286–298. [PubMed]
17. Wiencke JK, Zheng S, Morrison Z, Yeh RF. Differentially expressed genes are marked by histone 3 lysine 9 trimethylation in human cancer cells. Oncogene. 2008;27:2412–2421. [PubMed]
18. McGarvey KM, Van Neste L, Cope L, Ohm JE, Herman JG, Van Criekinge W, Schuebel KE, Baylin SB. Defining a chromatin pattern that characterizes DNA-hypermethylated genes in colon cancer cells. Cancer Res. 2008;68:5753–5759. [PMC free article] [PubMed]
19. Yuan Y, Li CT, Wilson R. Partial mixture model for tight clustering of gene expression time-course. BMC Bioinformatics. 2008;9:287–303. [PMC free article] [PubMed]
20. Ji L, Tan KL. Identifying time-lagged gene clusters using gene expression data. Bioinformatics. 2005;21:509–516. [PubMed]
21. Supper J, Strauch M, Wanke D, Harter K, Zell A. Edisa: extracting biclusters from multiple time-series of gene expression profiles. BMC Bioinformatics. 2007;8:334–347. [PMC free article] [PubMed]
22. Madeira SC, Teixeira MC, Sá-Correia I, Oliveira AL. IEEE/ACM Trans. Comput. Biol. Bioinform. 2008. [(12 October 2009, date last accessed)]. Identification of regulatory modules in time series gene expression data using a linear time biclustering algorithm. 99, http://doi.ieeecomputersociety.org/10.1109/TCBB.2008.34. [PubMed]
23. Madeira SC, Oliveira AL. A polynomial time biclustering algorithm for finding approximate expression patterns in gene expression time series. Alg. Mol. Biol. 2009;4:8–46. [PMC free article] [PubMed]
24. Li J, Liu G, Li H, Wong L. Maximal biclique subgraphs and closed pattern pairs of the adjacency matrix: a one-to-one correspondence and mining algorithms. IEEE Trans. Knowl. Data Eng. 2007;19:1625–1637.
25. Spellman PT, Sherlock G, Zhang MQ, Iyer VR, Anders K, Eisen MB, Brown PO, Botstein D, Futcher B. Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. Mol. Biol. Cell. 1998;9:3273–3297. [PMC free article] [PubMed]
26. Cherry JM, Adler C, Ball C, Chervitz SA, Dwight SS, Hester ET, Jia Y, Juvik G, Roe T, Schroeder M, et al. SGD: Saccharomyces Genome Database. Nucleic Acids Res. 1998;26:73–79. [PMC free article] [PubMed]
27. Prelic A, Bleuler S, Zimmermann P, Wille A, Buhlmann P, Gruissem W, Hennig L, Thiele L, Zitzler E. A systematic comparison and evaluation of biclustering methods for gene expression data. Bioinformatics. 2006;22:1122–1129. [PubMed]
28. Gonalves JP, Madeira SC, Oliveira AL. Biggests: integrated environment for biclustering analysis of time series gene expression data. BMC Res. Notes. 2009;2:124–134. [PMC free article] [PubMed]
29. Curtis RK, Oresic M, Vidal-Puig A. Pathways to the analysis of microarray data. Trends Biotechnol. 2005;23:429–435. [PubMed]
30. Monteiro PT, Mendes ND, Teixeira MC, d'O;rey S, Tenreiro S, Mira NP, Pais H, Francisco AP, Carvalho AM, Lourenco AB, et al. Yeastract-discoverer: new tools to improve the analysis of transcriptional regulatory associations in Saccharomyces cerevisiae. Nucleic Acids Res. 2008;36:D132–D136. [PMC free article] [PubMed]
31. Schlitt T, Brazma A. Current approaches to gene regulatory network modelling. BMC Bioinformatics. 2007;8(Suppl. 6):S9. [PMC free article] [PubMed]
32. Guelzim N, Bottani S, Bourgine P, Képès F. Topological and causal structure of the yeast transcriptional regulatory network. Nat. Genet. 2002;31:60–63. [PubMed]
33. Wu WS, Li WH. Systematic identification of yeast cell cycle transcription factors using multiple data sources. BMC Bioinformatics. 2008;9:522–534. [PMC free article] [PubMed]
34. Veen M, Stahl U, Lang C. Combined overexpression of genes of the ergosterol biosynthetic pathway leads to accumulation of sterols in Saccharomyces cerevisiae. FEMS Yeast Res. 2003;4:87–95. [PubMed]
35. Vandeputte P, Tronchin G, Larcher G, Ernoult E, Bergès T, Chabasse D, Bouchara JP. A nonsense mutation in the erg6 gene leads to reduced susceptibility to polyenes in a clinical isolate of candida glabrata. Antimicrob. Agents Chemother. 2008;52:3701–3709. [PMC free article] [PubMed]
36. Caraux G, Pinloche S. Permutmatrix: a graphical environment to arrange gene expression profiles in optimal linear order. Bioinformatics. 2005;21:1280–1281. [PubMed]
37. Ronen M, Botstein D. Transcriptional response of steady-state yeast cultures to transient perturbations in carbon source. Proc. Natl Acad. Sci. USA. 2006;103:389–394. [PubMed]
38. Nayak B, Xie P, Akagi S, Yang Q, Sun L, Wada J, Thakur A, Danesh FR, Chugh SS, Kanwar YS. Modulation of renal-specific oxidoreductase/myo-inositol oxygenase by high-glucose ambience. Proc. Natl Acad. Sci. USA. 2005;102:17952–17957. [PubMed]

Articles from Nucleic Acids Research are provided here courtesy of Oxford University Press