Recent high-throughput efforts such as ENCODE have generated a large body of genome-scale transcriptional data in multiple conditions (e. g. , cell-types and disease states). Leveraging these data is especially important for network-based approaches to human disease, for instance to identify coherent transcriptional modules (subnetworks) that can inform functional disease mechanisms and pathological pathways. Yet, genome-scale network analysis across conditions is significantly hampered by the paucity of robust and computationally-efficient methods. Building on the Higher-Order Generalized Singular Value Decomposition, we introduce a new algorithmic approach for efficient, parameter-free and reproducible identification of network-modules simultaneously across multiple conditions. Our method can accommodate weighted (and unweighted) networks of any size and can similarly use co-expression or raw gene expression input data, without hinging upon the definition and stability of the correlation used to assess gene co-expression. In simulation studies, we demonstrated distinctive advantages of our method over existing methods, which was able to recover accurately both common and condition-specific network-modules without entailing ad-hoc input parameters as required by other approaches. We applied our method to genome-scale and multi-tissue transcriptomic datasets from rats (microarray-based) and humans (mRNA-sequencing-based) and identified several common and tissue-specific subnetworks with functional significance, which were not detected by other methods. In humans we recapitulated the crosstalk between cell-cycle progression and cell-extracellular matrix interactions processes in ventricular zones during neocortex expansion and further, we uncovered pathways related to development of later cognitive functions in the cortical plate of the developing brain which were previously unappreciated. Analyses of seven rat tissues identified a multi-tissue subnetwork of co-expressed heat shock protein (Hsp) and cardiomyopathy genes (Bag3, Cryab, Kras, Emd, Plec), which was significantly replicated using separate failing heart and liver gene expression datasets in humans, thus revealing a conserved functional role for Hsp genes in cardiovascular disease.
The increasingly cheaper and rapid accumulation of large -omics datasets across several experimental conditions has prompted generation of a wealth of data on biological networks. This growth of network data now permits their large scale applications to biomedical research, including analysis of gene function, metabolic and signaling pathways, as well as disease-related or cell function-related networks [1], [2]. However, reconstructing and interpreting large biological networks, such as co-expression networks, protein-protein interaction networks or genetic networks, with different features (e. g. , sparse or densely interconnected, etc.) poses many challenges, advocating efficient and flexible methods for network inference and pattern discovery. An important level of complexity in current network analysis regards its extension to multiple conditions, for instance different species [3], cell-types [4] or disease states [5], [6]. For example, reconstruction of networks across multiple disease-states is becoming a useful approach for efficient drug-target discovery, as networks can inform the “biological context” (e. g. , pathways, cellular processes) where genes operate and therefore can help designing better therapeutic interventions [7]. In genetic studies of complex diseases researchers increasingly focus on groups of highly interconnected genes within larger networks (referred to as clusters, modules or subnetworks) to elucidate specific cellular and molecular processes that might represent functional disease mechanisms and pathological pathways [8]–[10]. While several computational tools for network analysis in single datasets or conditions are available, only few computationally efficient methods for genome-scale network analysis across multiple conditions have been developed to date. These methods can be broadly classified into two main categories: (i) methods to find the “difference” between networks across conditions or to pinpoint condition-specific networks [11]–[14], or (ii) methods to identify the common parts in networks across conditions [15]–[17]. More recently, tensor-based computational frameworks [15] or probabilistic Markov blanket search algorithms [18] have been proposed to learn network structures across conditions. However, these methods are either heavily influenced by the choice of input parameters (e. g. , number of clusters, number of nodes within a cluster, cluster interconnectivity) [15] or, being based on probabilistic graphical modelling, they become prohibitively slow for high number of conditions since they are trying to learn the structure of large graphs [18]. Complementary to the above approaches, spectral methods, such as Singular Value Decomposition (SVD), have been also proposed to investigate patterns of connectivity between nodes within a single network [19], [20] or for comparing two networks [21]. Generally, any network can be described as a graph, which is denoted as comprising a set of vertices or nodes together with a set of edges [22]. The graph may be represented by a square, symmetric, real-valued matrix of size whose entries denote the relationship between the corresponding nodes. In the affinity matrix, the element, called weight, represents the strength of connection between vertices and. For instance, in gene regulatory (or co-expression) networks, the nodes might represent genes (or mRNAs expression) and edges represent the strength of gene-gene interactions (or mRNAs co-expression). Generalized Singular Value Decomposition (GSVD) can be used to identify sub-network structures and for comparative analysis of genomic datasets across two conditions [11], [23]. Given two matrices and [24], [25], their GSVD is given by (1) where and have orthonormal columns, is invertible, with, with. The ratios are the generalized singular values of and. In this setup, the common factor is informative of the cluster structure shared across the two data matrices. Recently, a novel mathematical formulation, higher-order GSVD (HO GSVD), which is constructed for more than two data matrices has been proposed [26]. Under this framework, the matrices, each with full column rank (i. e. , the maximum number of linearly independent column vectors of is), are decomposed as (2) where is composed of normalized left basis vectors, with and the latent factor matrix is composed of normalized right basis vectors. The HO GSVD can be also derived in the special case of square, symmetric, full rank affinity matrices, , where each element represents the weight of the edge between node and in the th condition. It has been previously employed to compare multiple datasets with identical column size in order to detect their common substructures of columns (i. e. , observations) [26]. Yet, another useful application of the HO GSVD to genomics is to set it to discover gene networks across multiple conditions and pinpoint “common” and “differential” cluster structures. In this paper, we build on the flexible HO GSVD mathematical framework and propose a new, parameter-free computational algorithm (Cross-Conditions Cluster Detection or C3D) for automatic detection of both similarity and dissimilarity clustering patterns in large weighted (and unweighted) networks across several conditions (). The original HO GSVD model has been employed for analysis of datasets that had varying number of genes (), the same number of observations () (i. e. , arrays/time points in [26]) across conditions and with. As such, this illustrative application of the HO GSVD in genomics was aimed at the identification of common structures within the observations [26]. Here, we built on the initial HO GSVD to extract sub-structures (i. e. , common and differential clusters) from genes across multiple conditions () by applying the decomposition to the transposed expression matrix. We show how this enables a more general application of the HO GSVD framework to genome-scale network analysis of genomic data (e. g. , microarray, RNA-seq) in multiple conditions. Besides, a distinctive feature of our method is in its capability to take as an input either the raw expression matrices or co-expression matrices, allowing flexibility in the choice of the co-expression measures (e. g. , Spearman, Kendall, mutual information, etc.). Figure 1 illustrates the working principle of the C3D algorithm. The input data for C3D can be provided into different formats to be used by the HO GSVD: (i) the raw expression data matrices () or (ii) the co-expression data matrices (). In the former case, a first data initialization step is conducted where the input expression matrices, with the same number of genes are converted to co-expression matrices by scaling their variance to 1 and taking their quadratic form. In the second step (HO GSVD-based algorithm), an approximate HO GSVD is employed to identify a common basis, with representing the dimension of the GSVD common subspace, for the decomposition of the input datasets and identify the common and differential correlation structures. The HO GSVD-based algorithm computes a square matrix, which is built on the arithmetic mean of all pairwise quotients where denotes the Moore-Penrose inverse of the co-expression matrix [24] (see Methods section). The first eigenvectors of (according to the norm of the corresponding eigenvalues) are then used to identify an approximate decomposition of the input co-expression matrices and form the decomposition basis. Specifically, each selected column vector of is used to reorder the input data matrices such that candidate “common” (or “differential”) clusters can be identified. In the third step (cluster nodes selection and validation), we employ a mixture model approach to classify genes and assign them to each cluster based on a misclassification error rate (MER). Finally, we implemented an empirical cluster validation procedure to identify the conditions where clusters are present and assess the level of significance for clusters within each condition. To demonstrate the increased power and benefits of our HO GSVD-based algorithm, we carried out an extensive simulation study and benchmarked C3D against commonly used methods that were designed to detect either common (WGCNA [16], [17]) or differential network structures (DiffCoEx [13]) across multiple conditions. We show that our approach has higher power and stability in detecting both common and differential co-expression clusters across all simulated conditions, while being two to seven fold less computationally intensive than alternative methods. In contrast with alternative approaches that require specification of ad-hoc input parameters, the proposed method has the distinctive advantage of being parameter-free, which makes it a powerful tool for real data exploration and analysis. To substantiate this claim, we applied C3D to publicly available transcriptomic datasets in rats and humans and identified several multi-tissue gene co-expression networks that were associated with specific functional processes relevant to phenotypic variation and disease.
We carried out a simulation study to compare our method with commonly used approaches for identification of “common” or “differential” clusters across multiple networks: (1) WGCNA and (2) DiffCoEx. The WGCNA method for detection of common clusters across co-expression networks employs a “soft” threshold to assign a connection weight to each gene pair and extract densely connected gene clusters that are present in all conditions. The DiffCoEx method follows a strategy similar to WGCNA but, instead, it focuses on detecting the differences in co-expression patterns (“differential” clusters) between multiple conditions. Additional details on the specific parameterizations used in for WGCNA and DiffCoEx analyzes are reported in Text S1. To simulate a realistic example of gene expression data from multiple conditions that represent a typical “small large ” scenario, we draw inspiration from a publicly available multi-tissue microarray dataset consisting of genome-wide expression profiles from recombinant inbred rat strains in seven tissues [27]. We simulated different types of clusters that are either detected in all conditions (“common” clusters) or are specific to a subset of conditions (“differential” clusters), Figure 2. We considered dense clusters of variable sizes (100–500 nodes) where each node is connected with all other nodes in the cluster with a given weight (), which is defined as the Pearson correlation between expression profiles of genes and. We simulated clusters with varying cluster densities (0. 1,0. 3,0. 5,0. 7), which were defined as the average Pearson correlation between any pair of nodes within a cluster. In addition to the simple case of a cluster common to all conditions and with the same size (Cluster pattern 1), we set out to evaluate the sensitivity of our and alternative approaches to detect clusters which are present only in a subset of conditions and that overlap partially across conditions. This is more likely to be relevant for analysis of pathways and gene networks across tissues or during development, where varying gene-sets can exert their function only at specific developmental times or in specific cell-types. To account for these more complex scenarios, we simulated “nested” (Cluster pattern 2) and partially “overlapping” (Cluster pattern 3) cluster structures (Figure 2). Cluster pattern 2 and Cluster pattern 3 have an intersection part, defined by the nodes in common to all conditions, and a union part, defined by the nodes in common to all conditions plus the nodes present in individual conditions. In summary, for each of the four cluster densities considered one dataset consisted of a and matrix in conditions, where each cluster type (Clusters patterns 1–3) was simultaneously present in the data matrix. To assess reliability of the results, for each of these data we generated 20 independent replicates, yielding a total of 560 simulated datasets. Similarly, to evaluate how the number of available observations affects the methods' performance we simulated datasets consisting of a and matrix in conditions (20 replicates, 560 datasets in total). See Text S1 for additional details. The True Positive Rate (TPR) and the False Positive Rate (FPR) are widely used as evaluation metrics for a classification model and can be used to quantitatively assess (and compare) methods performance [28]. The TPR defines how many correct positive results (simulated clusters genes within the called cluster) occur among all results called positive in the analysis by a given method. FPR, on the other hand, defines how many incorrect positive results occur among all results called positives. Typically, a and the corresponding indicate a perfect classifier (or a perfect method). In our simulation study, the best cluster detection method would yield both high TPR and low FPR levels for different cluster types, sizes and densities. For each simulated cluster type, Figure 3 shows the TP/FP rates for C3D, WGCNA and DiffCoEx methods as a function of the simulated cluster densities. For C3D we controlled the (local) misclassification error (i. e. , the probability to assign wrongly a gene to a cluster) to be less than 0. 05 or less than 0. 2, and required that each cluster is detected with, whereas for WGCNA and DiffCoEx we used two (default) parameterizations chosen according to the software guidelines (see Methods section). The C3D method outperformed WGCNA in the identification of clusters present in all conditions (Cluster pattern 1, Figure 3), and showed to have consistently high TPR (and very low FPR,) irrespective of the simulated cluster density. WGCNA performance varied considerably as a function of the simulated cluster density and, depending on the adopted parameterization, FPR levels were (reaching 20% in one case), Figure 3. Furthermore, we observed large variations in WGCNA performance (mostly in the TPR), which are indicated by the large standard deviations in TPRs calculated from the 20 replicated datasets. For more complicated patterns (“nested” and “overlapping” clusters), we compared C3D with WGCNA to detect the intersection part (100 nodes) of common clusters. Since WGCNA is designed to detect only those clusters shared across all conditions, for clusters present in a subset of conditions, we run WGCNA only in the set of conditions where the simulated clusters were present. For Cluster patterns 2–3, C3D and WGCNA performances were similar, reaching high TPR for detection of the intersection part of clusters with simulated (Figure 3). However, C3D showed higher TPRs than WGCNA to detect clusters with low densities (0. 1–0. 3), while controlling the FPR at low levels (, Cluster pattern 2 intersection). In the case of partially overlapping clusters present in a subset of conditions (Cluster patterns 2–3) we compared C3D with DiffCoEx in respect of detecting the union part (500 nodes) of “differential” clusters, and calculated TPR and FPR for detection of this cluster (indicated with a black square at the top of Figure 3). We found that C3D outperformed DiffCoEx across the simulated scenarios. In the case of the “nested” cluster structures that are present in 5 out of 7 conditions, C3D had consistently higher TPR levels than DiffCoEx, which showed comparable TPR levels only for detection of highly-dense clusters (i. e. , , Cluster pattern 2 union, Figure 3). However, similarly to what observed for WGCNA method, in this case DiffCoEx showed large variability in its performance across the 20 replicated datasets. The difference in performance between C3D and DiffCoEx was observed also in the more complicated case of partially overlapping cluster structures (Cluster pattern 3). In this case, C3D showed consistently higher TPR than DiffCoEx that reached a maximum as compared with of C3D. Both methods showed comparably low FPR () for detection of the union part of Cluster patterns 2–3 (Figure 3). Similarly to what observed for the simulated data with observations, C3D performed better than (or as good as) both WGCNA and DiffCoEx when benchmarked on simulated datasets with only observations (Figure S1). As expected, all methods had lower TPRs associated with the detection of low-density clusters, however also with a small number of observations, C3D showed significantly better (and more stable) results than WGCNA and similar performance as compared with DiffCoEx. Notably, for detection of “common” clusters present in all conditions (Cluster pattern 1), CD3 held high TPR levels (and FPR) whereas WGCNA' s performance dropped significantly, reaching a maximum TPR (Figure S1). These data show that C3D on balance performed better than WGCNA and DiffCoEx across all simulated scenarios. We underline that while WGCNA and DiffCoEx methods are specifically designed to detect either common or differential clusters, respectively, here we showed that C3D was equally or more accurate than both methods in the detection of common and differential cluster structures. We also highlight how C3D ability to detect correctly the simulated clusters was highly consistent across all runs on the replicated datasets, as shown by the small standard deviations of the mean TP and FP estimates (Figure 3). In contrast, we observed that both WGCNA and DiffCoEx performances varied appreciably across the replicated simulations, often resulting in large standard deviations of the mean TP and FP estimates. To better assess the reliability of the different methods we calculated the relative standard deviation of the TPR measured in all analyzed datasets. In 560 simulated datasets of size, the C3D method had a median RSD of (range 113. 36) whereas WGCNA and DiffCoEx have median (range 447. 2) and median (range 133. 39), respectively. Similarly, in 560 datasets of size we estimated the following RSDs of TPR: 12. 43 (range 113. 38) for C3D, 57. 52 (range 161. 89) for WGCNA and 87. 96 (range 120. 59) for DiffCoEx. The large RSDs of TPR calculated from the WGCNA and DiffCoEx analyzes originated because these methods often detected the simulated cluster (s) only in small number of replicates (e. g. , 2 out of 20). Besides, in a few cases the TP/FP rates of WGCNA and DiffCoEx were influenced by the adopted parameterization (for instance, FPR in the WGCNA analysis of Cluster pattern 1, Figure 3), suggesting that different choices of the input parameters can affect the detection of clusters (see Text S1 for additional details). The C3D algorithm is built on the HO-GSVD framework and as such does not require the user to specify ad-hoc parameters to detect common or differential clusters. In our implementation of the C3D algorithm the user can control the MER at a specified level before the cluster genes are empirically validated using a permutation-based procedure (see Methods section). In these simulation studies, we have used two different MERs (5% and 20%) to inform a suitable choice of MER that maximizes true positive without inflating false positive rates. On average, we observed a increase in the TPR when was adopted as compared with. However, we found no significantly higher FPR, which were always across all simulated datasets, this suggesting that using the less stringent in real data analyzes is likely to increase the detection of true gene clusters, without increasing significantly false positives. Finally, we used a standard desktop computer (Mac Pro, GHz Quad-core Intel Xeon with 20 Gb RAM) to evaluate the computational time required by C3D and compare it with WGCNA and DiffCoEx to analyze the simulated datasets. While the run time of C3D scales exponentially with the number of genes in the input matrices or the number of conditions, our Matlab implementation of C3D is relatively fast and requires only 1,200s to analyze a gene co-expression matrix in conditions and 10s to analyze a gene co-expression matrix in conditions (Figure S2). When compared with competing approaches, we assessed that to process simulated datasets of 1,000 and 10,000 genes (with observations and conditions) C3D requires significantly smaller CPU time than DiffCoEx (up to 2. 3 fold more CPU time) and WGCNA (up to 8. 2 fold more CPU time), respectively (Figure S2). To show how C3D provides a powerful, practical framework for real genome-scale analyzes and yields new biological insights into pathways and molecular networks, we report an application to two large multi-tissue gene expression datasets in rats and humans. Transcriptional profiling was carried out by Affymetrix microarray in the rat and mRNA sequencing (RNA-seq) in humans, respectively. The microarray dataset consisted of genome-wide expression profiles (probe sets) that were measured in seven tissues (adrenal, aorta, fat, kidney, left ventricle, liver and skeletal muscle) in a panel of recombinant inbred rat strains [29], which is a well characterized model of hypertension, metabolic syndrome and cardiovascular disease [27], [30], [31]. The RNA-seq datasets consisted of genome-wide transcriptomic data of human fetal neocortex, which have been generated to investigate the molecular mechanisms underlying differences in germinal zones of the developing human brain. The human dataset consisted of expressed genes which were analyzed in four regions of the fetal neocortex (ventricular zone (VZ), inner subventricular zone (ISVZ), outer subventricular zone (OSVZ) and cortical plate (CP) ) from six 13–16 weeks postconception human fetuses [32]. In both rat and human analyzes, to identify common and differential clusters we extracted the top ten eigenvectors (based on the modulus of the eigenvalues of the decomposition of) as candidates which are then used as input for the cluster nodes selection and validation step of the C3D algorithm (see Methods).
Building on the HO GSVD framework, we have developed a new algorithm (C3D) for efficient, parameter-free and automatic detection of co-expression clusters and networks in multiple conditions. Our method is designed for analysis of weighted (and unweighted) networks (input matrices) across conditions, enabling applications to diverse data types and structures. Although the original HO GSVD algorithm assumes the non-singularity of the co-expression matrix, by using the Moore-Penrose pseudo-inverse, our C3D algorithm can be applied to the non-invertible case. We show that when an exact HO-GSVD of the input matrices exists (as defined in (4), see Methods), our HO GSVD is able to extract the right decomposition basis through the eigen-decomposition of, whereas it finds an approximate decomposition of the data in the absence of an exact solution (Figure S4). In particular, our empirical simulations and real-case applications reveal that our approximate decomposition is able to capture both common and differential co-expression structures for a wide range of noise levels, suggesting that our algorithm can be useful for practical applications to genomic data. Here, through the HO GSVD of large-scale genomic datasets we aimed to uncover the complex interactions between genes (networks) that can occur within or across multiple conditions. One distinctive feature of our computational method is in the flexible and simultaneous identification of both “common” and “differential” sub-network structures across several conditions. Selecting informative vectors of, we provide different orderings of to reveal candidate clusters that are important to all conditions or specific to a sub-set of conditions; then, we can distinguish the specific conditions where the clusters are present using a permutation-based approach. This procedure allows to pinpoint automatically the specific conditions where the sub-network structures are present and, at the same time, to provide an empirical estimate of the statistical significance (empirical P-value) for each cluster identified. In simulation studies, we demonstrated how C3D outperforms competing approaches in accuracy and reliability while being computationally less demanding. We highlight how our method allowed accurate detection of clusters within complex structures (i. e. , “common”, “nested” and “overlapping” networks) by specifying only the desired level of statistical significance: misclassification error rate to assign genes to clusters and empirical P-value for cluster detection. In contrast with other approaches, C3D does not need the user to specify ad-hoc parameters related to the expected number of clusters or cluster density [15] or necessary to determine the optimal height cut-off in the gene clustering tree [13], [16], [17]. Typically, these unknown parameters need to be “finely tuned” on each dataset in order to obtain the best compromise between TP and FP for each cluster (see Text S1 for additional details). We also showed that the results obtained by two competing and widely-used methods (WGCNA and DiffCoEx) were less stable than those provided by C3D. This was apparent in the significantly smaller relative standard deviations in TPR calculated across simulated datasets in the C3D analyzes as compared with WGCNA and DiffCoEx. Since C3D utilised raw gene expression data matrices as input, the higher stability of C3D might be due to the reduced influence of the small number of observations on the stability of co-expression estimates, which can result in extreme patterns of correlation changes, corresponding to stable and fragile co-expression, as previously shown [62]. The high stability in the results and the parameter-free “nature” of the HO GSVD approach make the C3D algorithm a powerful computational tool for real genomic data exploration and analysis. To demonstrate this point, we reported an application of C3D to two large transcriptional datasets: (i) microarray-based gene expression profiles in seven rat tissues and (ii) RNA-seq-based gene expression analysis of germinal zones from human fetal neocortex. In the rat analysis, we reported several functionally enriched co-expression clusters, including a previously identified inflammatory gene network driven by the IRF7 transcription factor that represents a gene expression signature of macrophages within complex tissues. While this co-expression network was experimentally validated [27] it was not recovered by WGCNA, that surprisingly placed the IRF7 transcription factor and many regulated target genes in the group of “non-clustered” genes. In addition, our C3D analyzes revealed novel gene co-expression networks in sub-sets of tissues. For instance, we identified a network comprising Hsp and known cardiomyopathy genes, which suggested coordinated regulation of heat shock proteins genes in multiple tissues, and their potential functional role in cardiovascular disease [50]. While this network was not recovered by either WGCNA or DiffCoEx analyzes, we were able to replicate this new finding using separate cardiac and liver gene expression datasets in humans (Figure 4). In the study of human fetal neocortex we demonstrated previously undescribed co-expression between cell cycle and ECM-receptor interaction pathways and support their role in the proliferation and self-renewal of neural progenitors. In addition, our analyzes highlighted that pathways central to later cognitive functions (e. g. , calcium signaling, long-term potentiation, axon guidance) are present at an early stage in the developing human brain [61], which was not previously appreciated. These studies illustrated how our method can be effectively applied to leverage the vast stream of genome-scale transcriptional data that has risen exponentially over the last years, promising to aid the fine-scale characterization of both context-specific and systems-level networks and pathways.
In this step we assume the input data are non-square matrices, where the rows represent the observations and the columns indicate genes. The number of genes must be the same across datasets while the number of observations can differ. We first log transform the data and subtract for each gene its average gene expression to avoid capturing differences in average gene expression across conditions. We then calculate the co-expression matrices corresponding to each condition. Each represents the covariance matrix of the data in condition. As in classic principal component analysis, the columns of can be scaled to unit variance to work on the correlation matrices rather than the covariance. Alternatively, our algorithm can directly take any co-expression matrix as input. This feature of our algorithm allows to extract common and differential clusters from matrices based on different co-expression measures, including robust correlation (e. g. Spearman, Kendall) and non linear metrics such as mutual information [63]. Similarly to classic SVD, each observation from the input data can be characterized by its expression profile and represented by a data point in a dimensional space. The observations from all datasets are contained in a subspace of dimension, which thereafter is referred to as the HO GSVD subspace. Here, we aim at finding directions in the HO GSVD subspace that either capture the variability in gene expression that is common to all conditions (common factors) or that is specific to a subset of conditions (differential factors). Inspired by [26] we developed a general algorithm that allows computation of an approximate solution to the HO GSVD problem in the non full column rank case. In the HO GSVD, are decomposed into where, is a diagonal matrix with elements for and contain the right basis vectors of the HO GSVD subspace where. The right basis vectors allow to identify set of genes (clusters) with similar co-expression patterns, that are either specific to a subset of conditions or common to all conditions. Here we explain the derivation of our HO GSVD-based algorithm in the general case of non-square matrices. The derivation and discussion of the special cases (square, symmetric matrices with full rank and square, symmetric matrices with full rank) is reported in Text S1. In the most general case, we define the right basis vectors as the solution of the eigen-decomposition problem of the matrix (3) where is the arithmetic mean of all the pairwise quotients and denotes the Moore-Penrose inverse of the co-expression matrix [24]. Here the Moore-Penrose inverse is used as a substitute of since the invertibility of is not guaranteed when, which is the typical scenario in genomics. We now assume there is an approximate HO GSVD where is composed of orthonormal left basis vectors and. In this case, for all we have (4) and its Moore-Penrose inverse is given by (5) Therefore we have (6) since is full row rank. Hence we can rewrite as follows (7) When there exists a common subspace of dimension, with basis vectors, for which the decomposition of the co-expression matrices (4) is exact, equation (7) becomes an equality and the eigenvectors of will lead to the exact basis of the common subspace. In HO GSVD applications to genomics data, can be as large as the total number of observations (i. e. ,), and an exact common decomposition of the co-expression matrices might not be possible. In this case the eigenvectors of do not provide an exact decomposition of the subspace. Moreover, is not guaranteed to be non-defective and have a full set of real eigenvalues and eigenvectors. However, even in the absence of an exact common decomposition, the real part of the complex eigenvectors can be used to derive a low rank approximation of the common subspace and extract common and differential covariance structures from the data. To test the ability of our HO GSVD based algorithm to capture these covariance structures in the data in the presence of a “noisy” HO GSVD decomposition we performed an empirical simulation study (see Text S1 for details). Our simulations suggest that if a common subspace of dimension with basis vectors explains a significant fraction of the variance in the original datasets, the approximation (4) holds and the first eigenvectors of the matrix (corresponding to the largest eigenvalues of) will provide a good approximation of the basis vectors of the HO GSVD subspace (Figure S4). We selected two large gene expression datasets from rats and humans, where genome-wide expression profiles were assessed in the same subject/animal across multiple tissues. The rat datasets consisted of microarray-based expression profiles for probe sets that were measured in adrenal, aorta, fat, kidney, left ventricle, liver and skeletal muscle tissues in a panel of recombinant inbred rat strains [29]. Microarray expression data were retrieved from ArrayExpress, http: //www. ebi. ac. uk/arrayexpress/, (skeletal muscle, E-TABM-458; aorta, E-MTAB-322; liver, E-MTAB-323, fat and kidney, E-AFMX-7; heart, MIMR-222; adrenal, E-TABM-457); gene expression summaries were derived using robust multichip average (RMA) algorithm [66] and normalized using Z-score transformation before analysis with C3D. The human data were retrieved from the Gene Expression Omnibus (GEO) database (www. ncbi. nlm. nih. gov/geo) under accession number GSE38805. Briefly, total RNA from the VZ, ISVZ, OSVZ, and CP of six 13–16 wk postconception human fetuses was isolated from laser-capture microdissected Nissl-stained cryosections of dorsolateral telencephalon (see [32] for additional details on experimental procedures). RNA-seq data were expressed as fragments per kilobase of exon per million fragments mapped (FPKM) values and normalized on log2 scale, yielding an expression matrix of in neocortex regions, which were analyzed by C3D. The Matlab implementation of the C3D algorithm, detailed instructions to run the code and an example of the simulated datasets used in these studies can be downloaded from http: //www. csc. mrc. ac. uk/Research/Groups/IB/IntegrativeGenomicsMedicine/ contact information: enrico. petrettocsc. mrc. ac. uk or xiaolin. xiaocsc. mrc. ac. uk
Complex biological interactions and processes can be modelled as networks, for instance metabolic pathways or protein-protein interactions. The growing availability of large high-throughput data in several experimental conditions now permits the full-scale analysis of biological interactions and processes. However, no reliable and computationally efficient methods for simultaneous analysis of multiple large-scale interaction datasets (networks) have been developed to date. To overcome this shortcoming, we have developed a new computational framework that is parameter-free, computationally efficient and highly reliable. We showed how these distinctive properties make it a useful tool for real genomic data exploration and analyses. Indeed, in extensive simulation studies and real-data analyses we have demonstrated that our method outperformed existing approaches in terms of efficiency and, most importantly, reproducibility of the results. Beyond the computational advantages, we illustrated how our method can be effectively applied to leverage the vast stream of genome-scale transcriptional data that has risen exponentially over the last years. In contrast with existing approaches, using our method we were able to identify and replicate multi-tissue gene co-expression networks that were associated with specific functional processes relevant to phenotypic variation and disease in rats and humans.
biotechnology
algorithms
systems biology
computer applications
computer science
mathematics
algebra
genetics
applied mathematics
computing methods
biology
genomics
computational biology
numerical analysis
Multi-tissue Analysis of Co-expression Networks by Higher-Order Generalized Singular Value Decomposition Identifies Functionally Coherent Transcriptional Modules