Functional associations of proteins in entire genomes by means of exhaustive detection of gene fusions

Recent progress in genome analysis has shown that it is possible to predict protein interactions or, more generally, functional associations of proteins using genome sequences alone [1,2,3]. These powerful methods rely on the observation that pairs of genes encoding proteins of known function (usually interacting or forming a complex) tend to be found in other species as a fused gene encoding a single multifunctional protein [4]. This type of event is known as gene fusion and is a well-known process in molecular evolution [5]. Many of these gene fusion events appear to be selectively advantageous by decreasing the regulational load in the cell for a particular process [1,3,5]. Therefore, the detection of gene fusions in one genome (defined as 'composite' proteins) allows the prediction of functional associations between homologous genes that remain separate in another genome (defined as 'component' proteins).

Although gene fusion events appear to be relatively rare, the accurate detection of a gene fusion event in one genome allows interactions to be predicted between many proteins in other genomes. It is this kind of one-to-many relationship that makes this method unique for discovering possible interactions or functional associations between proteins, even for those of unknown function. Unlike previous methods that rely on gene proximity to predict functional coupling [6], this robust method can also detect distal genes within a genome that may be involved in the same process. Furthermore, we have previously demonstrated [1] the high precision of our algorithm, which with an additional constraint of minimum alignment overlap has now increased to over 86% (see Materials and methods). This family of sequence-based methods is analogous with and complementary to the experimental approaches for the detection of protein interaction [7].

Our algorithm was applied individually for each of the 24 genomes, against the remaining 23 genomes which are used as references (see also Materials and methods). Paralogy in the query genome makes it difficult to determine precisely the actual number of possible associations. As we have previously pointed out, paralogy in the query genome increases uncertainty, while paralogy in the reference genome increases the fidelity of the predictions [1]. It is for this reason that detected component and composite proteins from all genomes are subsequently clustered according to sequence similarity [8]. Each cluster should therefore indicate a distinct family of component or composite proteins. The analysis of the distribution of these gene fusion classes among genomes allows us to investigate the dynamics and distribution of this evolutionary process and to assess the extent of the predictive power of the approach.

The detection of gene fusion events yielded 132,812 component and 66,406 composite proteins in an all-against-all genome comparison, but these values represent multiple occurrences of the same proteins across species. Of these, there are 7,224 component and 2,365 composite unique proteins across the 24 species (a 18- and 28-fold reduction respectively). The multiple detection of these cases within or across genomes signifies that the majority of components and composites are observed more than once and therefore represent genuine cases (as opposed to sequencing artifacts, which are usually isolated cases).

The high precision of the method allows the prediction of 39,730 unique pairwise functional associations of the components with reference to the composite protein set. Eighty-six percent of the 66,406 predicted associations obtained from the total number of composite proteins yield a Z-score value of less than 3 (Figure 1), previously shown to result in virtually no false-positive cases [1]. This increased precision is due to the introduction of an additional constraint that does not permit any overlap between the component proteins (see Materials and methods). All the above results are available on the worldwide web (see Materials and methods). Some of these associations are known, but we estimate that more than half of them are newly detected cases, testable by using techniques from functional genomics [7].

Currently, the only species for which predictions can be extensively validated is the yeast Saccharomyces cerevisiae, given the ongoing work on transcript profiling [9] and two-hybrid technology [10]. For yeast, there are 440 distinct component cases (predicted by all other genomes as reference, excluding some highly paralogous Drosophila melanogaster homologs) involved in 706 predicted interactions, most of which are detected by their homology to composite proteins from Caenorhabditis elegans and D. melanogaster. Two examples of predicted protein pairs that are known to interact are CPA1 (YOR303w) with CPA2 (YJR109c) [11] and MET3 (YJR010w) with MET14 (YKL001c) [12], both derived from C. elegans homologs.

