Determination of strongly overlapping signaling activity from microarray data

Background As numerous diseases involve errors in signal transduction, modern therapeutics often target proteins involved in cellular signaling. Interpretation of the activity of signaling pathways during disease development or therapeutic intervention would assist in drug development, design of therapy, and target identification. Microarrays provide a global measure of cellular response, however linking these responses to signaling pathways requires an analytic approach tuned to the underlying biology. An ongoing issue in pattern recognition in microarrays has been how to determine the number of patterns (or clusters) to use for data interpretation, and this is a critical issue as measures of statistical significance in gene ontology or pathways rely on proper separation of genes into groups. Results Here we introduce a method relying on gene annotation coupled to decompositional analysis of global gene expression data that allows us to estimate specific activity on strongly coupled signaling pathways and, in some cases, activity of specific signaling proteins. We demonstrate the technique using the Rosetta yeast deletion mutant data set, decompositional analysis by Bayesian Decomposition, and annotation analysis using ClutrFree. We determined from measurements of gene persistence in patterns across multiple potential dimensionalities that 15 basis vectors provides the correct dimensionality for interpreting the data. Using gene ontology and data on gene regulation in the Saccharomyces Genome Database, we identified the transcriptional signatures of several cellular processes in yeast, including cell wall creation, ribosomal disruption, chemical blocking of protein synthesis, and, criticially, individual signatures of the strongly coupled mating and filamentation pathways. Conclusion This works demonstrates that microarray data can provide downstream indicators of pathway activity either through use of gene ontology or transcription factor databases. This can be used to investigate the specificity and success of targeted therapeutics as well as to elucidate signaling activity in normal and disease processes.


Data analysis flowchart
Data analysis flowchart. The data was downloaded from Rosetta Inpharmatics and filtered to include only genes and experiments that showed significant variation. Bayesian Decomposition analysis generated patterns and associated gene lists for all dimensionalities between 3 and 25. ClutrFree was used to interpret these results, including use of the MIPS database of ontologies.

Background
Many diseases develop because of errors in signaling, and newer therapeutics specifically target proteins involved in cellular signaling [1,2]. However, these therapies are not always effective [3], and the reason for failure, whether inherent poor interaction or complex cellular response, is unknown. In order to understand the development of disease and drug resistance in these cases, the recovery of the process that led to the specific cellular malfunction must be identified. Such errors generally involve the cellular signaling networks that control cell growth, differentiation, apoptosis, and motility [4,5]. Because of the extreme underlying biological complexity of these pathways, diseases that involve errors in signaling processes arise from a myriad of different cellular malfunctions, for example in cancers [6,7] and diabetes [8,9]. It is from this complex background that functional genomics attempts to glean insight to improve our understanding of diseases.
One of the major uses of microarrays has been elucidation of gene expression in cancer, often focused on refining cancer identification using computational and statistical approaches [10][11][12]. In addition, the discovery of biomarkers in the form of differential levels of production of mRNA has been a focus in a number of studies [13][14][15]. The fact that determination of the mRNA levels of a single gene is easier than using an entire array has driven the shift to the use of arrays to generate potential biomarkers, so that the expression levels of these individual genes can be screened for in a more economical way (see, for example, [16]). For diabetes, microarrays have been used to elucidate gene expression in both type I and type II diseases, and customized chips targeting genes of interest have been developed [17].
While these techniques are useful, they have certain limitations as regards more advanced uses in the elucidation of mechanisms operating in diseased tissues. New therapeutics specifically target proteins involved in cellular signaling [1][2][3][38][39][40]. As noted above, these therapies are not always effective, and a method to understand the reason for their ineffectiveness is highly desirable. If the fail-ure modes for the targeted therapeutics are understood, new therapeutics can be designed or combination therapies undertaken. In addition, to design new therapies that work alone or in combination with other therapies, an understanding of signaling networks is required. Microarray measurements can provide insight into these issues.
Unfortunately, the recovery of pathway information from transcriptional data requires complex analysis, since signaling protein activity is not generally linked to the mRNA expression levels of genes encoding the signaling proteins themselves [41], nor are protein levels tightly coupled to transcript levels even in yeast [42,43]. This makes it impossible to directly link an increase in mRNA expression of the gene encoding a signaling protein, such as the therapeutic target, with activity of the protein and therefore of the signaling pathway. Instead, an analysis must treat changes in mRNA levels as downstream indicators of activity.
An important issue to resolve in order to correctly interpret patterns in microarray data is the underlying dimensionality of the data, since statistical analysis of genes in groups relies on correct separation. The dimensionality provides an estimate of the number of patterns required to explain the variation in the data not related to noise, which is equivalent to the number of basis vectors required mathematically to describe the data or the number of principal components required to span the data.
We present here a new application of Bayesian Decomposition [44][45][46][47][48] and ClutrFree [49] that estimates dimensionality by measuring the consistency of assignment of genes to patterns. With this approach, transcriptional signatures are linked to signaling activities through gene ontology [50] using the MIPS database [51] and through analysis of transcription factor activity [52]. We demonstrate this technique on the Rosetta deletion mutant dataset [53], which is a compendium of genome-wide transcription measured for 6300 genes across 300 conditions (mostly deletion mutants, but some chemical treatments). Figure 1 details the workflow of our analysis. Previous studies of the compendium were performed using hierarchical clustering [53], non-negative matrix factorization [54], and Bayesian Decomposition [44]. The dimensionality of the data was estimated in various ways in these studies leading to estimates from 7 to 50 dimensions.

