Evolutionary rates at codon sites may be used to align sequences and infer protein domain function
© Durand et al; licensee BioMed Central Ltd. 2010
Received: 15 June 2009
Accepted: 24 March 2010
Published: 24 March 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.
Data sets aligned
*metazoan MYB1 and protozoan MYB1
ω ≤ 0.2
*metazoan MYB1 and metazoan MYB2
ω ≤ 0.3
protozoan MYB1 and metazoan MYB2
ω ≤ 0.3
*metazoan GK and protozoan GK
ω ≤ 1.3
€κ light chain VR and κ light chain VR
ω ≤ 7.0
*κ light chain VR and λ light chain VR
ω ≤ 8.2
*metazoan MYB1 and metazoan p53
ω ≤ 1.3
metazoan MYB1 and metazoan GK
ω ≤ 1.3
κ light chain VR and metazoan p53
ω ≤ 7.0
*metazoan p53 and λ light chain VR
ω ≤ 8.2
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.
FIRE, T-Coffee, ClustalW and MAFFT performances
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.
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.
- 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.009View ArticlePubMedGoogle Scholar
- Abouheif E: Developmental genetics and homology: a hierachical approach. Trends Ecol Evol 1997, 12: 405–408. 10.1016/S0169-5347(97)01125-7View ArticlePubMedGoogle Scholar
- Li WH, Grauer D: Fundamentals of molecular evolution. Sunderland: Sinauer Press; 1991.Google Scholar
- Yang Z, Nielsen R: Estimating synonymous and nonsynonymous substitution rates under realistic evolutionary models. Mol Biol Evol 2000, 17(1):32–43.View ArticlePubMedGoogle Scholar
- Delport W, Scheffler K, Seoighe C: Models of coding sequence evolution. Brief Bioinform 2009, 10(1):97–109. 10.1093/bib/bbn049View ArticlePubMedPubMed CentralGoogle Scholar
- 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-4View ArticlePubMedGoogle Scholar
- Yang Z: PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol 2007, 24(8):1586–1591. 10.1093/molbev/msm088View ArticlePubMedGoogle Scholar
- Stegmaier P, Kel AE, Wingender E: Systematic DNA-binding domain classification of transcription factors. Genome Inform 2004, 15(2):276–286.PubMedGoogle Scholar
- 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-9View ArticlePubMedGoogle Scholar
- 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.002View ArticlePubMedPubMed CentralGoogle Scholar
- Pilstrom L: The mysterious immunoglobulin light chain. Dev Comp Immunol 2002, 26(2):207–215. 10.1016/S0145-305X(01)00066-0View ArticlePubMedGoogle Scholar
- 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.0003685View ArticlePubMedPubMed CentralGoogle Scholar
- Thompson JD, Gibson TJ, Higgins DG: Multiple sequence alignment using ClustalW and ClustalX. Curr Protoc Bioinformatics 2002, Chapter 2(Unit 2):3.PubMedGoogle Scholar
- 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/gki198View ArticlePubMedPubMed CentralGoogle Scholar
- 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.4042View ArticlePubMedGoogle Scholar
- Ye Y, Godzik A: Flexible structure alignment by chaining aligned fragment pairs allowing twists. Bioinformatics 2003, 19(Suppl 2):ii246–255.View ArticlePubMedGoogle Scholar
- 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/btn507View ArticlePubMedPubMed CentralGoogle Scholar
- Berman H, Henrick K, Nakamura H: Announcing the worldwide Protein Data Bank. Nat Struct Biol 2003, 10(12):980. 10.1038/nsb1203-980View ArticlePubMedGoogle Scholar
- 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.10163View ArticlePubMedGoogle Scholar
- Durand PM, Coetzer TL: Utility of computational methods to identify the apoptosis machinery in unicellular eukaryotes. Bioinform Biol Insights 2008, 2: 101–117.PubMedPubMed CentralGoogle Scholar
- 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/btm404View ArticlePubMedGoogle Scholar
- Swofford D: PAUP* Phylogenetic analysis using parsimony (*and Other Methods). Sunderland: Sinauer Associates; 2002.Google Scholar
- 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.4762View ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.