Independent component analysis: mining microarray data for fundamental human gene expression modules

The wide use of high-throughput DNA microarray technology promises to provide an increasingly detailed view of the human transcriptome in many different contexts. As experimentalists continue to sample a variety of biological conditions and cell types, public gene expression repositories grow exponentially: the Gene Expression Omnibus (GEO), the largest of these databases, now contains over 400,000 individual microarrays ( www.ncbi.nlm.nih.gov/geo/ ). Although biomedical researchers typically design microarray experiments to explore specific biological contexts, meta-analyses that integrate data from multiple experiments have the potential to reveal relationships that are not accessible through any individual dataset [ 1 – 10 ]. One critical question in the post-genomic era is the identification of sets of functionally related genes that correspond to a particular biological process. While curated databases such as the Gene Ontology (GO) group genes with related functions, these resources are necessarily incomplete and prone to human error [ 11 ]. Furthermore, GO categories do not necessarily correspond to co-regulated transcriptional units [ 12 ]. Computational methods based on co-expression, therefore, have the potential to enhance the identification of genes that contribute to the same biological processes [ 1 , 7 , 13 ]. Given the wide sampling of diverse cellular conditions, a meta-analysis of human gene expression data should yield information about many of the key pathways and processes in human biology. Methods that extract gene expression modules (subsets of genes that are co-regulated across subsets of conditions) provide a means for identifying functionally related genes from microarray data [ 13 – 21 ]. Popular methods for this analysis, such as hierarchical clustering, k-means clustering, and self-organizing maps (SOMs), identify disjoint sets of genes that are co-expressed throughout a particular dataset. Because these approaches require that genes belong to only one transcriptional module, however, they do not capture an important feature of gene regulation: it is common for the same proteins to serve different functions depending on biological context [ 22 ]. Furthermore, simple clustering approaches may be inappropriate for large-scale gene expression meta-analysis because they capture global signals in the compendium and exclude more specific signals [ 23 ]. These challenges have led to methods like biclustering that simultaneously cluster subsets of genes over subsets of conditions [ 23 – 25 ]. For a more biologically relevant model for gene expression, researchers often turn to matrix decomposition methods [ 19 , 25 – 32 ]. Gene expression data lends itself to representation as a matrix, with columns corresponding to genes and rows corresponding to experiments. Like biclustering, matrix decomposition offers an improvement over traditional clustering methods by modeling microarray data as a combination of biological modules that can share genes. Furthermore, matrix decomposition provides a generative model that allows for unified analysis of multiple datasets. One common method, principal component analysis (PCA), provides a representation of microarray data in terms of a set of linearly uncorrelated axes [ 26 , 27 ]. Each axis, or principal component, explains the maximal variance not represented by the previous, orthogonal components. While PCA can identify interesting biological information, its linear transformation involves only second order statistics; like clustering methods based on co-expression, PCA may miss more complex relationships between genes [ 31 ]. Furthermore, the application of PCA to biological data assumes that gene expression follows a multivariate normal distribution. Recent studies demonstrate that microarray gene expression measurements typically follow a super-Gaussian distribution [ 33 , 34 ], suggesting that PCA applies a model that does not accurately reflect our knowledge of biology. Independent component analysis (ICA), in contrast, provides a more biologically plausible model for gene expression data by assuming non-Gaussian data distributions. A blind source separation algorithm, ICA models observations as a linear combination of latent feature variables, or components, which are chosen to be as statistically independent as possible. For microarray data, observations consist of microarray gene expression measurements, and independent components are interpreted to be transcriptional modules that often correspond to specific biological processes [ 28 ]. ICA has proven successful in a variety of biological inquiries, including identifying oscillating regulatory modules in yeast cell cycle data [ 28 ], investigating tumor-related pathways [ 35 – 39 ], classifying disease datasets [ 40 , 41 ], characterizing transcriptional regulators [ 42 , 43 ], identifying disease-specific biomarkers [ 44 ], and examining response to bacterial infection [ 45 ]. Furthermore, ICA outperforms PCA and other unsupervised methods in identifying co-regulated and biologically relevant gene modules in diverse datasets [ 31 , 39 ]. Although ICA produces stochastic component estimates with different initial conditions, clustering and averaging the results from multiple runs can yield robust independent component estimates [ 38 , 45 ]. While previous applications of ICA to microarray data have used at most hundreds of experiments, no large-scale meta-analysis to identify fundamental gene modules has been attempted. Compared with resources such as the Gene Ontology, ICA provides a data-driven method for exploring functional relationships and grouping genes into transcriptional modules. We apply ICA to a large microarray compendium initially comprised of 9,395 microarrays representing a diverse set of experimental conditions and cell types. We identify 423 fundamental components (FCs) of human biology, and show that these components yield gene expression modules with coherent functions. To evaluate the biological relevance of our fundamental components, we develop a method to perform differential expression analysis in the feature-space described by our FCs. Using this technique, we assess the ability of our method to identify known mechanisms of parthenolide (PTL), a preclinical drug under investigation for its ability to selectively induce apoptosis in multiple cancer types [ 46 – 50 ]. Known PTL effects can be divided into two intracellular signals. First, PTL induces oxidative stress, evidenced by increased levels of reactive oxygen species [ 51 ] and activation of c-Jun N-terminal kinase (JNK) [ 52 ]. Second, PTL inhibits inflammatory responses via inhibition of STAT3 [ 53 ] and the transcription factor NF-κB [ 54 , 55 ]. PTL treatment leads to apoptosis in various cancer cell lines [ 46 , 51 ], and activation of p53 has been associated with the AML-specific apoptosis mechanism [ 46 ]. We show that independent components derived from a diverse compendium offer module-level insight into transcriptional response that cannot be gleaned from the PTL dataset alone.

