Learning a Prior on Regulatory Potential from eQTL Data

Method Overview

We briefly review the Lirnet method, referring to the Methods for a full description. Lirnet uses genotype and expression data of genetically diverse individuals (eQTL data) and aims to learn a regulatory prior concurrently with reconstructing a regulatory network. Building on earlier work [7],[14], Lirnet clusters genes into modules with the assumption that expression of the target genes in each module is governed by the same regulatory program. As with several other methods for the reconstruction of regulatory networks, Lirnet can accommodate two types of regulators: values of genotype markers (genotype regulators), representing genetic polymorphisms on chromosomal regions [1]–[3]; and expression levels of genes that are known to have regulatory roles (expression regulators), representing activity levels of genes that might regulate that module [7],[9],[14],[15].

Lirnet's regulatory programs are based on linear regression, a choice designed to allow for the incorporation and learning of regulatory potentials. For each module m, the expression levels of genes in the module (denoted by y_m,_j) are modeled as a linear regression of candidate regulators (denoted by x₁,…,x_n): y_m,_j∼w_m,1x₁+w_m,2x₂+…+w_m,nx_n, where all genes in the module share the same weights w_m,k A regulator r that has a zero weight w_m,_r has no effect on the expression of the targets in module m. A biologically plausible regulatory program should have a small number of regulators with a non-zero weight. To achieve this goal, we use the LASSO method [16] to select only the most significant regulators. In its simplest form, LASSO adds a fixed penalty term to the objective function that introduces a uniform bias towards sparsity in the weights w.

Lirnet, inspired by our recent work on feature selection [17], incorporates regulatory potential by allowing different regulators to have different sparsity biases; a regulator r whose regulatory potential is low for a module m will have a stronger bias towards w_m,r = 0, and thus a lower probability of being selected as an active regulator for that module. The regulatory potential C_r of regulator r is defined to be a function of its regulatory features f_r. The method flexibly accommodates any property or combinations of properties of a regulator (Figure 2, Tables S1, S2, S3) that might be indicative of its likelihood of having a causal effect on its targets. These features can include features of the regulator alone, such as the location and significance of sequence polymorphisms, the function of the gene (transcription factor, signaling protein, etc.), and conservation of the polymorphic site. They can also include features that involve both the regulator and the targets, such as the enrichment of the module with genes having known relationships to regulator (e.g. transcription factor targets). The regulatory prior β encodes the importance given to each regulatory feature (see Methods), and is automatically learned from the data, allowing less relevant regulatory features to be ignored and others to manifest their significance.

The learning algorithm of Lirnet jointly estimates w_m,r's and β by maximizing a joint objective that involves both. The algorithm iterates three steps until convergence (Figure 1): (1) learning the regulatory program for each module by estimating the weights w_m,r's; (2) learning the regulatory priors β that reflect the importance of each regulatory feature; (3) computing the regulatory potential for each candidate regulator, module pair, based on the current β, thereby biasing their selection in the next iteration's regulatory programs. The output of Lirnet is thus threefold. First, it constructs a set of learned regulatory programs for the modules used in the analysis; for module m, these are all the regulators r with a non-zero weight w_m,r. Second, it constructs a quantitative regulatory potential both for genes and for specific sequence polymorphisms within them, allowing us to rank candidates for the causal sequence variation and to prioritize hypotheses for further testing. Third, it produces a set of regulatory priors, which may provide insight on the properties of a polymorphism that tends to induce an effect on its downstream targets.

Statistical Evaluation on Yeast and Human Data

To test the versatility and generality of Lirnet, we applied it to two very different eQTL data sets. The first is the eQTL data set of Brem and Kruglyak [3], which measured the mRNA expression and genotype of 112 S. cerevisiae strains derived as the F2 progeny of a BY/RM cross. The second is the expression data measured in the lymphoblastoid cell lines of the 60 unrelated HapMap individuals in the CEU data set, using only the genotypes of a subset of markers (the 500 K tag SNPs on the Affymetrix chips) as regulators (see Methods).

Figure 2 shows the regulatory priors of the most significant regulatory features in each of these data sets, as identified by the Lirnet algorithm. Several aspects of the features automatically chosen as important by the algorithm are revealing. First the top learned regulatory potentials between two organisms as different as human and yeast are remarkably consistent, supporting the robustness of our approach. For example, aside from the feature indicating the stop codon, which affects only 43 genes, the strongest positive weight on the regulatory potential in both datasets is given to a feature denoting whether the sequence variation correlates with changes in the expression level of the closest gene (cis-eQTL). This gene-level feature indicates that the polymorphism is already causal towards a change in the cell (the expression level of the gene) and hence may have additional downstream effects. The second most significant regulatory feature, also in both data sets, is the conservation score, consistent with the hypothesis that changes in residues conserved across millions of evolutionary years are more likely to have a causative influence. Also with a significant weight are features that evaluate the functional relevance of the position of the polymorphism, and the significance of the actual change. These features include, for example, the presence and type of amino acid changes, and polymorphisms in the 5′ or 3′ UTR. Surprisingly, in both data sets, synonymous SNPs weigh more heavily than non-synonymous ones. This phenomenon might arise from the fact that synonymous SNPs can have an effect on translational efficiency or mRNA destabilization [18], consistent with recent findings that such SNPs are under significant purifying selection [19],[20]. The selection of this regulatory feature by our automated method lends support to this hypothesis, and is worthy of further investigation.

