 Research article
 Open Access
 Published:
FastMG: a simple, fast, and accurate maximum likelihood procedure to estimate amino acid replacement rate matrices from large data sets
BMC Bioinformatics volume 15, Article number: 341 (2014)
Abstract
Background
Amino acid replacement rate matrices are a crucial component of many protein analysis systems such as sequence similarity search, sequence alignment, and phylogenetic inference. Ideally, the rate matrix reflects the mutational behavior of the actual data under study; however, estimating amino acid replacement rate matrices requires large protein alignments and is computationally expensive and complex. As a compromise, suboptimal precalculated generic matrices are typically used for proteinbased phylogeny. Sequence availability has now grown to a point where problemspecific rate matrices can often be calculated if the computational cost can be controlled.
Results
The most time consuming step in estimating rate matrices by maximum likelihood is building maximum likelihood phylogenetic trees from protein alignments. We propose a new procedure, called FastMG, to overcome this obstacle. The key innovation is the alignmentsplitting algorithm that splits alignments with many sequences into nonoverlapping subalignments prior to estimating amino acid replacement rates. Experiments with different large data sets showed that the FastMG procedure was an order of magnitude faster than without splitting. Importantly, there was no apparent loss in matrix quality if an appropriate splitting procedure is used.
Conclusions
FastMG is a simple, fast and accurate procedure to estimate amino acid replacement rate matrices from large data sets. It enables researchers to study the evolutionary relationships for specific groups of proteins or taxa with optimized, dataspecific amino acid replacement rate matrices. The programs, data sets, and the new mammalian mitochondrial protein rate matrix are available at http://fastmg.codeplex.com.
Background
Amino acid replacement rate matrices represent the estimates of instantaneous substitution rates between amino acids. The rates simultaneously capture aspects of DNAlevel mutation, the genetic code and proteinlevel selection strength, which varies based on similarity in chemical and physical properties. For example, we usually observe a high substitution rate between lysine and arginine (both positively charged and polar) and a low substitution rate between lysine and cysteine (neutral and nonpolar). Ideally, the replacement rate matrix parameters are optimized against the data under study, but in practice information content in typical sequence alignments is insufficient to do so. Instead, a small number of generic matrices are made available to researchers.
Amino acid replacement rate matrices are essential for many protein analyses, including estimating pairwise distances between protein sequences, or reconstructing protein phylogenetic trees using maximum likelihood or Bayesian frameworks [1, 2] and references therein. Amino acid replacement rate matrices can also be converted into score matrices for protein sequence alignment. Roles and applications of amino acid replacement rate matrices were summarized by Thorne [3].
A number of methods have been proposed to estimate the matrices from protein alignments since the time of Dayhoff [4]. These methods belong to either counting or maximum likelihood approaches. The counting methods are fast, but they are limited to only pairwise protein alignments and closely related amino acid sequences [4, 5]. The maximum likelihood methods have been designed to fully utilize the information contained in multiple protein alignments and the corresponding phylogenetic trees which must be estimated from the data [6–8]. This assumes that the trees are correct, which can be ameliorated by a Bayesian analysis over a set of plausible trees but this would increase an already large computational burden.
With the rapid rise in whole genome sequencing it is now increasingly common to have access to both large numbers of taxa and long concatenated sequence alignments. This creates the opportunity to estimate dataspecific amino acid replacement matrices but also requires efficient computational methods because estimating amino acid replacement rate matrices from protein alignments by maximum likelihood methods is a complex and timeconsuming process [7, 9, 10] and references therein.
Recently, a fully automated maximum likelihood estimation procedure was proposed and used to estimate matrices from different data sets [8, 10, 11]. It consists of two main steps: building maximum likelihood phylogenetic trees and estimating parameters based on the information contained in multiple protein alignments and the corresponding phylogenetic trees. Building maximum likelihood trees is itself a difficult problem because the number of possible trees increases exponentially with the number of sequences in the alignment [1, 12]. A number of maximum likelihood tree search heuristics have been developed to reduce the computational burden [13–16]; however, building maximum likelihood trees is still the most time consuming step in the estimation process. For example, in this study it took 98% (319 out of 324 hours) when estimating the amino acid replacement rate matrix from 100 large alignments in the HSSP (homologyderived structures of proteins) database [17].
In this paper, we propose a new maximum likelihood estimation procedure, FastMG, to work with large data sets. The key idea is to split large alignments into multiple nonoverlapping subalignments with fewer sequences (each subalignment contains at most k sequences) in order to substantially reduce the computational burden of building maximum likelihood trees. The matrices are then estimated from the joint maximum likelihood analysis of the smaller subalignments instead of from the large original alignments.
A preliminary study showed that the splitting strategy significantly decreased the running time of the estimation procedure [9]. Here we demonstrate that a naïve random splitting method compromises the quality of the results. In contrast, our “treebased splitting method” selects subtrees that retain enough information to estimate accurate amino acid replacement rates. Experiments with different large data sets showed that the FastMG procedure yields similar quality matrices in much less time than the standard estimation procedure.
Results and discussion
We assessed the performance of the FastMG procedure on three large data sets: HSSP data set [17], Pfam data set [18], and our concatenated protein alignment of mitochondrial proteins from 299 mammalian species with 3062 amino acid sites. The FastMG procedure was examined with the random splitting algorithm, the treebased splitting algorithm, and different k values targeting subalignment sizes of 8, 16, 24, 32, and 64 sequences. All matrices were estimated on a personal computer (Intel 2.66 GHz, 8 GB RAM). The PhyML software version 3.0 [14] was used to build maximum likelihood trees from alignments with options: 4 gamma categories model, no invariant sites, no bootstrap, SPR tree improvement, and JTT model. The XRATE software version 0.2 [19] was used to estimate model parameters using information in the alignments and corresponding phylogenetic trees. Let us denote:

FastMG^{R}: The FastMG estimation procedure with the random splitting algorithm.

FastMG^{T}: The FastMG estimation procedure with the treebased splitting algorithm.

$${\mathrm{M}}_{k}^{\mathrm{R}}$$
: Replacement rate matrix estimated from data set M using the FastMG^{R} procedure and threshold k (e.g., ${\mathrm{HSSP}}_{8}^{\mathrm{R}}$ is the matrix estimated by the FastMG^{R} procedure with k = 8).

$${\mathrm{M}}_{k}^{\mathrm{T}}$$
: Replacement rate matrix estimated from data set M using the FastMG^{T} procedure and threshold k (e.g., ${\mathrm{HSSP}}_{8}^{\mathrm{T}}$ is the matrix estimated by the FastMG^{T} procedure with k = 8).

M^{s}: Replacement rate matrix estimated from data set M using the standard maximum likelihood estimation procedure (e.g., HSSP^{S} is the matrix estimated by the HSSP data set using the standard estimation procedure).
We compared ${\mathrm{M}}_{\mathrm{k}}^{\mathrm{R}},\phantom{\rule{0.5em}{0ex}}{\mathrm{M}}_{\mathrm{k}}^{\mathrm{T}},\phantom{\rule{0.25em}{0ex}}\mathrm{and}\phantom{\rule{0.25em}{0ex}}{\mathrm{M}}^{\mathrm{s}}$ matrices in terms of both quality and running time. To avoid model bias due to overfitting, each data set consisted of two alignment sets: the training alignment set and the testing alignment set. The matrices were estimated from alignments in the training set and subsequently used to build maximum likelihood trees for alignments in the testing set. Likelihood scores for test set alignments were used to compare the quality of different matrices as used in other studies [7, 8, 11]. Moreover, we used the KishinoHasegawa test [20] to assess the statistical significance of the difference between two matrices as used in previous studies [8, 11].
HSSP data set
We selected 400 alignments from the HSSP database to assess the performance of estimation procedures. The 100 alignments with the largest number of sequences were used as the training alignment set for estimating matrices while the other 300 alignments were used as the testing alignment set. The 100 training alignments contained between 140 and 285 sequences (the mean and max pairwise distances between sequences were 0.481 and 1.286 respectively), for a total of 18325 sequences. The 300 testing alignments contained between 10 and 100 sequences (the mean and max pairwise distances between sequences were 0.635 and 1.846, respectively), for a total of 12854 sequences.
We examined the correlation between matrices estimated from the standard and the FastMG procedures. Table 1 shows high correlations between the 190 exchangeability coefficients of these matrices. As expected, the correlations increase with the increase of k. The results also show that treebased splitting gives exchangeability parameters that are closer to the standard procedure. The opposite is true for the 20 frequency parameters of the matrices, likely because subalignments created by random splitting tend to represent the residue composition of the entire data set, whereas treebased splitting gives subalignments that reflect the residue composition of individual clades. Another reason is likely because the random splitting algorithm selects taxa that are more distantly related and for which there was more time for the substitution process to reach equilibrium.
The more important impact on tree likelihood values in pairwise comparisons between HSSP^{S} and other matrices estimated by FastMG^{R} and FastMG^{T} procedures are represented in Table 2. The matrices estimated from the FastMG^{R} procedure were not as good as the HSSP^{S} matrix for any k value. In contrast, the FastMG^{T} procedure generated high quality matrices and FastMG^{T} with k ≥ 16 was at least as good as the standard estimation procedure.
Table 3 shows the running time of the standard estimation procedure and the FastMG procedure with different splitting algorithms and k values for the HSSP data set. The estimation time of the FastMG procedure increased linearly with the increase of k (e.g., it took 11.0 hours and 22.9 hours to estimate ${\mathrm{HSSP}}_{8}^{\mathrm{T}}$ and ${\mathrm{HSSP}}_{16}^{\mathrm{T}}$, respectively). This was an order of magnitude faster than the 323.7 hours needed to estimate HSSP^{S}. Interestingly, FastMG^{T} was noticeably and consistently faster than FastMG^{R}. The difference is not as large but suggests faster convergence for treebased splitting.
Pfam data set
We also examined different estimation procedures on alignments from the Pfam database. The training alignment set contained the 100 largest alignments from the Pfam database with in total 7640 sequences and a range of 46 to 202 sequences per alignment (the mean and max pairwise distances between sequences were 1.341 and 62.428, respectively). The testing alignment set consisted of 480 other alignments with in total 5434 sequences and a range of 5 to 41 sequences per alignment (the mean and max pairwise distances between sequences were 1.174 and 63.096, respectively). Note that the Pfam alignments tended to contain fewer sequences than the HSSP alignments.
We observed similar values and trends as for the HSSP data set, with very high correlations between the Pfam^{S} matrix and matrices estimated from the FastMG procedure (see Table 4).
The pairwise comparisons of tree likelihood values between the Pfam^{S} matrix and matrices estimated from the FastMG procedure are represented in Table 5. Again, we observed similar trends as for the HSSP data set. The quality of matrices increased with an increase of k, however, FastMG^{R} matrices were never as good as the Pfam^{S} matrix. In contrast, FastMG^{T} produced matrices that were in a majority of cases slightly better than the standard estimation procedure for all k values. For example, the ${\mathrm{Pfam}}_{16}^{\mathrm{T}}$ matrix was better than the Pfam^{S} matrix on 290 out of 480 testing alignments. Moreover, the KishinoHasegawa test showed that the ${\mathrm{Pfam}}_{16}^{\mathrm{T}}$ matrix was significantly better than the Pfam^{S} matrix on 104 testing alignments. The opposite was true only 41 times.
The running time of the standard estimation procedure and the FastMG procedure with different splitting algorithms and k values for the Pfam data set are represented in Table 6. We observed similar running time patterns as for the HSSP data set, however, the magnitude of the speed gain was less in this case (e.g. ~4 times for ${\mathrm{Pfam}}_{16}^{\mathrm{T}}$ compared to Pfam^{S}) because of the typically lower number of sequences in the Pfam alignments.
Mammalian protein alignment
Our mammalian protein alignment (Mam) consists of the concatenated 12 mitochondrial proteins (all except ND6) of 299 mammalian species and 3602 amino acid sites. The mean and max pairwise distances between sequences were 0.256 and 0.441, respectively. Because it is a single alignment, we used 10fold cross validation to examine the performance of different estimation procedures [21]. In particular, the 3602 amino acid sites were randomly distributed across 10 nonoverlapping partitions P_{1}, …, P_{10} each consisting of 360 or 361 sites. The validation was repeated 10 times. Let V_{ t } (t = 1…10) denote the t^{th} validation (i.e. the part P_{ t } was used as the testing data and the other 9 parts were used as the training data).
Table 7 shows the likelihood comparisons between the Mam^{S} matrix and the matrices estimated from the FastMG procedure at the 10 validations. Again, the matrices estimated from the FastMG^{R} procedure were never as good as the Mam^{S} matrix, while the matrices estimated by FastMG^{T} with k ≥ 16 were similar or slightly better than the Mam^{S} matrix. For example, FastMG^{T} with k = 16 gave slightly better likelihood scores than the standard estimation procedure in all 10 validations and significantly better scores in 4 validations.
We also examined the performance of the original MtMam matrix [22] and matrices estimated from the FastMG procedure. Table 8 shows that the matrices estimated by FastMG^{T} were significantly better than the original MtMam for all 10 validations.
The estimation time of different matrices is presented in Table 9. For this alignment with 299 species the FastMG procedure was an order of magnitude faster than the standard estimation procedure (e.g. FastMG^{T} with k = 16 was about 24 times faster than the standard estimation procedure). It also repeats the trend that the treesplitting method is both faster and yields more accurate matrices.
Trends and special considerations
At the start of our studies we anticipated that alignment splitting would result in only minimal deterioration of matrix quality. Instead we observed, on balance, a small improvement in the performance of FastMG^{T} matrices compared to those obtained without splitting. Although the effect is small it is consistent and suggests that a systematic effect is at play. Although further studies are needed to better understand these effects, it is plausible that the specific elimination of deeper and often less well defined branches contributes to the effect. In addition, treebased splitting results in alignments of more closely related sequences with lower risk of sequence alignment errors.
There were a number of alignments where the matrices inferred from the standard and FastMG^{T} procedures resulted in different tree topologies. This occurred in about 5075% of cases for HSSP alignments and 1025% of cases for the Pfam alignment (Tables 2 and 5, columns 4 and 5). This different rate of occurrences is expected because the HSSP alignments have more sequences and are therefore more prone to topology differences. We considered the possibility that achieving a significantly better likelihood score was due to finding a different tree topology. However, in cases where the likelihood score is significantly better only a minority of alignments has a different tree topology.
Our three test cases all include a considerable amount of sequence divergence and all show that treesplitting is superior to random splitting. A preliminary study on closely related influenza sequences showed that treesplitting is still optimal [9]. However, in extreme cases the number of observed substitutions upon treesplitting may become too small to be informative. In such a case random splitting may be preferred.
Another consideration is that our current FastMG procedure uses neighborjoining to create the tree needed for the treebased splitting algorithm. This scales well for alignments with hundreds or even thousands of sequences, but becomes inefficient for extremely large alignments. The performance of the FastMG procedure for such huge alignments will likely benefit from faster alignment splitting methods [23] and faster tree building methods (e.g. FastTree [24]) and this deserves further study if such cases become more common.
Conclusions
Amino acid replacement matrices are essential for many statistical methods to analyze protein sequences. Maximum likelihood methods typically generate the best replacement matrices because they can fully utilize information in the multiple protein alignments. However, for this application maximum likelihood analysis is complex and computationally expensive. Here we propose a modified maximum likelihood procedure to estimate amino acid replacement rate matrices from large data sets. The key component of the modified estimation procedure is the splitting algorithm that divides large alignments into multiple subalignments with fewer sequences that are subsequently used to estimate matrices. The extensive experiments showed that the FastMG^{T} procedure performed well with large data sets and reduces the running time as function of the number of sequences from approximately quadratic to linear, as we expect based on the time complexity of tree inference (see Methods section). FastMG^{T} with k ≥ 16 was about an order of magnitude faster than the standard estimation procedure while it did not reduce the quality of estimated matrices.
Experiments strongly suggest that k = 16 is a good choice for the FastMG^{T} estimation procedure. Even the analysis of the 100 largest alignments of the HSSP database took less than a day with FastMG^{T} and k = 16 on a typical desktop computer. Thus, our method now enables researchers to avoid generic precalculated matrices and instead estimate optimal amino acid replacement matrices for their particular needs from large data sets on their personal computers.
Methods
Model
As usual, the amino acid substitution process is assumed to be independent among amino acid sites. Although the data might not be compositionally homogeneous across the sequences in the alignment (e.g., different amino acid compositions among the clades), we assume that the standard model for the amino acid substitution process over the tree is a Markov process with timehomogeneous, timecontinuous, and timereversible properties [6–8] and references therein. The model is represented by a 20 × 20 instantaneous substitution rate matrix Q = {q_{ xy }} where q_{ xy }(x ≠ y) represents the number of substitutions from amino acid x to amino acid y per time unit. The diagonal elements q_{ xx } are assigned such that the row sums are all zero. Since the model is time reversible, the matrix Q can be decomposed into a symmetric exchangeability rate matrix R = {r_{ xy }} and an amino acid equilibrium frequency vector Π = {π_{ x }} such that q_{ xy } = r_{ xy }π_{ y } and ${q}_{xx}={\displaystyle \sum _{x\ne y}}{q}_{xy}\text{.}$
Model estimation procedure
Given a set of c protein alignments D = {D_{1}, …, D_{ c }}, the substitution model Q can be estimated by maximizing the likelihood L(D) using equation 1 as follows
where L(T_{ i }, ρ_{ i }, Q; D_{ i }) is the likelihood of protein alignment D_{ i } given phylogenetic tree T_{ i }; the rate variation model ρ_{ i }; and substitution model Q[7, 8].
Optimizing L(D) is a difficult problem because we have to simultaneously optimize parameters of T, Q, and ρ. Previous studies indicated that a good approximation of Q (Q’) can be obtained with nearoptimal trees (T’) and rate variation model (ρ’) [7, 8, 10] and references therein. Because trees and rate variation models are computed a priori, the likelihood L(D) can thus be approximated by equation 2 as follows
A better model Q can be estimated from alignments of D using an iterative approach as detailed in the 4step standard estimation procedure as follows [8]:
Standard estimation procedure