Dimensionality of the data
We propose a value for the Rosetta dataset dimensionality based on the average persistence calculations at each tree level made with ClutrFree using multiple Bayesian Decomposition simulations. The dimensionality has been inferred from the average persistence defined in the Methods section. As the number of basis vectors (i.e., patterns) k is increased, the curve shows a dramatic drop for k > 15 (see Figure 2). This drop is due to the reorganization of the groups of mutants constituting the basis vectors for k > 15 patterns, leading to an overly low average persistence.
The freedom to move between branches leads also to some loss of consistency in the annotations as one moves down a branch. This contrasts with the behavior of basis vectors obtained for less than 15 patterns where biological functions split logically as the number of patterns increases. We also observed a reduction in the number of genes related to each basis vector for 16 patterns in comparison to 15. Also, the standard deviation across samples of the obtained vectors is significantly higher for 16 patterns (1.7 × 10 -3 ) than for 15 patterns (7.9 × l0 -4 ) indicating that the Markov chain sampling is not as tightly constrained by the probability distribution. The behavior observed occurs because of the potential to overfit the data with 16 basis vectors allowing the algorithm to find multiple configurations to explain the variation in the data.

Identifying patterns and functions
Bayesian Decomposition retrieves the two linked matrices: the P matrix (pattern matrix) groups mutants that share cellular functions, which can be deduced from the genes linked to each pattern contained in the A matrix.
Each mutant (a column of the P matrix) can belong to multiple patterns, which models the fact that each mutant will have many cellular functions active. Each gene (a row of the A matrix) can be assigned to multiple patterns, reflecting the fact that evolution has led to genes being involved in multiple cellular processes. Interpretation of the results involves identifying cellular processes from the genes that are significantly expressed in a pattern (i.e., within a column of A).
We proceed by using the dimensionality estimate of 15 patterns and exploring for each pattern the genes associated with that pattern. These genes are interpreted using the MIPS ontology for yeast [51] in order to predict the cellular processes associated with a pattern. In addition, for patterns that can be linked to signaling pathways, we discuss the use of data on genes regulated by specific transcription factors and validate the results by analysis of specific key deletion mutants. For each pattern that shows enhancement of ontological terms we provide the terms, the enhancement (as defined in the methods section), and the p value for a hypergeometric test on the term.
We summarize the results in terms of patterns previously identified in other studies using this data set, then we present the new results isolating signatures for activity of the mating and filametation pathways.

