CoXpress: differential co-expression in gene expression data

Microarrays have become a standard tool for the exploration of global gene expression changes at the cellular level [1]. Data analysis often includes the use of a statistical test, such as a t-test or analysis of variance, to find genes differentially expressed in one set of conditions when compared to another, or the use of clustering algorithms in order to find groups of genes which behave similarly over a number of experiments [2]. However, these techniques may not detect differential co-expression patterns that exist between two biological states.

Statistical tests, such as the t-test or ANOVA, identify genes that are differentially expressed under one or more conditions. The output of such tests is a simple list of genes, with an associated test statistic and p-value [3]. There is no indication of which genes may be interacting with one another. Alternatively, clustering algorithms are often used to find groups of genes which display similar expression profiles across a dataset, and these clusters are subsequently analyzed visually for patterns of interest [4, 5]. Eisen et al used hierarchical cluster analysis to determine groups of co-expressed genes, and found that genes within those groups were functionally related [4], and the use of hierarchical cluster analysis is now a standard technique for analysing microarray data [2, 6]. Yeung et al [7] assessed the use of hierarchical clustering to find groups of co-regulated genes. Various clustering algorithms were used on a number of datasets, and the results evaluated by determining those genes that share a common transcription factor. Of the algorithms tested, MCLUST [8] and two hierarchical methods (based on the pearson correlation coefficient) showed the highest coincidence of correlated and co-regulated genes.

However, genes which show highly correlated patterns of expression in one biological state, but not in another, may not be highly correlated across the entire dataset, and therefore would not be associated with one another if a clustering algorithm is used. Variation may exist in the expression of a gene in different groups of individuals due to the presence of sub-populations, and this may lead to that gene being grouped incorrectly. Furthermore, clustering algorithms do not provide methods to identify groups that are behaving differently in different biological conditions.

Recent work has concentrated on alternative approaches to the discovery of co-expressed genes. Li [9] describes a method whereby genes whose expression is associated with differential co-expression patterns in other pairs of genes may be discovered, and Lai et al [10] describe a conceptually similar method whereby pairs of genes that display differential co-expression patterns between the normal and cancerous state may be discovered. Other approaches have centred on the construction of large gene co-expression networks. Lee et al [11] analysed 60 human microarray data sets to construct gene co-expression networks conserved across multiple data sets, and Stuart et al [12] constructed a gene co-expression network across different organisms, indicating that such relationships are evolutionarily conserved. However, neither of the above attempted to find group of genes differentially co-expressed between different conditions. Choi et al [13] tackled this problem by constructing normal and tumour co-expression networks from a variety of public datasets, comparing the results to find differences in co-expression patterns associated with cancer. In all of these cases, the networks were built by comparing genes pairwise, using some variation of the pearson correlation coefficient, to determine if a co-expression link exists between the two genes. These links were then joined to form a co-expression network.

Cluster analysis and network construction can be thought of as alternative methods for finding co-expressed genes. Whereas networks concentrate on conserved, pairwise comparisons, there is no guarantee that genes that are close in the network, but are not directly linked, have correlated expression profiles. Alternatively, cluster analysis produces groups of genes that are correlated above a certain level, defined by where the tree is cut and the clustering algorithm, but there is no indication of which particular pairs of genes are interacting. Kostka and Spang [14] described the first method to investigate differentially co-expressed groups of genes, using an additive model for scoring gene-gene co-expression and then a stochastic search algorithm to find groups of genes showing differential co-expression patterns. Jen et al [15] have produced ACT, the Arabidopsis Co-expression Tool, which allows users to calculate co-expression across user-defined data sets and uses a correlation cut-off to define groups of genes.

Here we describe coXpress, a simple and easy to use package that allows users to explore differential co-expression in an intuitive way. The package is aimed at biologists who want to analyse differential co-expression in their data set, which can be achieved in just 4 simple commands once the data has been loaded. CoXpress uses hierarchical cluster analysis to explore the relationship between genes, cutting the tree to form groups of genes that are co-expressed. This is an intuitive approach that many biologists are familiar with. CoXpress then uses a resampling approach to find those groups that are co-expressed in one set of experiments and not in another. The package should be used as first step in the analysis of co-expression, and is designed to complement the approaches described above.