Step 0: Input a set of multiple alignments D and a starting matrix Q (typically input only exchangeability rate matrix R, the frequency vector Π is estimated from the data D).

Step 1: Build phylogenetic tree T_{ i } and rate across site model ρ_{ i } for each alignment D_{ i } using maximum likelihood tree construction programs such as PhyML [14].

Step 2: Estimate a new exchangeability matrix Q’ using the approach described by Le and Gascuel [8] and the XRate software [19].

Step 3: Compare Q’ and Q, if they are nearly identical, return Q’ as the optimal model. Otherwise, assign Q ← Q’ and goto Step 1.
Previous studies have showed that this estimation procedure usually stops after three iterations. The procedure is called “standard maximum likelihood estimation procedure”.
Building maximum likelihood trees in Step 1 is a computationally expensive problem [12]. Although a number of heuristics have been proposed for searching maximum likelihood trees [13–16], Step 1 is still the bottleneck of the estimation process. It is roughly estimated [25] that the computing time of PhyML (and of other similar maximum likelihood approaches) is in O(n^{2}s), where n is the number of taxa and s the number of sites. Experiments with a number of different data sets have confirmed this approximation. Thus, it is expected that splitting the original alignments into two equallysized subalignments divides the total computing time by a factor two. This property explains why the computing time to build trees displayed in Tables 3, 6 and 9 is approximately linear as a function of k (size of subalignments).
Alignmentsplitting algorithms
Consider a multiple alignment D_{ i } of n sequences (d_{ i }^{1}, …, d_{ i }^{n}), splitting alignment D_{ i } is a process to divide alignment D_{ i } into nonoverlapping smaller subalignments D_{ i }^{1}, …, D_{ i }^{m} such that each sequence d_{ i }^{j=1…n} ∊ D_{ i } belongs to exactly one subalignment D_{ i }^{t}(t = 1…m). For example, alignment D_{ i } of 8 sequences (d_{ i }^{1}, …, d_{ i }^{8}) can be split into three subalignments D_{ i }^{1} = (d_{ i }^{1}, d_{ i }^{2}, d_{ i }^{8}), D_{ i }^{2} = (d_{ i }^{3}, d_{ i }^{4}, d_{ i }^{5}), D_{ i }^{3} = (d_{ i }^{6}, d_{ i }^{7}). Let k be the maximum number of sequences in a subalignment; we seek a method of alignment splitting that allows us to minimize k, to maximize computational efficiency, while retaining adequate amounts of information in subalignments for estimating the amino acid replacement rates.
Random splitting algorithm
Given a multiple alignment D_{ i } of n sequences (d_{ i }^{1}, …, d_{ i }^{n}) and a threshold k, we describe a naïve splitting algorithm, called “Random splitting algorithm”. The general idea of the algorithm is to randomly split sequences (d_{ i }^{1}, …, d_{ i }^{n}) into subalignments such that each subalignment contains at most k sequences. To prevent creating too small subalignments that might not contain enough information for estimating the replacement rates, the algorithm also requires that each subalignment needs to contain at least k/2 sequences. Thus, each subalignment contains from k/2 to k sequences. The random splitting algorithm is fully described in Algorithm 1.
The random splitting algorithm is very simple; however, its main drawback is that sequences of the same subalignment might be very distantly related. This could compromise estimation of amino acid replacement rates.
Treebased splitting algorithm
We designed a new splitting algorithm, called treebased splitting algorithm, to maximize homology between sequences in each subalignment. Let T_{ i } denote the phylogenetic tree relating sequences (d_{ i }^{1}, …, d_{ i }^{n}) of D_{ i }. The key idea of the treebased splitting algorithm is that sequences in the same subtree of T_{ i } should be split into the same subalignment. Figure 1 shows an example of 9 sequences related by a tree. The sequences can be split into 2 subalignments (s_{1}, s_{2}, s_{3}, s_{8}) and (s_{4}, s_{5}, s_{6}, s_{7}, s_{9}) where (s_{1}, s_{2}, s_{3}, s_{8}) sequences belong to the left subtree while (s_{4}, s_{5}, s_{6}, s_{7}, s_{9}) sequences belong to the right subtree.
The treebased splitting algorithm follows the Neighborjoining algorithm schema [26] to stepbystep group sequences into subalignments. The algorithm also requires that each subalignment contains at least k/2 sequences and at most k sequences. The treebased splitting algorithm is fully described in Algorithm 2.
Note that the distance matrix between sequences used in the Neighborjoining algorithm is estimated by the maximum likelihood method using the LG matrix [8]. Technically, we used BIONJ [27] (an improved version of the Neighborjoining algorithm) to split large alignments.
Fast maximum likelihood estimation procedure (FastMG)
The FastMG procedure consists of two phases: first, the large original alignments are split into nonoverlapping subalignments by one of the alignment splitting algorithms; then the matrix is estimated by joint maximum likelihood analysis of the smaller subalignments instead of the large original alignments. The FastMG procedure is described by the following 5steps
Fast maximum likelihood estimation procedure (FastMG)