Patterns identified in previous studies
Examination of pattern 1 shows expression of the overall common minimal processes necessary for survival, with 386 annotated genes associated with this pattern at a 3σ level. Measure of enhancement, e, of cellular functions, reveals two highly represented functional groups: 1) groups related to protein synthesis and 2) groups related to DNA synthesis. Group  This pattern, which essentially includes genes necessary for viability, contains all the mutants of the dataset, although the Ssn6∆ mutant shows a lower level for this pattern than other mutants. As the Ssn6∆ mutant exhibits substantially greater overall expression than any other mutant (including the Tup1∆ mutant with the second highest level), this may reflect the high association of the Ssn6∆ mutant seen in almost all patterns, which will have some gene overlap with this pattern.   16 15 Dimensionality (K) Average Persistence P Pattern 5 contains 172 annotated genes. Highly enhanced ontologies include transport-related functions: Transported Compounds (Term 20.01, e = 2.02, p < 10 -4 ), Ccompound and Carbohydrate Transport (Term 20.01.03, e = 2.60, p < 3 × 10 -4 ), and Cellular Transport, Transport Facilitation and Cellular Routes (Term 20, e = 1.62, p < 6 × 10 -4 ), in addition to other transport terms at p < 10 -3 . The pattern contains the two deletion mutants, Ssn6∆ and Tup1∆, and represents the strong response seen in the original study [53]. Ssn6p and Tup1p form a system of transcriptional repression that appears to be highly conserved in eukaryotes [55]. In yeast, the complex acts as a global transcriptional repressor over a large number of genes (more that 150), coordinating several cellular systems, including haploid specific genes, glucose repressible genes, and oxygen utilization genes [56]. Turning off this repression leads to a large overall increase in gene expression (the overall expression in these two mutants is many fold higher than in other mutants). Pattern 7 is related to the lack of cell wall functions (Cell Wall, Term 42.01, e = 0.0), as 28 cell wall genes (of 32 total) are absent from this pattern, while the other four genes have multiple annotations suggesting they have roles unrelated to Cell Wall function. Enhancement is present for Protein Modification (Term 14.07, e = 4.7, p < 10 -4 ) and Fermentation (Term 02.16, e = 4.5, p < 0.01). This pattern contains the mutants Gas1∆ and Fks1∆, which impair cell-wall synthesis, as well as the mutant YER083c∆, annotated as disrupting the cell wall in the original study [53]. The pattern contains other mutants disrupting ergosterol biosynthesis as well, including Erg2∆, She4∆, as well as YER044c∆. In addition, the pattern includes yeast treated with the drugs that are known to disrupt the cell wall, such as Tunicamycin and Glucosamine.

Patterns related to cellular signalling pathways
The two patterns that represent new insights into this data are 13 and 15, which appear related to two strongly coupled developmental pathways in yeast. Previous studies [44,53] have identified the mating pathway transcriptional response, however this has included both the filamentation response and the mating response. It is difficult to separate these signatures, as the mating and filamentation pathways share many common elements in a MAPK cascade as shown in Figure 3 [57][58][59].
Gene ontology (GO) was used to determine the biological function described by each pattern, with a term added specifically for transposable elements, as these are known to play a role during filamentation [60,61]. The terms that showed enhancement are summarized in Table 1. The patterns show strong overlap, since many genes are shared between the mating and filamentation responses. However, the filamentation ontology term is significantly higher only in pattern 15, which also shows a strong signature of transposable element genes. Meanwhile, the GO terms for meiosis and morphogenesis (such as for budding in S. cerevisiae) are significantly enhanced only in pattern 13. This allows association of pattern 13 with activation of the mating pathway, and pattern 15 with activation of the filamentation pathway. The strongly linked MAPK signaling pathways for mating and filamentation are shown schematically with black arrows indicating mating pathway signaling and gray arrows showing filamentation pathway signaling. The mating pathway is initiated by binding to Ste2p or Ste3p receptors, while the causative molecular trigger for filamentation is unclear. The pathways share many components.

