Email updates

Keep up to date with the latest news and content from BMC Bioinformatics and BioMed Central.

Open Access Methodology article

Simple integrative preprocessing preserves what is shared in data sources

Abhishek Tripathi12*, Arto Klami23 and Samuel Kaski23

Author affiliations

1 Department of Computer Science, P.O. Box 68, FI-00014, University of Helsinki, Finland

2 Helsinki Institute for Information Technology, Finland

3 Department of Information and Computer Science, Helsinki University of Technology, P.O. Box 5400, FI-02015 HUT, Finland

For all author emails, please log on.

Citation and License

BMC Bioinformatics 2008, 9:111  doi:10.1186/1471-2105-9-111


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2105/9/111


Received:13 November 2007
Accepted:21 February 2008
Published:21 February 2008

© 2008 Tripathi et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Bioinformatics data analysis toolbox needs general-purpose, fast and easily interpretable preprocessing tools that perform data integration during exploratory data analysis. Our focus is on vector-valued data sources, each consisting of measurements of the same entity but on different variables, and on tasks where source-specific variation is considered noisy or not interesting. Principal components analysis of all sources combined together is an obvious choice if it is not important to distinguish between data source-specific and shared variation. Canonical Correlation Analysis (CCA) focuses on mutual dependencies and discards source-specific "noise" but it produces a separate set of components for each source.

Results

It turns out that components given by CCA can be combined easily to produce a linear and hence fast and easily interpretable feature extraction method. The method fuses together several sources, such that the properties they share are preserved. Source-specific variation is discarded as uninteresting. We give the details and implement them in a software tool. The method is demonstrated on gene expression measurements in three case studies: classification of cell cycle regulated genes in yeast, identification of differentially expressed genes in leukemia, and defining stress response in yeast. The software package is available at http://www.cis.hut.fi/projects/mi/software/drCCA/ webcite.

Conclusion

We introduced a method for the task of data fusion for exploratory data analysis, when statistical dependencies between the sources and not within a source are interesting. The method uses canonical correlation analysis in a new way for dimensionality reduction, and inherits its good properties of being simple, fast, and easily interpretable as a linear projection.

Background

Combining evidence from several heterogeneous data sources is a central operation in computational systems biology. We assume several vector-valued data sources, such that each source consists of measurements from the same object or entity, but on different variables.

In modeling in general, when it is possible to make sufficiently detailed modeling assumptions, data integration is in principle straightforward. Given a statistical model of how transcriptional regulation works, for instance, the Bayesian framework tells how to integrate gene expression data, prior knowledge, and transcription factor finding data. Lots of practical problems of course remain to be solved. Alternatively, in a classification task of proteins to ribosomal or membrane proteins, for instance, integration is likewise straightforward: do the integration such that the classification accuracy is maximized. This has been done effectively in semidefinite programming for kernel methods [1] and using Gaussian Process prior within the Bayesian framework [2].

In exploratory analysis, that is, when "looking at the data" to start data analysis while the hypotheses are still vague, it is not as straightforward to decide how data sources should be integrated. The task of exploring data is particularly important for the current high-throughput data sources, to be able to spot measurement errors and obvious deviations from what was expected of the data, and to construct hypotheses about the nature of the data. Nowadays in bioinformatics applications this stage is typically done using dimensionality reduction and information visualization methods, and clusterings. A good exploratory analysis method is (i) fast to apply interactively, (ii) easily interpretable by the analyst, and (iii) widely applicable. Linear projection methods, as such or as preprocessing for clusterings and other methods, fulfill all these criteria.

Fusing the sources is not trivial since we need to choose from three very different options. If all sources are equally important and there is not special reason to do otherwise, it makes sense to simply concatenate the variables from all sources together, and then continue with the resulting single source. The classical linear preprocessing method for this case is Principal Component Analysis (PCA). The second option is suitable when one of the sources, such as the class indicator in functional classification tasks, is known to be of the most interest. Then it is best to include only those variables or features within each source that are informative of the class variable. A classical linear method applicable in this case is linear discriminant analysis. This second option is supervised, and only applicable when the class information is available.