Step 0: Input a set of multiple alignments D; a starting matrix Q (typically input only exchangeability rate matrix R, the frequency vector Π is estimated from the data D); and a threshold k.

Step 1: For each alignment D_{ i } ∊ D, split D_{ i } into subalignment D_{ i }^{1}….D_{ i }^{m} using either the random splitting algorithm or the treebased splitting algorithm.

Step 2: Build phylogenetic tree T_{ i }^{j}and rate across site model ρ_{ i }^{j} for each subalignment D_{ i }^{j} using maximum likelihood tree construction programs such as PhyML [14].

Step 3: Estimate a new matrix Q’ from subalignments using the approach described by Le and Gascuel [8] and the XRate software [19].

Step 4: Compare Q’ and Q, if they are nearly identical, return Q’ as the optimal model. Otherwise, assign Q ← Q’ and go to Step 2.
References
 1.
Felsenstein J: Inferring Phylogenies. 2004, Sunderland, MA, USA: Sinauer Associates
 2.
Yang Z: Computational Molecular Evolution. 2006, Oxford, UK: Oxford University Press
 3.
Thorne JL: Models of protein sequence evolution and their applications. Curr Opin Genet Dev. 2000, 10: 602605. 10.1016/S0959437X(00)001428.
 4.
Dayhoff M, Schwartz R, Orcutt B: A model of evolutionary change in proteins. Atlas Protein Seq Struct. 1978, 5: 345351.
 5.
