A method for predicting quantitative epistatic interactions

To address the problem of inferring epistatic interactions between untested gene pairs, we propose a constrained nearest neighbor-based approach, which exploits the similarity between interaction profiles. Our goal can be most succinctly stated using terms from set theory: given two E-MAPs containing overlapping sets of genes (A and B), we wish to predict scores for interactions between those genes present in A but not in B (A \ B) and those genes present in B but not in A (B \ A) (Figure 1). We achieve this by identifying the nearest neighbors from the genes present in both datasets (A ∩ B). Predictions between genes are only made when one of the genes has a neighbor that is similar enough, i.e. above a specified similarity threshold (using Pearson's correlation as a similarity metric).

Our method is distinct from existing nearest neighbor imputation methods in three ways: 1) we do not attempt to impute missing or erroneous values, rather we impute values that were never measured in the original screen; 2) we infer the novel interactions by combining pairs of independent but overlapping datasets; and 3) we only provide scores for those interactions which we can estimate accurately.

Epistasis scores can be accurately predicted by combining datasets

In order to assess the effectiveness of our method, we first performed a 'leave-one-out' style validation procedure using the Chromosome, RNA and Signalling datasets. The values for interactions between genes in (A \ B) and genes in (A ∩ B) are removed one at a time, and an attempt is made to predict them. This is repeated for interactions between genes in (B \ A) and genes in (A ∩ B), allowing us to assess the effect of altering the similarity threshold on the accuracy of the resulting predictions. We used two measures of accuracy: the Pearson's correlation between the predicted and actual interactions, and the normalized root mean squared error (NRMSE, see Eqn. 1). An improvement in accuracy should result in a higher correlation, and a lower NRMSE value.

Both measures show a similar trend (Figure 2). Beyond a minimum similarity threshold of ~0.4, there appears to be an almost linear relationship between the threshold used and the measured accuracy. Scatter plots constructed to compare predicted and experimentally observed epistasis scores show that our predictions have similar variance to independent E-MAP experiments at a correlation threshold of 0.6. Importantly, at a threshold of 0.6 the number of gene pairs misclassified into incorrect epistasis categories (positive classed as negative etc.) is very low (~1%). Thus, for our analysis we used a threshold of 0.6, preferring a smaller number of more accurate predictions to a significantly larger number of less accurate predictions (~35,000 vs ~160,000). However, all predictions made with a threshold of 0.4 and above are given in additional files 1, 2 and 3.

Predicted epistatic interactions overlap with known interactions and pathways

It has been widely observed that epistatically interacting gene pairs are more likely to share annotated biological properties than randomly selected gene pairs [24]. For instance, gene pairs that show strong epistatic interactions are likely to be involved in the same biological pathways, and so are likely to share Gene Ontology [35] annotations, and to display similar phenotypes. If our predicted epistatic interactions are accurate, then we would expect that they would be similarly enriched for shared annotations and phenotypes. They are also more likely to have been previously identified in genetic interaction screening experiments than randomly selected genes. Furthermore, enrichment for synthetic lethal interactions is likely to be stronger for negative than positive interactions, because synthetic lethal screens primarily report on negative growth phenotypes.

We therefore sought to validate our predictions by comparison with a number of additional data sets. We selected the strongly epistatic pairs from our set of predictions using the thresholds identified in [33] (positive, S-score > 2.0; and negative, S-score < -2.5), and asked whether they are enriched in a variety of annotated properties associated with genetically interacting gene pairs obtained from the literature including the presence of synthetic lethal (genetic) interactions, positive genetic interactions, genes sharing experimental phenotypes, and genes sharing database annotations (Gene Ontology).

As expected, pairs of genes predicted to interact both positively and negatively tend to share GO Process and SGD Phenotype annotations (Table 3), confirming that these genes generally operate within similar biological processes in the cell. Furthermore, the predicted negatively interacting pairs were at least 20-fold more likely to have previously been identified as Synthetic Sick/Lethal in published experiments than random pairs, while pairs predicted to interact positively were at least 13 times as likely to have previously been labeled as such (Table 3). These observations hold for all three pairs of E-MAPs considered, and so are likely to be widely applicable. Overall, validation using both internal (leave-one-out analysis) and external (comparison with annotated biological features) measures supports the assertion that our method successfully generates reliable predictions of epistatic relationships of both positive and negative polarity. The ability of our method to accurately predict the polarity of an predicted epistatic interaction is important, because distinguishing between positive and negative epistasis is critical to mapping the high level relationships between biochemical processes and protein complexes in the cell.