The third option is to include only those aspects of each source that are mutually informative of each other. Those are the shared aspects, and this task can be motivated through two interrelated lines of thought. The first is noise reduction. If the sources are measurements of the same entity corrupted by independent noise, then discarding the source-specific aspects will discard the noise. The second line of motivation is more abstract, to analyze what is interesting in the data. Here the different measurement sources can convey very different kinds of information of the entities being studied. One example is copy number aberrations and expression measurements of the same genes in cancer studies [3], and another is the activation profiles of the same yeast gene in several stressful treatments in the task of defining yeast stress response [4]. In these examples it is what is in common in the sources that we are really interested in. Note that the "noise" may be very structured; its effective definition is that it is source-specific.

Commonalities in data sources have been studied by methods that search for statistical dependencies between them. The earliest method was classical linear Canonical Correlation Analysis (CCA) [5], which has later been extended to nonlinear variants and more general methods that maximize mutual information instead of correlation. Yet, being fast, simple and easily understandable, the linear CCA still has a special place in the data analysis toolbox, analogously to the linear Principal Component Analysis which is still being used heavily instead of all modern dimensionality reduction and factor analysis techniques.

CCA addresses the right problem, searching for commonalities in the data sources. Moreover, being based on eigenvalue analysis it is fast and its results are interpretable as linear correlated components. It is not directly usable as a data fusion tool, however, since it produces separate components and hence separate preprocessing for each source. If the separate outputs could be combined in a way that is both intuitively interpretable and rigorous, the resulting method could become a widely applicable dimensionality reduction tool, analogously to PCA for a single source. Performing dimensionality reduction helps in avoiding overfitting, focusing on the most important effects, and reduces computational cost of subsequent analysis tasks.

In this paper we turn CCA into a data fusion tool by showing that the justified way of combining the sources is simply to sum together the corresponding CCA components from each source. An alternative view to this procedure is that it is equivalent to whitening each data source separately, and then running standard PCA on their combination. This is one of the standard ways of computing CCA, but for CCA the eigenvectors are finally split into parts corresponding to the sources. So the connection to CCA is almost trivial and it is amazing that, as far as we know, it has not been utilized earlier in this way.

Our contribution in this paper is to point out that CCA can be used to build a general-purpose preprocessing or feature extraction method, which is fast, and easily interpretable. There are two alternative interpretations. The first is the connection to CCA discussed above. The second is that it extends the standard practice of standardizing the mean and variance of each variable separately before dimensionality reduction. Now each data source is standardized instead of each variable.

We have developed a practical software tool for R that incorporates the subtle but crucial choices that need to be made to choose the dimensionality of the solution. The method is demonstrated on three collections of gene expression measurements.

A kernelized version of CCA (KCCA) has been used in specific data fusion tasks (see e.g. [6,7]) and it could be easily extended to be used in the same way as the linear CCA here. We will focus on the linear mappings for two practical reasons: Computation of the linear version is fast and the components are more easily interpretable. In particular, the kernelized version does not reveal which of the original features cause the dependencies between sources.

Results and Discussion

Algorithm

In this section we first explain a simple two-step procedure, based on whitening and PCA, for finding the aspects shared by the sources, and then show how the same fusion solution can equivalently be derived from the result of applying a generalized CCA to the collection. The two-step procedure provides the intuition for the approach: First remove the within-data variation, and then capture all the variation that is still remaining. The connection to CCA then demonstrates how the procedure provides a solution to the issue of combining the separate components CCA gives.

Denote a collection of p data sets by {X1,...,Xp}, where each Xi is a m × ni matrix such that m N, and N = ∑ni. The rows of the matrices correspond to the same object in each set, while the columns correspond to features that need not be the same in the data sets. For example, in traditional expression analyses the rows would be genes and the columns would be conditions, treatments, time points, etc. For notational simplicity, we assume zero mean data.

In the first step, each data set is whitened to remove all within-data correlations, and the data are scaled so that all dimensions have equal variance. The whitened version <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M1">View MathML</a> of a data matrix Xi is given by <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M2">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M3">View MathML</a> is the whitening matrix. The whitening matrix is simply <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M4">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M5">View MathML</a> is the covariance matrix of Xi.

After each data set has been whitened, the next step is to find the shared variation in them. This is done by principal component analysis (PCA) on the columnwise concatenated whitened data sets. Since all the within-data structure PCA could extract has been removed, it can only find variation shared by at least two of the data sets, and the maximum variance directions it searches for correspond to the highest between-data correlations.

Formally, applying PCA to the columnwise concatenation of the whitened data sets <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M6">View MathML</a> yields the factorization

CZ = V Λ VT,(1)

