Identification of Significant Domain Genetic Interactions Using Maximum Likelihood Estimation

By assuming the independences among domain genetic interactions, we estimated the probabilities of domain genetic interactions using the Maximum Likelihood estimation (MLE) method. The probability of a domain pair indicates its propensity to genetically interact. The high probability of a domain pair imply that number of SLGIs including the domain interaction is high comparing to the number of protein pairs including the domain interaction that are not genetically interacted. Then, we calculated an evidence score E [19] for each domain genetic interaction to estimate its relative significance. The E scores are not the standard likelihood ratio test as they are calculated using only the positive data [19]. The E scores indicate the relative significances of the domain pairs in the SLGIs. The probability and the E score give different view about the significances of domain genetic interactions. The high probabilities do not always correlate to high E score. For example, there are a large number of protein pairs including a single domain pair. If very few of them are SLGIs, the probability of the domain pair will be low. However, if that domain pair is the only domain pair in those SLGIs, it will have high E score. On the other hand, in case there are a large number of protein pairs containing the same multiple domain pairs, if most of the protein pairs are SLGIs, the probabilities of those domain pairs will be high. However, the E scores of those domain pairs will be relatively low as removing any of those domain pairs will have limited effect on the probabilities of SLGIs. Therefore, we used both the probability and E score to select significant domain genetic interactions.

Domain Genetic Interactions Rarely Overlap Domain Physical Interactions

To investigate the relationship between domain genetic interactions and the domain physical interactions, we compared the 3848 significant domain genetic interactions with the domain physical interactions from the iPfam database [20]. There are 4030 domain physical interactions of 1867 domains in iPfam database (2008 version). Among them, 1556 domain physical interactions of 1048 domains exist in yeast proteins we studied. There are 422 domains in both the genetic and physical interactions. We applied the Fisher's exact test in R to examine if domains participating genetic interactions are independent from domains participating in physical interactions. We obtained a p-value of 0.00025, which indicated a significant difference between two sets of protein domains. Furthermore, there are only 70 domain pairs overlapped between the 1556 domain physical interactions and 3848 domain genetic interactions. The Fisher's exact test obtained a p-value of 0, which showed that the domain genetic interactions significantly differ from the domain physical interactions. Thus, the domain genetic interaction is a new type of relationship among domains and can provide a complementary view about the relationships between domains.

The Properties of Domain Genetic Interaction Network

To obtain an overview of the domain genetic interactions, we modeled the domain genetic interactions as a network, in which each node represents a domain and each link represents a genetic interaction between two domains. Then, we examined the properties of this domain genetic interaction network. The average connectivity of the network is 7.392. The domain PF00022 (Actin) has the highest connections of 186. The average shortest path among all nodes is 3.31. Comparing the size of nodes, this number is very small. Furthermore, the average node clustering coefficient of the network is 0.159. These properties indicated that the domain genetic interaction network has the small world property [22]. The analysis of the connectivity distribution of this network showed a power-law distribution with an exponent degree of 1.45 (Figure 1), which indicated that the domain genetic interaction network is a scale free network [23]. These results demonstrated that the domain genetic interaction network follows the common principles of biological networks [24].

Domain Genetic Interactions in the SLGIs between Multidomain Proteins

Most proteins are multidomain proteins, which are created as a result of different genetic events, such as insertions and duplications [25, 26]. Multidomain proteins may have different functionalities due to different domains. Our identification of domain genetic interactions helps understand how the domains contribute to functionality defined by the SLGIs, and then help elucidate the functional relationships between proteins at domain level from their genetic interactions.

