Methodology article | Open | Published:
Evolutionary rates at codon sites may be used to align sequences and infer protein domain function
BMC Bioinformaticsvolume 11, Article number: 151 (2010)
Sequence alignments form part of many investigations in molecular biology, including the determination of phylogenetic relationships, the prediction of protein structure and function, and the measurement of evolutionary rates. However, to obtain meaningful results, a significant degree of sequence similarity is required to ensure that the alignments are accurate and the inferences correct. Limitations arise when sequence similarity is low, which is particularly problematic when working with fast-evolving genes, evolutionary distant taxa, genomes with nucleotide biases, and cases of convergent evolution.
A novel approach was conceptualized to address the "low sequence similarity" alignment problem. We developed an alignment algorithm termed FIRE (F unctional I nference using the R ates of E volution), which aligns sequences using the evolutionary rate at codon sites, as measured by the dN/dS ratio, rather than nucleotide or amino acid residues. FIRE was used to test the hypotheses that evolutionary rates can be used to align sequences and that the alignments may be used to infer protein domain function. Using a range of test data, we found that aligning domains based on evolutionary rates was possible even when sequence similarity was very low (for example, antibody variable regions). Furthermore, the alignment has the potential to infer protein domain function, indicating that domains with similar functions are subject to similar evolutionary constraints. These data suggest that an evolutionary rate-based approach to sequence analysis (particularly when combined with structural data) may be used to study cases of convergent evolution or when sequences have very low similarity. However, when aligning homologous gene sets with sequence similarity, FIRE did not perform as well as the best traditional alignment algorithms indicating that the conventional approach of aligning residues as opposed to evolutionary rates remains the method of choice in these cases.
FIRE provides proof of concept that it is possible to align sequences and infer domain function by using evolutionary rates rather than residue similarity. This represents a new approach to sequence analysis with a wide range of potential applications in molecular biology.
Investigations in molecular biology frequently require the analysis of sequence alignments and several methods are available for this purpose. Once a correct alignment is obtained, inferences may be made concerning phylogenetic relationships and putative functions . A fundamental problem arises when accurate sequence alignments cannot be obtained due to poor similarity, which may occur with homologous or analogous genes . Homologous genes, comprising orthologs (arising from speciation events) and paralogs (arising from gene duplication events) share common ancestry; however, sequence similarity may be low when they are rapidly evolving, evolutionary distant, or the sequences have significant nucleotide biases. Analogous genes have similar functions, but arise from convergent evolution and the absence of shared ancestry means there is little or no sequence similarity .
To address the limitation of poor sequence similarity in homologous or analogous sequences, a novel alignment strategy was conceptualized and the FIRE (F unctional I nference using R ates of E volution) algorithm developed. This method uses the evolutionary rate at codon sites, rather than individual residues, to align sequences. Evolutionary pressures are inferred from the parameter ω (ratio of non-synonymous (dN) to synonymous (dS) substitutions, corrected for opportunity) , which is typically used to investigate Darwinian selection at the molecular level. A non-synonymous rate significantly greater that the synonymous rate, ω (dN/dS) > 1, reflects positive selection, while neutral and purifying selection are inferred when ω = 1, and ω < 1, respectively. The evolutionary rate may vary across whole coding sequences, at individual codons within a sequence or along branches within a phylogenetic tree and numerous evolutionary models and software statistical packages for performing the analyses are available. For a recent overview of the subject see . The method reported here makes use of the evolutionary rate at codon sites to align sequences and demonstrates the potential to infer protein domain function in sequences that are subject to similar evolutionary constraints.
Results and Discussion
The aim of this study was to address the limitation of poor similarity when performing sequence alignments. The traditional approach of using the positional homology of residues to align sequences was therefore abandoned and the parameter ω employed instead. The question we asked is: can the selective pressures acting at codon sites across coding sequences, and not residue positional homology, be used to perform alignments? To investigate this question, we aligned homologous domains (orthologous and paralogous data sets), which typically have similar functions, using ω values at codon sites across the sequences. Next, if sequences with similar functions can be aligned using evolutionary rates, we tested the hypothesis that this approach may be used to infer protein domain function in the absence of significant sequence similarity. Domains with poor sequence similarity but similar function (such as the antibody data sets) were employed for this purpose.
The FIRE algorithm was developed in order to perform a pairwise alignment using ω MLEs (maximum likelihood estimates). The ω MLEs at codon sites were obtained from multiple sequence alignments (MSAs) of closely related orthologous sequences (see Methods below for details) and FIRE is therefore, in essence, aligning two MSAs or clades. FIRE was modified from the Needleman-Wunsch algorithm  and finds the pairwise alignment using the codon alignment (based on ω MLEs) to maximize the similarity metric. A codon score, cs, measures the similarity between two aligned ω values in the range [0,1]. The maximum difference between two ω values is capped to ω max and is parameterized - we chose 1.5 as a default, since it is biologically more meaningful to identify sites under positive selection than to emphasize the absolute values of sites with ω>1. Thus, cs(ω1, ω2) = 0, if |ω1-ω2|>ω max and cs(ω1, ω2) = 1-(|ω1-ω2|)/ω max otherwise. The FIRE score is the sum of the cs scores over all aligned codon pairs, normalized for sequence length by dividing the FIRE score by the number of codons in the longer sequence. The opening and extension gap penalties are parameterized and the defaults of 0.5 and 0.05, respectively, were used for the analyses in this study. The FIRE algorithm produces a normalized score, percentage similarity plot, histogram listing the number of codons in the alignment with similar scores per decile, and an alignment of the amino acid sequences. The FIRE software and a User Information file providing further details are freely available at http://dept.ee.wits.ac.za/~scott/fire and are attached as additional files 1 and 2.
The Bayes Empirical Bayes (model M2) or Naïve Empirical Bayes (model M3) posterior mean ω MLEs at codon sites were obtained for 15 data sets using the PAML (Phylogenetic Analysis using Maximum Likelihood) v4.0 software  and a FIRE alignment of each data set with every other set was performed (225 alignments for each model). Data sets included the following domains: (i) a highly conserved transcription factor MYB1 DNA-binding domain (DBD) ; (ii) MYB2, a paralog of MYB1; (iii) a conserved tumor suppressor p53 DBD ; (iv) a metabolic enzyme glycerol kinase (GK) ; and (v) light chain antibody variable regions . Variations in the following parameters were present across data sets and did not adversely affect PAML or FIRE analyses: domain length (90-504 codons), ω MLE range (0-9), dN range (0.2-55.0), and paralogous sequences (metazoan MYB1 and MYB2). The number of sequences per data set ranged from 4-12. The data sets with low sequence number were used to examine the effect of this on FIRE outputs, and as expected, sets with fewer sequences produced less accurate ω MLEs and decreased FIRE reliability. Sequence divergence (dS value across the tree) varied from 2.1 to 48.1, which is within PAML suggested limits of <50 . An exception was the protozoan MYB1 set (dS = 216.1); however, it is unlikely this led to erroneous results since high dS values falsely elevate ω values and in this set all the ω MLEs were <0.3.
The results of salient examples are discussed. FIRE scores and similarity plots are provided in Table 1 and Figure 1, respectively, and the corresponding FIRE alignments are documented in additional file 3. The results presented were obtained for ω MLEs obtained from PAML analyses under model M3 (NSsites = 3). In general, FIRE results for ω MLEs with model M3 were the same or better than M2, the likely reason being that ω MLEs under M3 are less constrained.
Following normalization for sequence length, homologous domains with similar functions produced FIRE scores >0.60 and FIRE plots with the majority of codon similarities >60%. For example, alignments for two orthologous sets (metazoan MYB1/protozoan MYB1) and two paralogous sets (metazoan MYB1/MYB2) produced scores of 0.93 and 0.94, respectively, and similarity plots nearing 100% over most of the sequence. The metazoan GK/protozoan GK alignment provided a comparison of two orthologous sets with a greater range in ω MLEs, significantly different sequence lengths (>100 codons) and non-contiguous evolutionary conserved codons . The difference in sequence length is responsible for the relatively low score (0.62), which is reflected by the gaps in the plot (Figure 2C). Removing the gaps prior to performing a FIRE alignment increased the score to 0.89.
Conserved domains with dissimilar functions produced poor alignments, for example the metazoan MYB1/GK alignment (FIRE score = 0.09). The MYB1/p53 DBD alignment provided an interesting test case. Both are transcription factors, however, the domains are implicated in very different biological functions and, according to our hypothesis, this difference should result in a poor FIRE alignment. This was indeed the case (FIRE score = 0.45). We did, however, note that FIRE produced false positives when two unrelated highly conserved domains (ω MLEs <0.3 across the entire sequence) of similar lengths were aligned (data not shown). Including other computational methods such as structure determination would be valuable to identify these cases.
The effect of sites under strong positive selection on this approach was observed from alignments that included antibody sets. As a result of positive selection, antibody sequences of the variable region demonstrate poor sequence similarity. Despite this, the FIRE scores (>0.60), plots and alignments between the two κ sets, and the κ and λ sets suggested that these domains are under similar evolutionary pressures, which correlates with their similar functions. Alignments of the κ and λ antibody sets with any of the domains that are functionally unrelated, for example MYB1, GK or p53, produced poor FIRE results.
To compare the performance of FIRE against the more conventional ClustalW , MAFFT  and structure-based T-Coffee  algorithms for the data sets in Table 1, we employed sequence-structure alignments (based on known 3D molecular structures) as an independent standard of truth. FATCAT (F lexible structure A lignmen T by C haining A ligned fragment pairs allowing T wists)  and DALI (D istance matrix Ali gnment program)  algorithms, which differ in their treatments of flexible structures, were used to align data sets 1-6 (negative control data sets cannot be aligned due to different 3D structures) based on structures extracted from the worldwide Protein Data Bank . As a measure of performance, the proportions of correctly aligned residue pairs (including gaps) obtained by FIRE, T-Coffee, ClustalW and MAFFT alignments were determined using FATCAT or DALI as the reference alignment. For example, in data set 1 FATCAT aligned 54 pairs of residues, and using this as the reference alignment, the proportions correctly aligned by FIRE, T-Coffee, ClustalW and MAFFT were 0.87, 1.00, 0.99 and 1.00, respectively. Irrespective of whether FATCAT or DALI was used as the reference alignment, FIRE, T-Coffee, ClustalW and MAFFT produced similar performances. Results with FATCAT as the reference alignment are summarized in Table 2 and data set 6 is provided as an alignment example in Figure 2.
The T-Coffee and MAFFT algorithms performed the same or better than FIRE or ClustalW for all data sets. The T-Coffee performance is unsurprising since (i) the same structure files used by FATCAT and DALI were included in the T-Coffee algorithm, and (ii) it is well known that structure-based alignments or a combination of structural information with other approaches (such as T-Coffee which combines structural data with homology-based methods) produce the most accurate alignments when sequence similarity is low . Although T-Coffee and MAFFT are state-of-the art methods and known to perform well, it is worth noting that in data sets 5 and 6, the performances of the three homology-based algorithms (T-Coffee, MAFFT and ClustalW) may be inflated relative to FIRE due to the presence of short stretches of conserved residues involved in stabilizing the tertiary structures of the antibody light chains. It was also observed that FIRE performed better than ClustalW in these same data sets, demonstrating the value of aligning sequences based on evolutionary rates when sequence similarity is low. One shortcoming of the FIRE algorithm is that the performance may actually decrease when sequence similarity is high. The likely reason is that if the two sequences being aligned share a long stretch of very low ω MLEs, such as occurs in highly conserved domains like the MYB transcription factor in data set 2, there is a risk that a gap (possibly due to an insertion/deletion elsewhere in one of the sequences) is introduced causing a misalignment over the conserved region.
The main findings from these results indicate that it is conceptually and methodologically possible to align functionally similar domains accurately using the evolutionary rate at codon sites. In addition, good alignment scores were only obtained for sequences coding for similar functions, indicating that domains with similar functions are subject to similar evolutionary constraints. This suggests that the FIRE approach may be valuable for inferring domain function in situations such as convergent evolution or when sequences are highly divergent. However, it was noted that for homologous genes where there is some sequence similarity, FIRE was not as accurate as MAFFT or T-Coffee. In these cases, the conventional algorithms remain the method of choice.
The value of the FIRE approach currently lies in its ability to align sequences independent of residue similarity. This is helpful for analyzing sequences with poor similarity, which typically occurs with evolutionary distant genes, convergent evolution and sequences with extreme nucleotide biases . It is known that structure-based methods are also valuable in these circumstances and it is likely, therefore, that a combination of the two approaches will offer the best strategy. Structural information may also be valuable for eliminating false positives and negatives produced by FIRE. Furthermore, all components of the FIRE output: normalized scores, plots, alignments and histograms should be evaluated in their biological context. Further experimentation and subsequent refinements to the FIRE algorithm will lead to improvements in method sensitivity and specificity.
FIRE provides proof of concept that it is possible to align sequences and infer domain function by using evolutionary rates. It complements the arsenal of available computational methods and represents a new approach to sequence analyses with a wide range of potential applications in molecular biology.
MLEs of the ω parameter
Coding sequences were extracted from NCBI http://www.ncbi.nlm.nih.gov and PlasmoDB v5.5 http://www.plasmodb.org databases. MSAs were performed with MAFFT  and phylogenetic trees were constructed with ClustalW2  and PAUP* . MSAs and phylogenetic guide trees were processed with PAML 4 (codeml algorithm, F3 × 4 codon model, Model = 0, NSsites = 2 and 3)  to obtain ω MLEs at codon sites. Each data set (comprising a list of ω MLEs) was aligned with every other data set with the program FIRE, which was specifically developed for this purpose (see "algorithm" section above). A list of accession numbers, sequence details, and PAML sequence and tree files are available from the corresponding author.
FIRE uses the rst output files from a PAML analysis (using either NSsites 2 or 3) to extract the ω MLEs and perform an alignment. Two examples of the PAML rst and mlc raw output files of the protozoan GK and λ light chain variable regions data set are provided in additional file 4. The rst files provide ω MLEs and mlc files provide statistical details regarding the analysis, for example dN, dS and kappa (transition/transversion ratio) values. The FIRE results for all data sets aligned with each other are available from the corresponding author.
T-Coffee, ClustalW and MAFFT MSAs
T-Coffee, ClustalW and MAFFT analyses were performed using online servers at the Swiss Institute of Bioinformatics http://tcoffee.vital-it.ch/cgi-bin/Tcoffee/tcoffee_cgi/index.cgi, the European Bioinformatics Institute http://www.ebi.ac.uk/Tools/ and the MAFFT homepage http://align.bmr.kyushu-u.ac.jp/mafft/software, respectively. The T-Coffee advanced algorithm combined a ClustalW alignment with the available PDB structure files, which were selected based on the corresponding data set sequence files and downloaded from the World Wide Protein Databank http://www.wwpdb.org. The E. coli crystal structure PDB file (ID: 1BO5) was used for a T-Coffee analysis of the glycerol kinase set due to a lack of metazoan and protozoan structural data. Default settings were used for ClustalW2 alignments. For the MAFFT analysis, the E-INS-i algorithm was used.
FATCAT and DALI alignments
Sequence-structure alignments were performed with FATCAT http://fatcat.burnham.org and DALI http://ekhidna.biocenter.helsinki.fi/dali_server online servers. The two PDB structures used for each alignment are representative sequences taken from the two clades being aligned with FIRE. The same structures were included in T-Coffee alignments. Due to a lack of structural data for the GK sequence set (data set 4), the FUGUE threading algorithm  was used to generate a second structure for E. histolytica GK (XM_650121.1) with the E. coli crystal structure (PBD ID: 1BO5) as a template. FATCAT, DALI and FUGUE structure alignments were used as reference alignments to which FIRE, T-Coffee, ClustalW and MAFFT alignments were compared.
Note added in proof
Our results complement a recent publication by S.L. Kosakovsky Pond et al. ("Evolutionary Fingerprinting of Genes", Mol Biol Evol 2010, 27:520-536), which demonstrated that probability distributions of evolutionary rates in coding sequences may be used as identifiers of genes. Furthermore, using rapidly evolving RNA viruses as test data, they found that genes within the same functional group have similar evolutionary fingerprints. The findings presented by Kosakovsky Pond et al. and the data in this manuscript suggest that the molecular signatures left behind by evolution represent a tier of information that is untapped by current sequence analysis methods.
Kuzniar A, van Ham RC, Pongor S, Leunissen JA: The quest for orthologs: finding the corresponding gene across genomes. Trends Genet 2008, 24(11):539–551. 10.1016/j.tig.2008.08.009
Abouheif E: Developmental genetics and homology: a hierachical approach. Trends Ecol Evol 1997, 12: 405–408. 10.1016/S0169-5347(97)01125-7
Li WH, Grauer D: Fundamentals of molecular evolution. Sunderland: Sinauer Press; 1991.
Yang Z, Nielsen R: Estimating synonymous and nonsynonymous substitution rates under realistic evolutionary models. Mol Biol Evol 2000, 17(1):32–43.
Delport W, Scheffler K, Seoighe C: Models of coding sequence evolution. Brief Bioinform 2009, 10(1):97–109. 10.1093/bib/bbn049
Needleman SB, Wunsch CD: A general method applicable to the search for similarities in the amino acid sequence of two proteins. J Mol Biol 1970, 48(3):443–453. 10.1016/0022-2836(70)90057-4
Yang Z: PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol 2007, 24(8):1586–1591. 10.1093/molbev/msm088
Stegmaier P, Kel AE, Wingender E: Systematic DNA-binding domain classification of transcription factors. Genome Inform 2004, 15(2):276–286.
Nedelcu AM, Tan C: Early diversification and complex evolutionary history of the p53 tumor suppressor gene family. Dev Genes Evol 2007, 217(11–12):801–806. 10.1007/s00427-007-0185-9
Martinez Agosto JA, McCabe ER: Conserved family of glycerol kinase loci in Drosophila melanogaster . Mol Genet Metab 2006, 88(4):334–345. 10.1016/j.ymgme.2006.01.002
Pilstrom L: The mysterious immunoglobulin light chain. Dev Comp Immunol 2002, 26(2):207–215. 10.1016/S0145-305X(01)00066-0
Durand PM, Naidoo K, Coetzer TL: Evolutionary patterning: a novel approach to the identification of potential drug target sites in Plasmodium falciparum . PLoS ONE 2008, 3(11):e3685. 10.1371/journal.pone.0003685
Thompson JD, Gibson TJ, Higgins DG: Multiple sequence alignment using ClustalW and ClustalX. Curr Protoc Bioinformatics 2002, Chapter 2(Unit 2):3.
Katoh K, Kuma K, Toh H, Miyata T: MAFFT version 5: improvement in accuracy of multiple sequence alignment. Nucleic Acids Res 2005, 33(2):511–518. 10.1093/nar/gki198
Notredame C, Higgins DG, Heringa J: T-Coffee: A novel method for fast and accurate multiple sequence alignment. J Mol Biol 2000, 302(1):205–217. 10.1006/jmbi.2000.4042
Ye Y, Godzik A: Flexible structure alignment by chaining aligned fragment pairs allowing twists. Bioinformatics 2003, 19(Suppl 2):ii246–255.
Holm L, Kaariainen S, Rosenstrom P, Schenkel A: Searching protein structure databases with DaliLite v.3. Bioinformatics 2008, 24(23):2780–2781. 10.1093/bioinformatics/btn507
Berman H, Henrick K, Nakamura H: Announcing the worldwide Protein Data Bank. Nat Struct Biol 2003, 10(12):980. 10.1038/nsb1203-980
Marchler-Bauer A, Panchenko AR, Ariel N, Bryant SH: Comparison of sequence and structure alignments for protein domains. Proteins 2002, 48(3):439–446. 10.1002/prot.10163
Durand PM, Coetzer TL: Utility of computational methods to identify the apoptosis machinery in unicellular eukaryotes. Bioinform Biol Insights 2008, 2: 101–117.
Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, Valentin F, Wallace IM, Wilm A, Lopez R, et al.: Clustal W and Clustal X version 2.0. Bioinformatics 2007, 23(21):2947–2948. 10.1093/bioinformatics/btm404
Swofford D: PAUP* Phylogenetic analysis using parsimony (*and Other Methods). Sunderland: Sinauer Associates; 2002.
Shi J, Blundell TL, Mizuguchi K: FUGUE: sequence-structure homology recognition using environment-specific substitution tables and structure-dependent gap penalties. J Mol Biol 2001, 310: 243–257. 10.1006/jmbi.2001.4762
We thank two anonymous reviewers for valuable suggestions, which significantly improved the manuscript. We also thank Abdelkrim Rachedi for helpful advice in performing the structure-based alignments. This work was supported by: the University of the Witwatersrand Medical Faculty Research Endowment Fund, the National Health Laboratory Service and the National Bioinformatics Network.
PMD conceptualized the study. PMD and SH performed the experiments. SH scripted the FIRE algorithm. PMD, SH and TLC analyzed the results and wrote the manuscript. All authors read and approved the final manuscript.