We have attempted to test the validity of our predictions by comparing the set of components to a list of potentially interacting gene products, using results from a large-scale two-hybrid experiment [10]. However, there is only one case shared between the 1,004 proteins involved in 957 putative interactions detected by the two-hybrid system and the complete set of 706 pairs in this analysis: YIL033C (SRA1) and YKL166C (TPK3) matching the C. elegans protein C09G4.2 and D. melanogaster protein CT10911. This very low count of common pairs may be expected by the sampling biases of the two rather independent methodologies, given that each approach can only detect a very small subset of the total number of actual interacting pairs in yeast. Interestingly, based on a simple conditional probability calculation, an estimate for the total number of detectable interactions in the yeast cell may be of the order of 675,000.

Another validation procedure for the S. cerevisiae predictions was obtained by comparing all 706 component pairs against their expression profiles from publicly available gene expression data. We have found that at least 20% of our predictions exhibit very strong correlations across gene expression experiments. For each of the pairs, a profile from 87 experiments involving cell cycle [13], sporulation [14] and diauxic shift [9] was used to determine whether expression data corroborated our predictions for the association of the component proteins (see Materials and methods). The detected pairs of components from fusion analysis clearly exhibit similar patterns of expression for the above mentioned experiments (Figure 2, inset). Despite the noise levels for the gene expression data as a result of the limited number of experimental conditions available, with some random pairs exhibiting significant correlations, there are twice as many predicted associations than random, above the threshold of average correlation value of 0.5. With a higher threshold of 0.55, precision is increased, with four times as many predicted associations over the random background. Above this value, 92 predicted functional associations (20% of 536 available pairs, see Materials and methods) exhibit high correlation across all experiments (Figure 2). Below that threshold, it is very difficult to estimate the precision rate of our predictions, because of the high amount of noise and the rather limited number of publicly available gene expression data sets. This comparison between fusion detection and transcript profiling contrasts with previous approaches [2], where expression data was used as a filtering step for the detection of functional associations, and not as a validation criterion.

We have analyzed the S. cerevisiae predictions and detected many interesting cases, which appear to be hitherto undetected functional associations between yeast proteins. Two of these are discussed in some detail here. First, MXR1 (peptide methionine sulfoxide reductase, involved in anti-oxidative processes) [15] and YCL033C (function unknown) are predicted to be functionally associated by virtue of gene fusion in three other species - Helicobacter pylori (both strains), Haemophilus influenzae and Treponema pallidum. This observation is supported by experimental results [16]. MXR1 is 39% identical to the amino terminus of the H. pylori composite proteins and YCL033C is 38% identical to the carboxyl terminus of these proteins. It appears that YCL033C is a selenoprotein, also homologous to the human SelX protein, which may be involved in scavenging reactive oxygen species [16]. These two proteins may be associated to protect the yeast cell from oxidative damage.

Second, another interesting observation involves yeast proteins MSS4 (phosphatidylinositol 4-phosphate kinase), which is involved in a signaling pathway responsible for the cell-cycle-dependent organization of actin cytoskeleton [17], and CCT3 (cytoplasmic chaperonin subunit gamma) which is involved in microtubule and actin assembly [18]. A central domain of CCT3 is 25% identical to a large domain of C. elegans protein VF11C1L.1 and the carboxy-terminal domain of MSS4 is 29% identical to its carboxyl terminus. Thus, these two proteins are predicted to cooperate in cell-cycle-dependent cytoskeleton organization and assembly.

The distribution of components and composites differs dramatically between species. There are 7,224 component cases, with an average of 350 cases per genome, exhibiting significant variation (Figure 3a, blue bars). The query genome sequences detected 2,365 composite cases, with an average of 115 cases per genome (Figure 3b, blue bars). Interestingly, we have observed some relatively small genomes containing composite proteins, which may yield predictions for components of higher organisms. For instance, there are 71 proteins (forming 30 families) in the C. elegans genome that match a fused protein gene in Mycobacterium tuberculosis. Two such examples are the component pair T06C10.1/C49H3.7 matching composite Rv0957 and the component pair W04C9.1/Y65B4B_12.b matching both composites Rv1272c and Rv1273c. Another clear prediction is the S. cereuisiae component pair YER052c/YJRl39c (encoding HOM3/HOM6 respectively) matching composite MetL (3847.PRO) from Escherichia coli and other species.