where the orthonormal matrix V contains the eigenvectors, Λ is a diagonal matrix of projection variances, and CZ is the covariance matrix of Z.

Projecting Z onto the first d eigenvectors Vd corresponding to the d largest eigenvalues gives the d principal components, which are the optimal d-dimensional representation in terms of the shared variance. The whole data collection becomes integrated into

Pd = ZVd,(2)

where Pd is of size m × d and contains a d-dimensional feature vector for each of the analyzed objects. The idea is then simply to use this new representation for any further analysis, which can be made using any method that operates on vectorial data. The whole procedure can be seen as fusing the collection of data sets into a single set, while at the same time reducing the total dimensionality of the data to find the most reliable shared effects.

As mentioned in the Background section, the above two-step procedure is equivalent to running CCA on the collection and summing the separate components from each source. The connection is derived here for two data sets. The proof extends easily for several data sets, for one of the many alternative generalizations of CCA.

CCA is a method for finding linear projections of two sets of variables so that the correlation between the projections is maximal. CCA is often formulated as a generalized eigenvalue problem

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M7">View MathML</a>

(3)

where Cij denotes the (cross-)covariance of Xi and Xj. The eigenvalues λ of the solution appear as pairs 1 + ρ1, 1 -ρ1,...,1 + ρp, 1 - ρp, 1,...,1, where p = min(n1, n2), and (ρ1,...,ρp) are the canonical correlations. The canonical weights corresponding to the canonical correlations are <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M8">View MathML</a>, i = 1,...,p.

In conventional use of CCA we are usually interested in the correlations, the canonical weights ui, and the canonical scores, defined as projections of X1 and X2 on the corresponding canonical weights. Next we show how the combined data set (2) can be obtained from the canonical scores, thus providing a way of using CCA to find a single representation that captures the dependencies.

For a single component, (1) can be equivalently written as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M9">View MathML</a>

where α is the variance, v is the corresponding principal component, and <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M10">View MathML</a> denotes the (cross-)covariance of <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M1">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M11">View MathML</a>. Due to the whitening, the <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M12">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M13">View MathML</a> are identity matrices. We can alternatively write <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M14">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M15">View MathML</a>, leading to

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M16">View MathML</a>

where Iv has been subtracted from both sides. Equivalently,

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M17">View MathML</a>

(4)

Let us denote <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M18">View MathML</a> by diag [W1, W2], and multiply the right hand side of (4) by the identity matrix I = diag[W1, W2]-1diag [W1, W2]. On the right side of the equation we then have the term diag[W1T, W2T]-1 diag[W1, W2]-1 = diag [C11, C22] based on the definition of the whitening matrix, and thus (4) can be written as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M19">View MathML</a>

(5)

where <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M20">View MathML</a>. Adding <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M21">View MathML</a> to both sides gives equation structurally identical to (3). Both methods thus lead to the same eigenvalues, i.e. λ = α, and the eigenvectors are related by a linear transformation,

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M22">View MathML</a>

The combined representation (2) of d dimensions can be written in terms of canonical scores as Pd = ZVd = [X1, X2]diag[W1, W2]Vd = [X1, X2] <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M23">View MathML</a> = X1U1,d + X2U2,d, where U1,d and U2,d are the first d canonical directions of the two data sets.

CCA can be generalized to more than two data sets in several ways [8], and the two-step procedure described here is equivalent to the one formulated as solving a generalized eigenvector problem Cu = λDu, where C is the covariance matrix of the column-wise concatenation of the Xi and D is a block-diagonal matrix having the dataset-specific covariance matrices Cii on its diagonal. Here u is a row-wise concatenation of the canonical weights corresponding to the different data sets. The proof follows along the same lines as for two data sets, and again the combined data set for any <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M24">View MathML</a> dimensions can be written in terms of the generalized CCA results as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M25">View MathML</a>

where each Ui,d contains the d eigenvectors corresponding to the d largest eigenvalues.

In summary, the simple linear preprocessing method of whitening followed by PCA equals computing the generalized CCA on a collection of data sets and summing the canonical scores of the data sets. In practice it does not matter in which way the result is obtained, but the two-step procedure illustrates more clearly why this kind of approach is useful for data integration. Furthermore, it is not limited to linear projections, and the same motivation could be extended to different kind of models. In practice implementing the first step might, however, be difficult in more complex models.

Choice of dimensionality