Jones DT, Taylor WR, Thornton JM: The rapid generation of mutation data matrices from protein sequences. Comput Appl Biosci CABIOS. 1992, 8: 275282.
 6.
Adachi J, Hasegawa M: Model of amino acid substitution in proteins encoded by mitochondrial DNA. J Mol Evol. 1996, 42: 459468. 10.1007/BF02498640.
 7.
Whelan S, Goldman N: A general empirical model of protein evolution derived from multiple protein families using a maximumlikelihood approach. Mol Biol Evol. 2001, 18: 691699. 10.1093/oxfordjournals.molbev.a003851.
 8.
Le QS, Gascuel O: An improved general amino acid replacement matrix. Mol Biol Evol. 2008, 25: 13071320. 10.1093/molbev/msn067.
 9.
Le DV, Dang CC, Le QS, Le VS: A Fast and Efficient Method for Estimating Amino Acid Substitution Models. Proceedings of The Third International Conference on Knowledge and Systems Engineering. Edited by: Ho TB, McKay RI, Nguyen XH, Bui TD. 2011, New York, NY, USA: IEEE Publishing, 8591.
 10.
Dang CC, Lefort V, Le VS, Le QS, Gascuel O: ReplacementMatrix: a web server for maximumlikelihood estimation of amino acid replacement rate matrices. Bioinformatics. 2011, 27: 27582760. 10.1093/bioinformatics/btr435.
 11.