Improved mapping of epistatic relationships among complexes using inferred interactions

A prime motivation for identifying epistatically interacting genes is to improve our understanding of how the various activities of the cell are coordinated. Kelley and Ideker [21] combined synthetic lethal interactions with physical protein-protein interaction, protein-DNA interaction, and pathway information from the KEGG database [36] to construct 'between pathway' and 'within pathway' models for genetic interactions. This work found that 'between pathway' interactions were the more common class, perhaps because negative interactions typically reflect redundant behavior that is more likely to be found between parallel biochemical pathways than within a single pathway focusing on a given biochemical function. Segre and coworkers [23], used a flux balance model to predict both positive and negative epistatic interactions between gene pairs, and showed that sets of genes involved in coordinated activities tend to show aligned ('monochromatic') epistatic polarity. Bandyopadhyay and coworkers [20] extended these approaches by combining E-MAP data with protein-protein interaction data, to identify modules defined by physical interactions and the largely monochromatic interactions between them.

In order to show how our predicted epistatic interactions can supplement and extend these overviews of the cell, we created a combined E-MAP, consisting of the published RNA, Chromosome and Signalling E-MAPs, augmented with new predictions arising from this manuscript. We mapped the resulting combined E-MAP onto a recently produced high quality list of yeast protein complexes [37], using the method of Bandyopadhyay and coworkers [20, 34]. We identified pairs of complexes bridged by genetic interactions which were significantly more negative or positive than one would expect by chance (P < 0.001). Additionally we identified complexes whose internal interactions were similarly 'monochromatic'.

Our inferred interactions have two main uses in this context. First, they provide additional evidence for previously proposed connections among protein complexes, and second, they establish new connections. By comparing the resulting network of linked complexes before and after the addition of predicted interactions, we can see which links are a direct result of our predicted interactions. In total 105 'inter complex' links were significantly 'monochromatic' after the addition of our predictions, in other words a set of previously unknown inter-complex links identified with the help of inference. In contrast, the statistical significance of only one 'intra complex' link increased after including our predictions. This apparent discrepancy chiefly arises due to the composition of the E-MAPs published to date, where complexes tend to be represented in only a single E-MAP. Interactions between complexes therefore frequently correspond to links between E-MAPs (Table 5 and Figure 3).

Dataset	Intra Complex Links	Inter Complex Links
Combined without predictions	20	674
Combined with predictions	21	761

Counts of the significant monochromatic epistatic interactions that occur within and between physical complexes defined by Pu et al. Inter complex links are genetic interactions that bridge two physical complexes, while intra complex links are genetic interactions within a single complex. Interactions are counted here if they are significantly more negative or positive than one would expect at random (P < 0.001).

The largest connected component of these 105 novel 'inter-complex' links. Is shown in Figure 4A. This figure is effectively an overview of how cellular processes are organized into modular arrangements, while operating at different hierarchical levels. At the level of the protein complex, individual proteins physically interact, while at the pathway or function level, protein complexes communicate with each other through interactions that are reflected at the epistatic level.

For instance, the cytoplasmic small ribosomal subunit is central to a cluster of protein complexes linking protein translation (guanyl nucleotide exchange factor; prefoldin complex) to DNA biology(e.g. replication factor C, RecQ helicase-Topoisomerase III) and complexes mediating gene expression/chromatin biology (e.g. NuA3 HAT complex). This is consistent with a coordinating role for the ribosome in regulating different aspects of these processes during the cell cycle [38]. Interestingly, the DNA biology gene clusters interact negatively with the ribosome (suggestive of a supportive or co-operative function), while the prefoldin and guanyl nucleotide exchange factor genes interact positively (suggestive of a modifying or regulating relationship).

Insights into chromatin modifying machines

