BACKGROUND
the "in silico" functional annotation of proteins generated by large scale sequencing projects is an important challenge in biology.
here we propose a rigorous protocol for multigenic families' annotation.
1/importance of phylogenetic reconstruction
because important functions are conserved during evolution, the first step in analysis is to determine homologous sequences.
more specifically, orthologs are more likely to share the same function while paralogs can undergo functional shifts.
developments and improvements in database similarity search programs such as blast  <cit>  allowed the rapid identification of potential homologous sequences and therefore allowed functional prediction of several thousand of genes and proteins present in databases.
however, the closest blast is often not the nearest neighbor  <cit> .
indeed, similarity-based approaches do not consider all the information from comparative and evolutionary biology.
they do not differentiate between orthologs and paralogs among homologs.
so, phylogenetic approaches, taking into account duplication and speciation events are necessary to robustly produce functional annotation of new, uncharacterized proteins  <cit> .
2/necessity to enlarge databases
for proteins' phylogenetic reconstruction, protein databases containing the proteomes of completely sequenced species along with individually submitted protein sequences are usually used   <cit> .
yet, the vast majority of species are not fully sequenced and most of their protein sequences are still unknown.
however, a lot of transcriptional information is carried by growing gene expression databases, concerning normal or pathological tissues   <cit> .
these mrnas could be used for  reconstruction of unknown proteins in "not yet sequenced" species.
in parallel, translation of est contigs can be used to enlarge the spectrum of species containing homologs when one analyses a protein family.
3/importance of expression patterns' determination for an accurate annotation
it should be noted that phylogenetic analysis can only give information at the biochemical function level.
furthermore, while orthologs can have very similar "molecular function", they can exhibit different "macroscopic functions", due to a transcriptional shift for example.
to produce an accurate proteic functional annotation, one must have complete sequence information given by phylogenetic reconstruction, with expression patterns analysis.
this is the second reason why using data from expression databases is interesting.
analysis of expression divergence between paralogs and orthologs have been recently published in human and mouse.
it appears that gene expression profiles diverge between paralogs.
orthologs can diverge in their expression pattern too  <cit> .
moreover, orthologs that have undergone recent duplication have less strongly correlated expression profiles than those that have not  <cit> .
up to now, there is no study examining expression divergence between homologs that takes into account a broad spectrum of species and after a phylogenetic reconstruction.
4/our approach
we present here a new way to functionally annotate proteins in silico, taking into account all these concepts.
in a first step, we reconstructed the phylogeny of a protein family, using an enlarged database containing 1/full length proteins from ncbi nr protein database  <cit>  and 2/translation of est contigs from ncbi dbest database  <cit> .
we used a new software platform, figenix  <cit> , adapted to this kind of phylogenomic reconstruction.
in a second step, we created an automated pipeline to couple these phylogenetic reconstructions with expression pattern data.
we then compared the evolution rate of the different paralogous sequences with their patterns of expression.
we validated this approach with the peptidyl arginine deiminase multigenic family.
these genes encode peptidyl arginine deiminase proteins  that are implied in rheumatoid arthritis , a systemic autoimmune disease that primarily manifests as a chronic symmetric polyarthritis which gradually destroys articular cartilage and bone  <cit> .
pads transform arginine residues into citrulline in a calcium-dependant manner  <cit> .
it was recently shown by a case-control linkage disequilibrium study that pad type  <dig> is a susceptibility locus for rheumatoid arthritis in a japanese population  <cit> .
these data were not retrieved in a caucasian population  <cit> .
up to now, no phylogenomic annotation coupled with expression data of this protein family has been proposed.
we believe this analysis will shed new light on this incompletely characterized gene family.
RESULTS
1/phylogenetic reconstruction
first, we built a phylogenetic tree with full length pad proteins  .
we found five paralogs of pad proteins  encoded by five different genes as previously described  <cit> .
we found pad members in birds , amphibians  and teleosteans .
this tree supported by high bootstrap values suggested the following steps for evolution of pad genes in craniates: pads' common ancestor was present as far back as in the last common ancestor of teleostean and mammals.
pad- <dig> first diverged by duplication from the common ancestor before the radiation of mammals.
we next built a phylogeny with the same protein dataset completed by additional proteins built from translations of est contigs.
first, blast searches on est databases allowed us to identify other species in which pad proteins are likely present.
we found pads in additional teleosteans  and in cephalochordates .
because ests are not full length sequences, they can not be used in this raw state to reconstruct a phylogenetic tree.
to correct for this, we built larger sequences by merging ests relative to the same gene into contigs.
est contigs were translated and included in a local database and used for phylogenetic reconstruction along with protein from ncbi nr database.
because ests contigs often do not correspond to full length transcripts, we had to build several trees with alignment only on parts of the protein .
a tree built for each domain is available in additional file online .
some ests could not be integrated into a contig  and were used in this state as sequence queries to reconstruct a phylogeny.
an example of a tree built with a singleton  is shown in figure  <dig>
all the trees built with partial sequences were validated only if they were in agreement with the topology of the reference tree .
this allowed us to classify all the contigs and singletons into paralogy or orthology groups.
this analysis allowed us to propose a more detailed evolutionary story for the pad genes: the common ancestor of the pad gene was at least present in the common ancestor of the euchordates  .
it means that the pad gene present in amphioxus is co-orthologous to all pads .
analysis of trees' topologies allowed us to conclude that duplications of pads certainly occurred after the speciation of craniata or more probably after the speciation of sarcopterygian, but before the mammals' radiation.
so, phylogenetic trees including translated est contigs allowed us to enlarge the dataset of species and to classify these "new proteins" into a paralogy group.
when looking at the branches lengths, which are directly correlated with the sequence evolution rate, it appeared that some members of pad genes evolved faster than the others.
branches lengths were longer within the pad- <dig> paralogy group  whereas branches lengths were the shortest within the pad- <dig> paralogy group .
it simply means that the similarity between species is higher for pad- <dig> than it is for pad- <dig>  this indicates that pad-2s' group is the most highly conserved.
furthermore, the distance from the common ancestor is lower for pad-2s' group than it is for the others, suggesting that pad-2's group has been submitted to high negative selective pressure and has probably retained the ancestral biochemical function.
2/coupling expression pattern with phylogenetic classification
the phylogenetic reconstruction allowed us to classify all the contigs or singletons in phylogenetic groups for the pads phylogeny.
we then developed software agents to extract transcriptional data  from the ncbi dbest database  <cit> .
we first normalized expression data , because expression data are not normalized in ncbi dbest.
we then built a table classifying ests according to their position in the tree  and to the tissue in which they are expressed.
results are illustrated in a simplified figure  <dig>  raw data are in table  <dig> .
given the current est coverage, we focused on human and mouse.
the bias could occur when, say, the coverage of the zebrafish is much lower than human or mouse.
est data are all represented in the table, but statistical analysis took into account only human and mouse data.
for homo sapiens and mus musculus, data are congruent with expression patterns data available on unigene  <cit> .
coupling phylogenetic classification with expression data permitted us to perform updated footprints of the transcriptional pattern for each paralogs group in this multigenic family.
one group of paralogs, pad- <dig>  showed a broad tissular expression.
as the expression pattern was different across paralogous pad genes, we could only compare a few tissue types for differential levels of expression.
the statistical test described by audic and claverie  <cit>  was used.
significant differences  were shown in thymus, whole blood, lung, inner ear, mixed liver and spleen library and leucocytes .
3/coupling sequence evolution rate and expression pattern
when observing tissue distribution, we noticed that it was conversely correlated  to sequence evolution rate  .
this was remarkable when we compared branches lengths and tissue distribution between pad- <dig>  and pad- <dig>  .
discussion
1/using phylogenomics and data from est databases
using est data to complete data from protein databases is commonly performed.
yet, if one uses only the best blast hits approach to classify ests data, there is a large risk of misclassification into paralogy or orthology groups.
thus, expression data for a given paralogy group could be undermined by the wrong assignment of est data into a group.
in our example, we found additional expressed sequences from teleostean and cephalochordates in est databases.
when performing analysis with reciprocal best blast hit on nr protein database; est from amphioxus  and teleostean is classified as an ortholog of pad- <dig> only, but not as a co-ortholog of all pads as is shown by phylogenetic analysis.
this illustrates why phylogenetic reconstruction is essential for accurate est analysis and able to better avoid misclassification.
our approach enlarges the dataset of homologs-containing species and robustly and accurately classifies multigenic family members into paralogy and orthology group.
comparison of our results to those of chavanas et al.
<cit> , showed that our method allowed a more complete reconstruction of these genes' evolutionary histories, including a broader spectrum of species .
our evolutionary scenario of duplications which occurred after the radiation of mammals is consistent with the co-localization of pad genes in the same genomic region , which suggests cis-duplications.
notice that one third of reliably-inferred alternative mrna isoforms are candidates for nonsense-mediated mrna dacay , an mrna surveillance system  <cit> .
these transcripts, rather than being translated into proteins, are expected to be degraded and may be subject to regulated unproductive splicing and translation.
in our study, we did not take into account this nmd-type alternative splicing, and took all information  to reconstruct the most parsimonious story of a gene family.
however, when looking at the overall data, we could suppose that some est found in some tissus are not translated into protein but are regulatory elements.
phylogenetic trees show that some pad family members evolve faster than others, with the longest branches lengths seen in the pad- <dig> paralogy group and the shortest branches lengths in the pad- <dig> paralogy group.
these data suggested that pad-2s have probably kept the ancestral function and that other pads may have run a neofunctionalization, still unknown.
for example, the biochemical condition  for enzyme activation could be different according to each pads.
such an hypothesis has just been supported by yamamoto's group who showed small different enzymatic properties between pad- <dig> and pad- <dig>  <cit> .
any kind of expression data could be used to enhance the accuracy of the analysis.
the method used by xun gu et al.
and developed in genetics  <cit>  and pnas  <cit>  is particularly interesting because they used microarray data in yeast to develop a statistical framework for studying the evolution of genes after duplications.
on the basis of a bayesian-based method, they reconstructed the evolutional trace of expression diversity after gene duplication.
their conclusions are in agreement with ours, showing that the expression of duplicate genes tends to evolve asymmetrically, one copy maintaining the ancestral expression profile while the others evolve rapidly.
2/when analyzing transcriptional data , we confirmed that transcriptional shift occurred between paralogs as previously shown  <cit> .
pad- <dig> shows a widespread expression whereas pad- <dig> is expressed essentially during embryonic development.
pad- <dig>  - <dig> and - <dig> have a restricted expression pattern .
for orthologous genes, it is very difficult to compare different tissues in different species because 1/analysis is never done under the same experimental conditions and 2/one can never be sure that tissues evolved under the same physiological conditions in the different species.
3/our analysis allowed us to correlate  sequence evolution rate and tissular expression distribution of paralogs .
it is interesting to notice that shorter branches lengths are correlated with ubiquitous tissular expression.
one can hypothesize that one copy  has probably kept the ancestral function, with ubiquitous expression.
other copies may have slightly different biochemical function, and underwent shifts in their expression pattern which is different and more limited.
these data are consistent with studies showing that housekeeping genes evolve more slowly than tissue specific genes, and are under stronger specific constraints  <cit> .
our data suggest that the same phenomena could exist in multigenic families but, it must be confirmed in other examples.
CONCLUSIONS
we believe our method of adding ests and expression data to phylogenetic analysis provides a new way to annotate multigenic families.
more than classical phylogeny, it allows highlighting of transcriptional shifts between paralogs.
it shows that functional shifts can occur in differential tissue expression rather than in biochemical function of the protein.
it also shows a correlation between lower sequence evolution rate  and larger tissular expression distribution.
for drug development, it points out the fact that when one analyses an orthologous protein in a species phylogenetically close to humans, one should keep in mind that tissular distribution of a protein can be different between species before extrapolation of the function to human.
this type of analysis is in its infancy and must be extended to other multigenic families and to all kinds of expression databases, including database where expression data are normalized such as in unigene, sage and micro arrays.
in the future, the availability of more and more sequence information from different species will enable tracing a genome's evolutionary history, down to the expression data level, with more accuracy.
methods
1-phylogenetic reconstruction with full-length proteins, from ncbi nr database
pad- <dig> was chosen as the query sequence for the first phylogenetic analysis.
to build the phylogeny, the following steps, as described elsewhere  <cit>  were used:
blast
search blast hits on nr protein database  <cit>  with e-value filter set to 1e- <dig>
multiple alignment
elimination of sequences disturbing alignments and doubles.
calculation of an alignment with clustalw software  <cit> .
elimination of large gaps.
domains searching with hmmpfam software  <cit> .
composition test for residues on tree-puzzle  <cit> .
for each domain, elimination of non-monophyletic "repeats" in a tree built with neighbor joining algorithm on clustalw software  <cit> .
elimination of sequences with divergent residue composition by using the amino-acid composition test from tree-puzzle software  <cit> .
elimination of sites not under neutrality  by statistical methods for testing functional divergence after gene duplication  <cit> .
selection of congruent domains with hompart test from paup package  <cit> , and building of a new alignment by merging preserved parts of domains' alignments.
tree building
from this alignment, three phylogenetic trees were generated using the following methods: neighbor joining  <cit> , maximum parsimony in paup  <cit> , maximum likelihood in puzzle  <cit> .
by comparing topologies of these trees with pscore command  from paup package and kishino-hasegawa test  <cit>  from tree-puzzle package, production of a fusion of these trees in a unique consensus tree according to the results of these tests.
comparison of this fusion tree with the ncbi tree of life, , allowed deducing paralogous and orthologous proteins, using the algorithm described by zmasek et al.
<cit> .
trees were rooted with the midpoint method.
a fusion of the congruent three trees is noted npl-a at the first root .
respective bootstrap values are noted in front of each node.
analysis of the tree and rerooting
our phylogenies are built with the midpoint method rooting.
in order to avoid misclassification induced by fast evolving sites, sites which are not under neutrality are generally eliminated, using statistical methods for testing functional divergence after gene duplication  <cit> .
in the phylogenetic "tree of reference", fast evolving sites could not be removed by gu software because not enough sequences were present in the different paralogy groups.
so, the pad- <dig> paralogy group, which evolved faster than the others was automatically placed in out-group to equilibrate the tree.
this was not necessarily reflective of the true evolutionary history of the pad family.
indeed, when looking at the branches lengths between species belonging to the same paralogy group, particularly mus musculus, rattus norvegicus and homo sapiens, we can deduce that duplication of the pad ancestor occurred after the separation of sarcopterygii and actinopterygian.
all of the species external to these groups in the tree of life were placed in out-group.
indeed, if one placed pad- <dig> in the out-group, the conclusion concerning pad' evolutionary events would change.
one would conclude that pads duplicated before separation of euchordates, but not after the emergence of sarcopterigian.
in this hypothesis, cephalochordates and teleosteans pads reside in the pad- <dig> paralogous group.
this scheme is based on the hypothesis that pad- <dig>  - <dig>  - <dig>  - <dig> were all lost during evolution in teleosteans and cephalochordates.
the tree built with padi- <dig> in the out-group is farless parsimonious with regards to evolutionary events and is deducted from too many supposed gene losses.
the same conclusions could be arrived at with other pad groups placed in the out-group.
so, a second tree, more parsimonious regarding evolutionary events was built with an out-group species rooting .
all this phylogenetic reconstructions were performed with the figenix automated platform  <cit> .
2-building the tree including translations of ests contigs
in a second step, the tree was completed with sequences extracted from est database.
blast
we first searched by blast request in ncbi dbest database, sequences similar to pads .
creation of a local database
these short ests were put in a local database.
this local database was then completed by mrna  of completely sequenced proteins .
contigs generation
because ests are not full length mrnas of known proteins, they cannot be included in a phylogenetic building in this raw state.
we built a specific protein database using est clustering.
in order to reduce sequencing error, 10% of the est sequence length was eliminated from each tip and only the  <dig> core central nucleotides were kept.
each species' est group was analyzed with cap <dig>   <cit>  to produce rna contigs .
these rna sequences were translated into proteins .
phylogenetic reconstruction
phylogenetic reconstruction was performed with the same method as described above.
the phylogenetic analysis was at first built with the whole pad- <dig> protein as the query sequence.
some short sequences  were eliminated during the phylogenetic process.
for these sequences, phylogenetic trees were built using five arbitrary created domains of pad- <dig> , used as the sequence query .
using this method, most of the proteins translated from est contigs could be analyzed and classified in a group of paralogs.
however,  <dig> sequences were not overlapped by these arbitrary domains.
these singletons were analyzed as sequence query in separated phylogenies.
the trees were validated only if their topologies were similar to the topology of the reference tree .
as an example, we provide a tree built with a protein from an amphioxus pad est singleton .
est contigs incorporated in the trees are grouped with complete sequences .
five trees built with this method are available on the supplementary material on line.
transcriptional pattern analysis
phylogenetic analysis allowed classification of each est into a paralogy group.
information concerning tissue or organ provenance is given for each est in the ncbi dbest database  <cit> .
we developed software "agents" in figenix platform to automatically collect information corresponding to each est.
this expression information is: ests content of each contigs and tissue or organ expression corresponding to each est.
these data were pooled and classified in a summarizing table .
for an accurate analysis of expression, we normalized expression data relative to the number of clones in each library .
data were congruent with expression profile available for homo sapiens and mus musculus available in unigene  <cit> .
differential expression significance and data consistency
in order to compare the expression levels from est libraries between the different paralogous copies of pad genes, we used the statistical test described by audic and claverie  <cit> .
for each tissue type and cell category, we compared the est expression data  for a given species between the different paralogous copies of pad genes.
when a library was shared between two paralogous copies, we directly compared the two ratios.
the threshold used for differential expression was set to p > 95% either for under or over expression.
when several different libraries were hit by a pad gene for a given tissue or cell type, we first checked whether each library was consistent with all the other ones by applying the same statistical test  <cit> .
when a given library was different from more than 50% of other libraries , we removed the significantly different library.
only tissue and cells not associated with a pathological state were considered for comparison, as expression can be altered in cancers and other diseases.
similarly, genetically modified or selected mice variants were removed from the comparison as they may not reflect normal expression.
authors' contributions
nathalie balandraud initiated the study, collected the data and performed most steps of bioinformatic and statistical analysis presented in the paper.
she developed the initial framework for the analysis of the est genomic sequences.
she produced all the phylogenies, all the figures in the paper and table  <dig>  and wrote the manuscript.
philippe gouret designed and supervised all the bioinformatic study with the creation of «figenix» software.
etienne g.j.
danchin participated in creation of «figenix» software, in the annotation of the gene set, participated in result interpretation, designed the statistical tests applied in the study and in writing later versions of the manuscript.
mathieu blanc built the software agent to retrieve est tissue expression data.
daniel zinn co-participated in the creation of the software agent for ests.
jean roudier participated in the inititation of the study and in writing latest version of the manuscript.
pierre pontarotti supervised the study and co-wrote the manuscript.
supplementary material
additional file 1
supplementary trees built with est contigs.
a simplified tree of life.
table 2: normalized complete table
click here for file
acknowledgements
we thank the "société française de rhumatologie"  that provided funds.
we thank dr. nathalie lambert and dr. daniel gautheret for useful discussion and manuscript correction.
we thank nick boespflug for manuscript correction.

 the "in silico" functional annotation of proteins generated by large scale sequencing projects is an important challenge in biology. here we propose a rigorous protocol for multigenic families' annotation. 1/importance of phylogenetic reconstruction because important functions are conserved during evolution, the first step in analysis is to determine homologous sequences. more specifically, orthologs are more likely to share the same function while paralogs can undergo functional shifts. developments and improvements in database similarity search programs such as blast allowed the rapid identification of potential homologous sequences and therefore allowed functional prediction of several thousand of genes and proteins present in databases. however, the closest blast is often not the nearest neighbor. indeed, similarity-based approaches do not consider all the information from comparative and evolutionary biology. they do not differentiate between orthologs and paralogs among homologs. so, phylogenetic approaches, taking into account duplication and speciation events are necessary to robustly produce functional annotation of new, uncharacterized proteins. 2/necessity to enlarge databases for proteins' phylogenetic reconstruction, protein databases containing the proteomes of completely sequenced species along with individually submitted protein sequences are usually used. yet, the vast majority of species are not fully sequenced and most of their protein sequences are still unknown. however, a lot of transcriptional information is carried by growing gene expression databases, concerning normal or pathological tissues. these mrnas could be used for reconstruction of unknown proteins in "not yet sequenced" species. in parallel, translation of est contigs can be used to enlarge the spectrum of species containing homologs when one analyses a protein family. 3/importance of expression patterns' determination for an accurate annotation it should be noted that phylogenetic analysis can only give information at the biochemical function level. furthermore, while orthologs can have very similar "molecular function", they can exhibit different "macroscopic functions", due to a transcriptional shift for example. to produce an accurate proteic functional annotation, one must have complete sequence information given by phylogenetic reconstruction, with expression patterns analysis. this is the second reason why using data from expression databases is interesting. analysis of expression divergence between paralogs and orthologs have been recently published in human and mouse. it appears that gene expression profiles diverge between paralogs. orthologs can diverge in their expression pattern too. moreover, orthologs that have undergone recent duplication have less strongly correlated expression profiles than those that have not. up to now, there is no study examining expression divergence between homologs that takes into account a broad spectrum of species and after a phylogenetic reconstruction. 4/our approach we present here a new way to functionally annotate proteins in silico, taking into account all these concepts. in a first step, we reconstructed the phylogeny of a protein family, using an enlarged database containing 1/full length proteins from ncbi nr protein database and 2/translation of est contigs from ncbi dbest database. we used a new software platform, figenix, adapted to this kind of phylogenomic reconstruction. in a second step, we created an automated pipeline to couple these phylogenetic reconstructions with expression pattern data. we then compared the evolution rate of the different paralogous sequences with their patterns of expression. we validated this approach with the peptidyl arginine deiminase multigenic family. these genes encode peptidyl arginine deiminase proteins that are implied in rheumatoid arthritis, a systemic autoimmune disease that primarily manifests as a chronic symmetric polyarthritis which gradually destroys articular cartilage and bone. pads transform arginine residues into citrulline in a calcium-dependant manner. it was recently shown by a case-control linkage disequilibrium study that pad type is a susceptibility locus for rheumatoid arthritis in a japanese population. these data were not retrieved in a caucasian population. up to now, no phylogenomic annotation coupled with expression data of this protein family has been proposed. we believe this analysis will shed new light on this incompletely characterized gene family.  1/phylogenetic reconstruction first, we built a phylogenetic tree with full length pad proteins. we found five paralogs of pad proteins encoded by five different genes as previously described. we found pad members in birds, amphibians and teleosteans. this tree supported by high bootstrap values suggested the following steps for evolution of pad genes in craniates: pads' common ancestor was present as far back as in the last common ancestor of teleostean and mammals. pad- first diverged by duplication from the common ancestor before the radiation of mammals. we next built a phylogeny with the same protein dataset completed by additional proteins built from translations of est contigs. first, blast searches on est databases allowed us to identify other species in which pad proteins are likely present. we found pads in additional teleosteans and in cephalochordates. because ests are not full length sequences, they can not be used in this raw state to reconstruct a phylogenetic tree. to correct for this, we built larger sequences by merging ests relative to the same gene into contigs. est contigs were translated and included in a local database and used for phylogenetic reconstruction along with protein from ncbi nr database. because ests contigs often do not correspond to full length transcripts, we had to build several trees with alignment only on parts of the protein. a tree built for each domain is available in additional file online. some ests could not be integrated into a contig and were used in this state as sequence queries to reconstruct a phylogeny. an example of a tree built with a singleton is shown in figure all the trees built with partial sequences were validated only if they were in agreement with the topology of the reference tree. this allowed us to classify all the contigs and singletons into paralogy or orthology groups. this analysis allowed us to propose a more detailed evolutionary story for the pad genes: the common ancestor of the pad gene was at least present in the common ancestor of the euchordates. it means that the pad gene present in amphioxus is co-orthologous to all pads. analysis of trees' topologies allowed us to conclude that duplications of pads certainly occurred after the speciation of craniata or more probably after the speciation of sarcopterygian, but before the mammals' radiation. so, phylogenetic trees including translated est contigs allowed us to enlarge the dataset of species and to classify these "new proteins" into a paralogy group. when looking at the branches lengths, which are directly correlated with the sequence evolution rate, it appeared that some members of pad genes evolved faster than the others. branches lengths were longer within the pad- paralogy group whereas branches lengths were the shortest within the pad- paralogy group. it simply means that the similarity between species is higher for pad- than it is for pad- this indicates that pad-2s' group is the most highly conserved. furthermore, the distance from the common ancestor is lower for pad-2s' group than it is for the others, suggesting that pad-2's group has been submitted to high negative selective pressure and has probably retained the ancestral biochemical function. 2/coupling expression pattern with phylogenetic classification the phylogenetic reconstruction allowed us to classify all the contigs or singletons in phylogenetic groups for the pads phylogeny. we then developed software agents to extract transcriptional data from the ncbi dbest database. we first normalized expression data, because expression data are not normalized in ncbi dbest. we then built a table classifying ests according to their position in the tree and to the tissue in which they are expressed. results are illustrated in a simplified figure raw data are in table. given the current est coverage, we focused on human and mouse. the bias could occur when, say, the coverage of the zebrafish is much lower than human or mouse. est data are all represented in the table, but statistical analysis took into account only human and mouse data. for homo sapiens and mus musculus, data are congruent with expression patterns data available on unigene. coupling phylogenetic classification with expression data permitted us to perform updated footprints of the transcriptional pattern for each paralogs group in this multigenic family. one group of paralogs, pad- showed a broad tissular expression. as the expression pattern was different across paralogous pad genes, we could only compare a few tissue types for differential levels of expression. the statistical test described by audic and claverie was used. significant differences were shown in thymus, whole blood, lung, inner ear, mixed liver and spleen library and leucocytes. 3/coupling sequence evolution rate and expression pattern when observing tissue distribution, we noticed that it was conversely correlated to sequence evolution rate. this was remarkable when we compared branches lengths and tissue distribution between pad- and pad-. discussion 1/using phylogenomics and data from est databases using est data to complete data from protein databases is commonly performed. yet, if one uses only the best blast hits approach to classify ests data, there is a large risk of misclassification into paralogy or orthology groups. thus, expression data for a given paralogy group could be undermined by the wrong assignment of est data into a group. in our example, we found additional expressed sequences from teleostean and cephalochordates in est databases. when performing analysis with reciprocal best blast hit on nr protein database; est from amphioxus and teleostean is classified as an ortholog of pad- only, but not as a co-ortholog of all pads as is shown by phylogenetic analysis. this illustrates why phylogenetic reconstruction is essential for accurate est analysis and able to better avoid misclassification. our approach enlarges the dataset of homologs-containing species and robustly and accurately classifies multigenic family members into paralogy and orthology group. comparison of our results to those of chavanas et al., showed that our method allowed a more complete reconstruction of these genes' evolutionary histories, including a broader spectrum of species. our evolutionary scenario of duplications which occurred after the radiation of mammals is consistent with the co-localization of pad genes in the same genomic region, which suggests cis-duplications. notice that one third of reliably-inferred alternative mrna isoforms are candidates for nonsense-mediated mrna dacay, an mrna surveillance system. these transcripts, rather than being translated into proteins, are expected to be degraded and may be subject to regulated unproductive splicing and translation. in our study, we did not take into account this nmd-type alternative splicing, and took all information to reconstruct the most parsimonious story of a gene family. however, when looking at the overall data, we could suppose that some est found in some tissus are not translated into protein but are regulatory elements. phylogenetic trees show that some pad family members evolve faster than others, with the longest branches lengths seen in the pad- paralogy group and the shortest branches lengths in the pad- paralogy group. these data suggested that pad-2s have probably kept the ancestral function and that other pads may have run a neofunctionalization, still unknown. for example, the biochemical condition for enzyme activation could be different according to each pads. such an hypothesis has just been supported by yamamoto's group who showed small different enzymatic properties between pad- and pad-. any kind of expression data could be used to enhance the accuracy of the analysis. the method used by xun gu et al. and developed in genetics and pnas is particularly interesting because they used microarray data in yeast to develop a statistical framework for studying the evolution of genes after duplications. on the basis of a bayesian-based method, they reconstructed the evolutional trace of expression diversity after gene duplication. their conclusions are in agreement with ours, showing that the expression of duplicate genes tends to evolve asymmetrically, one copy maintaining the ancestral expression profile while the others evolve rapidly. 2/when analyzing transcriptional data, we confirmed that transcriptional shift occurred between paralogs as previously shown. pad- shows a widespread expression whereas pad- is expressed essentially during embryonic development. pad- - and - have a restricted expression pattern. for orthologous genes, it is very difficult to compare different tissues in different species because 1/analysis is never done under the same experimental conditions and 2/one can never be sure that tissues evolved under the same physiological conditions in the different species. 3/our analysis allowed us to correlate sequence evolution rate and tissular expression distribution of paralogs. it is interesting to notice that shorter branches lengths are correlated with ubiquitous tissular expression. one can hypothesize that one copy has probably kept the ancestral function, with ubiquitous expression. other copies may have slightly different biochemical function, and underwent shifts in their expression pattern which is different and more limited. these data are consistent with studies showing that housekeeping genes evolve more slowly than tissue specific genes, and are under stronger specific constraints. our data suggest that the same phenomena could exist in multigenic families but, it must be confirmed in other examples.  we believe our method of adding ests and expression data to phylogenetic analysis provides a new way to annotate multigenic families. more than classical phylogeny, it allows highlighting of transcriptional shifts between paralogs. it shows that functional shifts can occur in differential tissue expression rather than in biochemical function of the protein. it also shows a correlation between lower sequence evolution rate and larger tissular expression distribution. for drug development, it points out the fact that when one analyses an orthologous protein in a species phylogenetically close to humans, one should keep in mind that tissular distribution of a protein can be different between species before extrapolation of the function to human. this type of analysis is in its infancy and must be extended to other multigenic families and to all kinds of expression databases, including database where expression data are normalized such as in unigene, sage and micro arrays. in the future, the availability of more and more sequence information from different species will enable tracing a genome's evolutionary history, down to the expression data level, with more accuracy. methods 1-phylogenetic reconstruction with full-length proteins, from ncbi nr database pad- was chosen as the query sequence for the first phylogenetic analysis. to build the phylogeny, the following steps, as described elsewhere were used: blast search blast hits on nr protein database with e-value filter set to 1e- multiple alignment elimination of sequences disturbing alignments and doubles. calculation of an alignment with clustalw software. elimination of large gaps. domains searching with hmmpfam software. composition test for residues on tree-puzzle. for each domain, elimination of non-monophyletic "repeats" in a tree built with neighbor joining algorithm on clustalw software. elimination of sequences with divergent residue composition by using the amino-acid composition test from tree-puzzle software. elimination of sites not under neutrality by statistical methods for testing functional divergence after gene duplication. selection of congruent domains with hompart test from paup package, and building of a new alignment by merging preserved parts of domains' alignments. tree building from this alignment, three phylogenetic trees were generated using the following methods: neighbor joining, maximum parsimony in paup, maximum likelihood in puzzle. by comparing topologies of these trees with pscore command from paup package and kishino-hasegawa test from tree-puzzle package, production of a fusion of these trees in a unique consensus tree according to the results of these tests. comparison of this fusion tree with the ncbi tree of life,, allowed deducing paralogous and orthologous proteins, using the algorithm described by zmasek et al.. trees were rooted with the midpoint method. a fusion of the congruent three trees is noted npl-a at the first root. respective bootstrap values are noted in front of each node. analysis of the tree and rerooting our phylogenies are built with the midpoint method rooting. in order to avoid misclassification induced by fast evolving sites, sites which are not under neutrality are generally eliminated, using statistical methods for testing functional divergence after gene duplication. in the phylogenetic "tree of reference", fast evolving sites could not be removed by gu software because not enough sequences were present in the different paralogy groups. so, the pad- paralogy group, which evolved faster than the others was automatically placed in out-group to equilibrate the tree. this was not necessarily reflective of the true evolutionary history of the pad family. indeed, when looking at the branches lengths between species belonging to the same paralogy group, particularly mus musculus, rattus norvegicus and homo sapiens, we can deduce that duplication of the pad ancestor occurred after the separation of sarcopterygii and actinopterygian. all of the species external to these groups in the tree of life were placed in out-group. indeed, if one placed pad- in the out-group, the conclusion concerning pad' evolutionary events would change. one would conclude that pads duplicated before separation of euchordates, but not after the emergence of sarcopterigian. in this hypothesis, cephalochordates and teleosteans pads reside in the pad- paralogous group. this scheme is based on the hypothesis that pad- - - - were all lost during evolution in teleosteans and cephalochordates. the tree built with padi- in the out-group is farless parsimonious with regards to evolutionary events and is deducted from too many supposed gene losses. the same conclusions could be arrived at with other pad groups placed in the out-group. so, a second tree, more parsimonious regarding evolutionary events was built with an out-group species rooting. all this phylogenetic reconstructions were performed with the figenix automated platform. 2-building the tree including translations of ests contigs in a second step, the tree was completed with sequences extracted from est database. blast we first searched by blast request in ncbi dbest database, sequences similar to pads. creation of a local database these short ests were put in a local database. this local database was then completed by mrna of completely sequenced proteins. contigs generation because ests are not full length mrnas of known proteins, they cannot be included in a phylogenetic building in this raw state. we built a specific protein database using est clustering. in order to reduce sequencing error, 10% of the est sequence length was eliminated from each tip and only the core central nucleotides were kept. each species' est group was analyzed with cap to produce rna contigs. these rna sequences were translated into proteins. phylogenetic reconstruction phylogenetic reconstruction was performed with the same method as described above. the phylogenetic analysis was at first built with the whole pad- protein as the query sequence. some short sequences were eliminated during the phylogenetic process. for these sequences, phylogenetic trees were built using five arbitrary created domains of pad-, used as the sequence query. using this method, most of the proteins translated from est contigs could be analyzed and classified in a group of paralogs. however, sequences were not overlapped by these arbitrary domains. these singletons were analyzed as sequence query in separated phylogenies. the trees were validated only if their topologies were similar to the topology of the reference tree. as an example, we provide a tree built with a protein from an amphioxus pad est singleton. est contigs incorporated in the trees are grouped with complete sequences. five trees built with this method are available on the supplementary material on line. transcriptional pattern analysis phylogenetic analysis allowed classification of each est into a paralogy group. information concerning tissue or organ provenance is given for each est in the ncbi dbest database. we developed software "agents" in figenix platform to automatically collect information corresponding to each est. this expression information is: ests content of each contigs and tissue or organ expression corresponding to each est. these data were pooled and classified in a summarizing table. for an accurate analysis of expression, we normalized expression data relative to the number of clones in each library. data were congruent with expression profile available for homo sapiens and mus musculus available in unigene. differential expression significance and data consistency in order to compare the expression levels from est libraries between the different paralogous copies of pad genes, we used the statistical test described by audic and claverie. for each tissue type and cell category, we compared the est expression data for a given species between the different paralogous copies of pad genes. when a library was shared between two paralogous copies, we directly compared the two ratios. the threshold used for differential expression was set to p > 95% either for under or over expression. when several different libraries were hit by a pad gene for a given tissue or cell type, we first checked whether each library was consistent with all the other ones by applying the same statistical test. when a given library was different from more than 50% of other libraries, we removed the significantly different library. only tissue and cells not associated with a pathological state were considered for comparison, as expression can be altered in cancers and other diseases. similarly, genetically modified or selected mice variants were removed from the comparison as they may not reflect normal expression. authors' contributions nathalie balandraud initiated the study, collected the data and performed most steps of bioinformatic and statistical analysis presented in the paper. she developed the initial framework for the analysis of the est genomic sequences. she produced all the phylogenies, all the figures in the paper and table and wrote the manuscript. philippe gouret designed and supervised all the bioinformatic study with the creation of «figenix» software. etienne g.j. danchin participated in creation of «figenix» software, in the annotation of the gene set, participated in result interpretation, designed the statistical tests applied in the study and in writing later versions of the manuscript. mathieu blanc built the software agent to retrieve est tissue expression data. daniel zinn co-participated in the creation of the software agent for ests. jean roudier participated in the inititation of the study and in writing latest version of the manuscript. pierre pontarotti supervised the study and co-wrote the manuscript. supplementary material additional file 1 supplementary trees built with est contigs. a simplified tree of life. table 2: normalized complete table click here for file acknowledgements we thank the "société française de rhumatologie" that provided funds. we thank dr. nathalie lambert and dr. daniel gautheret for useful discussion and manuscript correction. we thank nick boespflug for manuscript correction.