CoXpress is released as a package for R. R is a freely available, open-source statistical package [16] that is widely used in the biological community. R has very powerful statistical and graphical capabilities, and many add-on packages are freely available. The bioconductor project [17, 18] provides a huge number of add-on packages for R, covering a wide range of biological data analysis applications, and the implementation of coXpress in R provides seamless integration with many of these packages. CoXpress is written in the native R language and has been fully tested on both windows and linux. R is available for windows, linux, unix and MacOS (including MacOS X).

The input for coXpress is a matrix of data, with rows representing genes and columns representing microarrays. The R data.frame object is most convenient, and can be created by reading in a text file (using the read.table function), an Excel spreadsheet (using the RODBC library) or from existing R objects, created by the packages from the bioconductor project such as affy [19], limma [20] or marray [18].

The genes are first clustered based on their expression values in a subset of experiments (termed subset 1), using the cluster.gene function. This function wraps the dist, cor and hclust functions that are built in to R, and thus provides a simple interface to hierarchical clustering. When a correlation coefficient is used as the distance measure, the distance measure is calculated as 1 - r, where r is the pearson correlation coefficient. The resulting tree is cut at a user-defined value, using the cutree function, to form groups of genes that are co-expressed in subset 1. These groups are then examined in both subset 1 and a second set of experiments, defined by the user, which we will term subset 2.

Groups of size 1 are ignored as there can be no co-expression. Groups of size two are handled by the cox.pairs function. The cox.pairs function uses the cor.test function in R to test if the genes are significantly correlated in subset 1 and subset 2. Thus, a pair of genes significantly correlated in subset 1 and not significantly correlated in subset 2 can be described as differentially co-expressed.

Groups with more than two members are handled by the coXpress function. The flow of analysis in coXpress is represented in figure 1. For each group of size n, where n ≥ 3, the pairwise correlation coefficients of the group in subset 1 are calculated. These are then summarised using the t-statistic, the use of which is discussed below. Then, m random groups of size n are created by randomly re-sampling the data matrix. For each of these random groups, the pairwise correlation coefficients of the group in subset 1 are calculated and again summarised using the t statistic. Thus, a distribution of t statistics is created, of size m, from randomly assigned groups of size n. The observed t statistic is then compared to the distribution of random t statistics. The proportion of random statistics greater than the observed is used as a "probability of randomness" for the group in subset 1. This process is then repeated for subset 2. A group which is found to be non-random in group 1 and random in group 2 is said to be differentially co-expressed. These groups will be highly correlated in subset 1 but show little or no correlation in subset 2. To find the reverse, the process must be repeated, but basing the original groups on a cluster analysis of the data based on subset 2.

Conceptual flow of analysis through the coXpress function. This figure shows the conceptual steps taken by the coXpress function in order to find differentially co-expressed groups of genes in microarray data. Steps should be followed in numerical order.

Full size image

The t-statistic is used here not as a test of significance, but as a means of summarising a set of pairwise correlation coefficients into a single value. Correlation coefficients are on the scale:

• 1 ≤ r ≤ +1

where 1 represents positive correlation, -1 represents negative correlation and 0 represents no relationship. The t-statistic is used here to summarise the "difference from zero" of a group of pairwise correlation coefficients. The exact formula for this is:

where x is the vector of unique, pairwise correlation coefficients, $\overline{x}$ is the mean of x and se(x) denotes the standard error of x. A group of highly correlated genes will have a mean correlation close to 1 and a small standard error, resulting in a large value for t. However, a group of uncorrelated genes will have a mean close to 0 and a relatively large standard error, resulting in a small value for t. The observed t statistic is compared against m random t statistics in order to calculate a probability of randomness.

The AML/ALL leukaemia dataset

The utility of coXpress is demonstrated using gene expression data from the leukaemia microarray study of Golub et al [21]. This dataset represents gene expression measurements from 38 tumour mRNA samples, 27 acute lymphoblastic leukaemia (ALL) cases and 11 acute myeloid leukaemia (AML) cases. The HU6800 Affymetrix array was used, which contains 6800 probesets. The dataset has been filtered such that genes with negative values in any sample have been removed, resulting in 2568 genes.

Using coXpress, the genes were first clustered according to their expression levels in the 27 ALL samples, using the cluster.gene function. The distance measure used was 1 - r, where r is the pearson correlation coefficient. The resulting tree was cut at a distance of 0.4, representing a correlation coefficient of 0.6, using the cutree function.