Dang CC, Le QS, Gascuel O, Le VS: FLU, an amino acid substitution model for influenza proteins. BMC Evol Biol. 2010, 10: 99110. 10.1186/147121481099.
 12.
Chor B, Tuller T: Maximum likelihood of evolutionary trees: hardness and approximation. Bioinformatics. 2005, 21: 97106. 10.1093/bioinformatics/bti1027.
 13.
Guindon S, Gascuel O: A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol. 2003, 52: 696704. 10.1080/10635150390235520.
 14.
Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O: New algorithms and methods to estimate maximumlikelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010, 59: 307321. 10.1093/sysbio/syq010.
 15.
Le VS, von Haeseler A: IQPNNI: moving fast through tree space and stopping in time. Mol Biol Evol. 2004, 21: 15651571. 10.1093/molbev/msh176.
 16.
Stamatakis A, Ludwig T, Meier H: RAxMLIII: a fast program for maximum likelihoodbased inference of large phylogenetic trees. Bioinformatics. 2005, 21: 456463. 10.1093/bioinformatics/bti191.
 17.
Schneider R, de Daruvar A, Sander C: The HSSP database of protein structuresequence alignments. Nucleic Acids Res. 1997, 25: 226230. 10.1093/nar/25.1.226.
 18.
Bateman A, Birney E, Cerruti L, Durbin R, Etwiller L, Eddy SR, GriffithsJones S, Howe KL, Marshall M, Sonnhammer ELL: The Pfam protein families database. Nucleic Acids Res. 2002, 30: 276280. 10.1093/nar/30.1.276.
 19.
