BACKGROUND
the availability of many completely sequenced genomes provides rich material for studying protein evolution.
gene duplications, gene losses, gene acquisitions and horizontal transfer of genes make it very difficult to reconstruct the exact evolutionary history of a protein family.
a widely used approach to study such history is to find orthologs and paralogs by comparing completely sequenced genomes.
orthologs are genes derived from a single ancestral gene as a result of the speciation event, while paralogs are genes that result from gene duplication events  <cit> .
how can we resolve these evolutionary relationships for hundreds of organisms having thousands of genes using just an ordinary personal computer?
one possible approach is to create orthologous groups for each node of the evolutionary tree  <cit>  and to compare orthologous groups belonging to different nodes of the tree.
since each orthologous group is represented by a multiple alignment of protein sequences, we need a very fast algorithm to compare multiple alignments.
profiles represent a very sensitive technique to represent a protein family and to compare protein families.
they were introduced by michael gribskov and co-workers  <cit>  and they have proved to be a valuable tool for finding weak similarities between distant proteins belonging to one family or superfamily and for improving sensitivity of database searches  <cit> .
two views have been established on the composition of profiles: continuous and discrete.
the continuous model of profiles imposes no restriction on amino acids and their counts present in an alignment column, while the discrete model assumes that only certain amino acid residues can be found at a certain position in the multiple alignment of a protein family.
the efforts of several authors were directed to improve profiles by sequence weighting  <cit>  and by introducing pseudocounts  <cit> .
as it was pointed by roman tatusov and co-workers  <cit> , the most efficient method in improving the quality of profiles by adding pseudocounts is the dirichlet mixture.
this mixture is a linear combination of dirichlet distributions.
although no restrictions are imposed on frequency profiles in this model, dirichlet distributions can be considered as discrete condensation points in 20-dimensional space.
the discrete view on profiles is exemplified by the prosite database  <cit> , which states what kind of amino acids can be present in a single position of a protein signature.
existing procedures for profile-profile comparisons  <cit>  use dynamic programming and complex scoring functions between profiles like dot product to score a column from one profile aligned to a column of the other profile.
these methods were successfully used to detect weak similarities between different protein families, to recognize folds and to predict 3-d structure.
a pairwise alignment of profiles using dynamic programming requires a2*m*n operations, where a is the alphabet size , m and n are protein lengths.
fast heuristic blast-like procedures require a discrete alphabet and it is not possible to use them to compare a pair of profiles.
if we are able somehow to convert profiles to discrete values, then we can substitute an alignment column with a symbol denoting the column type.
to this end, all alignment columns from the blocks database were converted to frequency distributions and a special clustering procedure was applied to them.
it appeared that more than 98% of alignment columns in multiple alignments from standard databases of protein multiple alignments belong to just  <dig> clusters, each with its own dominating amino acid residue .
therefore, it is possible to substitute an alignment column with this dominating residue that can be also thought as the "ancestral" amino acid residue from which the alignment column evolved.
in the case we cannot assign an alignment column to its "ancestral" residue, a special symbol "x" is used, denoting a "garbage" column.
when "ancestral" sequences are obtained it is possible to find a similarity score by applying the rigorous smith-waterman algorithm  <cit> .
but it is more desirable to develop an algorithm that might be somewhat not so rigorous but fast enough to allow a researcher to perform very computationally intensive tasks, for example, to compare orthologous groups of different taxons on a personal computer in reasonable time.
implementation
clustering procedure
each alignment column from a database of multiple alignments was converted to the frequency vector using a simple formula:
f→=fi=nin
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqadegabaaabagafmozaymbasaacqgh9aqpcqggoaakcqwgmbgzdawgaawcbagaegymaedabeaakiabcycasiabdagamnaabaaaleaacqaiyagmaeqaaogaeiilawiaesojgskaeiilawiaemozay2aasbaasqaaiabikdayiabicdawaqabagccqggpaqkaeaacqwgmbgzdawgaawcbagaemyaakgabeaakiabg2da9maalaaabagaemota40aasbaasqaaiabdmgapbqabaaakeaacqwgobgtaaaaaaaa@44d6@
where ni is the number of times amino acid i occurs in the column, n is the number of sequences in the multiple alignment.
thus, we obtained the data set d of all frequency vectors computed from all alignments columns.
following shmuel pietrokovski  <cit> , we define similarity between column frequency vectors as the pearson correlation coefficient r:
r=∑i=120⋅∑i=1202⋅∑i=1202
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgybgccqggoaakcuwgmbgzgawcgaqbaiabcycasiqbdagamzaalyaagagaeiykakiaeyypa0zaasaaaeaadaaewbqaaiabcicaoiqbdagamzaafawaasbaasqaaiabdmgapbqabagccqghsislcuwgmbgzgaqegaqbaiabcmcapiabgwsixlabcicaoiqbdagamzaagawaasbaasqaaiabdmgapbqabagccqghsislcuwgmbgzgaqegagbaiabcmcapawcbagaemyaakmaeyypa0jaegymaedabagaegomaijaegimaadaniabgghildaakeaadagcaaqaamaaqahabagaeiikagiafmozaymbauaadawgaawcbagaemyaakgabeaakiabgkhitiqbdagamzaaryaafagaeiykakyaawbaasqabeaacqaiyagmaaaabagaemyaakmaeyypa0jaegymaedabagaegomaijaegimaadaniabgghildgccqghfly1daaewbqaaiabcicaoiqbdagamzaagawaasbaasqaaiabdmgapbqabagccqghsislcuwgmbgzgaqegagbaiabcmcapmaacaaaleqabagaegomaidaaaqaaiabdmgapjabg2da9iabigdaxaqaaiabikdayiabicdawaqdcqghris5aawcbeaaaaaaaa@6e12@
we have chosen a variation of the classical "k-means" procedure  <cit>  due to the simplicity of its implementation and easy analysis how it depends on its main parameter.
since we do now know the number of clusters in the data set d, we can setup a correlation coefficient threshold rthresh and then merge clusters ci and cj if the distance between them is less than rthresh.
our modified "k-means" procedure will run as follows.
<dig>  for each frequency vector x find all frequency vectors y such that r<rthresh.
these frequency vectors will form an initial cluster.
if a frequency vector is already included in an initial cluster, it is not considered as a seed for a new cluster.
compute the initial means μ <dig> ...,μk.
<dig>  assign each frequency vector x in d to the cluster ci whose mean μi is the nearest to x.
<dig>  recompute the means of all the clusters.
μi=∑x∈cix
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacqwf8oqbdawgaawcbagaemyaakgabeaakiabg2da9iabcicaoiabigdaxiabc+caviabcyha8jabdoeadnaabaaaleaacqwgpbqaaeqaaogaeiifawnaeiykakyaaabuaeaaieqacqgf4baeasqaaiab+hha4jabgigiolabdoeadnaabaaameaacqwgpbqaaeqaaawcbeqdcqghris5aaaa@4363@
<dig>  merge clusters ci and cj if r<rthresh.
<dig>  repeat steps  <dig>   <dig> and  <dig> until the number of clusters does not change.
computing frequency column clusters from blocks, cog and pfam
to find frequency column clusters in the blocks database  <cit> , each block from blocks was processed to give a new multiple alignment.
this was done to balance the set of sequences in the alignment, so that closely homologous sequences were removed.
the processing procedure scanned the sequences in each block from the first sequence to the last.
the first sequence from the block was always included in the new multiple alignments.
the ith  sequence was included in the new alignment only if it was less than 65% identical to all previously included sequences.
thus, all sequences in the resulted multiple alignment were less than 65% identical with each other.
after the extraction procedure, only alignments having at least  <dig> sequences were considered for to further processing to have enough statistical material.
in total,  <dig> alignment columns were obtained.
the average number of sequences in a multiple alignment was  <dig>  the average protein identity in a multiple alignment was 27%.
fig.
<dig> shows the number of clusters as dependent on the correlation coefficient threshold r*.
the plot shows that when r*<  <dig>  the number of clusters is almost constant, whereas after  <dig>  the number of clusters sharply increases.
in the interval  <dig> <r*< <dig> , the number of clusters equals  <dig>  additional file  <dig> contains  <dig> average frequency vectors corresponding to each frequency cluster.
each cluster is dominated by a single amino acid, suggesting that columns belonging to one cluster are result from the evolution of just one ancestral amino acid.
to test that this dependence of the number of clusters on the correlation coefficient threshold is the general feature of multiple alignments we applied the same clustering procedure to the cog  <cit>  and pfam   <cit>  databases of multiple alignments of protein families.
in these cases we did not do any filtering of sequences in multiple alignments but we only weighted sequences using the position-based method of sequence weighing  <cit> .
we obtained the same dependence of the number of clusters on the correlation coefficient threshold with flat regions corresponding to  <dig> clusters when the correlation coefficient threshold varies from  <dig>  to  <dig>  and slow decrease when it is less than  <dig>  and abrupt rise when it is greater than  <dig>  .
cluster averages were very similar to those shown in additional file  <dig> with one dominating amino acid in each cluster, though diagonal values were different reflecting the average identity of protein sequences in multiple alignments of different databases .
any clustering procedure will also attempt to create a hierarchy of clusters.
it is interesting to know if there are any other subclusters inside the  <dig> clusters obtained by our procedure.
to answer this question, we applied our procedure to all columns in the blocks database that belong to the alanine cluster.
fig.
<dig> shows that the alanine cluster does not contain any subtypes, since when the correlation coefficient threshold is less than  <dig>  we obtain only one cluster, and when it exceeds  <dig>  we have a sharp rise in the number of clusters.
similar curves were obtained for all other  <dig> clusters .
to evaluate the size of clusters on the simplex  in the 20-dimensional space, we generated  <dig> random frequency vectors and counted how many of them fell inside our  <dig> reference clusters.
a frequency vector was considered to be inside a cluster if the correlation coefficient between it and the cluster average vector exceeded  <dig> .
only 4% of frequency vectors were inside our  <dig> clusters, proving the fact that the  <dig> clusters occupy a negligible part of the total 20-dimensional frequency vector space.
to evaluate how columns in standard multiple alignments are covered by clusters from additional file  <dig>  we analyzed all  <dig> columns available in the cog database by converting them to frequency vectors and finding the nearest cluster for each column.
if the correlation coefficient between the column frequency vector and average the nearest cluster exceeded  <dig> , the column was considered to belong to the cluster.
only  <dig> columns did not belong to any cluster, less than 2% of the total number of columns.
random test
to support the idea of an ancestral residue, the set of artificial alignments was generated in the following way.
parameters of the generation process were chosen to closely mimic parameters of the clustering procedure done on the blocks database.
<dig> random protein sequences each having the length of  <dig> amino acids were generated with the standard background frequencies of amino acids.
this produced  <dig> alignment columns, which is close to the number of alignment columns used in the clustering procedure done on blocks .
each generated sequence was considered to be an ancestral sequence, and the process of column evolution in each column of the generated sequence was modeled by random generation of  <dig> amino acid residues to obtain the same depth of the generated multiple alignment as the average depth of multiple alignments in our clustering procedure done on blocks.
each of these  <dig> amino acid residues was generated from its "ancestral" residue with probabilities determined by  <dig> transition probability matrices  <cit> , corresponding to blosum <dig> , blosum <dig>  blosum <dig> , blosum <dig> .
one of these  <dig> matrices was randomly assigned to a generated column with the following probabilities  <dig>  ,  <dig>  ,  <dig>   and  <dig>  .
these four probabilities were experimentally found values to simulate column evolution as closely as possible as it exists in blocks.
the average protein identity in thus generated protein sequences was 32%, which was close to the value obtained in the previous procedure .
the same clustering procedure was applied to randomly generated frequency vectors, and similar  <dig> clusters were obtained.
the pearson correlation coefficients between corresponding cluster averages from both groups of clusters is > <dig> .
this proves that the structure of clusters in protein multiple alignments and in the random test is basically the same.
fig.
<dig> shows the histogram of frequencies of alanine, serine and leucine in the alanine cluster from blocks and from our random simulation.
these histograms run pretty close one against each other proving that our random procedure can essentially reproduce the column evolution that exists in nature.
assignment of an alignment column to a frequency cluster
to assign columns to clusters, we weight sequences in a multiple alignment to eliminate any distortions in frequency vectors that could be caused by unequal representation of similar sequences in a multiple alignment.
then for all columns frequency vectors are calculated and the correlation coefficient is found between a frequency vector and the nearest cluster average.
if this correlation coefficient exceeds a threshold value we assign a column to its nearest cluster, otherwise it is replaced with the "x" symbol.
how phog-blast works
to find a similarity score between two multiple alignments, we convert both multiple alignments to "ancestral" sequences as it was described in the previous section.
then our method works in a very similar way as the blast  <cit>  and fasta  <cit>  algorithms.
<dig>  hash all l-tuples from the query "ancestral" sequence to the hash table.
<dig>  lookup all l-tuples from the "ancestral" sequence from the database in the hash table.
for each l-tuple hit increment the number of hits at a particular diagonal.
<dig>  leave only diagonals with nmin or more l-tuple hits for further processing.
if no such diagonals are found, declare both compared multiple alignments as non-similar and quit the procedure.
<dig>  at diagonals, extend the leftmost l-tuple hit leftward until the total score becomes negative.
memorize the maximum score achieved during this extension and the left end of the segment corresponding to this score.
in the same way, extend this l-tuple rightward.
the result of this extension is a maximal segment pair .
choose the l-tuple which is the left most to the resulting msp if such l-tuple can be found.
apply the same extension process to this l-tuple.
scan the diagonal from left to right until no l-tuple hits are found anymore.
thus, a diagonal with nmin or more l-tuple hits can give one or more msps.
for further processing, leave only msps exceeding some diagonal threshold score scorethresh.
if no such msps are found, declare both compared multiple alignments as non-similar and quit the procedure.
<dig>  order msps by their left ends belonging to one member of the pair.
<dig>  using the following recurrent relationship, obtain the score of the maximal scoring chain of msps:
score1=s1scorei+1=max⁡j,]:j≤i,lj≥ri+ <dig> mj≥pi+1)
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqaaegabaaabagaem4uamlaem4qamkaem4ta8kaemouailaemyrau0aasbaasqaaiabigdaxaqabagccqgh9aqpcqwgzbwcdawgaawcbagaegymaedabeaaaoqaaiabdofatjabdoeadjabd+eapjabdkfasjabdweafnaabaaaleaacqwgpbqacqghrawkcqaixaqmaeqaaogaeyypa0zaacbeaeaacyggtbqbcqgghbqycqgg4baeasqaaiabdqgaqbqabagcdaqadaqaaiabdohaznaabaaaleaacqwgpbqacqghrawkcqaixaqmaeqaaogaey4kasiaem4uamlaem4qamkaem4ta8kaemouailaemyraukaemoaaomaeyoei0iaem4rackaeyoei0iaem4zacmaeyyxictagiyba0maeiyyaemaeiieagnaei4waslaeiikagiaemibaw2aasbaasqaaiabdqgaqbqabagccqghsislcqwgybgcdawgaawcbagaemyaakmaey4kasiaegymaedabeaakiabcmcapiabcycasiabcicaoiabd2gatnaabaaaleaacqwgqbgaaeqaaogaeyoei0iaemicaa3aasbaasqaaiabdmgapjabgucariabigdaxaqabagccqggpaqkcqggdbqxcqgg6agocqwgqbgacqghkjyocqwgpbqacqggsaalcqwgsbabdawgaawcbagaemoaaogabeaakiabgwmizkabdkhaynaabaaaleaacqwgpbqacqghrawkcqaixaqmaeqaaogaeiilawiaemyba02aasbaasqaaiabdqgaqbqabagccqghljyscqwgwbacdawgaawcbagaemyaakmaey4kasiaegymaedabeaaaogaayjkaiaawmcaaaaaaaa@9248@
where nmsp is the total number of msps, si is the score of the ith msp, li and ri are msps' left and right ends in the first sequence, mi and pi are msps' left and right ends in the second sequence, g is the fixed gap penalty for starting a gap, g· max  is the variable gap extension penalty.
following the above description of the algorithm, we provide a summary of all parameters that can influence the phog-blast score between any two multiple alignments:
profile clusters.
substitution matrix between "ancestral" sequences.
correlation coefficient threshold r*.
tuple size l.
minimum number of tuples on one diagonal nmin.
msp score threshold scorethresh.
gap starting penalty g.
gap continuation penalty g.
selection of parameters
extensive computational experiments with phog-blast allowed us to determine parameter values that provide the best phog-blast performance.
for some parameters, we used traditional values.
for other parameters, we had to carry out a number of computational experiments to determine their best values.
to determine the best value for a parameter, we varied it over a certain interval, with all other parameter values remaining constant, and monitored the phog-blast performance.
profile clusters  were used when converting frequency vectors to "ancestral" sequences.
we used the well-known blosum <dig> substitution matrix to compare amino acid residues of "ancestral" sequences.
our experiments with phog-blast showed that the correlation coefficient threshold r* that determines whether an alignment column should belong to a cluster is  <dig> .
it is this correlation coefficient that determines the lower cutoff point on curves showing the dependence of the number clusters on the correlation coefficient threshold employed in the clustering procedure .
remarkably, it is the same for all three databases – blocks, pfam and cog, while the upper cutoff point varies in the range from  <dig>  to  <dig> .
with the correlation coefficient threshold  <dig> , only less than 2% of alignment columns in blocks, pfam and cog are out of clusters.
computational experiments with phog-blast have also allowed us to determine the best values for other parameters of the algorithm: nmin is from  <dig> to  <dig>  l =  <dig> and scorethresh is from  <dig> to  <dig>  3-tuples are probably not conserved in protein families from distant organisms, and for l =  <dig> there is a lot of non-homologous diagonals, so it is difficult to choose nmin.
the conservation of one to three 2-tuples on one diagonal in protein families is an empirical observation that we can use in building the efficient and fast algorithm.
for gap penalties, traditional values for them work very well: g =  <dig> and g =  <dig>
RESULTS
to test the phog-blast procedure, we compared it with psi-blast  <cit> , a well-known algorithm for fast profile search, and with compass  <cit> , one of the best rigorous profile-profile comparison methods.
the testing procedure was performed on the following databases: cog  <cit> , pfam  <cit>  and balibase <cit> .
we have randomly split each multiple alignment from each database into two subalignments and thus obtained two sets of multiple alignments .
then we applied the two similarity search algorithms to be tested and found bbhs .
ideally, both subalignments that were derived from one multiple alignment should find each other as the bbh.
the lower the number of bbhs that are not found, the better the scoring algorithm works.
we used psi-blast in the following way.
using the program formatdb, we made two databases of sequences extracting sequences from the multiple alignments from the databases databasea and databaseb.
then for each multiple alignment in databaseb we ran the program blastpgp with the option -b against the database of sequences made from databasea.
this will find the best sequence in databasea for the profile made from a multiple alignment in databaseb.
for this best sequence we found a multiple alignment in databasea, to which it belongs, and for this multiple alignment we ran the program blastpgp with the option -b against the database of sequences made from databaseb.
if the best sequence found in database b belongs to the multiple alignment, from which we started the matching procedure, then this multiple alignment belongs to the set of bbhs.
previously developed rigorous profile-profile comparison methods  <cit>  use very sophisticated models of profile representation and of course our algorithm can not compete with these methods if it is applied to detect remote homologies.
when we used phog-blast to find bbhs  while building the phog database  <cit> , bbhs with scores less than a given threshold  were discarded.
such bbhs introduce us into the area of remote homology, and in this area of research the model of profile clusters is too simple to represent the reality.
to test how phog-blast works in the twilight area, we compared phog-blast with compass.
we used compass as it was described in its documentation.
from the compass web site  <cit>  we downloaded  <dig> pfam alignments containing at least one sequence from the fssp database  <cit> .
for each such pfam alignment we found a fssp family as it was described in  <cit> .
we used compass and phog-blast to perform all-against-all comparisons between these  <dig> pfam alignments and to find bbhs.
a bbh was considered to be true positive if both pfam alignments belonged to one fssp family and false positive if both pfam alignments belonged to different fssp families.
fig.
<dig> shows the sensitivity curves for both methods.
of course, compass showed better performance in the twilight area.
however, it took compass  <dig> hours to perform the test, while phog-blast required only half an hour of computation time.
all tests were performed on  <dig>  ghz pentium  <dig> computer with  <dig> gb of ram running red hat linux  <dig> .
to test how compass and phog-blast work at smaller evolutionary distances , we used these algorithms to restore protein families in a similar fashion as we did it when compared phog-blast with psi-blast.
to do this test, we have randomly selected  <dig> pfam alignments having at least  <dig> sequences from the earlier mentioned downloaded  <dig> pfam alignments.
from each such pfam alignments, we randomly selected  <dig> sequences, and from these sequences we made two subalignments each having  <dig> sequences randomly selected from this group of  <dig> sequences.
we used compass and phog-blast to find bbhs between these two databases each having  <dig> alignments consisting of  <dig> sequences.
all tests were performed on  <dig>  ghz pentium  <dig> computer with  <dig> gb of ram running red hat linux  <dig> .
it took compass  <dig> minutes to perform the test, while phog-blast required only  <dig> minutes of computation time.
compass was able to correctly restore  <dig> bbhs with no false positives, while phog-blast restored  <dig> bbhs with no false positives.
this observation together with our experience of using phog-blast when we built the phog database  <cit>  demonstrates that phog-blast used in this setting does not create false positive bbhs and finds  <dig> % less bbhs than rigorous profile-profile methods.
discussion
we developed a new algorithm for finding similarity between protein multiple alignments by finding a way to convert multiple alignments into "ancestral" sequences.
the justification for this conversion was our clustering procedure done on frequency distributions and the random simulation of the column evolution.
our algorithm shows performance comparable with psi-blast, or better, but it runs much faster.
our results demonstrate that despite the seeming infiniteness of column space, it is possible to split it into  <dig> types and to indicate their evolutionary origin.
our "ancestral" sequences are similar to consensus sequences used for representing protein multiple alignments.
however, we think that our approach to convert a multiple alignment into a consensus sequence is more precise, since we have a clearer evolutionary base for such conversion.
this approach helped us to develop a very precise algorithm for finding similarity between multiple alignments.
in building frequency distributions we did not use pseudo-counts.
any inclusion of pseudo-counts in any form when building frequency distributions only deteriorated the performance of phog-blast.
since we found  <dig> condensation points in the column space, this contrasts with the approach taken in  <cit> , where an arbitrary decision was made that the dirichlet mixture consists of nine dirichlet densities.
this raises an intriguing question: what if we assume that each frequency distribution belongs to one of  <dig> column clusters modeled as closely as possible as dirichlet densities?
since mixture coefficients that link dirichlet densities together reflect the probabilities with which each dirichlet density occurs in the mixture, these coefficients can be put simply equal to the frequencies of  <dig> amino acids in protein databases.
our methodology of testing profile algorithms can also provide a framework for comparing different sequence weighting methods.
within our testing framework, the smaller the number of bbhs that are not found, the better a particular sequence weighting method works.
we have not found any subtypes for profile clusters.
however, alignment columns that are descendents of one ancestral amino acid residue differ in the speed of evolution.
when we convert columns to their clusters, fast evolving columns and slow evolving columns would acquire the same symbol.
we have not so far devised any way to include this additional information into our method.
CONCLUSIONS
since phog-blast deals with "ancestral" sequences, and not with profiles, it can be a useful tool for fast comparing multiple alignments of protein families.
this task frequently arises in different areas of comparative proteomics and protein evolution.
our accompanying paper  <cit>  shows how we used phog-blast in the development of the phog database and in the automatic reconstruction of orthologs and paralogs from protein complements of different species.
availability and requirements
phog-blast is freely available in java programming language from  and requires the java runtime environment.
it is also available as an additional file with this manuscript .
the phog-blast software is provided with no guarantee or warranty of any kind.
it may be distributed under the terms of the gnu general public license.
authors' contributions
im did all computation and wrote the manuscript.
am developed the approach for this research, provided overall guidance and revised this manuscript.
supplementary material
additional file 1
average frequencies of frequency column clusters obtained from blocks.
dominating amino acids are shown in bold face.
click here for file
additional file 2
phog-blast
click here for file
acknowledgements
we are grateful to s. sunyaev and m. roytberg for discussion regarding the comparison of different profile methods.
we would like to thank m. gelfand for critically reviewing this manuscript.
this study was partially supported by grants from the russian academy of sciences , the howard hughes medical institute , and the russian fund of basic research .

 the availability of many completely sequenced genomes provides rich material for studying protein evolution. gene duplications, gene losses, gene acquisitions and horizontal transfer of genes make it very difficult to reconstruct the exact evolutionary history of a protein family. a widely used approach to study such history is to find orthologs and paralogs by comparing completely sequenced genomes. orthologs are genes derived from a single ancestral gene as a result of the speciation event, while paralogs are genes that result from gene duplication events. how can we resolve these evolutionary relationships for hundreds of organisms having thousands of genes using just an ordinary personal computer? one possible approach is to create orthologous groups for each node of the evolutionary tree and to compare orthologous groups belonging to different nodes of the tree. since each orthologous group is represented by a multiple alignment of protein sequences, we need a very fast algorithm to compare multiple alignments. profiles represent a very sensitive technique to represent a protein family and to compare protein families. they were introduced by michael gribskov and co-workers and they have proved to be a valuable tool for finding weak similarities between distant proteins belonging to one family or superfamily and for improving sensitivity of database searches. two views have been established on the composition of profiles: continuous and discrete. the continuous model of profiles imposes no restriction on amino acids and their counts present in an alignment column, while the discrete model assumes that only certain amino acid residues can be found at a certain position in the multiple alignment of a protein family. the efforts of several authors were directed to improve profiles by sequence weighting and by introducing pseudocounts. as it was pointed by roman tatusov and co-workers, the most efficient method in improving the quality of profiles by adding pseudocounts is the dirichlet mixture. this mixture is a linear combination of dirichlet distributions. although no restrictions are imposed on frequency profiles in this model, dirichlet distributions can be considered as discrete condensation points in 20-dimensional space. the discrete view on profiles is exemplified by the prosite database, which states what kind of amino acids can be present in a single position of a protein signature. existing procedures for profile-profile comparisons use dynamic programming and complex scoring functions between profiles like dot product to score a column from one profile aligned to a column of the other profile. these methods were successfully used to detect weak similarities between different protein families, to recognize folds and to predict 3-d structure. a pairwise alignment of profiles using dynamic programming requires a2*m*n operations, where a is the alphabet size, m and n are protein lengths. fast heuristic blast-like procedures require a discrete alphabet and it is not possible to use them to compare a pair of profiles. if we are able somehow to convert profiles to discrete values, then we can substitute an alignment column with a symbol denoting the column type. to this end, all alignment columns from the blocks database were converted to frequency distributions and a special clustering procedure was applied to them. it appeared that more than 98% of alignment columns in multiple alignments from standard databases of protein multiple alignments belong to just clusters, each with its own dominating amino acid residue. therefore, it is possible to substitute an alignment column with this dominating residue that can be also thought as the "ancestral" amino acid residue from which the alignment column evolved. in the case we cannot assign an alignment column to its "ancestral" residue, a special symbol "x" is used, denoting a "garbage" column. when "ancestral" sequences are obtained it is possible to find a similarity score by applying the rigorous smith-waterman algorithm. but it is more desirable to develop an algorithm that might be somewhat not so rigorous but fast enough to allow a researcher to perform very computationally intensive tasks, for example, to compare orthologous groups of different taxons on a personal computer in reasonable time. implementation clustering procedure each alignment column from a database of multiple alignments was converted to the frequency vector using a simple formula: f→=fi=nin mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqadegabaaabagafmozaymbasaacqgh9aqpcqggoaakcqwgmbgzdawgaawcbagaegymaedabeaakiabcycasiabdagamnaabaaaleaacqaiyagmaeqaaogaeiilawiaesojgskaeiilawiaemozay2aasbaasqaaiabikdayiabicdawaqabagccqggpaqkaeaacqwgmbgzdawgaawcbagaemyaakgabeaakiabg2da9maalaaabagaemota40aasbaasqaaiabdmgapbqabaaakeaacqwgobgtaaaaaaaa@44d6@ where ni is the number of times amino acid i occurs in the column, n is the number of sequences in the multiple alignment. thus, we obtained the data set d of all frequency vectors computed from all alignments columns. following shmuel pietrokovski, we define similarity between column frequency vectors as the pearson correlation coefficient r: r=∑i=120⋅∑i=1202⋅∑i=1202 mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaacqwgybgccqggoaakcuwgmbgzgawcgaqbaiabcycasiqbdagamzaalyaagagaeiykakiaeyypa0zaasaaaeaadaaewbqaaiabcicaoiqbdagamzaafawaasbaasqaaiabdmgapbqabagccqghsislcuwgmbgzgaqegaqbaiabcmcapiabgwsixlabcicaoiqbdagamzaagawaasbaasqaaiabdmgapbqabagccqghsislcuwgmbgzgaqegagbaiabcmcapawcbagaemyaakmaeyypa0jaegymaedabagaegomaijaegimaadaniabgghildaakeaadagcaaqaamaaqahabagaeiikagiafmozaymbauaadawgaawcbagaemyaakgabeaakiabgkhitiqbdagamzaaryaafagaeiykakyaawbaasqabeaacqaiyagmaaaabagaemyaakmaeyypa0jaegymaedabagaegomaijaegimaadaniabgghildgccqghfly1daaewbqaaiabcicaoiqbdagamzaagawaasbaasqaaiabdmgapbqabagccqghsislcuwgmbgzgaqegagbaiabcmcapmaacaaaleqabagaegomaidaaaqaaiabdmgapjabg2da9iabigdaxaqaaiabikdayiabicdawaqdcqghris5aawcbeaaaaaaaa@6e12@ we have chosen a variation of the classical "k-means" procedure due to the simplicity of its implementation and easy analysis how it depends on its main parameter. since we do now know the number of clusters in the data set d, we can setup a correlation coefficient threshold rthresh and then merge clusters ci and cj if the distance between them is less than rthresh. our modified "k-means" procedure will run as follows. for each frequency vector x find all frequency vectors y such that r<rthresh. these frequency vectors will form an initial cluster. if a frequency vector is already included in an initial cluster, it is not considered as a seed for a new cluster. compute the initial means μ...,μk. assign each frequency vector x in d to the cluster ci whose mean μi is the nearest to x. recompute the means of all the clusters. μi=∑x∈cix mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaaiigacqwf8oqbdawgaawcbagaemyaakgabeaakiabg2da9iabcicaoiabigdaxiabc+caviabcyha8jabdoeadnaabaaaleaacqwgpbqaaeqaaogaeiifawnaeiykakyaaabuaeaaieqacqgf4baeasqaaiab+hha4jabgigiolabdoeadnaabaaameaacqwgpbqaaeqaaawcbeqdcqghris5aaaa@4363@ merge clusters ci and cj if r<rthresh. repeat steps and until the number of clusters does not change. computing frequency column clusters from blocks, cog and pfam to find frequency column clusters in the blocks database, each block from blocks was processed to give a new multiple alignment. this was done to balance the set of sequences in the alignment, so that closely homologous sequences were removed. the processing procedure scanned the sequences in each block from the first sequence to the last. the first sequence from the block was always included in the new multiple alignments. the ith sequence was included in the new alignment only if it was less than 65% identical to all previously included sequences. thus, all sequences in the resulted multiple alignment were less than 65% identical with each other. after the extraction procedure, only alignments having at least sequences were considered for to further processing to have enough statistical material. in total, alignment columns were obtained. the average number of sequences in a multiple alignment was the average protein identity in a multiple alignment was 27%. fig. shows the number of clusters as dependent on the correlation coefficient threshold r*. the plot shows that when r*< the number of clusters is almost constant, whereas after the number of clusters sharply increases. in the interval <r*<, the number of clusters equals additional file contains average frequency vectors corresponding to each frequency cluster. each cluster is dominated by a single amino acid, suggesting that columns belonging to one cluster are result from the evolution of just one ancestral amino acid. to test that this dependence of the number of clusters on the correlation coefficient threshold is the general feature of multiple alignments we applied the same clustering procedure to the cog and pfam databases of multiple alignments of protein families. in these cases we did not do any filtering of sequences in multiple alignments but we only weighted sequences using the position-based method of sequence weighing. we obtained the same dependence of the number of clusters on the correlation coefficient threshold with flat regions corresponding to clusters when the correlation coefficient threshold varies from to and slow decrease when it is less than and abrupt rise when it is greater than. cluster averages were very similar to those shown in additional file with one dominating amino acid in each cluster, though diagonal values were different reflecting the average identity of protein sequences in multiple alignments of different databases. any clustering procedure will also attempt to create a hierarchy of clusters. it is interesting to know if there are any other subclusters inside the clusters obtained by our procedure. to answer this question, we applied our procedure to all columns in the blocks database that belong to the alanine cluster. fig. shows that the alanine cluster does not contain any subtypes, since when the correlation coefficient threshold is less than we obtain only one cluster, and when it exceeds we have a sharp rise in the number of clusters. similar curves were obtained for all other clusters. to evaluate the size of clusters on the simplex in the 20-dimensional space, we generated random frequency vectors and counted how many of them fell inside our reference clusters. a frequency vector was considered to be inside a cluster if the correlation coefficient between it and the cluster average vector exceeded. only 4% of frequency vectors were inside our clusters, proving the fact that the clusters occupy a negligible part of the total 20-dimensional frequency vector space. to evaluate how columns in standard multiple alignments are covered by clusters from additional file we analyzed all columns available in the cog database by converting them to frequency vectors and finding the nearest cluster for each column. if the correlation coefficient between the column frequency vector and average the nearest cluster exceeded, the column was considered to belong to the cluster. only columns did not belong to any cluster, less than 2% of the total number of columns. random test to support the idea of an ancestral residue, the set of artificial alignments was generated in the following way. parameters of the generation process were chosen to closely mimic parameters of the clustering procedure done on the blocks database. random protein sequences each having the length of amino acids were generated with the standard background frequencies of amino acids. this produced alignment columns, which is close to the number of alignment columns used in the clustering procedure done on blocks. each generated sequence was considered to be an ancestral sequence, and the process of column evolution in each column of the generated sequence was modeled by random generation of amino acid residues to obtain the same depth of the generated multiple alignment as the average depth of multiple alignments in our clustering procedure done on blocks. each of these amino acid residues was generated from its "ancestral" residue with probabilities determined by transition probability matrices, corresponding to blosum, blosum blosum, blosum. one of these matrices was randomly assigned to a generated column with the following probabilities,, and. these four probabilities were experimentally found values to simulate column evolution as closely as possible as it exists in blocks. the average protein identity in thus generated protein sequences was 32%, which was close to the value obtained in the previous procedure. the same clustering procedure was applied to randomly generated frequency vectors, and similar clusters were obtained. the pearson correlation coefficients between corresponding cluster averages from both groups of clusters is >. this proves that the structure of clusters in protein multiple alignments and in the random test is basically the same. fig. shows the histogram of frequencies of alanine, serine and leucine in the alanine cluster from blocks and from our random simulation. these histograms run pretty close one against each other proving that our random procedure can essentially reproduce the column evolution that exists in nature. assignment of an alignment column to a frequency cluster to assign columns to clusters, we weight sequences in a multiple alignment to eliminate any distortions in frequency vectors that could be caused by unequal representation of similar sequences in a multiple alignment. then for all columns frequency vectors are calculated and the correlation coefficient is found between a frequency vector and the nearest cluster average. if this correlation coefficient exceeds a threshold value we assign a column to its nearest cluster, otherwise it is replaced with the "x" symbol. how phog-blast works to find a similarity score between two multiple alignments, we convert both multiple alignments to "ancestral" sequences as it was described in the previous section. then our method works in a very similar way as the blast and fasta algorithms. hash all l-tuples from the query "ancestral" sequence to the hash table. lookup all l-tuples from the "ancestral" sequence from the database in the hash table. for each l-tuple hit increment the number of hits at a particular diagonal. leave only diagonals with nmin or more l-tuple hits for further processing. if no such diagonals are found, declare both compared multiple alignments as non-similar and quit the procedure. at diagonals, extend the leftmost l-tuple hit leftward until the total score becomes negative. memorize the maximum score achieved during this extension and the left end of the segment corresponding to this score. in the same way, extend this l-tuple rightward. the result of this extension is a maximal segment pair. choose the l-tuple which is the left most to the resulting msp if such l-tuple can be found. apply the same extension process to this l-tuple. scan the diagonal from left to right until no l-tuple hits are found anymore. thus, a diagonal with nmin or more l-tuple hits can give one or more msps. for further processing, leave only msps exceeding some diagonal threshold score scorethresh. if no such msps are found, declare both compared multiple alignments as non-similar and quit the procedure. order msps by their left ends belonging to one member of the pair. using the following recurrent relationship, obtain the score of the maximal scoring chain of msps: score1=s1scorei+1=max⁡j,]:j≤i,lj≥ri+ mj≥pi+1) mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaafaqaaegabaaabagaem4uamlaem4qamkaem4ta8kaemouailaemyrau0aasbaasqaaiabigdaxaqabagccqgh9aqpcqwgzbwcdawgaawcbagaegymaedabeaaaoqaaiabdofatjabdoeadjabd+eapjabdkfasjabdweafnaabaaaleaacqwgpbqacqghrawkcqaixaqmaeqaaogaeyypa0zaacbeaeaacyggtbqbcqgghbqycqgg4baeasqaaiabdqgaqbqabagcdaqadaqaaiabdohaznaabaaaleaacqwgpbqacqghrawkcqaixaqmaeqaaogaey4kasiaem4uamlaem4qamkaem4ta8kaemouailaemyraukaemoaaomaeyoei0iaem4rackaeyoei0iaem4zacmaeyyxictagiyba0maeiyyaemaeiieagnaei4waslaeiikagiaemibaw2aasbaasqaaiabdqgaqbqabagccqghsislcqwgybgcdawgaawcbagaemyaakmaey4kasiaegymaedabeaakiabcmcapiabcycasiabcicaoiabd2gatnaabaaaleaacqwgqbgaaeqaaogaeyoei0iaemicaa3aasbaasqaaiabdmgapjabgucariabigdaxaqabagccqggpaqkcqggdbqxcqgg6agocqwgqbgacqghkjyocqwgpbqacqggsaalcqwgsbabdawgaawcbagaemoaaogabeaakiabgwmizkabdkhaynaabaaaleaacqwgpbqacqghrawkcqaixaqmaeqaaogaeiilawiaemyba02aasbaasqaaiabdqgaqbqabagccqghljyscqwgwbacdawgaawcbagaemyaakmaey4kasiaegymaedabeaaaogaayjkaiaawmcaaaaaaaa@9248@ where nmsp is the total number of msps, si is the score of the ith msp, li and ri are msps' left and right ends in the first sequence, mi and pi are msps' left and right ends in the second sequence, g is the fixed gap penalty for starting a gap, g· max is the variable gap extension penalty. following the above description of the algorithm, we provide a summary of all parameters that can influence the phog-blast score between any two multiple alignments: profile clusters. substitution matrix between "ancestral" sequences. correlation coefficient threshold r*. tuple size l. minimum number of tuples on one diagonal nmin. msp score threshold scorethresh. gap starting penalty g. gap continuation penalty g. selection of parameters extensive computational experiments with phog-blast allowed us to determine parameter values that provide the best phog-blast performance. for some parameters, we used traditional values. for other parameters, we had to carry out a number of computational experiments to determine their best values. to determine the best value for a parameter, we varied it over a certain interval, with all other parameter values remaining constant, and monitored the phog-blast performance. profile clusters were used when converting frequency vectors to "ancestral" sequences. we used the well-known blosum substitution matrix to compare amino acid residues of "ancestral" sequences. our experiments with phog-blast showed that the correlation coefficient threshold r* that determines whether an alignment column should belong to a cluster is. it is this correlation coefficient that determines the lower cutoff point on curves showing the dependence of the number clusters on the correlation coefficient threshold employed in the clustering procedure. remarkably, it is the same for all three databases – blocks, pfam and cog, while the upper cutoff point varies in the range from to. with the correlation coefficient threshold, only less than 2% of alignment columns in blocks, pfam and cog are out of clusters. computational experiments with phog-blast have also allowed us to determine the best values for other parameters of the algorithm: nmin is from to l = and scorethresh is from to 3-tuples are probably not conserved in protein families from distant organisms, and for l = there is a lot of non-homologous diagonals, so it is difficult to choose nmin. the conservation of one to three 2-tuples on one diagonal in protein families is an empirical observation that we can use in building the efficient and fast algorithm. for gap penalties, traditional values for them work very well: g = and g =  to test the phog-blast procedure, we compared it with psi-blast, a well-known algorithm for fast profile search, and with compass, one of the best rigorous profile-profile comparison methods. the testing procedure was performed on the following databases: cog, pfam and balibase. we have randomly split each multiple alignment from each database into two subalignments and thus obtained two sets of multiple alignments. then we applied the two similarity search algorithms to be tested and found bbhs. ideally, both subalignments that were derived from one multiple alignment should find each other as the bbh. the lower the number of bbhs that are not found, the better the scoring algorithm works. we used psi-blast in the following way. using the program formatdb, we made two databases of sequences extracting sequences from the multiple alignments from the databases databasea and databaseb. then for each multiple alignment in databaseb we ran the program blastpgp with the option -b against the database of sequences made from databasea. this will find the best sequence in databasea for the profile made from a multiple alignment in databaseb. for this best sequence we found a multiple alignment in databasea, to which it belongs, and for this multiple alignment we ran the program blastpgp with the option -b against the database of sequences made from databaseb. if the best sequence found in database b belongs to the multiple alignment, from which we started the matching procedure, then this multiple alignment belongs to the set of bbhs. previously developed rigorous profile-profile comparison methods use very sophisticated models of profile representation and of course our algorithm can not compete with these methods if it is applied to detect remote homologies. when we used phog-blast to find bbhs while building the phog database, bbhs with scores less than a given threshold were discarded. such bbhs introduce us into the area of remote homology, and in this area of research the model of profile clusters is too simple to represent the reality. to test how phog-blast works in the twilight area, we compared phog-blast with compass. we used compass as it was described in its documentation. from the compass web site we downloaded pfam alignments containing at least one sequence from the fssp database. for each such pfam alignment we found a fssp family as it was described in. we used compass and phog-blast to perform all-against-all comparisons between these pfam alignments and to find bbhs. a bbh was considered to be true positive if both pfam alignments belonged to one fssp family and false positive if both pfam alignments belonged to different fssp families. fig. shows the sensitivity curves for both methods. of course, compass showed better performance in the twilight area. however, it took compass hours to perform the test, while phog-blast required only half an hour of computation time. all tests were performed on ghz pentium computer with gb of ram running red hat linux. to test how compass and phog-blast work at smaller evolutionary distances, we used these algorithms to restore protein families in a similar fashion as we did it when compared phog-blast with psi-blast. to do this test, we have randomly selected pfam alignments having at least sequences from the earlier mentioned downloaded pfam alignments. from each such pfam alignments, we randomly selected sequences, and from these sequences we made two subalignments each having sequences randomly selected from this group of sequences. we used compass and phog-blast to find bbhs between these two databases each having alignments consisting of sequences. all tests were performed on ghz pentium computer with gb of ram running red hat linux. it took compass minutes to perform the test, while phog-blast required only minutes of computation time. compass was able to correctly restore bbhs with no false positives, while phog-blast restored bbhs with no false positives. this observation together with our experience of using phog-blast when we built the phog database demonstrates that phog-blast used in this setting does not create false positive bbhs and finds % less bbhs than rigorous profile-profile methods. discussion we developed a new algorithm for finding similarity between protein multiple alignments by finding a way to convert multiple alignments into "ancestral" sequences. the justification for this conversion was our clustering procedure done on frequency distributions and the random simulation of the column evolution. our algorithm shows performance comparable with psi-blast, or better, but it runs much faster. our results demonstrate that despite the seeming infiniteness of column space, it is possible to split it into types and to indicate their evolutionary origin. our "ancestral" sequences are similar to consensus sequences used for representing protein multiple alignments. however, we think that our approach to convert a multiple alignment into a consensus sequence is more precise, since we have a clearer evolutionary base for such conversion. this approach helped us to develop a very precise algorithm for finding similarity between multiple alignments. in building frequency distributions we did not use pseudo-counts. any inclusion of pseudo-counts in any form when building frequency distributions only deteriorated the performance of phog-blast. since we found condensation points in the column space, this contrasts with the approach taken in, where an arbitrary decision was made that the dirichlet mixture consists of nine dirichlet densities. this raises an intriguing question: what if we assume that each frequency distribution belongs to one of column clusters modeled as closely as possible as dirichlet densities? since mixture coefficients that link dirichlet densities together reflect the probabilities with which each dirichlet density occurs in the mixture, these coefficients can be put simply equal to the frequencies of amino acids in protein databases. our methodology of testing profile algorithms can also provide a framework for comparing different sequence weighting methods. within our testing framework, the smaller the number of bbhs that are not found, the better a particular sequence weighting method works. we have not found any subtypes for profile clusters. however, alignment columns that are descendents of one ancestral amino acid residue differ in the speed of evolution. when we convert columns to their clusters, fast evolving columns and slow evolving columns would acquire the same symbol. we have not so far devised any way to include this additional information into our method.  since phog-blast deals with "ancestral" sequences, and not with profiles, it can be a useful tool for fast comparing multiple alignments of protein families. this task frequently arises in different areas of comparative proteomics and protein evolution. our accompanying paper shows how we used phog-blast in the development of the phog database and in the automatic reconstruction of orthologs and paralogs from protein complements of different species. availability and requirements phog-blast is freely available in java programming language from and requires the java runtime environment. it is also available as an additional file with this manuscript. the phog-blast software is provided with no guarantee or warranty of any kind. it may be distributed under the terms of the gnu general public license. authors' contributions im did all computation and wrote the manuscript. am developed the approach for this research, provided overall guidance and revised this manuscript. supplementary material additional file 1 average frequencies of frequency column clusters obtained from blocks. dominating amino acids are shown in bold face. click here for file additional file 2 phog-blast click here for file acknowledgements we are grateful to s. sunyaev and m. roytberg for discussion regarding the comparison of different profile methods. we would like to thank m. gelfand for critically reviewing this manuscript. this study was partially supported by grants from the russian academy of sciences, the howard hughes medical institute, and the russian fund of basic research.