The dimensionality of the projection can be chosen to be fixed, such as two or three for visualization, or alternatively an "optimal" dimensionality can be sought. In this section we introduce our suggestion for optimizing the dimensionality. Intuitively, the dimensionality should be high enough to preserve most of the shared variation and yet low enough to avoid overfitting. The first few components contain most of the reliable shared variation among the data sets, while the last components may actually represent just noise, and thus dropping some of the dimensions makes the method more robust.

The maximum dimensionality is the sum of the dimensionalities of the data sets, but in practice already a considerably smaller dimensionality is often sufficient, and in fact leads to a better representation due to suppression of noise. Note also that in the case of two data sets the number of unique projections is only the minimum of the data dimensionalities.

In a nutshell, we increase the dimensionality one at a time, testing with a randomization test that the new dimension captures shared variation. To protect from overfitting, all estimates of captured variation will be computed using a validation set, i.e., for data that has not been used when computing the components (dimensions). The randomization test essentially compares the shared variance along the new dimension to the shared variance we would get under the null-hypothesis of mutual independency. When the shared variance does not differ significantly from the null-hypothesis, the final dimensionality has been reached.

To compute the shared variance of the original data, we divide the data into training, <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M26">View MathML</a> and validation data, <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M27">View MathML</a>. The two step procedure described in the Algorithm subsection is applied to the training data to compute the eigenvectors Vt and the whitening matrix Wt, where Wt is a block diagonal matrix containing the whitening matrices for each matrix in training data. The fused representation for the validation data is computed as <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M28">View MathML</a>, where Xv is the columnwise concatenation of the validation data matrices. Variance in the fused representation is now our estimate of shared variance. We average the estimate over 3 different splits into training and validation sets.

To compute the shared variance under the null hypothesis, random data sets are created from the multivariate normal distribution with a diagonal covariance matrix where the values in diagonal equal the columnwise variances of <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M26">View MathML</a>. The shared variance for the random data is computed in the same way as described above. We repeat the process for 50 randomly created data sets.

The shared variance in the original data is then compared to the distribution of shared variances under the null hypothesis, starting from the first dimension. When the dimensions no longer differ significantly (we used 2% confidence level), we have arrived at the "optimal" dimensionality and the rest of the dimensions are discarded.

Note that assuming normally distributed data in the null hypothesis is consistent with the assumptions implicitly made by CCA. The underlying task is to capture all statistical dependencies in the new representation, and finding correlations (as done by CCA) is equivalent to that only for data from the normal distribution. For considerably non-normal data the choice of dimensionality may not be optimal, but neither is the method itself. Therefore transforming the data so that it would roughly follow normal distribution (such as taking logarithm of gene expression values) would be advisable.

Implementation

We have implemented the method, including the choice of dimensionality and the validation measures presented in the section Validation measures, as an open-source package for R [See Additional file 1].

Additional file 1. A software package in R. An R implementation of the method including the source codes and documentation of the software.

Format: GZ Size: 193KB Download fileOpen Data

Experiments

Validation on gene expression data

We first validate the method on three gene expression data sets (described in Section Methods), by checking how well it preserves the shared variation in data sets and discards the data-specific variation.

In case of two data sets an estimate of mutual information can be computed directly from the canonical correlations as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M29">View MathML</a>

based on the assumption of normally distributed data. Consequently we started by confining to pairs of data sources. Figure 1 shows the results for one of the pairs in each collection; the rest are analogous. It is evident that the method retains the shared variation between data sets and the shared variation increases with increasing number of dimensions in the combined data.

thumbnailFigure 1. Mutual Information. Mutual information for two data sets as a function of the reduced dimensionality. Each subgraph represents mutual information curve for two data sets corresponding to each data collection. The curves for other pairs in each data collection show a similar pattern.

For more than two variables, the measures explained in the Methods Section are used. We compare the results with PCA of the concatenated data matrices. PCA is equally fast, linear, and unsupervised. Note that the proposed CCA-based method is also unsupervised as no class information is used. Furthermore, since both methods have a global optimum, differences in performance cannot be due to optimization issues. The only difference then is related to the main topic of this paper: whether to model all information in the whole data collection, as PCA does, or only the mutual dependencies.

Shared variance (6) and data-specific (7) variance captured by the fused data were computed for each of the three data collections. The presented results are averages over five-fold cross-validation, and the variances have always been computed for the left-out data. In addition to the PCA comparison, we provide baseline results obtained with random orthonormal projections that have uniform distribution on the unit sphere.