Klosterman PS, Uzilov AV, Bendaña YR, Bradley RK, Chao S, Kosiol C, Goldman N, Holmes I: XRate: a fast prototyping, training and annotation tool for phylogrammars. BMC Bioinformatics. 2006, 7: 428453. 10.1186/147121057428.
 20.
Kishino H, Hasegawa M: Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in hominoidea. J Mol Evol. 1989, 29: 170179. 10.1007/BF02100115.
 21.
Kohavi R: A Study of Crossvalidation and Bootstrap for Accuracy Estimation and Model Selection. Proceedings of the 14th International Joint Conferences on Artificial Intelligence. 1995, Montreal: Morgan Kaufmann Publishers Inc, 11371143.
 22.
Yang Z, Nielsen R, Hasegawa M: Models of amino acid substitution and applications to mitochondrial protein evolution. Mol Biol Evol. 1998, 15: 16001611. 10.1093/oxfordjournals.molbev.a025888.
 23.
Blackshields G, Larkin M, Wallace IM, Wilm A, Higgins DG: Fast embedding methods for clustering tens of thousands of sequences. Comput Biol Chem. 2008, 32 (4): 282286. 10.1016/j.compbiolchem.2008.03.005.
 24.
Price MN, Dehal PS, Arkin AP: FastTree 2 – approximately maximumlikelihood trees for large alignments. PLoS ONE. 2010, 5: e949010.1371/journal.pone.0009490.
 25.