We downloaded our microarray compendium from the Gene Expression Omnibus (GEO) [ 56 ], selecting all GEO Series (GSEs) run on the Affymetrix Human Genome U133 Plus 2.0 array (GEO accession GPL570 ) available on May 28, 2008. We filtered these arrays to remove samples that represented species other than human, and we removed GSEs with missing data so that imputation was not necessary. Normalized microarray data are often uploaded to GEO, but we downloaded only unprocessed CEL files in order to standardize the normalization procedures across all arrays in the compendium. After applying these filters, the resulting dataset consisted of 298 GEO Series comprised of 9,395 arrays. We applied a two-step normalization pipeline as previously described [ 10 ]. For each series, we aggregated and normalized probe level information using robust multi-array average (RMA) [ 57 ], transformed each expression value using log base 2, and removed technical bias resulting from variation in hybridization conditions and starting material using the R package bias 0.0.3 [ 58 ]. This within-series normalization identified probe or arrays outliers within each dataset. We mapped probes to genes using the Bioconductor annotation package hgu133plus2 1.16.0 [ 59 ], and calculated expression values for each gene by averaging the values of probesets measuring the same gene. Finally, we performed quantile normalization [ 60 ] on the entire compendium using the limma R package (version 2.18.2) [ 61 ] in order to reconcile broader differences between datatsets and ensure that all arrays were on the same scale prior to applying ICA. This produced a compendium comprising 20,099 genes and 9,395 arrays. To reduce the contributions of over-represented conditions, and to eliminate rare signals that the compendium did not sample sufficiently, we grouped similar arrays using hierarchical clustering. First, we centered and scaled the expression values within an array, and centered the expression values for each gene. We calculated pairwise correlations between arrays using Spearman’s rank correlation coefficient, r , a metric robust to outliers that performs well with microarray data [ 62 , 63 ]. We defined pairwise distances, d , by (1) d x , y = 1 − | r x , y | , where x and y are the gene expression measurements of individual microarrays. Using these distances, we applied agglomerative hierarchical clustering to the compendium using average linkage. To filter the resulting clusters, we required that nodes pass two cutoffs. First, in order to cap the heterogeneity of arrays within a cluster, we empirically limited the maximum pairwise intra-cluster distance, d max , to 0.3. Second, in order to eliminate conditions or signals that did not have sufficient representation within the compendium, we required a minimum cluster size, c min , of 5. We generated representative nodes for clusters that passed these cutoffs by using average linkage to calculate centroids. We called these representative cluster nodes “meta-samples,” and combined them to create a meta-compendium that reduced the dimensionality of the data. We examined the sensitivity of the meta-compendium and independent components extracted to the parameters d max and c min and found that our results were robust with different parameters (data not shown). ICA identifies latent variables, or independent components, that provide a new set of basis vectors for a dataset. In other words, ICA performs dimension reduction by describing the data in the feature-space identified by independent components. As opposed to gene modules, which typically consist of unordered subsets of genes, components contain scores for each gene in the domain. To extract a robust set of independent components from the meta-compendium, we applied iterated ICA and averaged the results. Here we briefly describe the ICA model; for a more thorough review of the procedural details of the application of ICA to microarray data, see Kong et al. [ 64 ]. Let X represent the m × n matrix of microarray data measuring the expression of n genes in m experiments. ICA models this expression matrix as a linear combination of independent biological processes by decomposing X into a k × n source matrix S and a m × k mixing matrix A such that X = AS , where k is a user-supplied parameter no greater than the minimum of m and n . For our analyses, we set k = m to extract the maximum number of components from the meta-compendium. These components, or rows of S , are independent in the sense that the gene weights in each component reflect samplings of independent random variables. In the context of gene expression, this suggests that the sets of genes comprising the groups strongly contributing to each component have independent compositions. We performed ICA using the R interface to the FastICA algorithm (version 1.1–11), which attempts to maximize the non-Gaussianity of the component distributions [ 65 ]. For microarray data, these distributions are typically super-Gaussian: only a small number of genes contribute heavily to a specific biological process, while the majority of genes do not contribute [ 28 ]. We used the contrast function (2) g ( u ) = log cosh u and allowed a convergence tolerance of 0.0001. To account for the stochastic nature of ICA, we ran the algorithm 20 times and pooled the independent component estimates in a method similar to the Icasso software [ 66 ]. We clustered component estimates with partitioning around medoids [ 67 ], a robust version of k -means, using the distance metric based on Spearman’s correlation coefficient in Eq. (1) . We defined the fundamental components (FCs) of our analysis as the k medoids that result from this clustering analysis. To visualize the independent component estimates, we used classical multi-dimensional scaling to represent the 20,099-dimension component matrix in two dimensions. To examine the biological processes associated with each component, we made use of two sources of information yielded by the ICA decomposition: the distribution of gene weights in a component (rows of S ), and the distribution of a component’s expression in arrays in the compendium (columns of A ). The rows of the source matrix, S , describe the contributions of individual genes to each component. Each component contains an expression value for each gene. ICA produces components with unit variance and zero mean, so the contributions of each gene are relative. We call the genes that significantly contribute to each component the active genes , and identified them by using a weight threshold [ 37 ] of three standard deviations from the mean. Similar to Lee and Batzoglou [ 31 ], we generated two sets of active genes for each component: one for genes with positive loadings and one for genes with negative loadings. Our initial observations, as well as those of Frigyei et al. [ 38 ], suggested that this division results in a greater number of significant annotations. For a given component, we called the module with more genes the dominant module . Since the signs of independent component expression values are arbitrary, we reoriented each component so that the loadings for the dominant module were positive. To investigate the biological functions of our fundamental components, we used Alexa et al.’s elim algorithm [ 68 ] to identify enriched Biological Process GO categories in the active gene modules, assessing significance with Fisher’s exact test. We annotated each gene cluster with GO categories with a Bonferroni-corrected p-value <0.05. To further characterize the biological processes associated with each component, we identified significant canonical pathways using Ingenuity Pathway Analysis (Ingenuity Systems, www.ingenuity.com ). Ingenuity uses Fisher’s exact test to calculate a p-value determining the probability that the association between genes in the dataset and the canonical pathway is explained by chance alone. We defined significant pathways as those with a p-value <0.05 after correcting for false discovery rate with the Benjamini-Hochberg method. Independent components provide an orthogonal basis with which to describe the original data, so we modeled each original experiment as a linear combination of FCs. To generate a mixing matrix that describes the original arrays, we projected the full 9,395-array compendium into the feature-space defined by the FCs using Eq. (3) . (3) A = X S T This resulted in an (array × component) mixing matrix, A , whose columns describe the contributions of a component to each array in the compendium. A high absolute weight in this matrix for component i and array j indicates that the gene expression pattern present in component i plays a significant role in defining the expression profile of array j . Because component expression forms complex distributions [ 37 ], we annotated each FC with the 100 individual arrays in the compendium that have the highest expression of that component. In other words, we found the 100 arrays with the highest scores in each column of A , which we called the active arrays . The active arrays define the context in which a gene module is expressed. We annotated each gene module with the GEO Series containing the active arrays, ordering the GSEs by the number of active arrays contributed by each series. From the full compendium, we selected and normalized GSE7538 , a dataset consisting of 12 pairs of PTL-treated or untreated primary acute myelogenous leukemia (AML) specimens [ 8 ]. We projected this dataset onto the basis vectors defined by our fundamental components ( Equation 3 ) to determine the levels of component expression in each array. Analogous to the use of differential expression analysis to identify differentially expressed (DE) genes of interest, we used the quantitative expression levels of components to identify DE components. Specifically, we applied the linear modeling and empirical Bayes methods of the limma 2.18.2 R package [ 61 ] to identify differentially expressed components between treated and untreated AML samples. To compare this approach with existing DE gene identification methods, we used the moderated t-statistic from limma to identify DE genes from the original dataset. We also applied ICA to the PTL dataset alone, using the methodology described in Section 2.2, to identify the maximum possible number of independent components (24). Using only the data from GSE7538 , we again used linear modeling to identify DE components. For each comparison, we adjusted the p-values for differential expression using the Benjamini-Hochberg correction. We examined genes with p-values below 0.05 and components with p-values below 0.001. We also restricted our search to components whose differential expression had the same sign for all 12 treated/untreated pairs; while this strategy might eliminate true signal, it also removes inconsistent patterns that may represent heterogeneity within the tumors. To identify novel regulators of parthenolide response, we used the dominant gene modules from DE components to create protein interactions networks with Ingenuity Pathway Analysis (Ingenuity Systems, www.ingenuity.com ). Ingenuity overlaid active genes onto a global molecular network developed from information contained in Ingenuity’s Knowledge Base, and it algorithmically generated networks based on gene connectivity. We ranked genes in order of connectivity, and examined the highly-connected hub genes. To generate a list of known PTL-associated genes, we used the Pharmspresso text-mining tool to search the biomedical literature [ 69 ]. Pharmspresso searched all 322 PubMed articles with “parthenolide” in the title or abstract looking for co-occurrences of the term “parthenolide” with gene names in the same sentence, and obtained a list of 363 known PTL-associated genes.