large scale and reliable proteins' functional annotation is a major challenge in modern biology. coupling phylogenetic reconstruction and expression data is a promising way of analysis that could be applied to all multigenic families to investigate the relationship between molecular and transcriptional evolution and to improve functional annotation. we then extracted all corresponding expression data contained in est database this analysis allowed us 1/to extend the spectrum of homologs-containing species and to improve the reconstruction of genes' evolutionary history. 3/to show a correlation between paralogous sequences' evolution rate and pattern of tissular expression. to illustrate our method, we used an example protein family, the peptidyl arginine deiminases, probably implied in rheumatoid arthritis. however, up to now, phylogenetic annotation did not take into account expression data. we developed an approach enabling the combination of expression data and phylogenetic analysis. we completed the phylogenetic reconstruction of pads using an enlarged sequence database containing translations of ests contigs. 2/to deduce an accurate gene expression pattern for each member of this protein family.
 large scale and reliable proteins' functional annotation is a major challenge in modern biology. phylogenetic analyses have been shown to be important for such tasks. however, up to now, phylogenetic annotation did not take into account expression data. therefore, integrating such data, like ests in phylogenetic annotation could be a major advance in post genomic analyses. we developed an approach enabling the combination of expression data and phylogenetic analysis. to illustrate our method, we used an example protein family, the peptidyl arginine deiminases, probably implied in rheumatoid arthritis.  the analysis was performed as follows: we built a phylogeny of pad proteins from the ncbi's nr protein database. we completed the phylogenetic reconstruction of pads using an enlarged sequence database containing translations of ests contigs. we then extracted all corresponding expression data contained in est database this analysis allowed us 1/to extend the spectrum of homologs-containing species and to improve the reconstruction of genes' evolutionary history. 2/to deduce an accurate gene expression pattern for each member of this protein family. 3/to show a correlation between paralogous sequences' evolution rate and pattern of tissular expression.  coupling phylogenetic reconstruction and expression data is a promising way of analysis that could be applied to all multigenic families to investigate the relationship between molecular and transcriptional evolution and to improve functional annotation. 