Figure 2 shows examples of domain genetic interactions in three SLGIs. Figure 2A shows the domain genetic interactions between SGS1 and TOP1. The SGS1 has three domains and TOP1 has two domains. Only one of the three domains of SGS1, PF00570 (HRDC domain), has high probabilities to interact with two domains of TOP1: PF01028 (Eukaryotic DNA topoisomerase I catalytic core domain) and PF02919 (Eukaryotic DNA topoisomerase I DNA binding domain). The other two domains: PF00270 (DEAD/DEAH domain) and PF00271 (Helicase conserved C-terminal domain) of SGS1 have no genetic interaction with two domains of TOP1 (very low probabilities). Figure 2B shows the domain genetic interactions between TOP3 and TOP1. Both domains of TOP3: PF01131 (DNA topoisomerase domain) and PF01751 (Toprim domain) show high probabilities of genetic interaction with two domains of TOP1. Figure 2C shows the domain genetic interactions between RAD5 and RAD50. The RAD5 has four domains and RAD50 has two domains. The genetic interaction between RAD5 and RAD50 is mainly due to a single domain genetic interaction between PF08797 (HIRAN domain) of RAD5 and PF04423 (zinc hook domain) of RAD50. The PF00097 (Zinc finger domain), PF00176 (SNF2 family N-terminal domain) and PF00271 (Helicase conserved C-terminal domain) of RAD5 and PF02463 (RecF/RecN/SMC N terminal domain) of RAD50 have low contribution to the genetic interaction between RAD5 and RAD50. Those examples imply different domain genetic interaction architectures in SLGIs between yeast multidomaon proteins. The domain genetic interactions may exist between all domains of two proteins (TOP3 and TOP1), or between part of domains of one protein and all domains of the other protein (SGS1 and TOP1), or between part of domains of one protein and part of domains of the other protein (RAD5 and RAD50). The domain genetic interactions are able to help understand functional relationships between multidomain proteins at domain level.

We also investigated the domain genetic interactions of SLGIs between SGS1 and other proteins. We found that only the PF00570 (HRDC domain) of SGS1 has significant genetic interactions with other domains. The results implied that certain functionality of SGS1 may be only due to its HRDC domain, rather than its DEAD/DEAH domain and helicase conserved C-terminal domain. Previous study showed that the HRDC domain of SGS1 is required for its cellular functions involving topoisomerases [27]. Thus, the domain genetic interaction analysis can help understand how domains contribute to the different functionalities of multidomain proteins.

Prediction and Validation of Genome-wide SLGIs Using Protein Domains

Having established that there is a strong correlation between domain genetic interactions and SLGIs, we explored to predict the probabilities of protein pairs to be SLGIs using probabilities of domain genetic interactions. We were able to assign 599752 protein pairs with probabilities greater than 0. Table S1 (Additional File 1) lists the number of SLGIs predicted by MLE approach at different probability cutoffs. All predicted SLGIs are hosted on our website for searching and downloading.

We then compared the correlation coefficients of gene expressions of predicted SLGIs to those of known SLGIs and those of all possible protein pairs. We used a yeast cell cycle gene expression data [28], which contains 77 data points. We calculated the T-score and P-value for the null hypothesis that there is no difference between the means of predicted SLGIs and the means of known SLGIs and the null hypothesis that there is no difference between the means of predicted SLGIs and the means of all pairs. The results are shown in Table S2 (Additional File 2). The correlation coefficients of gene expression of predicted SLGIs with threshold greater than 0.3 are similar to those of known SLGIs. The gene expression correlation coefficients of predicted SLGIs using different thresholds significantly differ from those of all pairs except for SLGIs with probabilities greater than 0.85, which have only small numbers of SLGIs. Those results indicated that the correlation coefficients of gene expressions of predicted SLGIs are similar to those of known SLGIs, rather than to those of random pairs. Recently, it was reported that the SLGIs are likely to have similar GO annotations [29]. We studied the distribution of similarities of Gene Ontology (GO) annotations between predicted SLGIs and also compared them to those of known SLGIs and all possible protein pairs. As shown in Table S3-5 (Additional File 3), the mean similarities of GO annotations of predicted SLGIs always significantly differ from those of all pairs. At certain probability thresholds, the mean similarities of GO annotations of predicted SLGIs show no significant differences from those of known SLGIs. Those thresholds are 0.35 for biological process and cellular component and 0.25 for molecular function. As the probability thresholds increase, the mean similarities of GO annotations also increase, which will make them differ from those of known SLGIs. The studies of GO annotations similarities and expression correlation coefficients showed that the predicted SLGIs at probability threshold around 0.3 are similar to experimentally obtained SLGIs.