These groups were then examined in both the ALL and AML cases using the coXpress function. The observed t statistics in all cases were compared with the t statistics generated by randomly resampling the dataset 10,000 times for each group size. The resulting table contains one row for each group.

To test the robustness of the method to outliers, a bootstrapping approach was used. Each group was re-tested 1000 times, each time randomly selecting 75% of the observations for each leukaemia subtype (20 AML cases and 8 AML cases). The number of times each group was found to be differentially co-expressed by the coXpress method was recorded.

Table 1 shows the results filtered for groups that are non-random in the ALL subset, random in the AML subset, and with more than 6 members. As can be seen, there are 10 groups, varying in size from 7 to 34 members. The mean pairwise correlations for the groups are all above 0.6 in the ALL cases, yet show little or no correlation in the AML cases, with mean values ranging from -0.093 to 0.144. The robustness resampling method provides evidence that the groups found are robust to outliers, with nine out of ten groups being found in over 90% of the resampled data sets, and the other being found in 76%.

Figure 2 demonstrates the method of coXpress. These graphs show data from the largest of the groups, group 3, which has 34 members. Fig. 2A compares the distribution of pairwise correlation coefficients in the ALL subset with two random distributions. The blue graph is the distribution of observed correlation coefficients in the ALL subset for group 3, the red graph is the distribution of pairwise correlation coefficients from data generated by the random uniform distribution, and the green graph is the distribution of pairwise correlation coefficients from a group of genes randomly selected from the dataset. As can be seen, the observed distribution for this group in the ALL subset is very different from the two random distributions. Fig. 2B is an identical graph for the group based on the AML subset. This time, the observed distribution shows no difference compared to the two random distributions. The t-statistics for each distribution are shown on these graphs. Fig. 2C shows the observed t-statistic for group 3 in the ALL subset compared to the distribution of 10,000 randomly generated t-statistics, and Fig. 2D is the equivalent graph for the AML subset. Again, it is clear that this group in the ALL subset is non-random, yet is no different to random in the AML subset.

A differentially co-expressed group from the Golub dataset compared to random distributions. Group 3 (n = 34) from the Golub [21] data compared to random distributions. A) the distribution of pairwise correlation coefficients for group 3 in the ALL subset (blue) compared to the distribution of pairwise correlation coefficients from a group of the same size generated by the random uniform distribution (red) and the distribution of pairwise correlation coefficients from a group of the same size randomly selected from the dataset (green). B) Equivalent graph to A for the AML subset. C) A comparison of the observed t-statistic for group 3 in the ALL subset with a distribution of 10,000 random t-statistics generated by randomly resampling the dataset. D) Equivalent graph to C for the AML subset. Distributions were smoothed and drawn using density function in R [16]. Note that graph C has a broken x-axis.

Full size image

Figure 3 shows the top 3 groups in table 1 graphically. Fig. 3A is the largest of the groups, with 34 members. These 34 genes have a mean pairwise correlation of 0.70 in the ALL subset, but only 0.003 in the AML subset. Fig. 3B shows a smaller group, with 7 members, with a mean pairwise correlation of 0.72 in the ALL subset and -0.09 in the AML subset. Finally, fig. 3C shows a group with 11 members, with a mean pairwise correlation coefficient of 0.679 in the ALL subset and only 0.086 in the AML subset. These graphs were produced using the plot.compare.group and plot.cluster.genes functions.

Expression profiles for three differentially co-expressed groups in the Golub dataset. Expression profiles for three groups of differentially co-expressed genes from the Golub dataset [21]. A) Group 3 (n = 34) in 27 ALL samples (left) and 11 AML samples (right). B) Group 62 (n = 7) in 27 ALL samples (left) and 11 AML samples (right). C) Group 121 (n = 12) in 27 ALL samples (left) and 11 AML samples (right). Expression levels have been scaled and centred.

Full size image

Figure 4 shows the same three groups in a different way. Here, each plot is a representation of the correlation matrix of the group of genes in either the ALL or the AML subsets. Each coefficient in the correlation matrix is represented as a square, with the colour of the square representing the amount of correlation. The colour scale used is green to red, with green representing -1 (negative correlation), red representing +1 (positive correlation) and black representing 0 (no correlation). In all three groups, the correlation matrices are red for the ALL subset, yet are a mixture of black, green and red in the AML subset. This view of the data is more useful than simply considering the average pairwise correlation, as it shows all of the values in an intuitive way. These graphs were produced using the show.cor.matrices function.