This has been a key observation that dictated the all-against-all genome comparison in this analysis. In other words, when species A is used as a query against species B, the resulting set of component and composite proteins is different from that with the reverse comparison, when species B is used as a query against species A. The three principal factors in gene fusion during evolution appear to be paralogy, genome size and phylogenetic distance. For instance, larger genomes have more composite, possibly paralogous, proteins. At the same time, closely related species evidently show similar patterns of gene fusion. The results below address each of these factors in turn and examine their relative contribution to the gene fusion process and their effects on the prediction of functional association of proteins.

For every genome, both sets of component and composite proteins were subsequently clustered [8], to detect the degree of paralogy for these proteins (Figure 3, green bars; see also sequence clustering in Materials and methods). There are 2,534 component families, with an average of 105 families per genome (Figure 3a, green bars) and 1,323 composite families, with an average of 55 families per genome (Figure 3b, green bars). Comparing these numbers with the number of unique cases, it is evident that there is a paralogy level of two- to threefold per genome for the composite and component proteins, respectively. As mentioned above, this effect contributes to the confidence of the predictions, depending on whether paralogy is observed in the query or the reference genome.

Another characteristic of this process is the redundancy of both sets of component and composite cases: the number of instances of these may be high but they are widely present across species, falling into well-defined protein families. When all components and composites are clustered as a single set (as opposed to within species, above), sequence clustering results in 1,287 single component families and 621 single composite families (as represented in the current analysis for the 24 species). Comparing these numbers with the number of families per species, it is apparent that there is a further twofold reduction for both sets. This result indicates that gene fusion is widespread in evolution but forms a finite set. Different species may contain a common core of composite families, but also provide new families that are used to predict functional association. For instance, D. melanogaster provides far more composite families (more than 200) compared to C. elegans (fewer than 100) (Figure 3b, green bars). Genomes with unique composite families, such as D. melanogaster, contribute strongly to the majority of predicted interactions. It may also be that only certain classes of proteins are involved in gene fusion and that there is an upper limit for the predictive power of this approach obtainable from (currently available) 621 families.

Evidently, the number of component and composite proteins detected in each species is also dependent on genome size (Figure 4). When the above numbers for unique cases and families of components (Figure 4a) and composites (Figure 4b) are normalized by the number of open reading frames (ORFs) for the species examined, the patterns of distribution are significantly altered. For instance, Aquifex aeolicus and Thermotoga maritima appear to have a large number of components involved in gene fusion (more than 12% of their genes are involved in this process) (Figure 4a), whereas the absolute numbers are low (Figure 3a). This is also the case for composites, where, for example, S. cerevisiae yields as many cases as D. melanogaster in relative terms (4% of the genome) (Figure 4b), while the absolute counts are dramatically different (Figure 3b).

Finally, when the factors of paralogy and genome size are removed by sequence clustering and normalization, respectively, the effect of phylogenetic distance between species can be detected. A distance measure based on shared composite families has been devised (see Materials and methods) and was used to identify relationships between the 24 species examined. The fact that the tree based on this distance measure (Figure 5) does not significantly contradict other trees based on sequence alignments is a strong indication that our hypotheses about the factors involved in gene fusion are valid. This result also indicates that certain types of fusion events appear to be confined to specific phylogenetic groups, such as the Archaea, various bacterial clades and the Eukarya (Figure 5).

The exhaustive detection of gene fusion events in entire genome sequences allows the prediction of functionally associated components based merely on genome structure. The all-against-all species comparison is a necessary step because we have repeatedly observed fused, composite proteins in taxonomically lower organisms. The landscape of gene fusions appears to be a complex one, affected by paralogy, genome size and phylogenetic distance.

Although gene fusion is widely present across various phylogenetic groups, it is a process that may involve only certain types of proteins. Yet, this approach for the prediction of functional associations of proteins results in robust predictions for physical interactions, pathway involvement, complex formation and other types of functional associations of protein molecules.

With the present analysis, we delineate the available universe of fusion events and detect a set of 621 composite protein families from which predictions may be obtained. This approach results in 39,730 pairs of functionally associated proteins across 24 species, with high precision and coverage. This novel set of predictions is made available to the scientific community for the first time, and we believe that many of these cases can be subsequently verified by experimental methods.