A good example of how prediction consolidates previously proposed links between complexes, while highlighting new ones, is the SWR-C complex (Figure 4B). This 13 subunit complex is responsible for deposition of the H2A histone variant Htz1 into chromatin in order to promote gene expression and inhibit silencing by heterochromatin [39, 40]. Unsurprisingly, SWR-C has been functionally linked to other chromatin modifying complexes including the NuA4 histone acetyltransfersase and the Ino80-C chromatin remodeling complexes [41]. In fact, several proteins are shared among these complexes, including Rvb1p, Arp4p and Yaf9p. Connections between SWR-C and NuA4 or Ino80-C components have previously been observed in E-MAP experiments, and also between SWR-C and another chromatin modifying complex, COMPASS. COMPASS houses a histone H3K4 methyltransferase activity that contributes to gene silencing near telomeres, and is linked to SWR-1 activity via three previously observed E-MAP pairs: SDC1 with SWR1, SWC5 and VPS71. Our predicted interactions however, link the COMPASS subunit SDC1 to eight additional SWR-C subunits: five via negative epistasis (SWC3, VPS72, ARP4, ARP6, YAF9) and three (RVB1, RVB2, ACT1) via weakly positive epistatic interactions (Figure 4 B). We tested the SDC1-SWC5 combination in the small-scale E-MAP (Table 4) and found a very strong negative epistasis (S-score = -15) that agreed with our imputation from the Chromosome Biology/RNA Processing pair of E-MAPs (S-score = -7).

We also predicted an interaction between components of the SWR-C complex (VPS72, SWC5, ARP6, YAF9, SWC3) and the nuclear pore complex (NPC) components Nup170 and Nup188, Figure 4(B). We experimentally tested and validated four of these newly predicted epistatic interactions between these two complexes in the small-scale E-MAP where the predicted S-score agreed with the measured interaction (all S < -3). Interactions between these two components and SWR-C had not been observed before E-MAP screening, but are consistent with observations in the RNA processing E-MAP that connected SWR-C components with the NPC components Nup120p and Mex67p, and a recent report linking Htz1p deposition and the tethering of genes of the nuclear pore [42].

We also experimentally confirmed negative interactions between SWC5 and VPS72 and the nuclear importin KAP122, Figure 4(B), whose gene product binds to the NPC components Nup1p and Nup2p (predicted and observed S-scores = -10/-8 for both interactions). Thus our predicted epistatic interactions suggest new links between SWR-C and components of both the NPC and soluble nuclear transport factors, consistent with a biological rationale whereby histone exchange within nucleosomes is likely to be coordinated with several aspects of the biology of the nucleus [43].

New connections were also made between SWR-C and protein complexes involved in different steps of gene expression. The snRNP Brr1p is involved in snRNP biogenesis pre-mRNA splicing and the cognate gene was predicted to strongly negatively interact with ARP6, VPS72, SWC3, SWC5, and YAF9 Figure 4(B). Two of these interactions were tested and confirmed in the small-scale E-MAP (predicted and observed S-score for SWC5-BRR1 was -8/-9, and for VPS72-BRR1 is -8/-6), while further connections had already been reported in the RNA E-MAP(SWR1-BRR1 and VPS71-BRR1). Similarly, five members of SWR-C were predicted to interact with LEA1, a gene encoding a U2 snRNP component involved in telomere maintenance [44]. Two interactions were tested and confirmed (S-score for VPS72-LEA1 predicated/observed was -7/-7; for SWC5-LEA1 was -7/-8), Figure 4(B). Four members of SWR-C, VPS72, APR6, SWC3 and SWC5, were predicted to negatively interact with the APT cleavage and polyadenylation factor subcomplex component SYC1, with both VPS72 and SWC5 tested and confirmed (predicted and observed S-scores -3/-2 in both cases), Figure 4(B). Thus our predicted interactions link SWR-C to several aspects of gene expression. Interestingly, mRNA splicing and poly(A) cleavage, previously believed to be independent steps of gene expression, are now considered to be linked through large protein complexes that mediate surveillance mechanisms [45].

Evidence for links between the RPD3L complex and ribosome maturation