Our starting compendium consisted of 20,099 genes measured by 9,395 heterogeneous arrays, representing 298 GEO Series (GSEs) that spanned a wide range of human biology. While a majority of these data series contributed fewer than fifteen arrays, we noticed that the five largest GSEs constituted over one third of the total arrays in the compendium. For instance, GSE2109 contributed 1,974 primary tumor samples from various sources and GSE8052 contributed 404 lymphoblastoid samples from an asthma study. A compendium-wide analysis using Spearman’s rank correlation as a similarity metric revealed strong correlations within datasets and between experiments performed with similar cell types and conditions ( Supplementary Figure 1 ). We reasoned that the disproportionate representation of some transcriptional patterns might cause certain biological conditions to dominate the ICA results. At the same time, we did not wish for rare gene expression patterns caused by extreme experimental perturbations or technical errors to skew the compendium. Therefore, to reduce the contributions of over-represented conditions and to eliminate cellular signatures that were not replicated within the compendium, we applied a filter based on hierarchical clustering to the original compendium ( Figure 1 ). The resulting meta-compendium consisted of 423 meta-samples that included contributions from 5,818 of the original arrays. We next applied iterated ICA to the meta-compendium to identify 423 fundamental components that span a wide range of human biology. Since ICA produces components in an arbitrary order, we sorted the components based on the amount of variance of the meta-compendium they explain ( Supplementary Figure 2 ). Unlike PCA, ICA does not use variance explained as a criterion to select components. Still, ordering components based on this metric may provide biologically relevant information [ 28 ]: fundamental components that explain a large percentage of the variance in a compendium represent processes that are active in many biological or experimental conditions, while those that explain a relatively small percentage of the variance represent relatively infrequent transcriptional patterns that are active only under certain conditions. In this ordering scheme, the first independent component explained 3.2% of the variance in the meta-compendium. Over multiple runs, ICA reliably reproduced the ten components that explain the most variance ( Figure 2a ). To verify that low-ranking components also clustered together, we examined a sampling of components that explained varying fractions of the variance in the data ( Figure 2b ). Signals appeared robust even for the less well-represented components. Henceforth, we refer to components by their rank in this ordering. That is, Fundamental Component 1 (FC-1) explains the most variance in the meta-compendium, while FC-423 explains the least. To investigate whether our method of meta-compendium creation succeeded in reducing the contributions of over-represented conditions, we compared the percentage of variance explained by each FC in the meta-compendium to that in the original compendium ( Supplementary Figure 3 ). Some FCs explain more variance in the meta-compendium than in the original compendium; these FCs may represent relatively rare but reproducible biology that was up-weighted in the meta-compendium. Many FCs, however, explain 1.5- to 3-fold less variance in the meta-compendium, and may represent biology that was over-sampled in the original compendium. We examined the active array clusters for two representative FCs, 62 and 268, selected for their high reduction in variance explained and their relative separation in the component ordering ( Supplementary Figure 3 ). FC-62 explains 1.26% of the variance in the full compendium but only 0.31% of the variance in the meta-compendium, and the 100 arrays with the highest expression of the component all are part of GSE8052 , the 404-array asthma study mentioned previously ( Supplementary Figure 1 ). Similarly, the 100 arrays with the highest expression for FC-268 consist entirely of brain samples from a single GEO series. Since FC-268 explains relatively little variance in either compendium, it likely represents an over-sampled gene expression signature that does not play a large role in most biological conditions. To examine the biological processes and functions associated with each component, we made use of two sources of information yielded by ICA decomposition: the distribution of gene weights in a component, and the distribution of component expression in arrays in the compendium. First, to identify genes that contribute significantly to each component, we used a loading cut-off of 3.0, corresponding to three standard deviations above the zero mean. For each component we created two gene modules, one for significant genes with positive loadings and one for negative loadings. The mean module size was 155.95, with a standard deviation of 73.57, although typically one module associated with a component was larger than the other; the mean difference between the sizes of a component’s modules was 110.00, with a standard deviation of 80.57. For a particular component, we refer to the module with more genes as the dominant module. 14,932 out of 20,009 genes belong to at least one module. We then calculated GO category enrichment in the active gene sets. 292 out of the 423 FCs had enriched GO categories in at least one of the two active gene modules. Supplementary Table 1 shows the enriched GO categories for the dominant modules of the first ten components. We observed a weak correlation between the amount of variance explained by a FC and the total number of enriched GO categories found for its active genes ( Figure 3 ). However, in accordance with Roden et al.’s observations with PCA [ 70 ], many FCs that explain less than 0.5% of variance in the compendium have strong GO annotation. This suggests that these FCs explain relatively infrequent but biologically coherent transcriptional patterns. To assess the utility of our fundamental components, we analyzed a dataset from the initial compendium (GEO accession GSE7538 ) using three methods: 1) standard differential gene expression analysis, 2) differential expression analysis based on independent components derived from the dataset alone, and 3) differential expression analysis using FCs derived from the full compendium of data. GSE7538 contains acute myelogenous leukemia (AML) CD34 + cells from twelve patients treated with parthenolide (PTL), a sesquiterpene lactone with antiinflammatory and anti-cancer activities [ 46 , 47 ]. Linear modeling with empirical Bayes methods in gene-space identified 2574 differentially expressed (DE) genes at a significance level <0.05 after correcting for multiple hypothesis testing, and we annotated this set of genes with significant GO categories and Ingenuity Canonical Pathways. We obtained more enriched terms with higher scores when separating the genes into up-regulated and down-regulated sets, so we present those results here. The significant pathways identified from the genes up-regulated in PTL-treated CD34 + AML cells included NRF2-mediated Oxidative Stress Response, Hypoxia Signaling in the Cardiovascular System, and Protein Ubiquitination Pathway ( Table 1 ), all of which are consistent with oxidative stress-induced apoptosis caused by parthenolide. Differential expression analysis in gene-space also marked Glucocorticoid Receptor Signaling, which is related to the NF-κB signaling pathway [ 71 , 72 ]. Next we applied iterative ICA to the 24-experiment data series, identifying 24 independent components. We projected the data into this feature space, and identified one significant DE component that explains 31.6% of the variability between treated and untreated cells. This method predicted that the Protein Ubiquitination Pathway and NRF2-mediated Oxidative Stress Response are up-regulated with increased significance compared to the gene-space analysis ( Table 1 ). GO enrichment analysis in this feature space did not yield as many enriched GO categories as the same analysis in gene-space ( Supplementary Tables 2,3 ). Finally we analyzed the differential expression in response to PTL treatment in terms of our fundamental components. We projected the original dataset into FC-space and identified 29 DE components with a p-value threshold of 0.001. To eliminate signals reflecting parthenolide response heterogeneity, we limited our results to the 19 components that were regulated in the same direction in all twelve paired samples ( Figure 4 ). Unlike DE analysis using independent components from GSE7538 alone, no FC explained a large fraction of differential expression as compared to the others; rather, complex manipulations of a wide range of gene modules contributed to the expression changes between PTL-treated and untreated cells. In addition to GO enrichment and Ingenuity Pathway Analysis, we identified for each DE component the experiments that contain arrays with the highest expression scores for that component. Although neither gene-space nor experiment-specific feature-space analysis identified NF-κB as a significant pathway in PTL response, FC-space analysis identified NF-κB signaling as the top pathway associated with FC-66, the fundamental component with the greatest differential expression ( Table 2 ). FC-66 was also enriched for JAK/Stat signaling and SAPK/JNK signaling ( Supplementary Table 4 ), two pathways with known associations with PTL response [ 52 , 53 ]. Listing the experiments that highly express this component revealed that FC-66 plays an important role in leukemia ( Table 2 ). In the second most over-expressed fundamental component, FC-362, NRF2-Mediated Oxidative Stress Response appeared at the top of the list of enriched Ingenuity pathways, echoing the findings in gene-space and experiment-specific component-space and matching previous experimental results [ 51 ]. Intriguingly, N-Glycan Biosynthesis appeared at the top of the list for the third most over-expressed fundamental component, FC-84. Subsets of glycan structures have specific functions that affect cell differentiation, growth and migration, and have been implicated with cancer development and metastasis [ 73 ]. Of the three most under-expressed FCs, only FC-62 offered descriptive GO categories and Ingenuity pathways. The dominant gene module for this component showed enrichment for EIF2 Signaling, and a closer examination of the gene list for this module revealed a preponderance of ribosomal proteins. This suggests that PTL modulates translational activity in the cell, perhaps as part of the stress response to the increase in reactive oxygen species. Unexpectedly, FC-66 was over-expressed in PTL-treated samples, seemingly countering current belief that PTL inhibits NF-κB activation. To determine the reason for this apparent reversal, we examined the expression of genes in the Ingenuity NF-κB Signaling pathway for both FC-66 ( Figure 5 ) and the original data from GSE7538 ( Supplementary Figure 4 ). This analysis revealed that genes in the NF-κB pathway exhibit a wide range of transcriptional changes in response to parthenolide, regardless of whether the genes play inhibitory or stimulatory roles. In the original data, IκBα (NIKBIA), an inhibitor of the NF-κB complex, is down-regulated in PTL-treated samples, while other components of inhibitory complexes are up-regulated (IκBβ, IκBε). Similarly, subunits of the NF-κB complex itself show varied transcriptional responses. However, the direct target of PTL, IKKβ (IKBKB), had a negative loading in FC-66 and was down-regulated in PTL-treated compared to untreated cells, suggesting that transcription-level regulation and protein-level activation are coordinated for this key gene in parthenolide response. To explore the biological validity of our method’s predictions on the gene-level as well as the module-level, we created a protein interaction network from the dominant module of each differentially expressed FC. For each FC, we ranked genes by connectivity in this network, since highly-connected proteins, or hub genes , are likely more essential to biological function and survival [ 74 – 76 ] and potentially represent therapeutic targets [ 77 ]. We defined our hub genes as the most highly connected 1% of the 3278 genes active in at least one of the 19 DE components, yielding a list of 32 highly-connected genes ( Table 3 ). To identify novel predictions, we compared this list of hub genes to our text-mined catalog of known PTL-associated genes. Fourteen of the 32 hub genes appear on this text-mined list, including key regulators of PTL response such as TNF, NFKB1, STAT4 and FOS. Among the genes not on the Pharmspresso list, several appear to represent known genes that the text-mining algorithm did not detect. C/EBP-β (CEBPB), for instance, is a transcription factor that interacts with a known PTL-associated gene, C/EBP-α (CEBPA), whose loss of function leads to the development of AML [ 78 ]. Kawasaki et al. found that regulatory targets of C/EBP-β are differentially expressed in PTL-treated prostate cancer cells [ 79 ]. HGF, although not previously directly linked to PTL, depends on the STAT pathway for HGF/SF-Met mediated tumorigenesis [ 80 ], and so may be linked to the PTL response mechanism. Other genes on this list may represent new genes and pathways involved in PTL response. CNBP and EIF3A both contribute to FC-62, whose dominant gene module is associated with protein translation. CNBP plays a role in controlling proliferation and cell survival, for instance by regulating the CT element of the c- myc protooncogene [ 81 ]. A group of highly-connected proteins – LCK, ZAP70, LAT, CD247 and CD3E – contribute to the T Cell Receptor Signaling pathway in FC-43 ( Supplementary Table 4 ). Some of these genes also appear in the NF-κB Signaling pathway ( Figure 5 ). Lewis and coworkers reported aberrant expression of T-cell markers in AML [ 82 ], and an independent microarray study implicated T-cell receptor signaling in AML stem cell function [ 83 ], suggesting that expression of this pathway may be relevant to AML response to PTL.