The results are presented for each of the data sets in Figures 2, 3, and 4. In all cases it is easily seen that the proposed method retains clearly less data-specific variation than PCA (bottom subfigures), regardless of the dimension. The CCA-based method still keeps more variation than random projections, indicating that it is not purposefully looking for projection directions that would lose more variation than necessary.

thumbnailFigure 2. Shared and Data-specific variation for leukemia data. Shared (top) and data-specific (bottom) variation retained with CCA (solid line) and PCA (dashed line) as a function of the reduced dimensionality for the leukemia data. The values obtained by random projections (dash-dotted line and dotted confidence intervals) have been included for reference. The suggested dimensionality for the CCA-projection is marked with a tick.

thumbnailFigure 3. Shared and Data-specific variation for cell-cycle data. Shared (top) and data-specific (bottom) variation retained with CCA (solid line) and PCA (dashed line) as a function of reduced dimensionality for the cell-cycle data. The values obtained by random projections (dash-dotted line and dotted confidence intervals) have been included for reference. The suggested dimensionality for the CCA-projection is marked with a tick.

thumbnailFigure 4. Shared and Data-specific variation for stress data. Shared (top) and data-specific (bottom) variation retained with CCA (solid line) and PCA (dashed line) as a function of the reduced dimensionality for the stress data. The values obtained by random projections (dash-dotted line and dotted confidence intervals) have been included for reference. The suggested dimensionality for the CCA-projection is marked with a tick.

At the same time the proposed method retains more between-data variation (top subfigures) for wide range of dimensionalities in all cases. The difference is particularly clear for the leukemia data (Fig. 2) where the CCA-based approach is considerably better than the PCA. In stress data (Fig. 4) the difference is also clear, but PCA is also very good in comparison to the random baseline. For cell-cycle data (Fig. 3) the differences are smaller, but for dimensionalities between 3 and 9 the CCA-based method is still clearly better.

It is striking that in all three cases the PCA, which simply aims to keep maximal variation, is the best also in terms of the shared variation for dimensionality of one. A one-dimensional projection, however, loses a lot of the variation and is not too interesting as a summary of several data sets. Hence, this finding does not have a lot of practical significance.

One notable observation is that especially for the leukemia data (Fig. 2) the between-data variance of the CCA-method is, for a wide range of dimensionalities, higher than the corresponding value for the original collection. This does not, however, seem to have clear operational meaning but is merely a side-effect of the heuristic measure.

The curves of extracted variance can be contrasted to the suggested dimensionalities (see Section Choice of dimensionality), marked with ticks in the plots. For two of the three data sets the suggested dimensionality is very close to the maximum point of the between-data variance curve, and when increasing the dimensionality the result remains relatively constant, or even decreases for the leukemia data. While the amount of data-specific variation still keeps increasing, there is no longer a significant amount of shared variation available, and the chosen dimensionality is thus good in terms of these two measures. For the third data collection, the cell-cycle, the suggested dimensionality is somewhat lower than what is needed for maximally capturing the between-data variation. However, as seen in the next section, the chosen dimensionality is still very good for a practical application, providing the best result in the actual case study.

Prototypical Applications

In this section we will discuss a few prototypical ways in which the method could be applied. The method is a general-purpose tool for integrating a collection of data sets in such a way that the effects common to several sets are enhanced. After the integration step any analysis method operating on vectorial data can be used. Here some simple methods are used for demonstrational purposes. The applications are demonstrated on the same data sets that were used in the technical validation.

Shared effects in leukemia subtypes

Pediatric acute lymphoblastic leukemia (ALL) is a heterogeneous disease with subtypes that differ markedly in their cellular and molecular characteristics as well as their response to therapy and subsequent risk of relapse [9]. Combining the expression measurements of the five different ALL subtypes gives a representation where the genes that have similar (or more exactly, statistically dependent) expression profile in several subtypes are similar. Here we are interested in the genes that are highly (over or under) expressed, and thus study the equivalent of differential expression in the combined data set.

The fusion method was applied to combine the five ALL data sets, resulting in a 11-dimensional representation. After this we can proceed as if we only had one data source. It has a 11-dimensional feature vector for each gene, and we separate the 1% of genes that have the highest distance from the origo, implying highest total contribution to the shared variation. This set of genes is compared to the corresponding set obtained from a 11-dimensional PCA projection of the whole collection. In addition, a baseline result computed from the full concatenation of the original data sets is included.