Another notable connection is that established between proteins associated with the Box C/D small nucleolar RNAs (snoRNAs) and the Rpd3L histone deacetylase. Genes encoding the Box C/D snRNP associated proteins Nop1 and Sik1 were previously linked to Rpd3L components Pho23, Sin3, and Rtx2, but this was extended in our study to four additional Rpd3L-C genes, Cti6, Dep1, Sds3, Rpd3 (Figure 4C). It is currently unclear how these two complexes functionally interact. The Box C/D snoRNP is responsible for 2'-O-methylation of pre-RNA during ribosome maturation [46], while Rpd3L-C is involved in the regulation of a wide variety of yeast genes [47–49]. Perhaps the extensive genetic interactions between genes encoding the complex subunits, all positive, reflect some communication between ribosome biogenesis and global gene expression changes during growth or during the cell cycle. A number of replication defects that map to the snoRNP component Sik1 (also known as Nop56) have been reported [50], while human snoRNP associated proteins have been co-purified with proteins involved in DNA replication and transcription [51]. Very recently, Rpd3L-C proteins were implicated in replication timing events in yeast [52], so a plausible explanation for epistasis between these complexes could be based on the coordination of DNA replication or the regulation of gene expression or both.

Materials

The three E-MAP datasets analyzed in this paper can be obtained from the supplementary materials of their corresponding papers [32–34].

Genetic interactions can be defined as the divergence(ε) in the observed fitness of strains with two disrupted genes(w_ab) from the expected fitness. The expected fitness is calculated using the fitness of strains with individual gene disruptions, typically using a multiplicative model(w_aw_b).

E-MAPs model this divergence using the S-score, a modified t-score, which takes into account both the magnitude of the divergence, and the variance of the measurements [19]. It thus represents the magnitude of the observed effect and the confidence that it is the result of a true genetic interaction. Each E-MAP consists of a matrix of these S-scores, indicating the type and strength of interaction between each pair of genes under consideration. Negative scores indicate negative epistasis, i.e. the yeast grew more slowly than expected, while positive scores indicate positive epistasis, i.e. more rapid growth was observed. Scores close to zero indicate the probable absence of an interaction between two genes - i.e. they function in independent pathways in the cell.

The GO Slim mapping at the Saccharomyces Genome Database (SGD)[55] was used as the source of gene ontology annotations. These are high-level terms, so annotations which contained more than 1000 genes were filtered out. Phenotype data was also taken from the Saccharomyces Genome Database. Phenotypes associated with more than 175 genes were filtered out, resulting in the removal of terms such as 'inviable', 'viable', and 'haploinsufficient'. Both annotation sets were downloaded on 1^st February 2010.

Additionally we investigated whether our predicted interactions were more likely than random to have been previously identified as genetically interacting. For this we used annotations from the Biogrid, version 2.0.61[56]. Our 'positive' test set was comprised of gene pairs annotated with 'Positive Genetic' or 'Synthetic Rescue', while our 'negative' test set was comprised of gene pairs annotated with 'Synthetic Lethality' or 'Synthetic Growth Defect'.

Creating the combined E-MAP

Our combined E-MAP contained all interactions present in the three individual E-MAPs and our predicted interactions. It was created as follows: In cases where the interaction was present in more than one E-MAP, an average of all available interactions was used. In cases where both a measured interaction and a predicted interaction were present, then the measured interaction was used. In cases where there were multiple predictions for a single interaction, the results were averaged.

Identifying Inter and Intra Complex links

Given the large number of predictions generated, a method to aid their visualization and interpretation was necessary. We exploited the modular architecture of the cell, and focused on interactions at the level of the protein complex. Our set of known protein complexes was taken from a recent manually curated set developed by Pu et al[37] as a more up to date alternative to the widely used MIPS dataset. This set consists of 408 complexes backed by evidence from small scale experiments. To identify strong monochromatic genetic interactions between complexes, we used a method developed by Bandyopadhyay et al [20]. The median genetic interaction between proteins from two different complexes was compared to the median of 10 ⁶ equal-sized random samples of genetic interactions. Interactions were considered significant at P < 0.001, the same threshold used in [34]. The same method and threshold were used to identify strong genetic interactions within individual complexes.

Small Scale E-MAP Validation

Validation was performed using a small-scale E-MAP, and is evaluated in Table 4. The double mutant strains were grown and scored using the standard protocols described in [19, 18]. We evaluated the performance of our prediction using terms borrowed from the information retrieval community: precision and recall. These are defined as follows:

The small-scale E-MAP was taken to be the 'gold-standard' for identifying strong positive and negative interactions, and precision and recall were evaluated for the overlap between this E-MAP and each of the datasets under consideration in Table 4.