Images of the correlation matrices for three differentially co-expressed groups in the Golub dataset. Images of the correlation matrices for three groups of differentially co-expressed genes from the Golub dataset [21]. A) Group 3 (n = 34) in 27 ALL samples (left) and 11 AML samples (right). B) Group 62 (n = 7) in 27 ALL samples (left) and 11 AML samples (right). C) Group 121 (n = 12) in 27 ALL samples (left) and 11 AML samples (right). Each coefficient in the correlation matrix is represented as a square, with the colour of the square representing the amount of correlation. The colour scale used is green to red, with green representing -1, red representing +1 and black representing 0.

Full size image

In each of the differentially co-expressed groups, not all pairwise correlation coefficients will have decreased or changed. To examine which pairs of genes have changed, the inspect.group function should be used. Table 2 shows the ten pairwise correlation coefficients that have changed the most between the ALL and AML subsets in group 62. As can be seen, these pairs of genes are all positively correlated in the ALL subset but are negatively correlated in the AML subset. Table 3 shows the ten pairwise correlation coefficients that have changed the least between the ALL and AML subsets in group 62. Many of these pairs of genes are still positively correlated in the AML subset, but not to the same extent. It is important that each differentially co-expressed group is examined in this way to determine which of the pairs of correlated genes have changed and which have not.

The GOHyperG function of the GOstats package [22] was used to find GO terms over-represented in the differentially co-expressed groups. Group 3, with 34 members, is enriched for GO terms for lymph node development, cell organisation and biogenesis, and protein biosynthesis and transport. Group 62, which has 7 members, is enriched for GO terms for methyltransferase activity, DNA modification, protein transport and DNA and protein methylation. Group 121, with 11 members, is enriched for GO terms for nucleotidase activity, and RNA splicing, processing and metabolism.

The ALL subtype dataset

This dataset is from the Acute Lymphoblastic Leukaemia study by Yeoh et al [23]. Six subtypes of ALL leukaemias are represented in 248 cases. The six subtypes are T-ALL, E2A-PBX1, BCR-ABL, TEL-AML1, MLL rearrangement, and hyperdiploid >50. The HG_U95Av2 Affymetrix microarray was used which contains 12,600 probesets. The dataset has been filtered such that genes with negative values in any sample have been removed, resulting in 1516 genes present in the dataset.

Using coXpress, the genes were first clustered according to their expression levels in the BCR-ABL samples, using the cluster.gene function. The distance measure used was 1 - r, where r is the pearson correlation coefficient. The resulting tree was cut at a distance of 0.5, representing a correlation coefficient of 0.5, using the cutree function. These groups were then examined in both the BCR-ABL and T-ALL subsets.

Those groups of size two were analysed using the cox.pairs function. Table 4 shows three pairs of genes that are significantly positively correlated in the BCR-ABL subset, and significantly negatively correlated in the T-ALL subset.

Groups of N ≥ 3 were analysed in the BCR-ABL and T-ALL subsets using the coXpress function. The observed t statistics in all cases were compared with the t statistics generated by randomly resampling the dataset 10,000 times for each group size. The resulting table contains one row for each group.

To test the robustness of the method to outliers, a bootstrapping approach was used. Each group was re-tested 1000 times, each time randomly selecting 75% of the observations for each leukaemia subtype. The number of times each group was found to be differentially co-expressed by the coXpress method was recorded.

Table 5 shows the results filtered for groups that are non-random in the BCR-ABL cases, random in the T-ALL cases, and with more than 10 members. Figure 5 shows the top 3 groups in table 1 graphically. Figure 5A shows a group of 16 genes that have a mean pairwise correlation coefficient of 0.669 in the BCR-ABL subset, yet only 0.06 in the T-ALL subset. Figure 5B shows a group of 10 genes that have a mean correlation of 0.65 in the BCR-ABL subset and only 0.08 in the T-ALL data. Finally, Figure 5C shows a group of 13 genes that have an average correlation of 0.64 in the BCR-ABL data, yet only 0.04 in the T-ALL data. The robustness resampling method provides evidence that the groups found are robust to outliers, with twelve out of thirteen groups being found in over 80% of the resampled data sets, and the other being found in 68.6%.