Novel SLGIs predicted by MLE Approach

The MLE approach was able to predict novel SLGIs. Table 2 lists 17 novel SLGIs (not included in our training data) with probability >0.9. We predicted the MYO4/DYN1 pair to be SLGI with the highest probability of 0.9895. The MYO4 is one of two type V myosin motors. The other one is MYO2, which is known to genetically interact with DYN1 [29]. Among 17 SLGIs, 12 SLGIs are between transcription initiation factor genes and genes from RNA polymerase complex. Previously, many SLGIs between transcription initiation factor genes and RNA polymerase genes have already been reported [30–32]. We expected our novel SLGIs to help further elucidate the transcription machine. We then investigated genes involved in cellular response to stresses caused by DNA damage. We downloaded a list of 116 DNA repair and recombination genes from Kyoto Encyclopedia of Genes and Genomes (KEGG) database [33]. Then, we extracted SLGIs in which at least one protein of SLGI pair are related to DNA repair. Figure 3 shows the DNA repair related SLGIs with probability > 0.7. Of a total of 133 SLGIs, 22 SLGIs are novel. Some of new predicted SLGIs are supported by previous studies. For example, the TOP3 and RAD1 double mutant has shown extreme synergistic growth defects in a previous study [34]. A recent study showed that the RTT109 and YKu70 double mutant exhibits synergistic defects under hydroxyurea treatment [35]. The PAP2 and POL2 were also shown to genetically interact at high temperatures [36].

Protein 1 Name	Protein 1 Annotation	Protein 2 Name	Protein 2 Annotation	Probability
MYO4	type V myosin motors	DYN1	Cytoplasmic heavy chain dynein	0.9895
RPO21	RNA polymerase II largest subunit	BRF1	TFIIIB B-related factor	0.9816
DST1	General transcription elongation factor TFIIS	RET1	Second-largest subunit of RNA polymerase III	0.9797
STO1	Large subunit of the nuclear mRNA cap-binding protein complex	RPA190	RNA polymerase I subunit	0.9586
STO1	Large subunit of the nuclear mRNA cap-binding protein complex	RPO21	RNA polymerase II largest subunit	0.9586
DST1	General transcription elongation factor TFIIS	RPA135	RNA polymerase I subunit A135	0.9570
BNR1	Formin, nucleates the formation of linear actin filaments	DYN1	Cytoplasmic heavy chain dynein	0.9441
CEG1	Alpha (guanylyltransferase) subunit of the mRNA capping enzyme	RPA190	RNA polymerase I subunit	0.9437
CEG1	Alpha (guanylyltransferase) subunit of the mRNA capping enzyme	RPO31	RNA polymerase III subunit	0.9437
RPB5	RNA polymerase subunit	RPA190	RNA polymerase I subunit	0.9437
RPB5	RNA polymerase subunit	RPO31	RNA polymerase III subunit	0.9437
DST1	General transcription elongation factor TFIIS	RPA190	RNA polymerase I subunit	0.9265
DST1	General transcription elongation factor TFIIS	RPO31	RNA polymerase III subunit	0.9265
ARO1	Pentafunctional arom protein	SPT16	Subunit of the heterodimeric FACT complex	0.9198
RET1	Second-largest subunit of RNA polymerase III	SUA7	Transcription factor TFIIB	0.9035
BRF1	TFIIIB B-related factor	RET1	Second-largest subunit of RNA polymerase III	0.9035
BRF1	TFIIIB B-related factor	RPB2	RNA polymerase II second largest subunit	0.9035

Compensatory Pathways from Predicted SLGIs