A functional annotation tool, DAVID (Database for Annotation, Visualization and Integrated Discovery) [10] was used to annotate the gene lists to find the gene ontology (GO) enrichments in the biological processes category. The most enriched GO-terms were the same for both CCA and PCA, which is understandable as we are using two linear projection methods on the same collection. We picked the GO-terms which have p-values (Bonferroni corrected) lower than 0.01, and present the counts of the genes from these categories in Table 1. The notable observation is that the CCA-based method has higher count in all but one category, in which the counts are tied. Both methods thus reveal the same kinds of processes, all related to immune response, but CCA is more accurate and is able to include more genes related to these biological processes in the top 1% genes.

Table 1. GO enrichment by CCA and PCA

Classification of cell cycle regulated genes in yeast

The second prototype application is about cell-cycle regulation using the gene expression of Saccharomyces cerevisiae [11]. Expression measurements from 5 different experiments are combined with the proposed method. The combined data is used for the classification of cell-cycle regulated genes.

As the new representation is simply a real-valued vector for each gene, several alternative classifiers are applicable; here K-nearest neighbor (KNN) classifier is selected for demonstrational purposes. We use the cell-cycle regulated genes reported by [11] as the class labels, giving a two-class classification problem: either a gene is or is not cell-cycle regulated.

The leave-one-out classification accuracy of CCA and PCA projections is shown in Figure 5, together with a baseline obtained by using the full concatenation of the original data sets. It is evident that the CCA-based method provides a considerably better representation for separating the cell-cycle regulated genes from the rest. Already the one-dimensional CCA-projection gives a higher accuracy than what is obtained with an eight-dimensional PCA-projection, and the maximal accuracy is clearly higher for CCA and obtained already with a three-dimensional representation. This is exactly the dimensionality suggested by the procedure explained in Section Choice of dimensionality.

thumbnailFigure 5. KNN classification for cell-cycle data. The classification accuracy obtained using the combined representation as a function of dimensionality. The CCA-based combination (solid line) is clearly superior to PCA-based approach (dashed line) for a wide range of dimensionalities and obtains higher maximal accuracy. As a baseline, the classification accuracy obtained by the concatenation of all original data sets (dotted line) is also included.

Defining the environmental stress in yeast

We also study yeast gene expression data from [12,13], consisting of time series of gene expression measurements of Saccharomyces cerevisiae in various stressful treatments. The data is combined to study the genes related to general environmental stress response (ESR).

In [12] the environmental stress response of yeast was studied based on a broad collection of different treatments, out of which 9 are used in our experiment. The original analysis relied primarily on a hierarchical clustering of the whole collection, and was thus based on the overall similarity of the expression patterns. While it is able to cluster the genes into sensible categories, it is ideologically comparable to the PCA approach for preprocessing: It does not take into account that not all variation is equally important.

We suggest that it might be a better idea to focus on the variation shared by the different data sources, instead of trying to characterize the similarity based on all variation. Treatment-specific effects would be specific stress responses and if the task is to find a general response, its fingerprint is in the shared variation. Thus the analysis of environmental stress response should start with a preprocessing step like the one suggested here. We demonstrate how the results of such approach differ from those obtained by [12].

We applied a KNN classifier to the combined data space to classify the genes to belong to the three categories labeled in [12] (a gene is either up- or down-regulated ESR gene, or is not coordinately regulated in stress). The accuracies of CCA and PCA approaches in this task are presented in Figure 6. Again a baseline obtained by using the full concatenation of the original data sets in included. Though the accuracies are similar for some initial dimensionalities, we notice that the accuracy after preprocessing by PCA is higher, by a margin of roughly 0.5% to 1%, for a wide range of dimensionalities including the suggested dimensionality of combined representation, 22, obtained with the method of Section Choice of dimensionality. Also, for the higher dimensionalities the baseline method which simply uses the original data is better. As argued above, this does not tell that CCA was the worse preprocessing method, but instead suggests that the original classes have indeed been constructed based on all variation in the data, including treatment-specific responses. This is not desirable since the definition of an ESR gene is that it would be responsive to stress in general. As the data set has slightly less than 6000 genes this corresponds to a difference of roughly 30 to 60 misclassifications. This characterizes the scale of the disagreement between the two fundamentally different approaches to the preprocessing phase.