we developed a new algorithm called phog-blast for fast similarity search of profiles. since phog-blast can compare multiple alignments of protein families, it can be used in different areas of comparative proteomics and protein evolution. for example, phog-blast helped to build the phog database of phylogenetic orthologous groups. an essential step in building this database was comparing protein complements of different species and orthologous groups of different taxons on a personal computer in reasonable time. the need to compare protein profiles frequently arises in various protein research areas: comparison of protein families, domain searches, resolution of orthology and paralogy. when it is applied to detect weak similarity between protein families, phog-blast is less precise than rigorous profile-profile comparison method, though it runs much faster and can be used as a hit pre-selecting tool. phog-blast required 10– times less computer memory and computation time than psi-blast. to determine the optimal alphabet, we analyzed columns in reliable multiple alignments and obtained column clusters in the 20-dimensional profile space by applying a special clustering procedure. with these clusters, only less than 2% of columns in multiple alignments are out of clusters. we tested the performance of phog-blast vs. psi-blast on three well-known databases of multiple alignments: cog, pfam and balibase. on the cog database both algorithms showed the same performance, on pfam and balibase phog-blast was much superior to psi-blast.
 the need to compare protein profiles frequently arises in various protein research areas: comparison of protein families, domain searches, resolution of orthology and paralogy. the existing fast algorithms can only compare a protein sequence with a protein sequence and a profile with a sequence. algorithms to compare profiles use dynamic programming and complex scoring functions.  we developed a new algorithm called phog-blast for fast similarity search of profiles. this algorithm uses profile discretization to convert a profile to a finite alphabet and utilizes hashing for fast search. to determine the optimal alphabet, we analyzed columns in reliable multiple alignments and obtained column clusters in the 20-dimensional profile space by applying a special clustering procedure. we show that the clustering procedure works best if its parameters are chosen so that profile clusters are obtained which can be interpreted as ancestral amino acid residues. with these clusters, only less than 2% of columns in multiple alignments are out of clusters. we tested the performance of phog-blast vs. psi-blast on three well-known databases of multiple alignments: cog, pfam and balibase. on the cog database both algorithms showed the same performance, on pfam and balibase phog-blast was much superior to psi-blast. phog-blast required 10– times less computer memory and computation time than psi-blast.  since phog-blast can compare multiple alignments of protein families, it can be used in different areas of comparative proteomics and protein evolution. for example, phog-blast helped to build the phog database of phylogenetic orthologous groups. an essential step in building this database was comparing protein complements of different species and orthologous groups of different taxons on a personal computer in reasonable time. when it is applied to detect weak similarity between protein families, phog-blast is less precise than rigorous profile-profile comparison method, though it runs much faster and can be used as a hit pre-selecting tool. 