Dereeper A, Guignon V, Blanc G, Audic S, Buffet S, Chevenet F, Dufayard JF, Guindon S, Lefort V, Lescot M, Claverie JM, Gascuel O: Phylogeny.fr: robust phylogenetic analysis for the nonspecialist. Nucleic Acids Res. 2008, 36 (suppl 2): W465W469.
 26.
Saitou N, Nei M: The neighborjoining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol. 1987, 4: 406425.
 27.
Gascuel O: BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data. Mol Biol Evol. 1997, 14: 685695. 10.1093/oxfordjournals.molbev.a025808.
Acknowledgements
This work is financially supported by Vietnam National Foundation for Science and Technology (102.012013.04). BH and OG were supported by LABEX NUMEV to work on this study.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
LSV, LSQ and OG discussed ideas. DCC and LSV designed algorithms and experiments. BH created the concatenated mammalian protein alignment. DCC implemented the algorithms, conducted experiments, and wrote the draft manuscript. All authors analyzed experiment results. BH and OG revised the manuscript. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Cite this article
Dang, C.C., Le, V.S., Gascuel, O. et al. FastMG: a simple, fast, and accurate maximum likelihood procedure to estimate amino acid replacement rate matrices from large data sets. BMC Bioinformatics 15, 341 (2014). https://doi.org/10.1186/1471210515341
Received:
Accepted:
Published:
Keywords
 Amino acid replacement rate matrices
 Maximum likelihood methods
 Phylogenetic trees
 Protein alignments
 Large data sets