thumbnailFigure 6. KNN classification for stress data. The classification accuracy obtained using the combined representation as a function of dimensionality. The CCA-based combination (solid line) is consistently worse than the PCA-based approach (dashed line), implying that the class labels might not correlate that well with the true shared response. As a baseline, the classification accuracy obtained by the concatenation of all original data sets (dotted line) is also included.

This result hints that the definitions created after CCA-based preprocessing would be mostly the same as the ones given in [12], but for some roughly 5 – 10% of genes the classification should be changed.

Conclusion

We studied the problem of data fusion for exploratory data analysis in a setting where the sensible fusion criterion is to look for statistical dependencies between data sets of co-occurring measurements. We showed how a simple summation of the results of a classical method of canonical correlation analysis gives a representation that captures the dependencies, leading to an efficient and robust linear method for the fusion task. It does not solve the data integration task in general, but it shows that the criterion in the data fusion task should not necessarily be to keep all the possible information present in the data collection. Instead, we may want to focus on the aspects shared by different views. We showed how that can be achieved with simple and easily applicable methods.

We demonstrated the validity of the method on three different real gene expression data sets using technical criteria. We further presented three examples on how the method could be used as the preprocessing step in different kinds of analysis tasks.

Methods

Data

Leukemia data

We used the data from [9]. A subset of the data for each leukemia subtype, BCR-ABL, E2A-PBX1, MLL, TEL-AML1 and T-ALL, was chosen such that the patients in each subtype are homogeneous. A set of hyperdiploid samples was used as a control.

We used RMA (Robust Multi-array Analysis) to preprocess the data, and subtracted the mean of the hyperdiploid samples. In total we analyzed 22,283 genes for 31 patients, divided into 5 data sets.

Cell-cycle data

We used the cell-cycle data from [11]. It includes 3 different time courses and 2 different induction experiments. In the original work the data was used to label a set of 800 genes as potentially cell-cycle regulated genes, based on similarity with 104 experimentally verified cell-cycle regulated genes.

We preprocessed the data by imputing missing values with the K-nearest neighbor method, using K = 10. After that the data was Fourier-transformed, and the power spectrum was used for the analysis. In the end we had 5,670 genes, including 724 out of 800 cell-cycle regulated genes defined in [11]. The total number of features in the 5 data sets was 38.

Yeast stress data

We used the yeast gene expression data under various stress conditions from [12,13]. We picked 15 different conditions, 9 from [12] and 6 from [13], resulting in 97 dimensions in total. We then combined them in order to study genes related to general environmental stress response (ESR).

We normalized all time series with their respective zero-points, and imputed missing values by gene-wise averages within each data set. After combining the genes from both sources we got 5,998 genes, out of which 868 were identified as ESR-genes by [12].

Validation measures

The method aims to keep all the variance that is shared among the data sets, while ignoring the variation that is specific to only one of them. In this section we introduce measures on how well this is achieved in real applications. Since there is not straightforward way of quantifying the degree of dependency for several high-dimensional data sources (correlation is only defined for two variables, and estimation of multivariate generalizations of mutual information is difficult), we used two partly heuristic variance-based criteria as comparison measures.

Both measures are based on examining reconstructions of the original data sets. If an integrated representation Pd of full dimensionality is used then it is naturally possible to create a perfect reconstruction, but lower dimensionality introduces errors. We want to measure to what degree the preserved information was shared and to what degree specific to individual data sets.

The reconstruction <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M30">View MathML</a> of the ith data set is obtained by extracting the columns corresponding to the ith data set from <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M31">View MathML</a> and multiplying that by the inverse of the whitening matrix Wi. Here <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M32">View MathML</a> is the generalized inverse of Vd, defined as <a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M33">View MathML</a>.

The first criterion measuring the data-specific variation after the dimensionality reduction to d dimensions is defined as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M34">View MathML</a>

(6)

Each term in the sum is simply the variance of a single reconstruction, and the sum matches the total variation in the collection of data sets. The measure is further normalized so that the value for d = N, the full dimensionality, is one.

For the shared variation we measure the pairwise variation between all pairs of data sets. The measure uses the same reconstructed data sets, and is defined as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/9/111/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/9/111/mathml/M35">View MathML</a>

(7)