Ste2/ Ste3
Ste20 Ras2 mating signaling filamentation signaling In addition, we analyzed the 10 genes whose expression is most strongly linked to each pattern. These are shown in Table 2, which summarizes which genes are known to be regulated by the transcriptional activators related to the mating (Ste12p) and filamentation (Ste12p-Tec1p complex) pathways. The results show that the top 10 genes related to pattern 13 have 9 genes of known function, with 8 related to the mating response, of which five are known to be regulated by Ste12p. For pattern 15, 7 of the top 10 genes are known to be transposable element genes, with three other genes having unknown functions. This again links pattern 13 to mating and pattern 15 to filamentation.
In order to validate that the patterns were actually measuring activity of the mating and filamentation pathways, we explored deletion mutants related to these pathways [61,62]. The mating response in S. cerevisiae is mediated via a MAPK signaling cascade initiated by binding to the Ste2p or Ste3p membrane receptors (Figure 3).  [61]. Finally, the signal is greatly reduced for the Ste20∆ deletion mutant in pattern 15 relative to pattern 13, which agrees with previous work suggesting that the filamentation pathway is more dependent on Ste20p signaling than is the mating pathway [62].

Discussion
Microarrays and GeneChips™ have become the tools of choice for the investigation of genome-wide transcription in most biological systems. The resulting data comprises noisy estimates of transcription levels for roughly 6,000 genes in yeast to more than 20,000 genes in typical mammalian studies. Numerous statistical and data mining methods have been applied to this data in order to identify individual genes showing differential expression, to identify patterns related to physiological states, and to identify groups of genes comprising biological processes. These studies generally have not focused on the estimation of cellular signaling from the data, despite the prevalence of cellular signaling in many diseases.
As noted in the introduction, the recovery of signaling pathway information from transcription data requires complex analysis, since protein levels do not correlate well with mRNA levels and signaling protein activity is not generally linked to the mRNA expression levels of genes encoding the signaling proteins themselves. As such, changes in mRNA levels are limited to being a downstream indicator of activity. If a complete model for the transcription of genes, including all known transcriptional regulators and biological processes regulating transcription, was available, the inference of activity would be straightforward. Unfortunately, the network models and even gene annotations are still far from being complete. In addition, the growing evidence supporting the important role of non-coding RNAs in regulation of gene expression (including antisense transcripts and micro-RNAs, see for example [63][64][65]) further undermines the potential of using mRNA species as markers for proteins and their activities [66].
In order to overcome this incompleteness, we have created the method described here. We couple identification of transcriptional signatures with our Bayesian Decomposition algorithm to a consistency analysis for gene assignment to patterns determined by comparison of different dimensionalities using ClutrFree. This allows the identification of the correct dimensionality to be applied to subsequent ontology and transcription factor analyses. ClutrFree is also used to determine the ontological terms enhanced within each pattern and to obtain a list of genes tied to this pattern, which can then be linked to specific transcription factors. In this way, the biological processes associated with conditions can be identified, and inferences can be made on the activity of specific transcription factors. This then allows inference on the activity of signaling pathways, which cannot be obtained with methods previously applied to microarray data. Overall, the method requires many separate steps, each modeling an aspect of the biological system, in order to make proper inferences on signaling from the data.
In the application to the Compendium data presented here, our analysis was able to extract the common features for a set of mutants that eliminated related pathways. As in previous studies, the global transcriptional repressor complex Ssn6-Tupl has been isolated in a single group. In addition, patterns for cell-wall synthesis, ribosomal function, and the global functions necessary for continued viability of yeast were isolated. In contrast to previous analyses of this data, two pathways related to the MAPK cascade were isolated, one related to mating and the other to filamentation. Once the correct dimensionality was determined, Bayesian Decomposition was able to identify transcriptional signatures unique for each pathway. The assignment was validated by an investigation of the deletion mutants known to adversely affect these pathways.