Expression profiles for three differentially co-expressed groups in the ALL subtype dataset. Expression profiles for three groups of differentially co-expressed genes from the Yeoh et al dataset [23]. A) Group 47 (n = 16) in 15 BCR-ABL samples (left) and 43 T-ALL1 samples (right). B) Group 31 (n = 10) in 15 BCR-ABL samples (left) and 43 T-ALL1 samples (right). C) Group 89 (n = 13) in 15 BCR-ABL samples (left) and 43 T-ALL1 samples (right). Expression levels have been scaled and centred.

Full size image

The GOHyperG function of the GOstats package [22] was used to find GO terms over-represented in the differentially co-expressed groups. Group 47 with 16 members, is enriched for GO terms for hormone catabolism, glucocorticoid receptor signalling and glucocorticoid catabolism. Group 31 with 10 members contains two probes for a gene in the RAS oncogene family, and is enriched for GO terms for oxidoreductase activity and ubiquitin activating enzyme activity. Finally, group 89 with 13 members contains genes annotated as B-cell lymphoma and cancer susceptibility genes, as well as genes enriched for GO terms for endothelial cell migration, regulation of cell motility and migration, angiostatin binding and regulation of blood vessel endothelial cell migration.

It is clear that coXpress is capable of finding differentially co-expressed groups of genes in both data sets. The groups presented above are extremely highly correlated in one subset of experiments, yet show little correlation in another subset. Furthermore, these patterns of correlation are shown to be non-random in the first subset, and no different from random in the second subset. The results show that it is the overall correlation structures of these groups that have changed significantly and some pairs of genes are still highly correlated in the second subset. It is important that each group is examined using the inspect.groups function in order to determine which of the pairs of genes are still correlated and which are not. The groups found by coXpress could also feed into the network construction technique described by Choi et al [13] to determine which pairwise relationships are conserved and which are not. One would expect the differences between ALL and AML leukaemia in the Golub dataset to be larger than those between different ALL subtypes in the Yeoh dataset, and the fact that coXpress can still find groups with such different correlation structures demonstrates the power of the method.

The use of hierarchical cluster analysis, followed by cutting the tree, is an intuitive approach and one that is familiar to biologists. However, it has limitations. For example, each gene may only be in one group, which does not ring true for biological systems, where many genes have multiple functions. Also, the choice of where to cut the tree is arbitrary. A high cut-off will produce many small groups of genes that are very highly correlated, whereas a lower cut-off will produce fewer groups, of larger size, which are not as highly correlated. In reality the user must use a range of different cut-offs to see which performs best with their dataset. Other clustering algorithms, such as MCLUST [7, 8], have been shown to out-perform hierarchical cluster analysis, however, there is no reason why these algorithms could not be used to define the groups of genes prior to running the coXpress function.

There are several directions in which the software can be developed. At present, the user defines which subsets of experiments are analysed, however it is possible that coXpress could suggest, or improve, these groupings using an approach such as random forests or genetic algorithms. This may allow researchers to discover sub-populations in the system under study. Integration of other clustering algorithms with coXpress, such as MCLUST, may also improve the performance of the software. In particular, clustering or grouping algorithms that allow genes to be present in more than one group may be advantageous. Finally, the integration of network construction algorithms would allow researchers to further analyse and visualise the differentially co-expressed groups discovered by coXpress.

We describe coXpress, an open-source R package that allows researchers to analyse differential co-expression patterns in DNA microarray data. CoXpress contains several methods for the discovery and visualisation of differentially co-expressed genes. We show how coXpress can be used to find groups of differentially co-expressed genes in two publicly available microarray datasets. The groups found are shown to be highly correlated in one subset of experiments, yet show little or no correlation in a second subset of experiments. A comparison against random distributions is used to obtain a p-value for the co-expression of the genes in different subsets.

• Project Name: coXpress

• Project Home Page: http://coxpress.sf.net

• Operating Systems: Windows, Linux

• Programming Language: R

• Other Requirements: R, gplots, gtools, gdata (for heatmaps), hu6800, hgu95av2, plotrix (for examples)

• License: GNU GPL