Evolutionary rates at codon sites may be used to align sequences and infer protein domain function
© Durand et al. 2010
Received: 15 June 2009
Accepted: 24 March 2010
Published: 24 March 2010
Skip to main content
© Durand et al. 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 (Functional Inference using the Rates of Evolution), 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 (Functional Inference using Rates of Evolution) 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.
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.
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 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.
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.
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.
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.