Protein pathways are a part of gene network in the cell that can accomplish certain functionality. The SLGIs have been proposed to have high probability of occurrence in compensatory pathways [37]. Thus, the SLGIs within the pathways are rare and the SLGIs between pathways are significantly abundant. Identification compensatory pathways from synthetic lethal genetic interactions can be a powerful way to understanding cellular functional relationships. We expected our new predictions will increase the ability of understanding compensatory pathways. We applied the algorithm of Ma [38] to identify compensatory pathways from 7583 predicted SLGIs with probability higher than 0.3. Among 7583 SLGIs, 4497 SLGIs are novel predictions. Although Ma et al. have shown that physical interactions are enriched in discovered pathways [38], there is no assumption that proteins in those pathways are physically interacting. Totally, we obtained 167 pairs of compensatory pathways, which include 638 proteins and 3535 SLGIs. Then, we examined the GO term co-occurrences in each pathway using the SGD GO Term Finder [39]. The GO Term Finder calculated the p-values that reflect the probability of observing the co-occurrence of proteins with a given GO term in a certain pathway by chance based on a binomial distribution. Among 167 compensatory pathway pairs, 153 pairs are found significantly enriched GO terms with p-value less than 0.05. All 167 pairs of compensatory pathways and their most significant GO terms are listed on our http://www.genenetworks.net website.

Figure 4A and 4B list two pairs of compensatory pathways related to DNA double strand break (DSB) repair. Among total 28 SLGIs between those two compensatory pathways, only 12 of them are known. The DSB is a kind of lethal DNA damage in which both strands of double helix are cleaved. The cell maintains multiple mechanisms to repair double strand breaks with the homologous recombination (HR) and non-homologous end joining (NHEJ) as two major mechanisms. Pathways involving those two mechanisms had shown to be compensatory to each other in Drosophila [40]. The pathways on the left of two compensatory pairs involve in NHEJ. The KU complex (YKU80 and YKU70) is the damage detector of NHEJ [41]. Mutation of nucleoporins NUP84 and NUP133 was reported to be hypersensitive to DNA damage [42]. Furthermore, the nucleoporins (NUP84 and NUP133) had been reported to colocalize and coimmunoprecipitate with Slx5/Slx8 [43], which regulate a DNA repair pathway that counteract Rad5-independent HR [44]. Our discovery suggested that the nucleoporins (NUP84 and NUP133) may regulate the NHEJ pathway through Ku complex for double strand breaks repair. The RAD24 is a DNA damage checkpoint protein [45]. Studies also had shown that the YKU80 and RAD24 are in the same NHEJ pathways [45] to repair irradiation and methylmethanesulphonate (MMS) damages. The CSM3 is a DNA replication checkpoint protein. Genetic interaction between RAD24 and CSM3 may indicate the pathway on the right side of Figure 4 A can actually become two pathways: one is YKU80, RAD24, NUP84; the other is YKU80, CSM3 and NUP84.

The pathways on the right of those two compensatory pairs of Figure 4 A and 4B involve in HR. On the right of Figure 4A, the pathway involves three proteins: RAD57, RAD51 and DMC1. The DMC1 and RAD51 are known to form a complex [46] and have roles in recombination [47]. RAD51 and RAD57 are in the same protein family and it has been shown that RAD57/RAD55 bind with RAD51 [48]. On the right of Figure 4B, the pathway involves four proteins: SGS1, RRP6, MRE11 and RAD52. The MRE11 and SGS1 are part of a two-step mechanism to initial HR [49]. The RAD52 plays a major role in the single strand annealing and strand exchange.

Figure 4C also shows two compensatory DNA repair pathways. Among 18 SLGIs between two pathways, 11 are novel predictions. Many previous studies have supported our prediction of compensatory functionalities of those two pathways. Guillet and Boiteus [50] reported that the APN2 and MUS81/MMS4 have overlapping function to repair 3'-blocked single strand breaks (SSBs). Vance and Wilson [51] showed that TDP1and RAD1 function as redundant pathways. SGS1/TOP3 had also been reported to have overlapping functionality with the MMS4/MUS81 [52]. Proteins in these compensatory pathways involve in many DNA-repair pathways, such as base excision repair (BER), nucleotide excision repair (NER) and HR. Studies already showed that BER, NER and HR pathways have overlapping specificities [53]. Our prediction of these compensatory pathways may help to understand the overlapping functionalities among BER, NER and HR pathways.