Conclusion
Microarray studies have been widespread in biological and medical research, often focusing on identification of genes significantly correlated with various disease states. However, many diseases arise from disruptions in cellular signaling, and in these cases gene expression only provides a downstream indicator of signaling activity. This greatly complicates the analysis. The new approach introduced here recovered signatures allowing us to make validated inferences on strongly overlapping signaling pathways.
The results demonstrate that for Saccharomyces cerevisiae, the mating and filamentation pathways can be distinguished from transcriptional signatures determined from analysis of microarray data, despite the intrinsic high noise, confounding transcriptional activity, and tightly coupled nature of the pathways. The next step will be to apply these methods to more complex signaling networks in worms, flies, and mammals.

The Rosetta deletion mutant data set
The Rosetta Compendium comprises 300 conditions, including 276 deletion mutants, 11 tetracycline regulated genes, and 13 drug treatments, in S. cerevisiae growing in rich medium [53]. The data were generated from a two color cDNA microarray hybridization assay [67], and transcriptional profiles were measured both with technical replication and biological replication (151 mutants).
In parallel with the 300 experiments, 63 controls of wildtype S. cerevisiae were grown in identical conditions and compared against each other, permitting creation of a gene-specific error model. The data was downloaded from Rosetta Inpharmatics.

Data preprocessing
The data was filtered to retain only conditions characterized by at least a variation of 3 fold in a minimum of 2 genes. Then all genes that did not vary by 3-fold in at least 2 conditions were also removed, leading to a data matrix comprising 764 genes and 228 conditions.
The data used by Bayesian Decomposition included both the mean log ratio for each data point and the uncertainty in this measurement determined by the Rosetta error model. Since Bayesian Decomposition as applied here requires positivity [45], the log ratios were converted to ratios and the uncertainties propagated to uncertainties on the ratios. Although analysis of residuals suggested that seven dimensions fit the data [44], the analysis presented here suggests that this is due to overestimation of uncertainty in the data. Bayesian Decomposition is not highly sensitive to minor misestimations of noise however, so that this should not be a problem for this analysis.

Analysis with Bayesian Decomposition
Bayesian Decomposition (BD) has been applied to multiple types of data: in vivo spectroscopic data [68], medical imaging [69], microarray data from single cell organisms [44,45], mammalian model organisms [47], humans [70], and on phylogenomic sequence data [71]. A detailed description of the algorithm [46] and a review of applications [48] have been published.
Briefly, BD models the microarray data, comprising a matrix of estimates of the ratio of expression between the experimental condition and a control, as the result of the multiplication of two matrices describing behaviors across conditions (the P or pattern matrix) and the distribution of genes within those behaviors (the A or amplitude matrix). Naturally, the data, D, includes noise, so that the full relationship is defined by where p(A,P) is the prior probability and p (D | A,P) is given by the likelihood such that The prior is used here to require positivity and to minimize structure in the estimates of A and P [46].
The analysis with BD is similar mathematically to an analysis with singular value decomposition (SVD) or with principal component analysis (PCA), since all methods estimate two matrices that together reconstruct the data.
In both PCA and SVD, orthogonality conditions force each row of P to be linearly independent, deriving P from either the data matrix (SVD) or the covariance matrix (PCA) using deterministic algorithms. Since patterns of expression related to biological processes will not generally be independent, BD uses the MCMC approach to avoid orthogonality. The resulting rows of P are usually easier to relate to biological processes than those from SVD or PCA.
BD was run at each posited dimensionality, K (as in equation 1), between 3 and 25, generating a mean estimate and an uncertainty (i.e., standard deviation from samples of the posterior probability distribution) for each element of A and P. The dimensionality is equivalent to the number of patterns, as the patterns can be viewed as nonorthogonal basis vectors for D.