again normalized so that the full dimensionality gives the value one. It is worth noticing that the sum of pairwise variations is not a perfect measure for the shared variation for collections with more than two data sets, but it is computationally simple and intuitive.

Availability and requirements

Project name: drCCA;

Project home page: http://www.cis.hut.fi/projects/mi/software/drCCA/ webcite;

Operating system(s): Platform independent;

Programming language: R

License: GNU LGPL;

Any restrictions to use by non-academics: Read GNU LGPL conditions

Authors' contributions

AT implemented the software and carried out the experiments. All authors participated in developing the algorithm, designing of the experiments, and writing of the manuscript. All authors read and approved the final manuscript.

Acknowledgements

The authors are with the Adaptive Informatics Research Centre. This work was supported in part by the Academy of Finland, decision number 207467, in part by the IST Programme of the European Community, under the PASCAL Network of Excellence, IST-2002-506778, and in part by a grant from the University of Helsinki's Research Funds. This publication only reflects the authors' views.

References

  1. Lanckriet GR, Bie TD, Cristianini N, Jordan MI, Noble WS: A statistical framework for genomic data fusion.

    Bioinformatics 2004, 20(16):2626-2635. PubMed Abstract | Publisher Full Text OpenURL

  2. Girolami M, Zhong M: Data Integration for Classification Problems Employing Gaussian Process Priors. In Advances in Neural Information Processing Systems 19. Edited by Schölkopf B, Platt J, Hoffman T. Cambridge, MA: MIT Press; 2007:465-472. OpenURL

  3. Berger JA, Hautaniemi S, Mitra SK, Astola J: Jointly analyzing gene expression and copy number data in breast cancer using data reduction models.

    IEEE/ACM Transactions on Computational Biology and Bioinformatics 2006, 3(1):2-16. PubMed Abstract | Publisher Full Text OpenURL

  4. Nikkilä J, Roos C, Savia E, Kaski S: Explorative modeling of yeast stress response and its regulation with gCCA and associative clustering.

    International Journal of Neural Systems 2005, 15(4):237-246. PubMed Abstract | Publisher Full Text OpenURL

  5. Hotelling H: Relations between two sets of variates.

    Biometrika 1936, 28:321-377. OpenURL

  6. Farquhar JDR, Hardoon DR, Meng H, Shawe-Taylor J, Szedmak S: Two view learning: SVM-2K, Theory and Practice. In Advances in Neural Information Processing Systems 18. Edited by Weiss Y, Schölkopf B, Platt J. Cambridge, MA: MIT Press; 2006:355-362. OpenURL

  7. Yamanishi Y, Vert JP, Nakaya A, Kanehisa M: Extraction of Correlated Gene Clusters from Multiple Genomic Data by Generalized Kernel Canonical Correlation Analysis.

    Bioinformatics 2003, 19:i323-i330. PubMed Abstract | Publisher Full Text OpenURL

  8. Kettenring J: Canonical analysis of several sets of variables.

    Biometrika 1971, 58(3):433-451. Publisher Full Text OpenURL

  9. Ross ME, Zhou X, Song G, Shurtleff S, Girtman K, Williams W, Liu HC, Mahfouz R, Raimondi S, Lenny N, Patel A, Downing J: Classification of pediatric acute lymphoblastic leukemia by gene expression profiling. [http://www.bloodjournal.org/cgi/content/abstract/102/8/2951] webcite

    Blood 2003, 102(8):2951-2959. PubMed Abstract | Publisher Full Text OpenURL

  10. Dennis G, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, Lempicki RA: DAVID: Database for Annotation, Visualization, and Integrated Discovery.

    Genome Biology 2003, 4(5):P3.

    Epub 2003 Apr 3.

    PubMed Abstract | BioMed Central Full Text OpenURL

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

    Molecular Biology of the Cell 1998, 9:3273-97. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  12. Gasch AP, Spellman PT, Kao CM, Carmel-Harel O, Eisen MB, Storz G, Botstein D, Brown PO: Genomic Expression Programs in the Response of Yeast Cells to Environmental Changes.

    Molecular Biology of the Cell 2000, 11:4241-4257. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  13. Causton HC, Ren B, Koh SS, Harbison CT, Kanin E, Jennings EG, Lee TI, True HL, Lander ES, Young RA: Remodeling of Yeast Genome Expression in Response to Environmental Changes.

    Molecular Biology of the Cell 2001, 12:323-337. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL