BACKGROUND
a major problem in biology is how to convert the data provided by dna or protein sequences into functional information.
for this reason, a significant fraction of molecular evolution studies are focused on the statistical characterization of the patterns of change of dna or protein sequences.
they are based on the general idea that modifications in the function of a gene are often related to changes in the selective regime acting on it, in such a way that a characteristic imprint is generated in its sequence.
for example, if a gene is modifying its function by a process that involves positive selection, we would expect to find very rapid changes in the amino acid sequence of the encoded protein, at a rate much higher than expected under neutral evolution.
on the contrary, if regions of a gene become functionally constrained, under negative selection, change in those regions will be very slow  <cit> .
different methods have been proposed to determine how selection acts on biological sequences .
several of them have been devised to compare the synonymous  and non-synonymous  rates of change in coding regions, often with the purpose of testing whether positive selection has been acting upon those regions  <cit> .
however, relative values of ka/ks >  <dig>  that are strong evidence for positive selection – being higher than the expected rates under the null hypothesis of neutral evolution  – are rarely found.
most frequently, it is determined that sequences are under strong negative selection.
for example, in a recent work in which  <dig> orthologous genes were compared in human and mouse, it was estimated that, in average, ka/ks =  <dig>   <cit> .
other methods are focused on non-synonymous changes alone, and try to model whether amino acid changes occur homogeneously or, alternatively, are concentrated on particular positions or regions of the proteins  <cit> .
several other methods have been proposed to explore whether the rates of change of genes vary among different lineages, especially among species or groups of species  <cit> .
specific tools devised to analyze functional divergence between duplicate genes have been also developed  <cit> .
this is due to the fact that functional shifts occur especially often in association with gene duplication events.
either by acquisition of new functions by one of the duplicates  or by division of the functions of the original gene between the two paralogs , duplication is often expected to radically alter the selective forces acting on the duplicate genes  <cit> .
however, most of the approaches developed so far have some serious limitations.
first, many of them are computationally cumbersome.
authors often have to choose between using the simplest tools available, that are relatively insensitive, or performing more precise analyses, but with a limited number of genes.
a second problem is that many of the methods devised are based on complex mathematical models of how dna sequences change, and it is often unknown how deviations from the assumptions implicit in those models may affect the conclusions obtained.
in fact, there are considerable discussions in the literature on whether certain types of analyses tend to generate spurious significant results .
therefore, there is a general need for tools that combine the features of being intuitive , fast and also sufficiently sensitive.
in this study, we describe a new bioinformatic tool, called uvpar, which may be used to detect regional changes in constraints in the protein sequences of duplicate genes.
the uvpar algorithm is a substantial refinement of an analytical strategy devised before by one of us and already successfully used for characterizing ancient functional changes in a family of atpases/atp synthases  <cit> .
the basis of that strategy is to determine, combining sliding-window analyses of the degree of amino acidic conservation and permutation tests, the regions of duplicate genes that have evolved at different rates in two species.
in its current implementation, uvpar allows for the fully automatic analysis of a large number of protein sequences in a short period of time.
as examples of its potential, we show the results of two analyses.
first, we performed uvpar comparisons for hox <dig> duplicate genes in six vertebrate species.
we demonstrated that some hox <dig> genes suffered in the past positively selected changes in their sequences associated to functional shifts after gene duplication events  <cit> .
here, we compare these previous results with the ones obtained with uvpar in order to establish how our novel method relates to other approaches.
as a second example, we generated a comprehensive study of the rbr gene family  <cit>  that involves several hundreds of analyses, to show how the program can be used at a large scale.
in summary, uvpar is a useful novel tool for studies focused on the characterization of functional shifts in proteins, most especially in cases in which the genes of interest are part of complex gene families.
implementation
let us consider the case in which the phylogenetic analyses of the sequences of certain genes has established that a duplication occurred, generating two paralogous genes, before two lineages of organisms split .
after both lineages became separated, the two genes accumulated differences until the present day .
the question that we want to tackle is whether sequence changes have accumulated differently in those genes in both lineages since their separation.
this often can be visualized as a difference in rates in one or several of the genes in a phylogenetic tree .
uvpar is a program specialized in analyzing these situations.
it takes into account the previous knowledge of orthology/paralogy generated by these phylogenetic analyses.
moreover, uvpar analyses are based on the multiple-sequence alignments used to generate the phylogenetic trees.
uvpar is written in c and we have compiled versions for windows and linux operating systems.
in its simplest implementation, the program uses as input a text file containing the sequences of the proteins encoded by four genes – two duplicate genes in two different species, such as a <dig>  b <dig>  a <dig> and b <dig> in figure  <dig> – in fasta format.
however, the program also accepts larger datasets, in which pairs of duplicates of multiple species are included in a single text file.
in this case, the program automatically performs analyses for all the possible combinations of species in the dataset.
the program interprets the sequences in the text file as the output of a multiple-sequence alignment .
once imported, the dataset is filtered, in such a way that positions in which gaps are present are eliminated from all the sequences.
the program then detects the sequence with smallest size  and from then on, all the analyses are performed with the first n amino acidic positions of each sequence.
after these two filters, uvpar progressively reads the sequences of the first two proteins  and calculates a value of similarity for each amino acid position of the two sequences, blk <dig>  according to a given amino acid substitution matrix.
the user may choose among the following matrices: pam <dig>  pam <dig>  pam <dig>  blosum <dig>  blosum <dig> and blosum <dig>  blosum <dig> is the one that we use by default.
then it does the same for the two duplicates of the second species and calculates a similar value, blk <dig>  finally, the program establishes the difference between the two bl values, that we named constraint value, c <dig> = blk <dig> - blk <dig>  <cit> .
c <dig> is thus a value obtained for "gene quartets", defined as pairs of duplicates in two different species.
if both duplicates are equally different in the two species, the constraint value should be about zero for all the positions.
the key of the strategy implemented in uvpar is to analyze the distribution of c <dig> values observed to determine whether some positive or negative values are significantly clustered together.
to do so, the vector of c <dig> values is randomly shuffled a number t of times  and the shuffled vectors are compared with the one determined for the real paralogs.
the comparison is performed for windows of increasing size , from w =  <dig> to w = n -  <dig>  for each window size, the maximum and minimum sums of c <dig> values, that we call s, are determined for both the original and the shuffled sequences.
then, uvpar compares the maximum and minimum sum in the original sequence with those in the set of shuffled sequences.
this set provides a distribution of maximum  values against which the values obtained in the real sequence is contrasted.
two hypotheses must be tested for each window size , and therefore it is convenient to use bonferroni's correction.
thus, when the value in the original sequence is found beyond the top or bottom  <dig>  % of the values of the simulated sequences, it is considered significant and uvpar selects the corresponding window for further analyses.
the second part of the analysis performed by uvpar is the comparison of all significant windows for a given gene quartet.
uvpar often detects multiple windows of different sizes that are significant.
in many cases those windows are nested and therefore they refer in part to the same positions in the analyzed sequences.
the biological interpretation of the significant results requires establishing the windows that best explain the detected constraint changes.
this problem can be tackled in different ways.
in our previous study  <cit> , we solved the problem of nested windows in the simplest way, just choosing the window with the largest  value of s and discarding the rest.
this approach was appropriate for the particular cases analyzed in that study, but it has been determined to be insufficient for complex cases found in our subsequent analyses.
this is due to two reasons.
first, we have found that ties in the s values often appear in significant windows of different sizes.
second, we have found significant partially overlapping windows, a case that was never detected before.
therefore, in uvpar we have implemented a much more precise way of dealing with multiple significant windows of different sizes.
this is a major improvement on the strategy described on our previously published work
the method is as follows: the program separates the significant windows into two groups, corresponding to those significant in the analyses involving minimum values and those derived from the analyses of maximum values.
within each group, windows are ordered according to their sizes.
the program then compares, for one of those two groups, all the significant windows, starting with the smallest one and moving progressively to larger windows.
every time that a window contained in a larger one is found that has an absolute s value equal or higher than that in the larger window, the latter is eliminated.
when this process has finished, the program performs the same search for all the significant windows in the other group.
the goal of this iterative screening is to eliminate all the large windows that are significant simply because they contain a small highly deviant region.
for example, we found cases in which short regions had such positive c <dig> values that even when negative values were added, the larger windows that contains both the highly positive "seed" and the additional negative values were still significantly positive.
those larger windows are eliminated after the comparisons described.
at the end of this first screening, we are left with all the nested significant windows in which an increase of length is associated with an increase of absolute s value plus all non-nested windows.
then, uvpar performs a second screening with the remaining windows, this time starting with the largest significant window and comparing it with all the smaller ones.
now, all windows contained in a larger one and at the same time with absolute s value lower than that found in the larger window are eliminated.
the program keeps repeating this search until all windows of both groups  have been analyzed.
this second screening allows the elimination of all windows that can be extended to larger ones.
logically, the larger windows have to be preferred if their s values are higher.
after these two screenings, we are left with sets of non-nested significant windows, that may or may not partially overlap.
as a final step, uvpar analyzes all the remaining windows and eliminates those with an absolute average c value – obtained by dividing the absolute s value by the window size – lower than  <dig> .
this last step is devised to eliminate a perverse effect that occurs in rare cases in which a few highly positive or negative c <dig> values are clustered together in one of the extremes of the vector.
then, it occasionally occurs that this region not only produces, as expected, a small highly significant window of a particular sign, but also generates a very large complementary window of the opposite sign that often spans the rest of the sequences.
this happens even when the c <dig> values outside of the short significant region are randomly distributed around zero or, more often, are almost all of them zero.
the few cases that we have found in which this effect was present were eliminated when the  <dig>  cutoff value was used, and therefore we have considered it to be a logical last filter to be implemented by default in the uvpar algorithm.
after all the searches are finished, the program generates an output file that details all the relevant parameters: 1) the number of times  that the sequences have been shuffled to generate the distribution of s values; 2) the sequences analyzed; 3) the detailed amino acidic alignment, with all the corresponding c <dig> values; 4) all the significant windows, before the filtering process; and, 5) all the significant windows, after the filtering process is completed.
thus, the algorithm can be summarized quite simply as follows:
import file with the number of pairs of duplicates to be analyzed and corresponding sequences in fasta format
read t value
filter sequences to eliminate gaps, estimate final size n
repeat
calculate c <dig> values
if  <dig> ==  <dig> for all k) then skip gene quartet
else:
for  from w =  <dig> to w = n -1
estimate maximum, minimum s values for the original sequences
shuffle c <dig> values, t times
generate distribution of maximum, minimum s values from the shuffled sequences
compare maximum, minimum s values for the original and the shuffled sequences, select windows according to a significance level p <  <dig>  %
filter significant windows, increasing size
filter significant windows, decreasing size
filter significant windows for an average c <dig> value > | <dig> |
generate output file, including: t values, sequences, alignments, c <dig> values, significant windows before and after filtering
RESULTS
speed and performance
the software that we used in our original study was very slow, what precluded both to use large t values and the analysis of large datasets.
due to several highly significant algorithmic improvements, uvpar is about  <dig> times faster.
this makes possible to perform the same searches several times and with different t values to determine the degree of error associated to each number of permutation tests.
when we checked for the impact of t values, we determined that to obtain a reasonably precise estimation of the probability associated to a particular value in the distribution of shuffled s values, it is convenient to choose values of t larger than  <dig>  which we previously used due to computational limitations.
we have fully repeated the analyses presented in  <cit>  using uvpar and, although all the conclusions of the work are correct, we have found that the error associated to the estimation of p values with t =  <dig> was considerable .
uvpar analyses are fast enough as to easily handle values of t =  <dig>  which we consider now to be the minimum number of permutation tests to be used in this type of studies.
if the sequences are short enough,t may be increased to  <dig> or even  <dig>  as examples, using t =  <dig> and a standard pc computer , sequence quartets of  <dig> amino acids can be analyzed in about  <dig> seconds, quartets of  <dig> amino acids in about  <dig> minutes and quartets of long sequences of  <dig> amino acids in about  <dig> minutes.
time only vary slightly depending on the number of significant windows found.
these results mean that uvpar can be used at a large scale.
as an example, the analysis of a whole family of proteins that we will detail in section  <dig> , that involved more than  <dig> individual uvpar searches  required a computation time of only  <dig> hours, about  <dig>  minutes per analysis.
simulations to determine the reliability and sensitivity of uvpar searches
to check for the number of false positives and false negative results generated by uvpar under different conditions, we explored, using simulations, a tree containing four species, with two paralogs each .
we used covtree  <cit>  to apply different changes in the rates of a particular region of one of the paralogs in the branch previous to the split between two of the species .
then, we used uvpar to search for coherent significant windows, i. e. windows that included at least part of the region in which the rate shifts were applied and that were present in all the four comparisons in which we would expect to find significant results after the rate shifts .
simulations were performed using sequences with different degrees of heterogeneity of evolutionary rates among sites.
heterogeneity was modeled as a gamma distribution with a variable α parameter.
values of α from  <dig>   to α =  <dig>  were used, to cover the range of results detected for real proteins .
results are shown in table  <dig> for windows of size w =  <dig> .
we found only  <dig> – 3% of cases in which coherent windows are present in the control simulations  no matter the value of α. when the rate of evolution for the region of interest is increased, more cases are detected, and an increase in the average c <dig> is also observed.
if rates are increased ten times and α values are relatively large , coherent windows are detected in 35% to 55% of the cases .
negative results in these simulations are caused by random changes that obscure the effect of the rate increase.
very low α values  in general precluded the detection of regional shifts by uvpar.
these results suggest that, in cases in which multiple sequences can be analyzed, systematic positive uvpar results most likely will correspond to real regional changes.
they also indicate that uvpar losses sensitivity when α values are very small and that the program does not detect small shifts in evolutionary rates; changes must be strong to generate a significant signal for uvpar analyses.
functional shifts in hox <dig> genes
hox genes are critical in many basic developmental processes, most especially in the determination of cell fates along the anterior-posterior body axis  <cit> .
vertebrates have a large number of hox genes, generated by gene and genome duplications  <cit> .
in a previous work, we showed that episodes of positive selection affected hox <dig> genes in vertebrate species  <cit> .
more precisely, we determined that both the ancestral duplication that gave rise to the paralogous hox-a <dig> and hox-b <dig> genes and the tetraploidization that occurred in the amphibian lineage that includes the model species xenopus laevis radically altered the selective regime acting on those genes.
in the first case, the duplication was followed by a period of time in which hox-a <dig> genes diversified under positive selection.
in the second case, one of the hox-b <dig> genes that were produced in the xenopus ancestor by the genome duplication process  also changed under positive selection.
the ratios of non-synonymous vs. synonymous substitutions in the branches affected by those two selective shifts were estimated to be about eight to ten times larger that the average ratio for the rest of branches of the tree  <cit> .
three distinct types of analyses suggested that positive selection acted on the regulatory n-terminal region of the hox <dig> proteins, while the highly conserved homeodomain was not affected  <cit> .
the large shift in evolutionary rates and the relative high values of α for these sequences  mean that this is a favorable case to test whether uvpar is able to generate results comparable to those found with more complex methods.
figure  <dig> shows the results for uvpar analyses  for the six vertebrates in which both duplicates, hox-a <dig> and hox-b <dig>  were analyzed in our previous work .
alignments were the same used in  <cit> .
uvpar results are congruent with two of the main conclusions of our previous work.
first, comparisons involving the xenopus hoxb7-b gene showed significant shifts when compared with hox <dig> genes of other vertebrates in  <dig> out of  <dig> comparisons .
the exception, the comparison involving mus musculus, is mainly due to a few homoplastic residues in the mus musculus hoxb <dig> sequence which preclude the overall values to become significant.
second, all significant windows were found outside of the homeodomain .
it must be noted that the third main result – the rapid differentiation after the hox-a7/hox-b <dig> duplication – could not be possibly detected here because that process occurred before all lineages that can be analyzed with uvpar became separate.
in summary, these uvpar results agree with our previous findings and show that uvpar is able to provide sensitive evidence for functional shifts comparable to those obtained by more complex and time-consuming mathematical approaches.
uvpar analysis of rbr ubiquitin ligases
we decided to demonstrate the ability of the program to handle large datasets by analyzing proteins of a family, called rbr, that we defined some time ago  <cit> .
rbr proteins are ubiquitin ligases with important roles in the control of protein degradation and several of them are known to be involved in human diseases .
rbr proteins are characterized by having the rbr supradomain, that is composed by two ring fingers  that are separated by an ibr domain  <cit> .
we checked for modifications in the functional constraints by analyzing the rbr supradomain of multiple rbr proteins in several model organisms.
we took advantage of the recent completion by our group of a comprehensive analysis in which we generated protein alignments for  <dig> rbr proteins and phylogenetic trees that allowed us to determine the orthology/paralogy relationships for all members of the family  <cit> .
from that analysis, we selected all the available sequences for  <dig> rbr genes  in seven model organisms .
the estimated average value of α for these sequences is  <dig> .
therefore sites evolve homogeneously enough for uvpar analyses.
we took into account that not all the genes are present in all the species analyzed , and some of the genes had species-specific duplications.
once these peculiarities were sorted out, we performed all the possible combinations of gene quartets for the seven species.
a total of  <dig> analyses were generated.
in  <dig> cases , we found significant results.
they are summarized in tables  <dig> and  <dig>  according to, respectively, the species and genes involved .
table  <dig> shows that a single species, the fish danio rerio, present a number of significant results that are well above the average found for the set of species.
we think that, similarly to what occurred for the xenopus hox <dig> genes shown above, this may be related to the genome duplication that happened in the evolutionary lineage in which danio emerged  <cit> .
if we compare the data for the different rbr genes , we observe that the significant cases appear quite homogeneously.
ari <dig> is the only gene with a number of significant results that is in average higher than that of the rest of rbrs.
this is in part explained by the large number of analyses  in which the ari1/ari <dig> comparisons were significant.
figure  <dig> shows the results for those comparisons.
it becomes evident examining this figure that many of the positive results are correlated.
in drosophila melanogaster there are two ari <dig> genes, ari-1a and ari-1b .
significant shifts in the ibr domain are observed between the pair ari-1a/ari- <dig> of drosophila melanogaster and the ari1/ari <dig> genes of all vertebrate species.
also, we see important changes along the whole rbr supradomain when genes of the urochordate ciona intestinalis are compared with vertebrate genes.
therefore all these results can be explained by a few major changes in particular evolutionary lineages.
in favorable cases, we can establish in which lineage the changes must have occurred.
for example, considering that d. melanogaster/vertebrate comparisons involving the d. melanogaster ari-1b gene are in general non-significant, we can conclude that a functional shift in the ibr domain of the d. melanogaster ari-1a gene occurred since the ari1-a/ari1-b duplication that has made ari-1a different from their vertebrate orthologs.
in summary, the results shown in figure  <dig> suggests that functional diversification between the rbr supradomains of the ari <dig> and ari <dig> genes has occurred just a few times and that those changes most frequently affected the ring <dig> and ibr domains.
this result points towards an important role of the n-terminal domains of the rbr supradomain in acquisition of novel functional properties in some species, perhaps to act on new ubiquitination substrates, while the non-canonical ring <dig> would be less significant.
in figure  <dig>  we show the distribution of maximum and minimum values in the shuffled sequences and the relative position in this distribution for one of the significant values observed.
bell shaped, quasi-symmetrical distributions as those shown in that figure are common.
they correspond to the typical extreme value distributions that are expected for large t and w values .
we have however observed in some cases highly asymmetrical distributions, often with multiple peaks, when window sizes are small.
discussion
we have shown that the uvpar program to may be used to quickly detect modifications in the functional constraints of duplicate genes.
favorable circumstances for the program are a reasonably homogenous rate of amino acid substitution among sites and an intense change in the selective regime acting on one of the paralogs.
considering that under positive selection ka/ks >  <dig> while in real proteins, the average value is much lower , we may expect uvpar to be particularly useful to detect positive selection, most especially associated to functional shifts after gene duplication.
in respect to other related methods, that often rely on debatable a priori assumptions and difficult mathematical calculations, uvpar results are easy to obtain and highly intuitive.
uvpar analyses depend on just two a priori suppositions.
first, that protein changes can be measured with a matrix such as blosum <dig>  second, that the shifts affect contiguous amino acids, i. e. regions of the proteins, and therefore they can be detected by examining their primary sequences, using a sliding-window approach.
there is strong evidence that regional shifts in proteins sequences often occur, in genome-scale analyses  <cit> .
uvpar advantages are most obvious when the goal is to quickly explore, across many species, complex gene families that include several duplicate genes.
this is clearly shown by the analysis of the rbr family described above.
we have also demonstrated that general conclusions about what genes or species are more prone to suffer functional shifts can be deduced from the comprehensive study of gene families .
the fact that many genes can be examined in a very short time opens up the possibility of performing genome-scale studies, impossible to generate with any other related method.
an obvious statistical problem is that, when a large number of quartets are examined with uvpar tests, the quite permissive significance level used  is expected to generate a number of false positive results.
our experience, as well as the simulations shown above, suggests examining three aspects of the data to obtain conclusions not affected by that problem.
these three aspect are: 1) the significance level of each positive result and the average value of c <dig> for the significant windows.
in our simulations, false positive windows have low average c <dig> values .
if it is deemed necessary, the analyses can be made more strict by either changing the significance level or by increasing the conventional cutoff value c <dig> > | <dig> |; 2) the logic in structural or functional terms of the results obtained ; and, 3) the systematic consistency of the results for multiple species.
in our opinion, internal consistency of the data is a particularly useful way to filter the data for real positives.
this is clearly observed in our simulations, in which we checked for the consistency of uvpar results for multiple analyses, finding that false positive coherent results are rare .
thus, when it is found, as we have shown for the ari1/ari <dig> data , that the genes of multiple species of a lineage  generate the same results when compared with genes of an outgroup species , it is good evidence for a real shift to have occurred.
this means that uvpar finds several times the same imprint even although the analyzed sequences are considerably distinct.
after all, orthologous genes of fishes, birds and mammals are separated by hundreds of millions of years of independent evolution.
the combination of uvpar results for multiple species may be often used to trace with precision the moment and the lineage in which a change in functional constraints occurred.
in this study, we have shown that the duplication of hox <dig> genes in xenopus laevis and of ari- <dig> genes in drosophila melanogaster were associated with modifications in the constraints of at least one of the duplicates.
these conclusions can be easily deduced when the data are presented in a favorable format, as that shown in figures  <dig> and  <dig>  our approach may thus complement other analytical methods devised to detect selective processes acting on particular branches of a phylogenetic tree  <cit> .
CONCLUSIONS
uvpar generates very fast and still sensitive information about modifications of the selective constraints in duplicate genes.
it can be used at a large scale, to analyze complex gene families in multiple species.
its speed, simplicity of use and also the fact that its mathematical assumptions are very intuitive makes uvpar an excellent tool for all groups interested in analyzing the impact of natural selection on gene sequences, especially at a genomic scale.
availability
uvpar is written in c. windows and linux versions of uvpar are available, free for academic users, at the following web page: .
in that page, an example of how to use uvpar can also be found.
non-academic users may obtain a license to use or distribute the program by contacting the corresponding author.
authors' contributions
va wrote the uvpar code and contributed to the improvements in the analytical strategy implemented in the program.
mg and il performed the hox and rbr analyses as well as the simulations presented above.
mg also contributed to the development of the analytical strategy and tested different versions of the program.
the basic strategy implemented in uvpar was developed by im, who also coordinated the research, contributed to the improvements in the strategy shown here and wrote the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
click here for file
acknowledgements
we thank antonio marco for critical reading of the manuscript.
our group is supported by grants gen2001-4851-c06- <dig> and saf2003- <dig>  and grant gv04b- <dig> .
j. i. l. is the recipient of a predoctoral fellowship .
figures and tables
μ =  <dig>
μ =  <dig>
μ =  <dig>
μ =  <dig>
<dig>  ±  <dig>
<dig>  ±  <dig>
<dig>  ±  <dig>
p

 a major problem in biology is how to convert the data provided by dna or protein sequences into functional information. for this reason, a significant fraction of molecular evolution studies are focused on the statistical characterization of the patterns of change of dna or protein sequences. they are based on the general idea that modifications in the function of a gene are often related to changes in the selective regime acting on it, in such a way that a characteristic imprint is generated in its sequence. for example, if a gene is modifying its function by a process that involves positive selection, we would expect to find very rapid changes in the amino acid sequence of the encoded protein, at a rate much higher than expected under neutral evolution. on the contrary, if regions of a gene become functionally constrained, under negative selection, change in those regions will be very slow. different methods have been proposed to determine how selection acts on biological sequences. several of them have been devised to compare the synonymous and non-synonymous rates of change in coding regions, often with the purpose of testing whether positive selection has been acting upon those regions. however, relative values of ka/ks > that are strong evidence for positive selection – being higher than the expected rates under the null hypothesis of neutral evolution – are rarely found. most frequently, it is determined that sequences are under strong negative selection. for example, in a recent work in which orthologous genes were compared in human and mouse, it was estimated that, in average, ka/ks =. other methods are focused on non-synonymous changes alone, and try to model whether amino acid changes occur homogeneously or, alternatively, are concentrated on particular positions or regions of the proteins. several other methods have been proposed to explore whether the rates of change of genes vary among different lineages, especially among species or groups of species. specific tools devised to analyze functional divergence between duplicate genes have been also developed. this is due to the fact that functional shifts occur especially often in association with gene duplication events. either by acquisition of new functions by one of the duplicates or by division of the functions of the original gene between the two paralogs, duplication is often expected to radically alter the selective forces acting on the duplicate genes. however, most of the approaches developed so far have some serious limitations. first, many of them are computationally cumbersome. authors often have to choose between using the simplest tools available, that are relatively insensitive, or performing more precise analyses, but with a limited number of genes. a second problem is that many of the methods devised are based on complex mathematical models of how dna sequences change, and it is often unknown how deviations from the assumptions implicit in those models may affect the conclusions obtained. in fact, there are considerable discussions in the literature on whether certain types of analyses tend to generate spurious significant results. therefore, there is a general need for tools that combine the features of being intuitive, fast and also sufficiently sensitive. in this study, we describe a new bioinformatic tool, called uvpar, which may be used to detect regional changes in constraints in the protein sequences of duplicate genes. the uvpar algorithm is a substantial refinement of an analytical strategy devised before by one of us and already successfully used for characterizing ancient functional changes in a family of atpases/atp synthases. the basis of that strategy is to determine, combining sliding-window analyses of the degree of amino acidic conservation and permutation tests, the regions of duplicate genes that have evolved at different rates in two species. in its current implementation, uvpar allows for the fully automatic analysis of a large number of protein sequences in a short period of time. as examples of its potential, we show the results of two analyses. first, we performed uvpar comparisons for hox duplicate genes in six vertebrate species. we demonstrated that some hox genes suffered in the past positively selected changes in their sequences associated to functional shifts after gene duplication events. here, we compare these previous results with the ones obtained with uvpar in order to establish how our novel method relates to other approaches. as a second example, we generated a comprehensive study of the rbr gene family that involves several hundreds of analyses, to show how the program can be used at a large scale. in summary, uvpar is a useful novel tool for studies focused on the characterization of functional shifts in proteins, most especially in cases in which the genes of interest are part of complex gene families. implementation let us consider the case in which the phylogenetic analyses of the sequences of certain genes has established that a duplication occurred, generating two paralogous genes, before two lineages of organisms split. after both lineages became separated, the two genes accumulated differences until the present day. the question that we want to tackle is whether sequence changes have accumulated differently in those genes in both lineages since their separation. this often can be visualized as a difference in rates in one or several of the genes in a phylogenetic tree. uvpar is a program specialized in analyzing these situations. it takes into account the previous knowledge of orthology/paralogy generated by these phylogenetic analyses. moreover, uvpar analyses are based on the multiple-sequence alignments used to generate the phylogenetic trees. uvpar is written in c and we have compiled versions for windows and linux operating systems. in its simplest implementation, the program uses as input a text file containing the sequences of the proteins encoded by four genes – two duplicate genes in two different species, such as a b a and b in figure – in fasta format. however, the program also accepts larger datasets, in which pairs of duplicates of multiple species are included in a single text file. in this case, the program automatically performs analyses for all the possible combinations of species in the dataset. the program interprets the sequences in the text file as the output of a multiple-sequence alignment. once imported, the dataset is filtered, in such a way that positions in which gaps are present are eliminated from all the sequences. the program then detects the sequence with smallest size and from then on, all the analyses are performed with the first n amino acidic positions of each sequence. after these two filters, uvpar progressively reads the sequences of the first two proteins and calculates a value of similarity for each amino acid position of the two sequences, blk according to a given amino acid substitution matrix. the user may choose among the following matrices: pam pam pam blosum blosum and blosum blosum is the one that we use by default. then it does the same for the two duplicates of the second species and calculates a similar value, blk finally, the program establishes the difference between the two bl values, that we named constraint value, c = blk - blk. c is thus a value obtained for "gene quartets", defined as pairs of duplicates in two different species. if both duplicates are equally different in the two species, the constraint value should be about zero for all the positions. the key of the strategy implemented in uvpar is to analyze the distribution of c values observed to determine whether some positive or negative values are significantly clustered together. to do so, the vector of c values is randomly shuffled a number t of times and the shuffled vectors are compared with the one determined for the real paralogs. the comparison is performed for windows of increasing size, from w = to w = n - for each window size, the maximum and minimum sums of c values, that we call s, are determined for both the original and the shuffled sequences. then, uvpar compares the maximum and minimum sum in the original sequence with those in the set of shuffled sequences. this set provides a distribution of maximum values against which the values obtained in the real sequence is contrasted. two hypotheses must be tested for each window size, and therefore it is convenient to use bonferroni's correction. thus, when the value in the original sequence is found beyond the top or bottom % of the values of the simulated sequences, it is considered significant and uvpar selects the corresponding window for further analyses. the second part of the analysis performed by uvpar is the comparison of all significant windows for a given gene quartet. uvpar often detects multiple windows of different sizes that are significant. in many cases those windows are nested and therefore they refer in part to the same positions in the analyzed sequences. the biological interpretation of the significant results requires establishing the windows that best explain the detected constraint changes. this problem can be tackled in different ways. in our previous study, we solved the problem of nested windows in the simplest way, just choosing the window with the largest value of s and discarding the rest. this approach was appropriate for the particular cases analyzed in that study, but it has been determined to be insufficient for complex cases found in our subsequent analyses. this is due to two reasons. first, we have found that ties in the s values often appear in significant windows of different sizes. second, we have found significant partially overlapping windows, a case that was never detected before. therefore, in uvpar we have implemented a much more precise way of dealing with multiple significant windows of different sizes. this is a major improvement on the strategy described on our previously published work the method is as follows: the program separates the significant windows into two groups, corresponding to those significant in the analyses involving minimum values and those derived from the analyses of maximum values. within each group, windows are ordered according to their sizes. the program then compares, for one of those two groups, all the significant windows, starting with the smallest one and moving progressively to larger windows. every time that a window contained in a larger one is found that has an absolute s value equal or higher than that in the larger window, the latter is eliminated. when this process has finished, the program performs the same search for all the significant windows in the other group. the goal of this iterative screening is to eliminate all the large windows that are significant simply because they contain a small highly deviant region. for example, we found cases in which short regions had such positive c values that even when negative values were added, the larger windows that contains both the highly positive "seed" and the additional negative values were still significantly positive. those larger windows are eliminated after the comparisons described. at the end of this first screening, we are left with all the nested significant windows in which an increase of length is associated with an increase of absolute s value plus all non-nested windows. then, uvpar performs a second screening with the remaining windows, this time starting with the largest significant window and comparing it with all the smaller ones. now, all windows contained in a larger one and at the same time with absolute s value lower than that found in the larger window are eliminated. the program keeps repeating this search until all windows of both groups have been analyzed. this second screening allows the elimination of all windows that can be extended to larger ones. logically, the larger windows have to be preferred if their s values are higher. after these two screenings, we are left with sets of non-nested significant windows, that may or may not partially overlap. as a final step, uvpar analyzes all the remaining windows and eliminates those with an absolute average c value – obtained by dividing the absolute s value by the window size – lower than. this last step is devised to eliminate a perverse effect that occurs in rare cases in which a few highly positive or negative c values are clustered together in one of the extremes of the vector. then, it occasionally occurs that this region not only produces, as expected, a small highly significant window of a particular sign, but also generates a very large complementary window of the opposite sign that often spans the rest of the sequences. this happens even when the c values outside of the short significant region are randomly distributed around zero or, more often, are almost all of them zero. the few cases that we have found in which this effect was present were eliminated when the cutoff value was used, and therefore we have considered it to be a logical last filter to be implemented by default in the uvpar algorithm. after all the searches are finished, the program generates an output file that details all the relevant parameters: 1) the number of times that the sequences have been shuffled to generate the distribution of s values; 2) the sequences analyzed; 3) the detailed amino acidic alignment, with all the corresponding c values; 4) all the significant windows, before the filtering process; and, 5) all the significant windows, after the filtering process is completed. thus, the algorithm can be summarized quite simply as follows: import file with the number of pairs of duplicates to be analyzed and corresponding sequences in fasta format read t value filter sequences to eliminate gaps, estimate final size n repeat           calculate c values           if == for all k) then skip gene quartet else:           for from w = to w = n -1           estimate maximum, minimum s values for the original sequences           shuffle c values, t times           generate distribution of maximum, minimum s values from the shuffled sequences           compare maximum, minimum s values for the original and the shuffled sequences, select windows according to a significance level p < %           filter significant windows, increasing size           filter significant windows, decreasing size           filter significant windows for an average c value > | |           generate output file, including: t values, sequences, alignments, c values, significant windows before and after filtering  speed and performance the software that we used in our original study was very slow, what precluded both to use large t values and the analysis of large datasets. due to several highly significant algorithmic improvements, uvpar is about times faster. this makes possible to perform the same searches several times and with different t values to determine the degree of error associated to each number of permutation tests. when we checked for the impact of t values, we determined that to obtain a reasonably precise estimation of the probability associated to a particular value in the distribution of shuffled s values, it is convenient to choose values of t larger than which we previously used due to computational limitations. we have fully repeated the analyses presented in using uvpar and, although all the conclusions of the work are correct, we have found that the error associated to the estimation of p values with t = was considerable. uvpar analyses are fast enough as to easily handle values of t = which we consider now to be the minimum number of permutation tests to be used in this type of studies. if the sequences are short enough,t may be increased to or even as examples, using t = and a standard pc computer, sequence quartets of amino acids can be analyzed in about seconds, quartets of amino acids in about minutes and quartets of long sequences of amino acids in about minutes. time only vary slightly depending on the number of significant windows found. these results mean that uvpar can be used at a large scale. as an example, the analysis of a whole family of proteins that we will detail in section, that involved more than individual uvpar searches required a computation time of only hours, about minutes per analysis. simulations to determine the reliability and sensitivity of uvpar searches to check for the number of false positives and false negative results generated by uvpar under different conditions, we explored, using simulations, a tree containing four species, with two paralogs each. we used covtree to apply different changes in the rates of a particular region of one of the paralogs in the branch previous to the split between two of the species. then, we used uvpar to search for coherent significant windows, i. e. windows that included at least part of the region in which the rate shifts were applied and that were present in all the four comparisons in which we would expect to find significant results after the rate shifts. simulations were performed using sequences with different degrees of heterogeneity of evolutionary rates among sites. heterogeneity was modeled as a gamma distribution with a variable α parameter. values of α from to α = were used, to cover the range of results detected for real proteins. results are shown in table for windows of size w =. we found only – 3% of cases in which coherent windows are present in the control simulations no matter the value of α. when the rate of evolution for the region of interest is increased, more cases are detected, and an increase in the average c is also observed. if rates are increased ten times and α values are relatively large, coherent windows are detected in 35% to 55% of the cases. negative results in these simulations are caused by random changes that obscure the effect of the rate increase. very low α values in general precluded the detection of regional shifts by uvpar. these results suggest that, in cases in which multiple sequences can be analyzed, systematic positive uvpar results most likely will correspond to real regional changes. they also indicate that uvpar losses sensitivity when α values are very small and that the program does not detect small shifts in evolutionary rates; changes must be strong to generate a significant signal for uvpar analyses. functional shifts in hox genes hox genes are critical in many basic developmental processes, most especially in the determination of cell fates along the anterior-posterior body axis. vertebrates have a large number of hox genes, generated by gene and genome duplications. in a previous work, we showed that episodes of positive selection affected hox genes in vertebrate species. more precisely, we determined that both the ancestral duplication that gave rise to the paralogous hox-a and hox-b genes and the tetraploidization that occurred in the amphibian lineage that includes the model species xenopus laevis radically altered the selective regime acting on those genes. in the first case, the duplication was followed by a period of time in which hox-a genes diversified under positive selection. in the second case, one of the hox-b genes that were produced in the xenopus ancestor by the genome duplication process also changed under positive selection. the ratios of non-synonymous vs. synonymous substitutions in the branches affected by those two selective shifts were estimated to be about eight to ten times larger that the average ratio for the rest of branches of the tree. three distinct types of analyses suggested that positive selection acted on the regulatory n-terminal region of the hox proteins, while the highly conserved homeodomain was not affected. the large shift in evolutionary rates and the relative high values of α for these sequences mean that this is a favorable case to test whether uvpar is able to generate results comparable to those found with more complex methods. figure shows the results for uvpar analyses for the six vertebrates in which both duplicates, hox-a and hox-b were analyzed in our previous work. alignments were the same used in. uvpar results are congruent with two of the main conclusions of our previous work. first, comparisons involving the xenopus hoxb7-b gene showed significant shifts when compared with hox genes of other vertebrates in out of comparisons. the exception, the comparison involving mus musculus, is mainly due to a few homoplastic residues in the mus musculus hoxb sequence which preclude the overall values to become significant. second, all significant windows were found outside of the homeodomain. it must be noted that the third main result – the rapid differentiation after the hox-a7/hox-b duplication – could not be possibly detected here because that process occurred before all lineages that can be analyzed with uvpar became separate. in summary, these uvpar results agree with our previous findings and show that uvpar is able to provide sensitive evidence for functional shifts comparable to those obtained by more complex and time-consuming mathematical approaches. uvpar analysis of rbr ubiquitin ligases we decided to demonstrate the ability of the program to handle large datasets by analyzing proteins of a family, called rbr, that we defined some time ago. rbr proteins are ubiquitin ligases with important roles in the control of protein degradation and several of them are known to be involved in human diseases. rbr proteins are characterized by having the rbr supradomain, that is composed by two ring fingers that are separated by an ibr domain. we checked for modifications in the functional constraints by analyzing the rbr supradomain of multiple rbr proteins in several model organisms. we took advantage of the recent completion by our group of a comprehensive analysis in which we generated protein alignments for rbr proteins and phylogenetic trees that allowed us to determine the orthology/paralogy relationships for all members of the family. from that analysis, we selected all the available sequences for rbr genes in seven model organisms. the estimated average value of α for these sequences is. therefore sites evolve homogeneously enough for uvpar analyses. we took into account that not all the genes are present in all the species analyzed, and some of the genes had species-specific duplications. once these peculiarities were sorted out, we performed all the possible combinations of gene quartets for the seven species. a total of analyses were generated. in cases, we found significant results. they are summarized in tables and according to, respectively, the species and genes involved. table shows that a single species, the fish danio rerio, present a number of significant results that are well above the average found for the set of species. we think that, similarly to what occurred for the xenopus hox genes shown above, this may be related to the genome duplication that happened in the evolutionary lineage in which danio emerged. if we compare the data for the different rbr genes, we observe that the significant cases appear quite homogeneously. ari is the only gene with a number of significant results that is in average higher than that of the rest of rbrs. this is in part explained by the large number of analyses in which the ari1/ari comparisons were significant. figure shows the results for those comparisons. it becomes evident examining this figure that many of the positive results are correlated. in drosophila melanogaster there are two ari genes, ari-1a and ari-1b. significant shifts in the ibr domain are observed between the pair ari-1a/ari- of drosophila melanogaster and the ari1/ari genes of all vertebrate species. also, we see important changes along the whole rbr supradomain when genes of the urochordate ciona intestinalis are compared with vertebrate genes. therefore all these results can be explained by a few major changes in particular evolutionary lineages. in favorable cases, we can establish in which lineage the changes must have occurred. for example, considering that d. melanogaster/vertebrate comparisons involving the d. melanogaster ari-1b gene are in general non-significant, we can conclude that a functional shift in the ibr domain of the d. melanogaster ari-1a gene occurred since the ari1-a/ari1-b duplication that has made ari-1a different from their vertebrate orthologs. in summary, the results shown in figure suggests that functional diversification between the rbr supradomains of the ari and ari genes has occurred just a few times and that those changes most frequently affected the ring and ibr domains. this result points towards an important role of the n-terminal domains of the rbr supradomain in acquisition of novel functional properties in some species, perhaps to act on new ubiquitination substrates, while the non-canonical ring would be less significant. in figure we show the distribution of maximum and minimum values in the shuffled sequences and the relative position in this distribution for one of the significant values observed. bell shaped, quasi-symmetrical distributions as those shown in that figure are common. they correspond to the typical extreme value distributions that are expected for large t and w values. we have however observed in some cases highly asymmetrical distributions, often with multiple peaks, when window sizes are small. discussion we have shown that the uvpar program to may be used to quickly detect modifications in the functional constraints of duplicate genes. favorable circumstances for the program are a reasonably homogenous rate of amino acid substitution among sites and an intense change in the selective regime acting on one of the paralogs. considering that under positive selection ka/ks > while in real proteins, the average value is much lower, we may expect uvpar to be particularly useful to detect positive selection, most especially associated to functional shifts after gene duplication. in respect to other related methods, that often rely on debatable a priori assumptions and difficult mathematical calculations, uvpar results are easy to obtain and highly intuitive. uvpar analyses depend on just two a priori suppositions. first, that protein changes can be measured with a matrix such as blosum second, that the shifts affect contiguous amino acids, i. e. regions of the proteins, and therefore they can be detected by examining their primary sequences, using a sliding-window approach. there is strong evidence that regional shifts in proteins sequences often occur, in genome-scale analyses. uvpar advantages are most obvious when the goal is to quickly explore, across many species, complex gene families that include several duplicate genes. this is clearly shown by the analysis of the rbr family described above. we have also demonstrated that general conclusions about what genes or species are more prone to suffer functional shifts can be deduced from the comprehensive study of gene families. the fact that many genes can be examined in a very short time opens up the possibility of performing genome-scale studies, impossible to generate with any other related method. an obvious statistical problem is that, when a large number of quartets are examined with uvpar tests, the quite permissive significance level used is expected to generate a number of false positive results. our experience, as well as the simulations shown above, suggests examining three aspects of the data to obtain conclusions not affected by that problem. these three aspect are: 1) the significance level of each positive result and the average value of c for the significant windows. in our simulations, false positive windows have low average c values. if it is deemed necessary, the analyses can be made more strict by either changing the significance level or by increasing the conventional cutoff value c > | |; 2) the logic in structural or functional terms of the results obtained; and, 3) the systematic consistency of the results for multiple species. in our opinion, internal consistency of the data is a particularly useful way to filter the data for real positives. this is clearly observed in our simulations, in which we checked for the consistency of uvpar results for multiple analyses, finding that false positive coherent results are rare. thus, when it is found, as we have shown for the ari1/ari data, that the genes of multiple species of a lineage generate the same results when compared with genes of an outgroup species, it is good evidence for a real shift to have occurred. this means that uvpar finds several times the same imprint even although the analyzed sequences are considerably distinct. after all, orthologous genes of fishes, birds and mammals are separated by hundreds of millions of years of independent evolution. the combination of uvpar results for multiple species may be often used to trace with precision the moment and the lineage in which a change in functional constraints occurred. in this study, we have shown that the duplication of hox genes in xenopus laevis and of ari- genes in drosophila melanogaster were associated with modifications in the constraints of at least one of the duplicates. these conclusions can be easily deduced when the data are presented in a favorable format, as that shown in figures and our approach may thus complement other analytical methods devised to detect selective processes acting on particular branches of a phylogenetic tree.  uvpar generates very fast and still sensitive information about modifications of the selective constraints in duplicate genes. it can be used at a large scale, to analyze complex gene families in multiple species. its speed, simplicity of use and also the fact that its mathematical assumptions are very intuitive makes uvpar an excellent tool for all groups interested in analyzing the impact of natural selection on gene sequences, especially at a genomic scale. availability uvpar is written in c. windows and linux versions of uvpar are available, free for academic users, at the following web page:. in that page, an example of how to use uvpar can also be found. non-academic users may obtain a license to use or distribute the program by contacting the corresponding author. authors' contributions va wrote the uvpar code and contributed to the improvements in the analytical strategy implemented in the program. mg and il performed the hox and rbr analyses as well as the simulations presented above. mg also contributed to the development of the analytical strategy and tested different versions of the program. the basic strategy implemented in uvpar was developed by im, who also coordinated the research, contributed to the improvements in the strategy shown here and wrote the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 click here for file acknowledgements we thank antonio marco for critical reading of the manuscript. our group is supported by grants gen2001-4851-c06- and saf2003- and grant gv04b-. j. i. l. is the recipient of a predoctoral fellowship. figures and tables μ = μ = μ = μ = ± ± ± p
we describe uvpar, a program designed to quickly test for changes in the functional constraints of duplicate genes. starting with alignments of the proteins encoded by couples of duplicate genes in two different species, uvpar detects the regions in which modifications of the functional constraints in the paralogs occurred since both species diverged. the development of fast and effective methods to detect modifications in the selective constraints of genes is therefore of great interest. the imprint of natural selection on gene sequences is often difficult to detect. sequences can be analyzed with uvpar in just a few minutes on a standard pc computer. however, many of those methods depend heavily on underlying assumptions regarding the mode of change of dna sequences and often require sophisticated mathematical treatments that made them computationally slow. we then describe a comprehensive study of the rbr family of ubiquitin ligases in which we have performed analyses involving duplicate genes in seven model species. to demonstrate the power of the program, we first show how the results obtained with uvpar compare to those based on other approaches, using data for vertebrate hox genes. these results show that uvpar can be used to generate sensitive analyses to detect changes in the selection constraints acting on paralogs. a significant increase in the number of functional shifts was observed for the species danio rerio and for the gene ariadne- 
 the imprint of natural selection on gene sequences is often difficult to detect. a plethora of methods have been devised to detect genetic changes due to selective processes. however, many of those methods depend heavily on underlying assumptions regarding the mode of change of dna sequences and often require sophisticated mathematical treatments that made them computationally slow. the development of fast and effective methods to detect modifications in the selective constraints of genes is therefore of great interest.  we describe uvpar, a program designed to quickly test for changes in the functional constraints of duplicate genes. starting with alignments of the proteins encoded by couples of duplicate genes in two different species, uvpar detects the regions in which modifications of the functional constraints in the paralogs occurred since both species diverged. sequences can be analyzed with uvpar in just a few minutes on a standard pc computer. to demonstrate the power of the program, we first show how the results obtained with uvpar compare to those based on other approaches, using data for vertebrate hox genes. we then describe a comprehensive study of the rbr family of ubiquitin ligases in which we have performed analyses involving duplicate genes in seven model species. a significant increase in the number of functional shifts was observed for the species danio rerio and for the gene ariadne-  these results show that uvpar can be used to generate sensitive analyses to detect changes in the selection constraints acting on paralogs. the high speed of the program allows its application to genome-scale analyses. 