BACKGROUND
sequence alignment is a fundamental tool in contemporary biology.
most algorithmic formulations of the problem seek an alignment maximizing a function known as the objective score.
objective scores are usually defined as a sum of terms for matching pairs of letters  and penalties for gaps.
while the effects of different substitution matrices and gap parameters have been extensively studied in the context of local alignment and homolog recognition , their effects on alignment accuracy, especially for global alignment, are not so well understood.
several heuristics are in common use, for example clustalw's choice of low-identity matrices for aligning low-identity sequences  <cit> , which have not to the best of my knowledge been empirically tested.
one factor impeding such testing is the lack of effective automated methods for optimizing parameters for a given objective function.
previous work in this area has included unsupervised expectation maximization  <cit> , discriminative training  <cit> , inverse parametric alignment  <cit>  and maximum-margin structured learning  <cit> .
katoh et al.
<cit>  reported using golden section search to optimize gap penalties, but this generally assumes a unimodal function which does not hold in this case.
exhaustive search has been attempted several times, for examples see  <cit> .
in this work i describe a new parameter optimization procedure, pop, compare it with the ipa method of  <cit> , and use it to investigate a number of questions related to global protein alignment accuracy, including: which is the best substitution matrix, do the best choices of matrix and gap penalties vary with sequence identity, should terminal gaps be penalized differently from internal gaps, and does the loss of precision due to integer rounding in substitution matrices degrade alignment accuracy?
ipa was chosen for comparison because it is a recently published method, code was readily available, and because it optimizes parameters exactly equivalent to those studied in this work.
to the best of my knowledge there has been no prior comparative study of alternative gap optimization methods, so ipa's effectiveness relative to other approaches is not known.
in the tests performed here, pop was shown to achieve from  <dig> % to  <dig> % higher training-set accuracy than ipa, and evidence is presented that pop finds parameters within about  <dig> % of the true optimum, i.e.
an improvement in training set error of up to an order of magnitude.
differences of a fraction of a percent will be shown to be informative in assessing differences between substitution matrices, for which a method with the improved precision of pop is required.
pop, like ipa, is designed to maximize training set accuracy, which may be achieved at some cost in generalization error, i.e.
worse performance on novel data.
in the case of pop this is by design: the problem of minimizing training set error is simpler than minimizing generalization error and the experiments reported here were not designed to compute the most biologically appropriate parameters for a given algorithm but rather to compare models with small and equal numbers of free parameters.
it is then reasonable to assume that training set accuracy is a good measure of the relative performance of those models.
let w = wi, i =  <dig> ..n be the parameters of interest , and q be the function to be optimized .
the goal is to find values wopt = argmax q that maximize q. this optimization problem is challenging for several reasons.
sufficiently small changes in w will leave all alignments unchanged and hence most q functions of interest have zero partial derivatives almost everywhere.
also, q is typically expensive to compute, requiring seconds to hours to evaluate at a single point, and is non-convex with many local maxima .
methods
ipa
kim and kececioglu  <cit>  described an inverse parametric alignment method that can be trained on so-called "partial examples" such as those found in balibase  <cit>  in which only a subset of columns are annotated as reliably aligned.
eagu kim  kindly provided a software implementation of this algorithm .
this implementation extended the method to allow separate terminal and internal gap penalties to be learned.
pop
the central idea in pop is to search for changes in parameters that result in "goldilocks" changes in q--not too big, not too small, but changes in just of the right size to have a good chance of indicating coarse trends rather than small-scale noise.
this approach can be regarded as a generalization of line search optimization .
as theoretical and empirical methods for dealing with non-convex, non-continuous q functions are currently limited, practical experience, especially examination of plots such as fig.
<dig>  is the best guide to determining the size of a "just right" change.
in the experiments reported here, values around  <dig> % were most effective.
pop proceeds in three stages.
in the first stage, an n-dimensional hypercuboid is explored by evaluating q at each point in a regularly spaced grid.
local maxima are identified as points at which q is greater than all neighboring points, and the best of these are used as starting points for the second stage.
the second and third stages use a hill-climbing strategy to approach a local maximum given a starting point.
the second stage uses a faster but less accurate variant of the hill-climbing method than the third stage, which starts from the best local optima found in the second stage.
the final result is the best maximum found by the third stage.
to save computation time, the three stages typically use increasingly large subsets of the training samples, with the first two stages using randomly chosen subsets and the third stage the entire training set.
hill-climbing
hill climbing repeats the following procedure until no improvement is found.
starting from a point w, each axis i is explored in both the positive and negative direction; i.e., all parameters are held fixed but one .
let δ be a proposed change in wi,  = wi + δ, w' = {wj, j ≠ i, }, q' = q, i.e.
q' is the value of q after a proposed move along one axis, and Δ = q' - q. in each direction a δ is found that gives Δ >  <dig> , or, failing that, a δ that reduces q by an amount that is small, but not too small, say in the range  <dig> % to 1%.
very small changes are undesirable because they will tend to be dominated by "noise"  rather than indicating a systematic trend.
large reductions in q will tend to take wi out of the range where it might make an improvement when combined with a change along another axis .
a heuristic designed to minimize these problems introduces a maximum  and minimum  reduction in q .
after computing the proposed changes, a new w is determined as described in making a move below.
proposed changes
a proposed change δ is determined as follows.
in the first iteration δ is set to a small fraction  of the absolute value of wi or a small value  if wi is zero, or a "hint" value set by the user.
in subsequent iterations the initial value of δ is its final value from the previous iteration.
having set the initial value, the following procedure is then repeated.
the value of Δ is computed.
if Δ =  <dig>  δ is multiplied by  <dig>  if Δ <  <dig> and -Δ <λ, i.e.
q is reduced by too small an amount, δ is doubled.
if Δ <  <dig> and -Δ >μ, i.e.
q is reduced by too much, δ is halved.
otherwise, Δ >  <dig>  and further exploration is tried.
first, smaller values of δ  are evaluated in an attempt to discover an intermediate maximum; if one is found this is repeated to investigate increasingly small changes.
if smaller values give no further improvement, larger values are tried  in a similar way.
a maximum improvement in q, denoted Γ, is imposed in order to avoid extreme changes in a single parameter when far from an optimum .
making a move
the search for proposed changes yields a set of 2n values {δi+, δi-, i =  <dig> ... n} where the + and - superscripts indicate moves in the positive and negative directions respectively.
all moves are considered that apply the proposed changes on from k =  <dig>   <dig> ... k axes simultaneously, where k is a user-settable parameter may be  <dig>   <dig> ... n. when k >  <dig>  new evaluations of q are required.
for example, k =  <dig> considers eight moves for each pair of parameters: , , , , , , , and  where + or -indicates an increase or decrease in the parameter and  <dig> indicates no change.
the move giving the best improvement in q is accepted, otherwise the routine terminates.
naively, one might expect that consideration of single-axis moves only  would suffice, but in practice it turns out that allowing moves on two or more axes sometimes gives a significantly better final result, with k =  <dig> sufficient in most cases.
it is desirable to keep k small as the number d of evaluated moves per iteration grows very rapidly with k:
it also turns out that allowing moves along axes that reduce q can, when combined with moves on other axes, improve q and that allowing this possibility also gives significantly improved optimization in some cases.
for example, a small increase in the gap open penalty and a small decrease in the gap extend penalty might each reduce q, but increase q when both changes are made simultaneously.
fast hill-climbing
a "fast" variant of the hill-climbing procedure sets k =  <dig>  and immediately applies any proposed move that is found to improve q. speed is also improved by increasing μ, λ and Γ. these modifications reduce the number of times the q function is invoked, saving execution time but sometimes yielding significantly inferior parameters.
directions to try, and the sign  first tried, may be selected at random or cycled to minimize systematic bias.
gap models and substitution matrices
let a model be the set of parameters associated with a given objective function.
in this work the substitution matrix is regarded as fixed; only gap penalty parameters are included.
a gap is a series of indel symbols; formally, a maximal consecutive sequence of indels.
if a gap includes the first or last column of an alignment it is described as terminal; all other gaps are internal.
while pop has been implemented for a wide variety of models, most of these are works in progress and i will therefore report results for just two models: g <dig>  in which the same affine penalties are applied to all gaps, and g <dig>  in which internal and terminal gaps have different open and extend penalties.
the parameters are g, e, g and e; the penalty for an internal gap of length l is g +  e and g +  e for a terminal gap.
thus g <dig> is a special case of g <dig> in which g = g and e = e. the objective score is then the sum of substitution scores minus gap penalties; a maximum-scoring global alignment is found using standard dynamic programming techniques.
the following substitution matrix types are considered: blosum  <cit> , rblosum  <cit> , pam  <cit> , jtt  <cit>  and vt/vtml  <cit> .
the rblosum matrices were constructed using a bug-fixed version of the program used to compute blosum matrices from the blocks database  <cit> .
surprisingly, the corrected blosum <dig> matrix was found to slightly degrade performance in homolog recognition  <cit> .
each matrix family is a series with members defined by a measure of evolutionary distance: percent identity cutoff in the case of blosum and rblosum, pam distance for the rest.
conventionally, integer valued matrices are used in which log-odds scores in fractional bits have been rounded to one or two significant figures.
presumably this is for historical reasons: in older computer processors integer arithmetic was faster than floating point.
this is no longer the case for many general-purpose processors, and regardless it is of interest to investigate whether the loss of precision due to integer rounding has an effect on alignment accuracy.
unless otherwise stated, full precision, one bit unit matrices were used.
benchmark data
reference data was taken from three protein alignment benchmarks: prefab  <cit>  version  <dig> and balibase versions  <dig> and  <dig>  there are  <dig>  pair-wise reference alignments in prefab v <dig>  extracting all pairs from the multiple alignments in balibase versions  <dig> and  <dig> gives  <dig>  and  <dig>  pairs respectively.
the large number of pairs in version  <dig> motivated the use of the more tractable version  <dig> for all-pairs tests.
subsets  <dig>  and  <dig>  of version  <dig> were used to investigate optimizing terminal gap penalties.
accuracy measure
the accuracy measure  for a single pair is the number of correctly aligned residue pairs divided by the number of residue pairs in the reference alignment.
the total accuracy score  is the weighted average of q over all pairs, where the weighting is uniform in the case of prefab and the inverse of the number of pairs in the original multiple alignment in the case of balibase.
it would have been desirable to estimate error bars using a method such as the bayesian bootstrap  <cit> ; however this proved to be infeasible due to limitations in available computer time .
for brevity, the reference sets will be denoted bali , prefab , termgaps , and intgaps .
to investigate the effects of evolutionary distance three subsets of balibase v <dig> pairs were constructed:  <dig>  randomly selected pairs with identities in the range 0-33% , 33-66%  and 66-99% , respectively.
these were selected from the full-length rather than domain-trimmed sequences.
RESULTS
comparison of pop and ipa
i compared pop and ipa using the termgaps set, where the biggest difference between g <dig> and g <dig> might be expected, and on bali.
the blosum <dig> matrix in 1/ <dig> bit units was used as this was hardcoded into ipa.
results are shown in table  <dig>  pop was found to be from  <dig> % to  <dig> % more accurate than ipa; these improvements are typical .
training set accuracy  for the g <dig> and g <dig> models on bali and termgaps using blosum <dig> 1/3-bit units.
in the numbers in the pop and ipa columns are the q values achieved by the gap para-meters reported by each method.
in all cases, pop achieves higher accuracy, the >acc column shows the amount by which the pop accuracy is better than ipa.
substitution matrix family
fig.
<dig> shows the results of optimizing model g <dig> on the bali and prefab sets.
the results are qualitatively similar on the two sets, giving confidence that they indicate general trends rather than artifacts of benchmark construction, of overtraining or of significantly suboptimal local maxima.
this is further confirmed by cross-training , which again gives similar, albeit noisier, results, as would be expected.
the results show vtml to be the best matrix series, with vtml > vt > blosum > jtt > pam holding for most members, though the differences between vtml, vt and blosum are small except at the extreme high and low-distance ends of the series.
the pam series is significantly inferior to vtml, giving accuracy scores around 2% lower.
matrix selection by evolutionary distance
it has been suggested that different substitution matrices might be more effective at different evolutionary distances.
for example, clustalw uses more distant matrices to align more distant sequences, and lassmann and sonhammer  <cit>  report that "softer", i.e.
more distant, matrices are better at aligning more more distant sequences while are equally good with more closely related sequences.
to investigate this, i optimized g <dig> on the id0_ <dig>  id33_ <dig> and id66_ <dig> sets, with the results shown in fig,  <dig>  interestingly, the plots are qualitatively similar for the three sets despite increasing pair-wise identity and the increasingly narrow variation in optimal q. remarkably, in the case of id66_ <dig> accuracies are all above  <dig> % and the difference between the best matrix  and worst  is only  <dig> %, yet trends observed at lower identities are still clearly discernible.
the peak in each curve that identifies the best matrix in each family is at approximately the same evolutionary distance on each set, showing that the best choice of matrix is almost independent of sequence divergence.
similar results are found with g <dig> , suggesting that the best matrix is also approximately independent of the gap model, as might be expected.
recommended matrix choice
in all the tests reported here, and others not shown, vtml <dig> is the best or close to best choice and is recommended as the general-purpose choice.
in the vt and blosum series, vt <dig> and blosum <dig> respectively are recommended.
in the clearly inferior pam series, pam <dig> to pam <dig> appear to be best.
improving clustalw accuracy
the above results lead to the conclusion that clustalw's heuristic use of distance-dependent matrices is probably not effective and suggest that a superior strategy would be to use a single high-accuracy matrix.
i tested this in clustalw v <dig> .
<dig>  <cit>  by setting the matrix to vtml <dig>  this was necessarily an integer-rounded version as floating-point matrices are not supported by the software; i chose to use 1/3-bit units.
the gap-open penalty was set to  <dig>  and gap-extend to  <dig> .
these penalties were chosen after trying a few reasonable values; no attempt was made at systematic optimization.
results are shown in table  <dig>  an 8% improvement is seen in the column score on balibase v <dig>  this suggests that a small reduction in pair-wise alignment errors due to a better choice of substitution matrix can have a cumulative effect in a multiple alignment and yield a substantial improvement.
however, it is not sufficient to raise the accuracy of clustalw to those of more recent methods such as muscle  <cit>  or probcons  <cit> , as shown in table  <dig>
improvement in accuracy of clustalw v <dig> .
<dig> by using the vtml <dig> matrix.
b <dig> is balibase v <dig>  p <dig> is prefab v <dig>  tc and sps are the average fraction of correctly aligned columns and residue pairs, respectively.
in prefab reference alignments are pair-wise so sps and tc are equivalent.
defaults is the accuracies using default parameters, vtml <dig> with that matrix and gap parameters given in the main text, and >acc is the improvement due to vtml <dig>  accuracies of the more recent methods muscle v <dig>  and probcons v <dig>  are given for comparison.
terminal gaps
the g <dig> model allows terminal gaps to be assigned different open and extend penalties versus internal gaps .
on bali and prefab optimized terminal gap open penalties are roughly half those of internal penalties, validating the scheme used in mafft  <cit>  and muscle in which terminal gaps have the same extension penalty as internal gaps and half the open penalty.
however, these results should be interpreted with caution since sequences in both bali and prefab are trimmed to domain boundaries.
multi-domain proteins are common, and different parameters and indeed different alignment algorithms may be appropriate when sequences have different domain organizations and are thus not globally alignable  <cit> .
optimized gap parameters reported by pop for the bali, prefab, termgaps and intgaps sets.
models are g <dig>  and g <dig> .
penalties are g , e , g  and e .
q is the accuracy score.
the matrix was the full-precision, one-bit unit version of vtml <dig>  one of the best matrices according to the tests reported here.
random element
pop has a random element due to the selection of random subsets of the training data in its first two stages.
it is therefore of interest to investigate how results vary for the same training data when different random number seeds are used.
i chose to do this using the termgaps set as practical experience shows it to have the most challenging parameter landscape of the sets considered here.
i ran pop ten times for the g <dig> model on this set and found the following characteristics of the resulting q values: mean  <dig> , standard deviation  <dig>  × 10- <dig>  maximum  <dig> , minimum  <dig> .
these results suggest a low sensitivity to subset selection and are consistent with finding the global optimum to within approximately three significant figures.
integer rounding
integer rounding causes a small, but consistent, degradation in accuracy of around  <dig> % to  <dig> % in vtml and  <dig> % to  <dig> % in blosum .
corrected blosum matrices
rblosum gave a small but again consistent improvement in accuracy over blosum of around  <dig> %.
this result is surprising considering that the opposite effect was found when testing homolog recognition  <cit> .
discussion
pop is an ad hoc algorithm without a strong theoretical foundation and has more heuristic parameters than one would ideally like.
it does not scale well to larger numbers of parameters.
implementing and using pop requires an understanding of the input data and some trial and error.
however, consistency of trends across different benchmarks, even when differences are very small, and consistency when run with different random number seeds combine to suggest that pop may find a global optimum to within approximately three significant figures on the pair-wise global alignment tests considered here, and regardless provide strong evidence that pop provides a sensitive test of differences between models with the same number of parameters.
it is also readily parallelized and, unlike some other approaches, can easily be applied to multiple alignment.
the hill-climbing in pop can start from parameters proposed by some other method, such as ipa, providing a lower bound on the training set error.
it should be noted that the reduced training set error achieved by pop may come at some cost in generalization error; this is a topic for further study.
compute resources required by pop are relatively modest , needing memory well within the range of current commodity computers and times  ranging from minutes to hours on the tests considered here.
the modern vtml matrices proved to give the best accuracy, in agreement with studies of homolog recognition  <cit> .
the vt and blosum series were not far behind, except for low-identity blosums which performed relatively poorly on all tests, even when aligning low-identity sequences.
the rblosum series  was marginally better, but still inferior to vtml.
full-precision matrices were also marginally better than the integer-rounded versions in common use.
it is natural to expect that improved pair-wise alignment accuracy will lead to improved multiple alignments, and this is a direction that deserves further exploration.
the accuracy of clustalw was significantly improved by using vtml <dig>  it is therefore of interest to review the matrices employed by other multiple aligners.
muscle uses vtml <dig>  which is close to vtml <dig> on most tests and therefore appears to have been a fortuitous choice.
mafft v <dig> used pam <dig>  a clearly suboptimal choice; possibly this partly explains the lower accuracy of mafft v <dig> relative to the similar algorithm of muscle v <dig>  it is not clear to me which matrices are used in more recent versions of mafft, though the authors report testing members of the blosum and jtt series  <cit> .
align-m  <cit>  uses blosum <dig>  given the rapid drop in accuracy observed between blosum <dig> and blosum <dig> on all benchmark tests, it seems likely that a significant improvement would result from using a different matrix.
kalign  <cit>  uses gonnet <dig>  <cit>  which i found to perform comparably with vtml <dig> .
mummals  <cit>  and probcons use blosum matrices, suggesting the possibility of a small improvement in those programs by using vtml instead.
this following table shows the cpu time  and memory  used to compute the results shown in table  <dig>
an examination of fig.
<dig> shows that blosum60- <dig> achieves almost the same accuracy as vtml150- <dig> on balibase but is roughly 1% worse on prefab.
this suggests that balibase may be biased towards blosum <dig> due to the use of sequence methods in alignment construction .
the source code for an example implementation of pop is available from the author's web site at http://www.drive <dig> com/pop.
this is designed for optimization of the g <dig> and g <dig> models described here, but can be modified relatively easily for other models by replacing the appropriate q and alignment functions.
as i can attest from an abundance of personal experience, manual optimization of alignment parameters is a tedious process that rarely leaves the practitioner feeling confident in the results.
having an automated method at one's experimentation with new and modified algorithms by enabling a relatively trustworthy and painless evaluation of their relative effectiveness.
i plan to use pop to explore ideas for improved pair-wise and multiple global alignment algorithms.
CONCLUSIONS
on the basis of this analysis, the vtml <dig> matrix is recommended as the most appropriate in general when global alignment accuracy is desired.
my results suggest a bias in balibase towards the blosum series of matrices of around 1% in accuracy.
while the effect is small, a bias of this magnitude could be significant in validations of multiple alignment methods because differences between the better methods on balibase are of comparable size.
bias towards substitution matrices or gap penalty functions is not unexpected as only 13% of balibase sequences have solved structures, and its alignments were therefore constructed mostly by the use of sequence rather than structure methods.
future studies of alignment accuracy should use data derived independently of sequence in order to avoid such biases.

 sequence alignment is a fundamental tool in contemporary biology. most algorithmic formulations of the problem seek an alignment maximizing a function known as the objective score. objective scores are usually defined as a sum of terms for matching pairs of letters and penalties for gaps. while the effects of different substitution matrices and gap parameters have been extensively studied in the context of local alignment and homolog recognition, their effects on alignment accuracy, especially for global alignment, are not so well understood. several heuristics are in common use, for example clustalw's choice of low-identity matrices for aligning low-identity sequences, which have not to the best of my knowledge been empirically tested. one factor impeding such testing is the lack of effective automated methods for optimizing parameters for a given objective function. previous work in this area has included unsupervised expectation maximization, discriminative training, inverse parametric alignment and maximum-margin structured learning. katoh et al. reported using golden section search to optimize gap penalties, but this generally assumes a unimodal function which does not hold in this case. exhaustive search has been attempted several times, for examples see. in this work i describe a new parameter optimization procedure, pop, compare it with the ipa method of, and use it to investigate a number of questions related to global protein alignment accuracy, including: which is the best substitution matrix, do the best choices of matrix and gap penalties vary with sequence identity, should terminal gaps be penalized differently from internal gaps, and does the loss of precision due to integer rounding in substitution matrices degrade alignment accuracy? ipa was chosen for comparison because it is a recently published method, code was readily available, and because it optimizes parameters exactly equivalent to those studied in this work. to the best of my knowledge there has been no prior comparative study of alternative gap optimization methods, so ipa's effectiveness relative to other approaches is not known. in the tests performed here, pop was shown to achieve from % to % higher training-set accuracy than ipa, and evidence is presented that pop finds parameters within about % of the true optimum, i.e. an improvement in training set error of up to an order of magnitude. differences of a fraction of a percent will be shown to be informative in assessing differences between substitution matrices, for which a method with the improved precision of pop is required. pop, like ipa, is designed to maximize training set accuracy, which may be achieved at some cost in generalization error, i.e. worse performance on novel data. in the case of pop this is by design: the problem of minimizing training set error is simpler than minimizing generalization error and the experiments reported here were not designed to compute the most biologically appropriate parameters for a given algorithm but rather to compare models with small and equal numbers of free parameters. it is then reasonable to assume that training set accuracy is a good measure of the relative performance of those models. let w = wi, i =..n be the parameters of interest, and q be the function to be optimized. the goal is to find values wopt = argmax q that maximize q. this optimization problem is challenging for several reasons. sufficiently small changes in w will leave all alignments unchanged and hence most q functions of interest have zero partial derivatives almost everywhere. also, q is typically expensive to compute, requiring seconds to hours to evaluate at a single point, and is non-convex with many local maxima. methods ipa kim and kececioglu described an inverse parametric alignment method that can be trained on so-called "partial examples" such as those found in balibase in which only a subset of columns are annotated as reliably aligned. eagu kim kindly provided a software implementation of this algorithm. this implementation extended the method to allow separate terminal and internal gap penalties to be learned. pop the central idea in pop is to search for changes in parameters that result in "goldilocks" changes in q--not too big, not too small, but changes in just of the right size to have a good chance of indicating coarse trends rather than small-scale noise. this approach can be regarded as a generalization of line search optimization. as theoretical and empirical methods for dealing with non-convex, non-continuous q functions are currently limited, practical experience, especially examination of plots such as fig. is the best guide to determining the size of a "just right" change. in the experiments reported here, values around % were most effective. pop proceeds in three stages. in the first stage, an n-dimensional hypercuboid is explored by evaluating q at each point in a regularly spaced grid. local maxima are identified as points at which q is greater than all neighboring points, and the best of these are used as starting points for the second stage. the second and third stages use a hill-climbing strategy to approach a local maximum given a starting point. the second stage uses a faster but less accurate variant of the hill-climbing method than the third stage, which starts from the best local optima found in the second stage. the final result is the best maximum found by the third stage. to save computation time, the three stages typically use increasingly large subsets of the training samples, with the first two stages using randomly chosen subsets and the third stage the entire training set. hill-climbing hill climbing repeats the following procedure until no improvement is found. starting from a point w, each axis i is explored in both the positive and negative direction; i.e., all parameters are held fixed but one. let δ be a proposed change in wi, = wi + δ, w' = {wj, j ≠ i, }, q' = q, i.e. q' is the value of q after a proposed move along one axis, and Δ = q' - q. in each direction a δ is found that gives Δ >, or, failing that, a δ that reduces q by an amount that is small, but not too small, say in the range % to 1%. very small changes are undesirable because they will tend to be dominated by "noise" rather than indicating a systematic trend. large reductions in q will tend to take wi out of the range where it might make an improvement when combined with a change along another axis. a heuristic designed to minimize these problems introduces a maximum and minimum reduction in q. after computing the proposed changes, a new w is determined as described in making a move below. proposed changes a proposed change δ is determined as follows. in the first iteration δ is set to a small fraction of the absolute value of wi or a small value if wi is zero, or a "hint" value set by the user. in subsequent iterations the initial value of δ is its final value from the previous iteration. having set the initial value, the following procedure is then repeated. the value of Δ is computed. if Δ = δ is multiplied by if Δ < and -Δ <λ, i.e. q is reduced by too small an amount, δ is doubled. if Δ < and -Δ >μ, i.e. q is reduced by too much, δ is halved. otherwise, Δ > and further exploration is tried. first, smaller values of δ are evaluated in an attempt to discover an intermediate maximum; if one is found this is repeated to investigate increasingly small changes. if smaller values give no further improvement, larger values are tried in a similar way. a maximum improvement in q, denoted Γ, is imposed in order to avoid extreme changes in a single parameter when far from an optimum. making a move the search for proposed changes yields a set of 2n values {δi+, δi-, i =... n} where the + and - superscripts indicate moves in the positive and negative directions respectively. all moves are considered that apply the proposed changes on from k =... k axes simultaneously, where k is a user-settable parameter may be... n. when k > new evaluations of q are required. for example, k = considers eight moves for each pair of parameters:,,,,,,, and where + or -indicates an increase or decrease in the parameter and indicates no change. the move giving the best improvement in q is accepted, otherwise the routine terminates. naively, one might expect that consideration of single-axis moves only would suffice, but in practice it turns out that allowing moves on two or more axes sometimes gives a significantly better final result, with k = sufficient in most cases. it is desirable to keep k small as the number d of evaluated moves per iteration grows very rapidly with k: it also turns out that allowing moves along axes that reduce q can, when combined with moves on other axes, improve q and that allowing this possibility also gives significantly improved optimization in some cases. for example, a small increase in the gap open penalty and a small decrease in the gap extend penalty might each reduce q, but increase q when both changes are made simultaneously. fast hill-climbing a "fast" variant of the hill-climbing procedure sets k = and immediately applies any proposed move that is found to improve q. speed is also improved by increasing μ, λ and Γ. these modifications reduce the number of times the q function is invoked, saving execution time but sometimes yielding significantly inferior parameters. directions to try, and the sign first tried, may be selected at random or cycled to minimize systematic bias. gap models and substitution matrices let a model be the set of parameters associated with a given objective function. in this work the substitution matrix is regarded as fixed; only gap penalty parameters are included. a gap is a series of indel symbols; formally, a maximal consecutive sequence of indels. if a gap includes the first or last column of an alignment it is described as terminal; all other gaps are internal. while pop has been implemented for a wide variety of models, most of these are works in progress and i will therefore report results for just two models: g in which the same affine penalties are applied to all gaps, and g in which internal and terminal gaps have different open and extend penalties. the parameters are g, e, g and e; the penalty for an internal gap of length l is g + e and g + e for a terminal gap. thus g is a special case of g in which g = g and e = e. the objective score is then the sum of substitution scores minus gap penalties; a maximum-scoring global alignment is found using standard dynamic programming techniques. the following substitution matrix types are considered: blosum, rblosum, pam, jtt and vt/vtml. the rblosum matrices were constructed using a bug-fixed version of the program used to compute blosum matrices from the blocks database. surprisingly, the corrected blosum matrix was found to slightly degrade performance in homolog recognition. each matrix family is a series with members defined by a measure of evolutionary distance: percent identity cutoff in the case of blosum and rblosum, pam distance for the rest. conventionally, integer valued matrices are used in which log-odds scores in fractional bits have been rounded to one or two significant figures. presumably this is for historical reasons: in older computer processors integer arithmetic was faster than floating point. this is no longer the case for many general-purpose processors, and regardless it is of interest to investigate whether the loss of precision due to integer rounding has an effect on alignment accuracy. unless otherwise stated, full precision, one bit unit matrices were used. benchmark data reference data was taken from three protein alignment benchmarks: prefab version and balibase versions and there are pair-wise reference alignments in prefab v extracting all pairs from the multiple alignments in balibase versions and gives and pairs respectively. the large number of pairs in version motivated the use of the more tractable version for all-pairs tests. subsets and of version were used to investigate optimizing terminal gap penalties. accuracy measure the accuracy measure for a single pair is the number of correctly aligned residue pairs divided by the number of residue pairs in the reference alignment. the total accuracy score is the weighted average of q over all pairs, where the weighting is uniform in the case of prefab and the inverse of the number of pairs in the original multiple alignment in the case of balibase. it would have been desirable to estimate error bars using a method such as the bayesian bootstrap; however this proved to be infeasible due to limitations in available computer time. for brevity, the reference sets will be denoted bali, prefab, termgaps, and intgaps. to investigate the effects of evolutionary distance three subsets of balibase v pairs were constructed: randomly selected pairs with identities in the range 0-33%, 33-66% and 66-99%, respectively. these were selected from the full-length rather than domain-trimmed sequences.  comparison of pop and ipa i compared pop and ipa using the termgaps set, where the biggest difference between g and g might be expected, and on bali. the blosum matrix in 1/ bit units was used as this was hardcoded into ipa. results are shown in table pop was found to be from % to % more accurate than ipa; these improvements are typical. training set accuracy for the g and g models on bali and termgaps using blosum 1/3-bit units. in the numbers in the pop and ipa columns are the q values achieved by the gap para-meters reported by each method. in all cases, pop achieves higher accuracy, the >acc column shows the amount by which the pop accuracy is better than ipa. substitution matrix family fig. shows the results of optimizing model g on the bali and prefab sets. the results are qualitatively similar on the two sets, giving confidence that they indicate general trends rather than artifacts of benchmark construction, of overtraining or of significantly suboptimal local maxima. this is further confirmed by cross-training, which again gives similar, albeit noisier, results, as would be expected. the results show vtml to be the best matrix series, with vtml > vt > blosum > jtt > pam holding for most members, though the differences between vtml, vt and blosum are small except at the extreme high and low-distance ends of the series. the pam series is significantly inferior to vtml, giving accuracy scores around 2% lower. matrix selection by evolutionary distance it has been suggested that different substitution matrices might be more effective at different evolutionary distances. for example, clustalw uses more distant matrices to align more distant sequences, and lassmann and sonhammer report that "softer", i.e. more distant, matrices are better at aligning more more distant sequences while are equally good with more closely related sequences. to investigate this, i optimized g on the id0_ id33_ and id66_ sets, with the results shown in fig, interestingly, the plots are qualitatively similar for the three sets despite increasing pair-wise identity and the increasingly narrow variation in optimal q. remarkably, in the case of id66_ accuracies are all above % and the difference between the best matrix and worst is only %, yet trends observed at lower identities are still clearly discernible. the peak in each curve that identifies the best matrix in each family is at approximately the same evolutionary distance on each set, showing that the best choice of matrix is almost independent of sequence divergence. similar results are found with g, suggesting that the best matrix is also approximately independent of the gap model, as might be expected. recommended matrix choice in all the tests reported here, and others not shown, vtml is the best or close to best choice and is recommended as the general-purpose choice. in the vt and blosum series, vt and blosum respectively are recommended. in the clearly inferior pam series, pam to pam appear to be best. improving clustalw accuracy the above results lead to the conclusion that clustalw's heuristic use of distance-dependent matrices is probably not effective and suggest that a superior strategy would be to use a single high-accuracy matrix. i tested this in clustalw v. by setting the matrix to vtml this was necessarily an integer-rounded version as floating-point matrices are not supported by the software; i chose to use 1/3-bit units. the gap-open penalty was set to and gap-extend to. these penalties were chosen after trying a few reasonable values; no attempt was made at systematic optimization. results are shown in table an 8% improvement is seen in the column score on balibase v this suggests that a small reduction in pair-wise alignment errors due to a better choice of substitution matrix can have a cumulative effect in a multiple alignment and yield a substantial improvement. however, it is not sufficient to raise the accuracy of clustalw to those of more recent methods such as muscle or probcons, as shown in table improvement in accuracy of clustalw v. by using the vtml matrix. b is balibase v p is prefab v tc and sps are the average fraction of correctly aligned columns and residue pairs, respectively. in prefab reference alignments are pair-wise so sps and tc are equivalent. defaults is the accuracies using default parameters, vtml with that matrix and gap parameters given in the main text, and >acc is the improvement due to vtml accuracies of the more recent methods muscle v and probcons v are given for comparison. terminal gaps the g model allows terminal gaps to be assigned different open and extend penalties versus internal gaps. on bali and prefab optimized terminal gap open penalties are roughly half those of internal penalties, validating the scheme used in mafft and muscle in which terminal gaps have the same extension penalty as internal gaps and half the open penalty. however, these results should be interpreted with caution since sequences in both bali and prefab are trimmed to domain boundaries. multi-domain proteins are common, and different parameters and indeed different alignment algorithms may be appropriate when sequences have different domain organizations and are thus not globally alignable. optimized gap parameters reported by pop for the bali, prefab, termgaps and intgaps sets. models are g and g. penalties are g, e, g and e. q is the accuracy score. the matrix was the full-precision, one-bit unit version of vtml one of the best matrices according to the tests reported here. random element pop has a random element due to the selection of random subsets of the training data in its first two stages. it is therefore of interest to investigate how results vary for the same training data when different random number seeds are used. i chose to do this using the termgaps set as practical experience shows it to have the most challenging parameter landscape of the sets considered here. i ran pop ten times for the g model on this set and found the following characteristics of the resulting q values: mean, standard deviation × 10- maximum, minimum. these results suggest a low sensitivity to subset selection and are consistent with finding the global optimum to within approximately three significant figures. integer rounding integer rounding causes a small, but consistent, degradation in accuracy of around % to % in vtml and % to % in blosum. corrected blosum matrices rblosum gave a small but again consistent improvement in accuracy over blosum of around %. this result is surprising considering that the opposite effect was found when testing homolog recognition. discussion pop is an ad hoc algorithm without a strong theoretical foundation and has more heuristic parameters than one would ideally like. it does not scale well to larger numbers of parameters. implementing and using pop requires an understanding of the input data and some trial and error. however, consistency of trends across different benchmarks, even when differences are very small, and consistency when run with different random number seeds combine to suggest that pop may find a global optimum to within approximately three significant figures on the pair-wise global alignment tests considered here, and regardless provide strong evidence that pop provides a sensitive test of differences between models with the same number of parameters. it is also readily parallelized and, unlike some other approaches, can easily be applied to multiple alignment. the hill-climbing in pop can start from parameters proposed by some other method, such as ipa, providing a lower bound on the training set error. it should be noted that the reduced training set error achieved by pop may come at some cost in generalization error; this is a topic for further study. compute resources required by pop are relatively modest, needing memory well within the range of current commodity computers and times ranging from minutes to hours on the tests considered here. the modern vtml matrices proved to give the best accuracy, in agreement with studies of homolog recognition. the vt and blosum series were not far behind, except for low-identity blosums which performed relatively poorly on all tests, even when aligning low-identity sequences. the rblosum series was marginally better, but still inferior to vtml. full-precision matrices were also marginally better than the integer-rounded versions in common use. it is natural to expect that improved pair-wise alignment accuracy will lead to improved multiple alignments, and this is a direction that deserves further exploration. the accuracy of clustalw was significantly improved by using vtml it is therefore of interest to review the matrices employed by other multiple aligners. muscle uses vtml which is close to vtml on most tests and therefore appears to have been a fortuitous choice. mafft v used pam a clearly suboptimal choice; possibly this partly explains the lower accuracy of mafft v relative to the similar algorithm of muscle v it is not clear to me which matrices are used in more recent versions of mafft, though the authors report testing members of the blosum and jtt series. align-m uses blosum given the rapid drop in accuracy observed between blosum and blosum on all benchmark tests, it seems likely that a significant improvement would result from using a different matrix. kalign uses gonnet which i found to perform comparably with vtml. mummals and probcons use blosum matrices, suggesting the possibility of a small improvement in those programs by using vtml instead. this following table shows the cpu time and memory used to compute the results shown in table an examination of fig. shows that blosum60- achieves almost the same accuracy as vtml150- on balibase but is roughly 1% worse on prefab. this suggests that balibase may be biased towards blosum due to the use of sequence methods in alignment construction. the source code for an example implementation of pop is available from the author's web site at http://www.drive com/pop. this is designed for optimization of the g and g models described here, but can be modified relatively easily for other models by replacing the appropriate q and alignment functions. as i can attest from an abundance of personal experience, manual optimization of alignment parameters is a tedious process that rarely leaves the practitioner feeling confident in the results. having an automated method at one's experimentation with new and modified algorithms by enabling a relatively trustworthy and painless evaluation of their relative effectiveness. i plan to use pop to explore ideas for improved pair-wise and multiple global alignment algorithms.  on the basis of this analysis, the vtml matrix is recommended as the most appropriate in general when global alignment accuracy is desired. my results suggest a bias in balibase towards the blosum series of matrices of around 1% in accuracy. while the effect is small, a bias of this magnitude could be significant in validations of multiple alignment methods because differences between the better methods on balibase are of comparable size. bias towards substitution matrices or gap penalty functions is not unexpected as only 13% of balibase sequences have solved structures, and its alignments were therefore constructed mostly by the use of sequence rather than structure methods. future studies of alignment accuracy should use data derived independently of sequence in order to avoid such biases. 
source code and test data is freely available from the author's web site at http://www.drive com/pop. here a new parameter optimization procedure, pop, is described and applied to the problems of optimizing gap penalties and selecting substitution matrices for pair-wise global protein alignments. the pam series is significantly worse, giving accuracies typically 2% less than vtml. while substitution matrices can readily be computed from reference alignments, it is challenging to compute optimal or approximately optimal gap penalties. it is also not well understood which substitution matrices are the most effective when alignment accuracy is the goal rather than homolog recognition. blosum matrices are found to be slightly inferior, even with the marginal improvements in the bug-fixed rblosum series. the hypothesis that more accurate alignments of distantly related sequences may be achieved using low-identity matrices is shown to be false for commonly used matrix types.
 while substitution matrices can readily be computed from reference alignments, it is challenging to compute optimal or approximately optimal gap penalties. it is also not well understood which substitution matrices are the most effective when alignment accuracy is the goal rather than homolog recognition. here a new parameter optimization procedure, pop, is described and applied to the problems of optimizing gap penalties and selecting substitution matrices for pair-wise global protein alignments.  pop is compared to a recent method due to kim and kececioglu and found to achieve from % to % higher accuracies on pair-wise benchmarks extracted from balibase. the vtml matrix series is shown to be the most accurate on several global pair-wise alignment benchmarks, with vtml giving best or close to the best performance in all tests. blosum matrices are found to be slightly inferior, even with the marginal improvements in the bug-fixed rblosum series. the pam series is significantly worse, giving accuracies typically 2% less than vtml. integer rounding is found to cause slight degradations in accuracy. no evidence is found that selecting a matrix based on sequence divergence improves accuracy, suggesting that the use of this heuristic in clustalw may be ineffective. using vtml is found to improve the accuracy of clustalw by 8% on balibase and 5% on prefab.  the hypothesis that more accurate alignments of distantly related sequences may be achieved using low-identity matrices is shown to be false for commonly used matrix types. source code and test data is freely available from the author's web site at http://www.drive com/pop. 