Persistence and dimensionality
The results of Bayesian Decomposition across the different estimated dimensionalities, K, were compared with ClutrFree in order to visualize stable basis vectors and a persistence measurement on them [49]. Each independent BD analysis provides a tree level, with each pattern represented by a node (see Figure 4). The analysis with the fewest patterns is placed at the top of the tree, and additional levels of analyses are added creating a tree from fewer to larger numbers of patterns. Connections in the tree are created in a greedy way. Level N+1 (e.g., 4) is connected to level N (e.g., 3) by finding the node in level N+1 with the highest correlation to a node in level N. The correlation is given by the Pearson correlation for the strength of the assignment of mutants to a pattern (i.e., between P kj for the nodes). These nodes are connected and removed. From the remaining nodes, the maximum correlation between a node at level N+1 and one at N is found again, and this process is repeated until only a single node remains at level N+1. This node is then connected to the node at level N that yields the highest correlation.
Following our previous work [49], we use a measure of persistence to quantify the robustness of a pattern across the variation of the number of patterns. An example calculation is shown in Figure 5. The assignment of a mutant to a pattern (i.e., a node) is binarized based on the mean and uncertainty of the assignment of the mutant to the pattern determined by the MCMC sampling, using a requirement that the mutant be assigned to the pattern at greater than 3σ above zero. In Figure 5, it is assumed there are four mutants in a pattern, thus there are four binarized values. Then at each node, each mutant is compared for consistency in presence of the mutant in linked nodes within the tree. For example, for the highlighted node in Figure 5, the first mutant is present in the node above and the node below, so it is present in all 3 connected nodes. For the second mutant it is present in 2, and the third and fourth mutants are absent. The average persistence for the node is therefore (3+2+0+0)/4 = 1.25 as noted in Figure 5.

Relationship of patterns across dimensionalities
We assessed the dimensionality of the data using measurements of persistence. The persistence was measured for analyses from 3 -25 patterns and the dimension chosen where a significant drop occurred in an otherwise slow monotonic decline, which was expected due to the branching nature of the tree. Figure 2 shows the significant drop between 15 and 16 dimensions, so 15 patterns were chosen for further analysis.

Ontology and function
To assign ontological terms to the genes contained in basis vectors, we annotated our data using the gene ontologies from the Comprehensive Yeast Genome Database (CYGD) hosted at the Munich Information center for Protein Sequences (MIPS) [51,72]. The analysis here utilizes the Functional Catalog (FunCat) format that describes each gene using a hierarchical ontological model.
Similar to persistence, we defined enhancement as a measure of the over-representation, or under representation, of a gene function in a subset of the data [47]. It is the ratio of the frequency of occurrence of genes annotated by a particular ontological term in the pattern to the frequency of occurrence of the same term in the whole dataset, with g p being the number of genes annotated with the term t in pattern p, n p the total number of genes in the pat-tern p, G the number of genes annotated by the term t in the data, and N being the total number of genes in the dataset. In addition, we apply a hypergeometric test to estimate a p-value for each term. Function was then determined by inspection of enhanced ontological terms.

Transcription factor analysis
The genes were also analyzed for the patterns determined to be related to signaling pathways by exploration of the ten genes most strongly associated with the pattern. Each gene was analyzed using the Saccharomyces Genome Database [73] to determine whether it was known to be associated with mating or filamentation processes and to determine if it was directly regulated by the Ste12p transcription factor. A sample calculation of the average persistence for a single node Figure 5 A sample calculation of the average persistence for a single node. The average persistence is calculated by comparing the persistence at each node in the tree given in Figure 4. Each assignment of each mutant (4 are shown here) to a pattern is binarized as described in the text, then the average persistence for a node is calculated by checking on the number of times the mutant assigned to the pattern occurs in the connected nodes. The mutant can occur in any branch below the node of interest to be considered as present. If it occurs in multiple child nodes at a single level, that is still treated as a single occurrence for that level. The average for a dimension is then the average of the persistence of all nodes at that level.