Aside from sequence-based features, we also provided the method a rough categorization of gene function, allowing it to learn which types of genes are likely to play a regulatory role. Despite receiving no prior knowledge about the relative importance of the various functional categories, the method automatically assigns high weight to functional categories with a regulatory role: In the human data, the highest weights among those features are given to genes involved in cell death, transport, cell growth, signal transduction, transcription, and cell communication. In the yeast data, ‘transcription regulator activity’, ‘telomere organization and biogenetics’, ‘protein folding’, ‘glucose metabolic processes’, and ‘RNA modification’ are chosen to be important. Notably, other work supports the differences between BY and RM in many of these processes, including glucose processing (BY and RM demonstrate dramatically different growth rates on a number of carbon sources including glucose; D. Pe'er, unpublished data) and telomere organization [7], showing the value of allowing the regulatory priors to be tailored to particular data sets and organisms.

Lirnet can also take advantage of other functional data, when available. For example, in the yeast data, the pairwise feature derived from ChIP-chip binding between the regulatory gene and targets in the module received a relatively high weight. We note, however, that the method is also effective when such data are not available for a given transcription factor and set of target genes (see Oaf1 example below) or when the features themselves are not available, as in the case of the human data.

We next tested whether the learned regulatory potentials improved the quality of the learned regulatory program, by computing the proportion of genetic variance (PGV) explained by the learned program. We compared Lirnet with a uniform regulatory potential (hereafter “flat” Lirnet), Lirnet with a learned regulatory potential, a standard single-marker linkage (as in [1]–[3]) and Geronemo [7]. The results (Figure 3) demonstrate that Lirnet explains a dramatically larger fraction of the variance for a much larger set of genes than all the other methods. For example, in the yeast data, Lirnet with learned regulatory potentials explains over 50% of the PGV for 1,644 genes, compared to 1,457 genes for flat Lirnet, 828 genes for Geronemo and 230 genes for the method of Brem & Kruglyak. The same advantage translates across the spectrum of PGV values, and is arguably even greater at the tail, where many genes that are very poorly explained by other methods have a considerable fraction of PGV explained by Lirnet. A more refined PGV analysis, with an independent test set, shows that this dramatic improvement does not arise from overfitting to the test data (Figure S1). In fact, the Lirnet model has a comparable number of parameters to Geronemo and fewer parameters than the method of Brem and Kruglyak, due to the use of modules.

Notably, even flat Lirnet considerably outperformed both the single-marker linkage and Geronemo methods, suggesting that Lirnet would still perform well in cases where the learned regulatory potential was less informed. The gap between Lirnet and flat Lirnet appears to increase when using more markers (see Figure S1D for human results on 100 K markers), and is larger in the YRI than in the CEU data. These results are consistent with a hypothesis that the benefits of a non-uniform regulatory prior are more pronounced when regulatory potentials are aggregated over smaller regions of linkage disequilibrium. Thus, one can expect the benefits of Lirnet's learned regulatory priors to grow as we move to denser genotyping.

Recovery of Known Regulatory Pathways

We evaluated how well the learned regulatory program recovers known regulatory interactions. We begin with a comprehensive analysis of the quality of the learned regulatory programs, demonstrating that they are consistent with other sources of data indicating regulatory interactions that were not provided to the Lirnet algorithm. We then provide a comparison to the state-of-the-art method of [9], which was applied to the same data.

As a gold-standard set of regulator-target relationships is not available, we constructed a comparison test set from various datasets: deletion and over-expression microarrays [21],[22]; chromatin immune-precipitation (ChIP-chip) binding experiments [23]; mRNA binding pull-down experiments [31]; transcription factor binding sites [65]; and a literature-curated set of signaling interactions from the Proteome database (http://www.proteome.com/). Although each of these data sets has its own limitations in terms both of false negatives and of false positives, agreement with these orthogonal data sources is a reasonable metric for evaluating the quality of a method's predictions. For a prediction that a regulator r regulates a module m, we defined it to be validated if there was significant overlap (hypergeometric p<0.01) between the members of m and the putative targets of r, suggested by one of the above datasets. We note that none of these datasets was used for constructing the regulatory features for Lirnet: Lirnet used only the ChIP-chip data set of [24], and all regulator-target pairs that appeared in these data were removed from the evaluation data [23].

Most of these data sets focus on regulatory relationships where r is a transcription factor, whereas Lirnet and the other methods we evaluate are also capable of identifying cases where r plays a different regulatory role, such as signaling, chromatin modification, or RNA degradation. Therefore, to increase the coverage of our validation effort, we also considered indirect regulatory relationships (Two-Step Cascade in Table 1), where a method predicted a regulator r that has some close relationship with a transcription factor t, and t is confirmed in the above data sets to regulate m. We considered cases where t and r have a reliable protein-protein interaction (PPI) (Xenarios et al. 2000); and cases where r phosphorylates t in the Proteome data set.

Table 1 summarizes the number of validated regulators for various models, applied to the same set of modules: Lirnet with a uniform regulator potential, Lirnet with the learned regulatory potential, Geronemo, and a random model. (See also Table S5 for a full list of Lirnet predictions and their support.) Overall, Lirnet recovers a larger fraction of the known regulatory interactions than the other methods. We note that the reference set supports only a subset of the predicted regulatory interactions. This fact is not surprising, as the data sources used for constructing the reference set focus on transcriptional regulation, whereas Lirnet and Geronemo cover a much larger set of regulatory relationships. Although we have made some attempt to expand our reference set to cover signaling interactions, the data set of literature-curated signaling interactions is only a small fraction of the total set of signaling interactions that presumably hold in yeast. Moreover, as shown in our previous study [7], a large part of the regulatory interactions in these data sets represent chromatin modification and post-transcriptional regulation, which are not represented in our reference set.

The BY/RM cross also exhibits a large amount of cis-regulation, which we did not explicitly model in the Lirnet analysis. Nevertheless, of 492 cis-regulated genes – those whose nearby marker is significantly predictive of its expression level (t-test p-value<1e-5; 12.8% of the 3152 genes used in our analysis), 307 genes (76.4%) are assigned to modules with cis-regulatory programs. More specifically, 149 genes are assigned to modules that have the genes' nearby markers as genetic marker regulators; 158 other genes are assigned to modules that have their nearby expression regulators having their markers in the regulators' regulatory programs, suggesting an indirect cis-effect from a locus to an expression level of a regulator in that region to a target.

We also compared our results with those of the recent Bayesian network method of [9]. As this method infers regulatory interactions for individual genes rather than modules, we used the number of regulator-target pairs validated in the reference set as our evaluation metric. For a fair comparison, we removed from the reference set any data sets that are used for learning either model, leaving only the deletion and over-expression microarray data [21],[22]. The results (Figure 4A) show that, for various levels of expression change in the deletion or over-expression data sets, the regulatory interactions inferred by Lirnet are more consistent with previously known regulatory relationships. In Figure 4A, we see that Lirnet recovers many more supported regulators than the Zhu et al. method. This large discrepancy is partially due to the fact that their method largely focuses on transcription factors, and hence is incapable of picking up many of the regulatory relationships that are uncovered by Lirnet. However, even if we focus attention only on the relationships that could be compared to the deletion/over-expression data (Figure 4B), and evaluate the fraction that were validated in these data, we see that Lirnet significantly outperforms the Zhu et al. method.

Predicting the Causative Sequence Polymorphism in a Large Chromosomal Region

One of Lirnet's key features is its ability to identify a specific causative regulator in a linked chromosomal region, an ability also presented in several other recent methods [9],[10]. Since other methods for comparison focused on identifying a causal gene not a specific SNP, we also prioritized genes for a direct comparison. We first compute the regulatory potential of all SNPs in the region of interest, relative to each module (to account for pairwise regulatory features). We then rank each gene based on the highest-scoring SNP associated with it. Figure S2 shows the overall distribution of the regulatory potentials for both SNPs and genes. We can see that the vast majority of SNPs and of genes have a fairly low regulatory potential. The distribution begins to tail off at a regulatory potential of about 0.694; only 0.22% of SNPs and 2% of genes have a regulatory potential that exceeds this value.

We first compare to the recent work of Zhu et al., who focus on 13 “hot spots” – chromosomal regions identified to regulate expression levels of a number of genes in a previous study [1]. Previous work has identified regulators for several of these hot spots, and Zhu et al. provide new experimental validation for a number of others. To compare to these results, we considered the genes linked to each hot spot as a “module”, and applied the learned regulatory prior (Figure 2) to individual SNPs in each region to compute each gene's regulatory potential (see Methods). We sorted the genes in each hot spot based on their regulatory potential, and listed the top 3 genes for each hot spot. Figure 4B compares the result of the suggested causative genes in each region between Lirnet and the method of Zhu et al. (2008). We see that, of the top Lirnet regulators, 14 regulators, spanning 11 hot spots, have experimental support (see Methods), in comparison to 8 regulators (7 hot spots) in the analysis of Zhu et al. Even if we consider only Lirnet's top regulator for each region, there is experimental support for 10 regulators.

We also compare to the recent method of Suthram et al. [10], which improves on earlier work of Tu et al. [11]. These methods consider a gene and a chromosomal region to which it is linked, and analyze the flow in a protein-protein interaction network to select a particular causal regulator within the region. Suthram et al. validate their results relative to a pre-defined set of 548 regulatory relationships, extracted from gene knockout or overexpression microarray studies [21],[25], similarly to our analysis above. The predicted network of Suthram et al. was not available, so we evaluated Lirnet using their protocol and the reference set, to allow for a direct comparison. For each target gene and linked region, we selected, as the Lirnet predicted regulator, the gene whose regulatory potential in the region was highest. We then evaluated these predictions using the reference set of Suthram et al. The results, shown in Table S6, show that Lirnet significantly outperforms both the method of Suthram et al. and the previous method of Tu et al. [11], according to this evaluation metric.

Biological Results

We also performed an in-depth analysis of some of the specific regulatory modules produced by Lirnet for the yeast data set, and evaluated its ability to identify both the correct regulators and the specific polymorphisms that gave rise to the expression change in the targets.

Transcriptional Regulators

One example of the predictive power of assigning regulatory potential to individual SNPs within a large chromosomal region is the Zap1 module (Figure 5A). The module contains ten target genes and two major regulators, the gene expression pattern of ZAP1, which encodes a transcription factor (TF) that activates genes in response to Zinc [26],[27], and a genetic region on chromosome 10 that contains ZAP1. Of the ten target genes in the module, six were among 40 probable Zap1 targets based on the presence of a consensus ZRE element and RNA expression patterns in zinc and in the absence of Zap1 (p = 5.7×10⁻¹⁰) [27]. While the causative role of Zap1 in this data has previously been affirmed a number of times [3], [7]–[9], Lirnet automatically identified polymorphisms within Zap1 as the ones most likely, within the linked region, to play a causal role (Table S7). The regulatory potential of the identified SNP is the highest over all yeast SNPs (Figure S2). The most significant regulatory feature by far in this identification was the known binding relationship between Zap1 and two of its target genes, but other features also played a role (Figure 5B). Thus, the method has identified a TF-target relationship for which there is significant biological support.

Importantly, however, Lirnet is also able to predict such relationships when relevant functional data such as binding assays are not available. One such example is the peroxisome module (Figure 6A), containing ten genes that are enriched for processes related to fatty acid metabolism (hypergeometric p<4.8×10⁻⁶) and peroxisome organization and biogenesis (hypergeometric p<5.5×10⁻⁶), nine of which we considered for further analysis (see Methods). Lirnet suggests two regulators: expression level of PIP2 (alias: OAF2), a gene that encodes a Zn(2)-Cys(6) TF that heterodimerizes with Oaf1 to regulate genes involved in peroxisomal functions via an ORE element [28], and a genetic region between nucleotides 51,324 and 52,943 on chromosome 1. Of the 11 genes in this region of chromosome 1 (Figure 6B; Table S8), Lirnet selected polymorphisms within OAF1 as having the highest regulatory potential; the regulatory potential value of OAF1 is within the top 1% over all genes (Figure S2). Several forms of data support the role of the Oaf1/Pip2 heterodimer in regulating this module. Of the nine target genes analyzed in the module, six contained the canonical ORE motif (p = 1.8×10⁻⁶). Moreover, five were in the top 1% of most significantly down-regulated genes in a microarray experiment that compared RNA expression levels of an oaf1Δ versus a wild-type (BY) strain under inducing conditions [29] (p = 8.0×10⁻⁹). We note that the PIP2 promoter itself is Oaf1-dependent and contains an ORE element [29]. Thus, differences in PIP2 expression patterns across the 112 segregants are also likely to be partly due to polymorphisms in OAF1. However, the fact that PIP2 expression was selected in addition to the OAF1 genotype demonstrates Lirnet's ability to identify multiple relevant regulators, even when they are correlated. We note that OAF1 was not identified by previous methods analyzing this data set.

The BY allele of OAF1 contains two non-conservative coding polymorphisms at conserved positions that are likely to alter its function: an R70W polymorphism in the DNA binding domain (the highest-scoring SNP) and a Q447P polymorphism in the ligand binding domain [30]. OAF1's high regulatory score (Figure 6B) is a combination of the correlation between its expression and genotype (cis-regulation), shared GO process annotations with the target genes, the presence of non-synonymous coding mutations and their effects on protein properties (e.g., pKa and pI), and (to a lesser extent) its function as a transcriptional regulator. Importantly, the ChIP-chip data set [24] used in our analysis did not contain Oaf1 binding information and therefore did not influence the choice of the regulator, demonstrating Lirnet's effectiveness even when binding data are not available. Moreover, another gene in the same region of chromosome 1 (PEX22) shares common functional annotations with many of the target genes, yet received a significantly lower score. This highlights the fact that functional annotations, although a useful source of prior knowledge, are not the primary cue used by Lirnet. Both of these characteristics are likely to play an important role in the analysis of data from other organisms, where binding data and functional annotations are both limited.

Post-Translational Network

Lirnet also suggested an intriguing hypothesis regarding a cascaded pathway involving two modules. The first (hereafter the “Puf3 module”) contains 153 co-expressed genes, highly enriched (p<10⁻⁹¹) for nuclear genes with mitochondrial functions, and very highly enriched (p<10⁻¹³⁰, Figure S3A) for genes whose mRNA transcripts are bound by the sequence-specific RNA-binding protein, Puf3 [31]. This enrichment is specific to Puf3 binding and not just coincident with the preponderance of mitochondrial genes (Figure S3B). Also, the Puf3 targets that are in the module show higher Puf3 motif score than those not (Figure S3C; Text S1). Indeed, RNA expression levels in a puf3Δ mutant in rich medium (YPD) showed a significant up-regulation (p<10⁻³⁷) mRNA levels of the module genes (Figure S3D). While these results suggested that the highly coherent expression profile of this module was due, at least in part, to regulation of RNA stability via Puf3, we found that neither PUF3 mRNA expression nor its genotype is correlated with expression of the module genes (Figures 7A, S4), suggesting that Puf3 itself was not the regulator driving the observed variability across the strains.

The Lirnet analysis identified several genes as being involved in the regulatory program of the Puf3 module (Figure 7A). Towards the top of the list, we find DHH1 (ranked 1^st) and KEM1 (ranked 4^th), two components of the dynamic cellular structures called cytoplasmic processing bodies (P-bodies) [32]–[35]. P-bodies are sites of RNA storage [32] that can modulate mRNA translation or degradation: RNA transcripts are translationally silenced while stored in the P-body [36] and can be subsequently degraded or released back into the translating pool [37]. P-bodies contain the catalytic subunits of the mRNA de-capping enzyme Dcp1/Dcp2 [38],[39], whose activity is regulated by Dhh1. However, the signals for determining mRNA localization to P-bodies and subsequent degradation or release have not been identified [32]. Thus, Lirnet suggested an intriguing regulatory connection between the Puf3-bound transcripts and a known posttranscriptional regulatory complex (P-bodies).

If Puf3 serves as a regulatory signal in one or more of the processes associated with P-bodies (RNA targeting, degradation, or release back into the translating pool), we would expect Puf3 protein to be localized to P-bodies. We therefore used fluorescence microscopy to test the subcellular localization of Puf3 in wild-type BY cells. Indeed, under certain conditions (see Methods), a Puf3-GFP fusion protein formed bright punctuate spots in the cytoplasm which co-localized with those of known P-body components, Dhh1 and Edc3 (Figure 7B). These results are consistent with those of a previous study [31] that reported punctuate cytoplasmic Puf3-GFP fluorescence in the BY strain background, but did not test for co-localization with P-bodies.

This finding demonstrates the role of p-bodies in the regulation of the Puf3 module genes, but does not elucidate the causal SNP responsible for the difference between strains. To identify this SNP, we explore the Lirnet predictions for the regulatory program determining P-body expression. Dhh1 and Kem1 are themselves members of another module that we call the post-transcriptional regulatory (PTR) module (Figure 8A). This module also contains other regulators of the Puf3 module, including GCN1 and GCN20, two members of a complex that regulates translational repression in response to nutrient starvation. Many other module members are associated, directly or indirectly, with post-transcriptional regulation (Figure S6). The sole regulator of the PTR module is a genotype marker located on Chromosome XIV (Figure 8A), the same region that had previously shown as linked to several of these targets and members of the Puf3 module [1],[40]. The smallest region linked to the PTR module spans more than 30 genes (Figure 8B; Table S9), making a systematic evaluation of the candidates a significant effort, even in a genetically tractable organism like yeast. We therefore used Lirnet to evaluate the genes in the region based on their regulatory potential for the PTR module. As shown in Figure 8B, MKT1 is ranked as the highest scoring gene in terms of the learned regulatory potential. The regulatory potential of MKT1 is within the top 1% over all genes (Figure S2).

Mkt1 interacts with the Poly(A)-binding protein associated factor, Pbp1, and is present at the 3′ end of RNA transcripts during translation [13]. Mkt1 contains a highly-conserved nuclease domain, homologous to the human XPG endonuclease [13] and is required for the translational regulation of an Mpt5/Puf5-dependent transcript (HO) [13]. In the BY strain background, MKT1 harbors two non-synonymous coding mutations (Figure S7). The first variation lies in its putative nuclease domain and makes a non-conservative amino acid change (G30D) in a residue highly conserved in yeast. Both this and the more conservative mutation (R453K) are outside a Pbp1 interacting domain, which maps to the C-terminus of the protein (residues 601–760) [13]. MKT1 has thus far been implicated in four diverse processes, including propagation of a double stranded RNA virus [41], growth at high temperature [42], efficiency of sporulation [43], and gene-specific translational regulation [13]. In many of these cases the BY variant, and, where examined, the G30D mutation specifically, appears to introduce a loss of function mutation [41],[43],[44]. Both these SNPs contributed significantly to MKT1's high regulatory potential, with the G30D SNP scoring somewhat higher. The properties that contributed the most to this selection are conservation, cis-regulation, the amino-acid properties of the coding SNP, and the common GO process with the targets. These properties are precisely the ones that an expert biologist would look for in a manual scan of the region; but Lirnet automatically learned the significance of these features and their relative importance, allowing it to correctly identify the correct polymorphism using a fully automated approach. Moreover, we note that the region contains a number of other plausible candidate genes, including transcription factors and a number of mitochondrial genes; nevertheless, Lirnet identified Mtk1 as top ranked.

To test the effect of the loss of Mkt1 function on the PTR and Puf3 module genes, we deleted the MKT1 open reading frame in the RM background and measured genome-wide RNA expression by DNA microarray analysis. Consistent with Lirnet's predictions, we observed a modest but highly significant down-regulation (KS p-value<10⁻²³) of the Puf3 module in the RM mkt1Δ strain (Figure 8C). While previous work demonstrated Mkt1's role in repressing the translation of an Mpt5/Puf5-dependent transcript [13], our results suggest that Mkt1 plays a role in the RNA stability of Puf3-dependent transcripts. In the PTR module, 16 of 40 genes were among the top 5% most up-regulated genes (hypergeometric p<10⁻¹⁰; Figure 8A), a set which includes the Puf3 module regulators DHH1, KEM1, GCN1, and GCN20, with expression changes (1.36 fold increase) similar to the difference between RM and BY. A similar response of the Puf3 and PTR module genes was observed in an RM strain harboring the BY allele of MKT1 (Figure S8), further supporting the role of MKT1 in the RNA expression differences seen in the original population. Taken together, these results provide strong evidence that MKT1 contains a causative variation for these modules, and further demonstrate Lirnet's ability to identify the correct causal regulator even in a large linked region

Thus, the Lirnet procedure automatically uncovered a comprehensive 3 tiered regulatory cascade in which MKT1 regulates P-body abundance, that consequently regulate Puf3 target abundance, providing significant detail and insight into the mechanism through which the Puf3 module is regulated (Figure S9). Other methods [8],[11] recently applied to these data produced no hypothesis regarding this pathway. The analysis of Brem and Kruglyak [3] linked these genes and many others to a region on Chromosome XIV (Figure S10), but no causal mutation was identified. Geronemo picked up Dhh1 as a key regulator but failed to identify the causal SNP or gene involved (Figure S11). The recent analysis of Zhu et al. [9], performed in parallel to our work, also identified the same genetic region, and provided experimental support for its causal role with respect to a large group of genes. However, their identification of MKT1 as the quantitative trait gene was manual, and based on biological intuition, rather than an automated method. Moreover, their method failed to elucidate the mechanism through which MKT1 regulates the gene expression, missing the role of both Puf3 and the P-bodies. Thus Lirnet is unique in its ability to automatically generate a comprehensive regulatory pathway from causal SNP via the intermediate regulator through which this SNP acts upon the linked genes (Figure S9).

Dataset and Regulators

We applied our analysis to the expression and the genotype data generated from 112 meiotic recombinant progeny of two yeast strains: BY4716 (BY; a laboratory strain) and RM11-1a (RM; a natural isolate) [3]. Our expression data and genotype data were selected as previously described [7]. Our 305 candidate expression regulators were selected using the process previously described [7], and are listed in our accompanying website (http://dags.stanford.edu/lirnet/).

We applied Lirnet to the eQTL dataset [4],[56] of human HapMap individuals −60 European (CEU) and 60 African (YRI) individuals. Among 47,297 probes in the expression data [4], we picked 7,324 whose standard deviation is greater than 0.03 and used them for our analysis. The phase II HapMap data [56] contain genotypes for nearly 4 million SNPs. To perform our experiment in a setting that is closer to that of a real association study, we selected only the SNPs that are on a commercial genotyping chip, namely Affymetrix GeneChip 100 k & 500 k, and used only their genotypes in our analysis.

Identifying Single Nucleotide Polymorphisms (SNPs) between BY and RM

We first identified orthologous genes between BY and RM. We downloaded the genome sequences of S288C (isogenic to BY) and RM from the Saccharomyces Genome Database (http://www.yeastgenome.org/) and Broad Institute of Fungal Genome Initiative (http://www.broad.mit.edu), respectively (sequences were retrieved on 12 January 2005). In order to define orthologous genes between BY and RM, we used reciprocal best BLAST hit [57] (protein sequences of S288C were downloaded from SGD on 12 January 2005). Out of 6,683 genes in 16 nuclear chromosomes, 6,292 (94.1493%) have reciprocal best matches between the two strains. We also retrieved the genomic sequences, 500 bp upstream/downstream of each orthologous pair. We aligned the DNA sequences of the ortholog pair by using LAGAN [58], and retrieved SNPs between the orthologs.

Regulatory Features for SNPs

We constructed a set of regulatory features that describes each single nucleotide variation (SNP) in terms of various intrinsic characteristics. We identified orthologs between BY and RM, and constructed a list of SNPs, as described above. For human regulatory features, we downloaded data from dbSNP containing a list of human SNPs and their various properties. For each SNP, we defined six kinds of features that can determine its regulatory potential. First, we characterized each SNP in terms of its location relative to genes, resulting in seven regulatory features (1 & 12–17 in Table S1). Each non-synonymous coding SNP can change various properties of the corresponding amino acid (AA), which can affect the regulatory role of its gene. Therefore, we described each non-synonymous coding SNP in 10 ways in terms of changes in various properties caused by the corresponding AA change based on various data sources [59],[60] (2–11 in Table S1). A sequence polymorphism on the genomic site that is strongly conserved is more likely to affect the regulatory network. Thus, we characterized each SNP in terms of the conservation score on its genomic site (18 in Table S1). The conservation score was computed based on comparison of protein sequences across a large group of species. For yeast data, we downloaded the aligned sequences from Wapinski et al. [61]. For human data, we downloaded the conservation scores from the UCSC human genome browser (“Most Conserved” track). We also incorporated the feature indicating whether the SNP is likely to regulate the expression of the gene in which it resides (19 in Table S1).

Because regulatory potential of a SNP is likely to be affected by the function of the gene where it resides, we defined a set of regulatory features that indicate whether the gene belongs to each of 87 Gene Ontology (GO) categories related to regulatory roles (20 in Table S1). For human data, we used 48 GO Slim biological process categories. Finally, a SNP might have different regulatory potential over different modules. We defined three ‘pair-wise features’ that describe how likely a SNP is to regulate a particular module (21–23 in Table S1). For each module, we picked GO categories – biological process and molecular function – that are significantly enriched in the module genes; and transcription factors (TFs) whose putative targets appear significantly in the module, based on the ChIP-chip binding data [24]. More precisely, we picked the GO categories and TFs whose hypergeometric p-value is below 10⁻³ after a false discovery rate (FDR) correction. Then, for a combination of a SNP and a module, we constructed three features based on whether: (1) the gene containing the SNP belongs to the module's GO process; (2) the gene containing the SNP belongs to the module's GO function; and (3) the SNP resides in the module's TF. In all cases, we took -log(p-value) to be the value of the regulatory feature, so that a regulator-module pair where the enrichment is highly significantly will have a higher-valued regulatory feature.

Overall, for each SNP n, this process results in a set of 22 regulatory features and 87 (for yeast)/ 48 (for human) features based on the gene function, listed in Table S1.

Regulatory Potentials for Regions (G-Regulators)

Based on the regulatory features of each individual SNP, we modeled the regulatory potential of each genetic marker, representing how likely sequence variations on the marker's chromosomal region regulate expression levels of genes. We also defined a regulatory potential for each e-regulator, representing how likely the regulator's expression is to regulate other genes' expression.

These potentials are based on the regulatory potential of individual SNPs. We model the probability that each SNP n causes expression variation (regulatory potential of n) as:

(1)

where β_k is the parameter called regulatory prior that determines the impact of each regulatory feature on the regulatory potential: higher values of β_k encode the fact that the presence of the feature f_nk increases the probability of having a regulatory effect. The learning algorithm of Lirnet automatically estimates the value of the β parameters from data. In our analysis, we focus only on regulatory features that are likely to increase the regulatory potential, and hence restrict β_k to be non-negative; this assumption can easily be relaxed in the context of other feature sets. A SNP with many important regulatory features (with high β's) will have a higher regulatory prior, but the sigmoid function introduces a saturation effect, preventing the regulatory potential from increasing unboundedly and swamping the data.

Due to linkage disequilibrium, each marker i can represent genotypes of the chromosomal region where it resides. We therefore define the regulatory potential of each genetic marker as an aggregate of the regulatory potentials of the individual SNPs in the corresponding chromosomal region. We assigned each SNP to the region associated with its nearest genotyped marker (or tag SNP). Then, for each region r, we aggregated the contributions of all SNPs (in the region), each modeled based on (1), by summing them up and taking a sigmoid function:

(2)

Therefore, a region that contains a number of SNPs with high regulatory potentials is likely to have a high regulatory potential, but the outer-most sigmoid function again prevents it from increasing unboundedly. We note that a region that contains a large number of moderately relevant SNPs can also achieve a high regulatory potential. This method of aggregation tends to prefer regions with more SNPs, which is arguably justified, as they are also more likely to contain a causal polymorphism. However, other methods of aggregation are also plausible. We experimented with several other approaches; the one selected achieved the highest performance in prediction of expression profiles in test data not used for training the model.

Regulatory Potentials for Expression Regulators

We also model the regulatory potential of candidate expression regulators based on their regulatory features. We used the regulatory features of SNPs (Table S1) for constructing those of an expression regulator. The regulatory features consist of five categories: (1) 7 features each representing the number of SNPs in the gene region having one of the features 1 & 12–17 in Table S1; (2) 1 feature representing the conservation score of the gene region (analogous to 18 in Table S1); (3) 1 binary feature indicating whether the gene is cis-regulated (analogous to 19 in Table S1); (4) 87 (for yeast)/ 48 (for human) binary features indicating whether the gene belongs to each of the GO categories listed in Table S12 (analogous to 20 in Table S1); and (5) three pairwise binary features indicating whether the gene belongs to a GO process category enriched in the module, whether the gene belongs to a GO function category enriched in the module and whether the gene is the TF whose putative binding targets are enriched in the module (analogous to 21–23 in Table S1).

We define the regulatory potential of r to be the probability that each candidate e-regulator r causes expression variation, which we model as follows:

(3)

where g_rk represents the k'th regulatory feature of e-regulator r (explained above) and α_k is the weight assigned to the k'th regulatory feature.

Learning Regulatory Programs using the Lirnet Algorithm

Lirnet attempts to reconstruct regulatory programs that define the regulatory interactions between each group of co-regulated genes (called a module) and its regulatory factors (regulators). Candidate regulators of a module consist of binary genotype values of genetic markers and expression levels of genes that are not in the module. We model the expression level of each gene g in a module m (denoted by y_m,g) as a linear combination of the potential regulators (denoted by x₁,…,x_n):

(4)

where ε represents a zero mean Gaussian noise, and x and y are standardized.

Our objective is to estimate the weights (w_m,1,…,w_m,n) for each module m, from the data that best reflect the regulatory relationship between x's and y. More precisely, given x and y, we aim to construct the network by maximizing the joint log-likelihood Log P( w , y | x ) = Log P( y | x , w )+Log P( w ), where for each module and its regulators P( y | x , w )∼Ν(Σ_iw_ix_i,σ²) and P( w ) represents the prior probability distribution of w. We model the prior probabilities on the weights based on the regulator's regulatory potential: For a regulator r, which can be either a region or an e-regulator, the prior probability is modeled as:

(5)

The regulatory potentials, Pr(Regulator r is causal), are defined in (2) and (3) as functions of β and α, respectively. C₀ and C₁ represent the maximum and minimum regularization parameters C_i, respectively. The prior on the weight w_r is an L₁ prior, which tends to move the weights of less relevant coefficients to 0 [16]; the larger C_r, the stronger the bias towards 0. As the regulatory potential increases, we have that C_r decreases, reducing the tendency of the learning algorithm to set the regulator's coefficient to 0.

To avoid the singularity problem (when the number of parameters is greater than the number of data instances) and a degeneracy problem that occurs when some of the x's are correlated, we introduced an additional L₂ regularization term (also called a ridge term) with a regularization parameter D. The Lirnet algorithm estimates α, β, and w by solving the following optimization problem:

(6)

where C_r's are defined in (5), θ = {α}∪{β}.

Unlike previous approaches that use tree regression [7],[62] or stepwise linear regression [48], this approach deals well with correlated regulators. In particular, when several regulators are highly correlated, both tree regression and stepwise linear regression will select one representative within the set, often arbitrarily; with that regulator selected, the correlated regulators have little explanatory power, and will not be added to the regulator set. This approach is susceptible to making arbitrary decisions based on random fluctuations in the data.

Lirnet uses an iterative coordinate descent algorithm to minimize the above objective function (Figure 1), iterating over two steps, where in one step we optimize over w's given the current α and β, and in the other step we optimize over α and β given the current w's. To learn θ, we solved the optimization problem of Eq. (6) with the current weights w's. To estimate w's, we modified least angle regression (LARS) [63] to allow it to handle the L₂ regularization term (the third term in (6)); LARS is known to be one of the most efficient algorithms for solving this type of regularized regression problem. The regularization parameters C₀, C₁, D and E were determined through a 10 fold cross validation procedure: The arrays were divided into 10 groups; in each run, we train on 9/10 of the arrays, and compute predictive accuracy on the held out 1/10; the parameters were selected to maximize the average performance over the 10 runs. The final set of regulators was determined by using the chosen D and E over the entire set of arrays.

Lirnet can be used both to construct a regulatory program for a pre-determined set of modules, or, as a step in an iterative procedure whereby modules are adapted dynamically to match the learned regulatory programs. In this iterative process, first developed in the module networks algorithm [14], one starts with an initial assignment of genes to modules, for which regulatory programs are subsequently learned. Each gene is then reassigned to the module whose current regulatory program best explains its expression profile. We ran Lirnet both on a set of Geronemo modules, and using this iterative process initialized from these modules. The PGV analysis (see below) showed essentially no difference between these two runs (data not shown). To facilitate comparison to the previous results, we therefore used the Lirnet program for the Geronemo modules.

Experiments on the Human Data

We applied Lirnet to the human eQTL dataset: genotype data from Phase II HapMap Project −60 European (CEU) and 60 African (YRI) individuals [56] – and expression profiles from the same individuals [4]. We treated SNPs on the Affymetrix GeneChip Human Mapping 100 k/ 500 k Array sets as genetic markers (i.e. tagging SNPs), and assumed we observed the genotype of only those SNPs. The regulatory features were constructed for 6,515,224 SNPs downloaded from NCBI dbSNP database (http://www.ncbi.nlm.nih.gov/SNP), based on various data sources such as dbSNP, UCSC genome browser (http://www.genome.ucsc.edu/), and gene ontology [64]. The list of regulatory features for the human data can be found in Table S3. We divided the SNPs into regions corresponding to each of the tag SNP, assigning each gene to its closest tag SNP, and defined the regulatory potential of each individual SNP and each region, according to Eq (1) & (2) in Methods, similarly to the experiments on the yeast data. The learned regulatory prior is listed in Figure 2 and Table S3.

Implementation of a Classical Method for Comparison

As a baseline method for comparison, we used the standard single-marker linkage model (as in [1]–[3]). For each (gene, marker) pair, we performed a linear regression using the gene's expression level as a response variable and the marker as a predictor, and chose the marker that achieves the best fit. For the comparison on the yeast data, we used the published results [1]. We compared those results to those of our implementation of single-marker linkage model (explained above), and the results were almost identical (data not shown). For the human data, we used our implementation.

Calculation of Percentage of Genetic Variance (PGV)

We estimated the percentage of genetic variance (PGV) explained by the identified genetic regulators, following the procedure of Brem & Kruglyak, as also used for Geronemo [7]. In brief, we randomly divided the data of 112 segregants into a detection set (training data) and an estimation set (test data). We used the method on the detection set to learn the regulation programs and (where relevant) the modules, and used the estimation set to calculate the PGV for these regulation programs. The PGV formula uses a corrected ANOVA, which automatically accounts for model complexity determined by the number of predictors. We repeated this process 10 times with different random splits of data, and estimated PGV of each gene by taking the average of its PGV over 10 runs.

We note that, in the protocol of Brem & Kruglyak, the set of regulators is chosen on the detection set, but the actual parameters are estimated using ANOVA on the estimation set. Thus, there is a risk that more complex regulatory program will be able to overfit the training data, producing misleadingly good results. Although the number of parameters in our model is not larger than the number in the Geronemo model, we wanted to demonstrate directly that overfitting is not a factor in these results. We therefore also used an alternative PGV protocol, where the entire regulatory program – both the choice of regulators and the parameters – are derived from the detection set alone, and then the resulting model is estimated on the test set. In the results from this protocol (Figure S1), the proposed models also considerably outperformed Geronemo.

Biological Evaluation of Learned Regulatory Programs

We constructed a set of putative regulator-target pairs for the biological evaluation of the methods. We used five kinds of datasets: (1) deletion and over-expression microarrays [21],[22]; (2) chromatin immunce-precipitation (ChIP-chip) binding experiments [23]; (3) mRNA binding pull-down experiments [31]; (4) transcription factor binding sites [65]; and (5) a literature-curated set of signaling interactions from the Proteome database (http://www.proteome.com). For (1), we considered as targets the genes whose expression changes are within the top 10% of all genes in terms of the magnitude of the expression change. For (2), for each transcription factor, we picked the genes with a p-value of p<0.01. For (3), (4) and (5), we downloaded the lists of putative targets suggested by the corresponding papers.

Biological Evaluation of Predicted Causative SNPs in a Large Chromosomal Region

For each of 13 regions that are identified to contain many candidate regulators of expression variation [1], we constructed a list of genes that have experimental supports of their regulatory role on the targets linked to that region, based on microarray data from deletion experiments [1], [9], [21]–[23] and ChIP-chip binding data [23] (see Figure 4B legend). For each region, we checked for each candidate regulator whether the candidate is “supported” by any of these data: (1) For the deletion microarray dataset, we considered the candidate regulator to be supported if there is a significant overlap (p<0.01; hypergeometric distribution) between its putative target based on the deletion data (within top 20% of differentially expressed genes) and the targets defined by the eQTL data [1]; (2) For the ChIP-chip binding data, we considered a regulator to be supported if there is a significant overlap between the putative targets with a binding significance of p<0.01 and the linked targets.

Ranking Genes in a Chromosomal Region

To identify the causal SNP in a chromosomal region chosen by Lirnet as a g-regulator, we ranked each SNP using its regulatory potential, computed from its regulatory features and the learned coefficients (as in Figure 2). We then ranked the genes according to the SNP of highest regulatory potential in the gene region (coding region, 500 bp upstream, 100 bp downstream).

Yeast Strains and Media

Unless stated, all S. cerevisiae strains used in this study are isogenic with a GAL2 ⁺ derivative of S288c [66] and are also isogenic to the BY parental strain [3]. Strains used to test the effect of deleting MKT1 in the RM background we constructed by transforming a URA3-marked deletion allele into RM11-1a. The “allele swap” strains (RM mkt1-by), replacing the RM allele of MKT1 (G30, R453) with the BY allele (D30, K453), were constructed by standard methods of PCR and transformation into the RM mkt1Δ::URA3 deletion strain. The presence of the appropriate variations and absence of any secondary mutations in the substituted region was confirmed by DNA sequence analysis. All strains were constructed by standard methods of PCR amplification and yeast transformation; details are available upon request. Microarray expression analysis of puf3Δ and mkt1Δ in the BY background used strains from the homozygous yeast deletion collection [67] (Open Biosystems) with the BY4743 isogenic parental strain as a control. Microarray expression analysis of mkt1Δ and mkt1-by in the RM strain background used the isogenic parental strain RM11-1a as a control. Strains containing Puf3, Dhh1, and Edc3 GFP protein fusions were taken from the collection described by Huh et al. [68]. Strains containing inframe protein fusions to the Red fluorescent protein, tdimer2, were constructed using pKT176 [69] by standard methods of yeast methods of PCR amplification and yeast transformation into either the strains containing the GFP tagged protein or the isogenic wild-type strain BY4741. Unless stated, all strains were grown in YPD medium and harvested in mid-log phase.

Microarray Analysis for the mkt1Δ Experiment

Total yeast RNA was isolated by hot phenol method [70]. For both standard and tiling array analysis (see Text S1), total RNA was converted to cDNA and labeled with Alexa 647 and Alexa 555 (Molecular Probes) using the Atlas PowerScript Fluorescent Labeling Kit (Clontech) and an oligo(dT) primer as described by the manufacturer. Labeled cDNA samples were hybridized to either a stock yeast expression array (Agilent-011445 Yeast Oligo Microarray G4140A) or a custom yeast tiling array (described below) and processed according to manufacturer's instructions (Agilent Technologies). Arrays were scanned using a ScanArray Express HT (Perkin Elmer) at a constant laser power of 90% and various photomultiplier tube gains as described in Dudley et al. (2002). Signal and background intensities were measured using GENEPIX image analysis software (Axon Instruments) and data from multiple intensity scans were combined onto a common scale using the MASLINER linear regression method [71]. The log₂ ratio of intensities of the signal and the background was calculated for each array element, and the standard normalization techniques described in Yang et al. [72] were applied to the log₂ ratio values. We used global normalization and intensity-dependent normalization by using LOWESS (locally weighted scatter-plot smoothing) [73] with parameters relevant to our experimental setting, single slides and single print tips.

Fluorescence Microscopy

Live yeast cells containing GFP and tdimer2 fluorescently tagged proteins were visualized with a Nikon Eclipse TE2000-E inverted microscope under 100× objective with oil. GFP was detected using a FITC filter, and tdimer2 was detected using HCRed1. Images were captured using a Hamamatsu Orca-ER CCD digital camera. Image capture and analysis used Metamorph 6.3R5 and Adobe Photoshop software. P-bodies, observed as bright punctate spots in the cytoplasm of cells containing a fluorescently labeled P-body protein, form in live cells after approximately 10 minutes in water or minimal medium lacking glucose under a microscope coverslip. Under the same conditions after approximately 12 minutes, a Puf3-GFP fusion protein formed similar fluorescent spots (Figure S5), most of which overlapped the P-bodies (Table S10). When present in the same cells, punctate spots of Puf3-GFP fluorescence overlap with the punctate pattern formed by known P-body components (69/75 = 92% of P-body spots are also Puf3 spots), showing localization of Puf3 to P-bodies (Table S10).

Analysis of the Oaf1 Module

Of the ten target genes in the peroxisome module, the proximity and orientation of one (YAL049C) suggested that its co-expression could be the result of cross hybridization to the OAF1 probe in the original microarray data; thus, it was removed from further consideration. To evaluate the dependence of the remaining target genes on Oaf1, we examined a published microarray dataset [29] comparing RNA expression oaf1Δ to a wild-type (BY) strain in the presence of oleate (an inducing condition). This dataset also included an estimate of the likelihood of differential expression [74]. We sorted RNA expression levels by the log10 ratios and filtered for λ values greater than 36.23 to arrive at the top 1% (63) most significantly down regulated genes (Table S11).

Accession Numbers

The revised RM11-1a PUF3 DNA sequence determined by this study will be deposited in GenBank (NCBI) prior to publication. All Microarray datasets will be deposited in the GEO database prior to publication.