To develop a data-driven view of the human transcriptome, the biological community needs analysis tools that leverage the growing repositories of gene expression data. We assembled a large cross-study compendium of human microarray data and extracted fundamental components of human biology using independent component analysis. These components described gene modules that have biologically coherent functions, and we used them to improve the differential expression analysis of a microarray dataset involving the preclinical drug parthenolide (PTL). Because we identified gene modules through a meta-analysis of expression data from GEO, our analysis depends fundamentally on the quality of the primary microarray data. Variation in hybridization conditions and starting mRNA quantity and quality could introduce spurious signals whose origins are technical rather than biological. To address this issue, we applied the technical bias correction proposed by Eklund and Szallasi [ 58 ] to all GEO data. Our initial compendium also reflected a non-uniform sampling of human biology due to biases in the biological conditions commonly investigated ( Supplementary Figure 1 ). In a study to mine functional gene relationships from a large microarray compendium, Huttenhower et al. applied Bayesian normalization to account for the over-representation of certain biological or experimental conditions [ 84 ]. Similarly, we wished to derive a broad but even sampling of human biology; we did not want our fundamental components to describe a disproportional number of pathways unique to cancer, for example. We applied a grouping method based on hierarchical clustering to create a meta-compendium that down-weighted the contributions from over-represented conditions and eliminated expression profiles that did not appear reproducibly in the compendium. In the process, we averaged microarray profiles within the same cluster, possibly removing information about some of the less well-represented transcriptional subunits present within datasets. However, we wanted to focus on modules active in multiple experiments throughout our collection, so we chose to retain well-represented transcriptional processes at the expense of less well-represented ones. Additionally, this procedure for normalizing the weights of different biological conditions resulted in the removal of over 3,000 arrays. Of these arrays, approximately 1,900 of them consisted of tumor samples derived from a single compilation dataset (GEO accession GSE2109 ). We note, however, that even with this normalization, the scope of our results depends on the sampling of biology available. In our analysis, we used ICA to identify gene modules that vary together across our gene expression compendium. While Fehrmann et al. used principal component analysis (PCA) for a similar analysis [ 21 ], ICA provides a complementary method for analyzing gene expression data [ 31 ]. At the same time, ICA does not rely on prior knowledge for identification of gene modules; although prior knowledge may enhance the analysis of well-studied biological processes [ 16 , 29 ], we chose to use gene expression data alone to guide the formation of our gene modules. In this context, ICA has two major procedural limitations. First, the user must specify the number of components to extract. Like other studies using ICA, we searched for the maximum possible number of components (423, the number of samples in the meta-compendium), assuming that more biological processes exist than we can describe with this number of sources [ 31 ]. Second, ICA is a stochastic algorithm. To extract robust biological signals from the meta-compendium, we repeated ICA twenty times. In contrast to approaches that mitigate unstable signals by eliminating them entirely [ 37 ], we clustered the results to obtain more reliable component estimates [ 66 ] using partitioning around medoids. To maximize our ability to comprehensively model microarray data, we included all 423 components in downstream analyses. One final challenge associated with matrix decomposition methods involved identifying active genes. Although alternative methods exist [ 38 ], we followed the approach of Liebermeister [ 28 ] and used a gene loading threshold. To confirm the biological relevance of these gene modules, we devised a method to analyze differential expression in the feature-space defined by our fundamental components. This follows naturally from the typical assumption in work with ICA that the active genes in each component are differentially expressed [ 28 , 35 ]. Although several groups have used principal component features to enhance the detection of individual differentially expressed genes [ 85 ] or gene pathways from GO and KEGG [ 86 ], the opportunity to assess the differential expression of entire gene modules has not been adequately explored. We exploited the component loadings in the mixing matrix of the ICA decomposition and used the moderated t test [ 61 ] to assess the differential expression of globally-derived gene modules within a single dataset. Analysis in feature-space reduced the number of statistical tests performed, decreasing the need for multiple hypothesis testing correction and increasing the sensitivity of biological signal detection [ 87 ]. We empirically selected different p-value cutoffs for DE gene analysis (p=0.05) and DE component analysis (p=0.001) to optimize the performance of each method. DE analysis in feature-space outperformed similar analysis in gene-space, identifying more known PTL-associated genes as significant ( Figure 6 ). Additionally, the pathways enriched in differentially expressed FCs proved to be highly consistent with previously characterized PTL effects on cellular processes, identifying both NF-κB signaling and oxidative stress as major mechanisms involved in response to PTL treatment. Because we modeled a gene expression profile as a linear combination of components, however, we must interpret pathway predictions with care. An active gene in a DE component may not necessarily be a DE gene in gene-space analysis. In our investigation of the mechanism of parthenolide, FC-66 loadings of genes in the NF-κB pathway did not necessarily mirror the direction of regulation in differential gene expression analysis ( Supplementary Figure 4 ). Complex gene networks such as the NF-κB signaling pathway contain proteins with both inhibitory and stimulatory effects, so the finding that FC-66 was over-expressed in PTL-treated cells did not necessarily indicate that the entire NF-κB pathway was both transcriptionally and functionally up-regulated. Rather, different genes responded to this stimulus according to their role in the pathway ( Figure 5 ). As a further demonstration of this transcriptional complexity, we note that our method identified the hub gene Cox-2 (PTGS2) as differentially up-regulated in FC-406 and down-regulated in FC-101 (data not shown). While we were unable to determine the net direction in which PTL response influences Cox-2 expression, our result corroborates previous data suggesting that Cox-2 protein activity changes in response to PTL treatment [ 88 , 89 ]. This demonstrates the enhanced sensitivity of our approach, as DE analysis in gene-space did not identify Cox-2 as an important gene in PTL response. After deriving gene modules using a data-driven approach, we used ontological annotations to provide a preliminary indication of the functional relationships between active genes. We found that for some of these gene modules, Ingenuity Pathway Analysis provided more descriptive relationships between genes, especially in the realm of signaling pathways (e.g. FC-66, Table 2 ). Our analysis of the mechanism of PTL, a drug that affects the NF-κB signaling pathway, particularly benefited from this form of annotation. In other cases, GO Biological Process categories labeled components where Ingenuity did not (e.g. FC-101, Table 2 ). Still, annotation of our fundamental components with GO categories left 131 components unlabeled, suggesting that GO does not completely describe some functional gene relationships. Indeed, some of the DE fundamental components identified in our analysis of PTL had neither descriptive GO categories nor significant Ingenuity Pathways. FC-208, for instance, was highly expressed in several cancer-related experiments, but was down-regulated in response to PTL ( Table 2 ). This gene module may represent a cancer-specific transcriptional subunit that is suppressed by PTL. Similarly, FC-101 had vague GO annotations including “defense response,” but was highly expressed in three independent AML experiments ( Table 2 ). The three genes with the highest weights in FC-101 consisted of DNTT, whose expression in AML cells correlates with outcome [ 90 , 91 ]; SOCS2, which is part of a gene signature prognostic for an unfavorable AML subtype [ 92 ]; and DDIT4, whose mRNA levels are induced in AML cell lines and normal CD34 + cells during differentiation [ 93 ]. Thus our method detected a transcriptional subunit active in AML that is not adequately characterized by traditional functional analysis methods. Since ICA defines components in a manner that allows a gene to participate in more than one transcriptional module, our approach also provides a mechanism for modeling context specificity [ 13 ]. As we highlighted using the example of Cox-2, our method defines the functional role of a gene in the context of other genes. Given the independence properties of components derived by ICA, we interpret these results as defining several specific contexts in which each gene operates. Furthermore, since we can identify experimental conditions that lead to the co-expression of sets of genes, we can match gene modules with their biological contexts. For example, the set of genes in the dominant module of FC-101 might be labeled as “up-regulated in AML” and “differentially expressed in PTL response.” This context awareness will prove invaluable as we continue to explore the functional relationships between genes at the genomic scale.

As opposed to a top-down resource such as the Gene Ontology (GO), we provide an entirely data-driven way to identify gene modules in a meta-analysis of expression data using ICA. These modules contain functionally coherent sets of genes and increase the descriptive power of differential expression analysis. The results of our analysis of a parthenolide-related dataset suggests that the fundamental components we have identified may serve as useful, reusable features; although we applied these modules to examine a microarray experiment from our original compendium, we can use the same components to investigate new expression data. In addition to differential expression analysis, these gene modules could serve as data-driven gene sets for Gene Set Enrichment Analysis [ 18 ] or as features for clustering and classifying diverse microarray expression data. Given the wide range of biological pathways represented in the compendium, these components could be applied to query a microarray database for experiments that modulate similar biological processes. Finally, we can use this data-driven framework to annotate genes with context-specific functional relationships.