Another interesting compensatory pathway pair is related to hydroperoxides response in the cell. As shown in Figure 4D, 11 of 15 SLGIs connecting two pathways are novel predictions. The three proteins (TSA1, TSA2 and PRX1) on the right are all thioredoxin peroxidase. They play a role in reduction of hydroperocxides in cytoplasm (TSA1 and TSA2) and mitochondrion (PRX1). Meanwhile, the five proteins (RAD52, RAD5, MDM31, MDM32, and MRE11) in the left pathway are related to DNA repair. RAD52, RAD5 and MRE11 are known related to double strand break repair [54–56]. The MDM31 and MDM32 relate the stability of mitochondrial DNA [57]. Those five proteins may relate to repair the damage created by the hydroperoxides. This compensatory pathway pair implies that removing hydroperoxides or repairing damage is alternative strategies for cells to survive from hydropeoxides.

Source of Data

Estimation of probabilities and significances of domain genetic interactions

We treated the protein genetic interactions P _m,n and domain genetic interaction D _i,j as random variables. The P _m,n = 1 if two proteins i and j genetically interact and P _m,n = 0 otherwise. The D _i,j = 1 if two domain i and j genetically interact and D _i,j = 0 otherwise. We estimated the probabilities of potential domain interactions Pr(D _i,j = 1) by maximizing the likelihood of observed genetic interactions using the Expectation-Maximization (EM) algorithm [17–19]. The EM algorithm iteratively estimates the maximum likelihood of the 'complete data' that combine the observed data and unobserved data. Here, the protein genetic interactions and the domain information of proteins are our observed data and the domain genetic interactions are our unobserved data. By assuming the independences among domain genetic interactions, the likelihood of observed protein genetic interactions based on domain genetic interactions can be obtained as:

(1)

where M_i,j be the number of genetic interacting pairs between domain i and j in all protein genetic interactions; N_i,j be the number of non genetic interacting domain pairs between domain i and j in protein genetic interactions; and K_i,j be the number of non genetic interacting protein pairs including domain i in one protein and j in the other one. The value of K_i,j is counted from all possible protein pairs with domain i in one protein and j in the other one with excluding the known genetic interacting protein pairs. The K_i,j will keep unchanged during EM computation. The a and b are pseudo counts to avoid the Pr(D_i,j = 1) or Pr(D_i,j = 0) to be zero when instances of domains i and j are rare. We set both a and b to 1 in our calculation. Initially, the M_i,j was set to the number of genetic interactions between domain i and j in experimental genetic interactions; N_i,j is set to 0. And Pr(D_i,j = 1) was initialized as following:

(2)

In each Expectation step of EM algorithm, we first estimated the expected values of E[M_i,j] and E[N_i,j] [19] using the current Pr(D_i,j = 1):

(3)

(4)

Then, we calculate the Pr(D_i,j = 1) using the E[M_i,j] and E[N_i,j] as following (Maximization step):

(5)

The EM algorithm was iterated the above Expectation and Maximization steps till the change of likelihood L is less than a pre-defined small value.

The evidence score E_i,j [19] of domain pair i and j is defined as the ratio between the probability that a pair of proteins, m and n, genetically interact given that the pair of domains, i and j, genetically interact and the probability that a pair of proteins, m and n, genetically interact given that the pair of domains, i and j, do not genetically interact:

(6)

The denotes the probability of genetic interacting between domains k and l given that the domains i and j do not genetically interact.

Prediction of SLGIs using probabilities of domain genetic interactions

We assumed that two proteins genetically interact (P _m,n = 1) if and only if at least one domain pair from the two proteins genetically interact (D _i,j = 1). Then, we calculated the probability of two proteins genetically interacting Pr(P _m,n = 1) as following:

(7)

A pair of proteins was predicted to be SLGI only if its probability is higher than a predefined threshold.