Friday, October 10, 2014

A "parsimonious" Bayesian supertree model for estimating species trees

When we have sequence alignments regarding several genes from a group of taxa, we usually want to extract the phylogenetic information common to all of them. However, in many cases such phylogenomic analyses depend on selecting one sequence from each species per gene family (=alignment), or excluding paralogs, or partitioning these paralogous sequences into loci, or utilizing only gene families without apparent paralogs. If we want to analyse all our data at once, without excluding sequences or whole alignments, we are left with few options.

We just published such an alternative, which is based on the idea that we can measure the disagreement between the phylogenetic trees representing each gene and a putative tree representing the species. Therefore, by using disagreement measures that allow for arbitrary mappings between the trees, we can handle gene trees with paralogs, multiple individuals from the same population or missing data. These measures we call "distances" [1], and we developed a probability distribution describing how these distances can work as penalties against very dissimilar gene and species tree pairs.

We can use any combination of the reconciliation distances, the recently developed mulRF distance, and (very experimentally) an approximate SPR distance to include into our multivariate penalty distribution. We are also experimenting with other distances, as we implement them. This penalty distribution is then incorporated into a hierarchical Bayesian model, which I call "parsimonious" since it doesn't use a fully probabilistic model to describe the coalescent processes, or the birth and death of new loci. It assumes instead that only the most parsimonious reconciliations are relevant to the model. (I was advised, however, that calling it a "parsimonious Bayesian model" could be confusing...)

The distance supertree model is based on several distance measures d(G,S) between each gene family tree G and the species tree S. A species trees that is more similar to all gene trees is more likely than a more distant one. Notice that d(G,S) can in fact be a vector with several distances.

We implemented this model into the software guenomu, which is available under a GPL license at http://bitbucket.org/leomrtns/guenomu. The input to the software is a set of files with the distribution of gene trees as estimated for each gene family, independently, and the output will be the posterior distribution of these gene family trees together with the distribution of species trees.  We tested our model on many data sets simulated with the SimPhy software -- which is able to simulate the evolution of gene families with duplications, losses, and the multispecies coalescent fully probabilistically -- followed by a quick-and-dirty emulation of a Bayesian phylogenetic inference [2].

The difference between the input and output (posterior) distribution of trees for each gene family is that the input trees were estimated independently -- let's say, by running MrBayes for each alignment representing a gene family -- while the posterior takes into account the other gene families through their common species tree. Therefore the posterior distribution is a re-sampled version of the input, and as we see in the figure below it improves the gene tree estimation.

Input and posterior distributions if gene trees across many simulations (average values over gene families, per simulated data set). The simulations were pooled by species tree size, where we can see that guenomu can reduce the uncertainty of the gene trees. The accuracy is the fraction of splits (=branches) successfully reconstructed. Figure adapted from doi:10.1093/sysbio/syu082
Our model was successful in reconstructing the species tree even for high levels of incomplete lineage sorting (short species tree branches, in coalescent units) coupled with duplications and losses. It also fared a bit better than iGTP, and much better than our implementation of distance matrix-based species tree inference methods [3]. Notice that only software that accepts gene trees with several tips from the same species can be compared. We were gladly surprised to see that iGTP under the duplication-loss cost also performed well, provided we use the gene tree frequencies as weights.


Violin plots showing the distribution of accuracies in species tree estimation, over all simulations. The two red distributions are for the consensus and MAP tree estimates using guenomu, while the brown and blue plots are for other reconstruction methods. The dendrogram at the top classifies the methods by accuracy. Figure adapted from doi:10.1093/sysbio/syu082


Notes:
[1] They are not proper metrics since they are not symmetric, for instance.

[2] Since our simulated gene families have hundreds of tips, simulating the alignments and then sampling the gene tree distributions with MrBayes or friends would take too long (we did this for smaller data sets only). We therefore created a program (available with guenomu) that would copy many times each tree, replacing randomly short branches by one of its alternative bipartitions.

[3] We must take into account that these matrix-based methods (like GLASS, SD, etc.) assume that all disagreement is due to the coalescent, which is not true under our simulations. Furthermore our implementation may not be as good as some established software. Therefore our results are not evidence against these methods. (I particularly love their idea of being able to work with the distance matrices.)

Reference:
de Oliveira Martins L., Mallo D. & Posada D. (2014). A Bayesian Supertree Model for Genome-Wide Species Tree Reconstruction, Systematic Biology, DOI: http://dx.doi.org/10.1093/sysbio/syu082

(The supplementary material is not available yet at DataDryad, as of today Oct 10, 2014. I assume it will go online soon, but if you want it please drop me a line)

Tuesday, December 31, 2013

If your code is not open, you're not doing science

This is a list of articles explicitly suggesting or tangentially indorsing Open Source Software in the context of academic publications. The list is far from complete and in no particular order, but it creates a strong defense of source code as an integral part of the scientific process (code is data after all!). 




Friday, February 8, 2013

The difference between the RF and the NNI distance

Just to complement my answer to a blog post, where I maintain that the Nearest-Neighbor Interchange (NNI) distance is not equivalent to the Robinson-Foulds (RF) distance, a simple example:



Where we can see that trees T1 and T2 differ only in the location of nodes A and B -- on these trees, we can naturally think of the nodes A, B, 1,..., 6 as representing leaves, but they might also be large subtrees.

The RF distance is the number of edges (=branches) that are unique to each tree (that's why it's also called the symmetric difference), and it may be normalized to one. If we highlight the unique edges on trees T1 and T2


We see that the (unnormalized) RF distance is 10. For dichotomic trees, the number of unique edges is the same on both trees.

The NNI distance is the minimum number of NNIs that must be applied to one tree such that it becomes equal to the other. One NNI branch swap will change exactly one edge, thus is very tempting to assume that the NNI distance can be found by looking at the distinct edges.

But the problem is when the same branch is involved in more than one path of the "NNI walk". The RF distance (divided by two, for fully resolved trees) is then a lower bound on the minimum number of NNIs. In our example:


The NNI distance between T1 and T2 is 6, one more than the RF distance since the edge splitting (1,2,3) and (4,5,6) is used twice in the NNI computation. The problem, as explained by Liam, is that simulating trees with a specified distance is hard, and the solution of using very large trees masks the cases where the distances disagree...

Reference:
Bryant D. (2004). The Splits in the Neighborhood of a Tree, Annals of Combinatorics, 8 (1) 1-11. DOI:

(Crossposted from Bioinformatics News and Reviews, my personal blog)

Friday, June 15, 2012

Postdoctoral positions to work with Phylogenomics in Vigo

[crossposted from BioinfoNews]

There are two postdoc positions available in our lab -- the Phylogenomics group of the University of Vigo, in Spain. You can read the announcement in evoldir or on the link below, from the European Commission. Basically we are looking for someone with a good theoretical/empirical background in phylogenomics (NGS, phylogenetics, HPC, gene-species trees discordance, etc.) The applicant should have a good combination of computational, statistical and practical analysis skills. In any case, please get in contact with David and keep in mind that a honest evaluation of your skills will be greatly appreciated.

EURAXESS - Postdoctoral positions in Phylogenomics -- University of Vigo, Spain
Two postdoctoral positions are available to work on phylogenomics methodology and applications under a project funded by the European Research Council (ERC) in David Posada’s lab at the University of Vigo, Spain (http://darwin.uvigo.es).

DESCRIPTION: Two postdoctoral positions are available to work on phylogenomics methodology and applications under a project funded by the European Research Council (ERC) in David Posada’s lab at the University of Vigo, Spain (http://darwin.uvigo.es).


Initial appointments will be made for one year, with possible extension. Gross annual salary including benefits will be around 24-30 kEuros, commensurate with experience. Starting date should be as soon as possible .

REQUIREMENTS: Candidates should have a doctoral degree in Biology, Computer Science, Statistics or Mathematics. Excellent bioinformatic skills and prior experience with phylogenomics are essential. Strong communication and teamwork abilities are fundamental. Familiarity with parallel and distributed computing environments is convenient. Advanced statistical skills would be a plus.

APPLICATION: Please send a letter of interest, C.V., and the names and contact details of two referees to David Posada at dposada@uvigo.es, indicating “postdoc phylogenomics” in the subject of the email. Questions and requests for more information should be directed at the same address. Review of applications will begin immediately, and continue until the positions are filled.
Random picture of the Vigo bay

PS: In September we will also have a Computational and Statistical Phylogenomics Meeting in town. Besides being a great opportunity for knowing Vigo, the invited talks represent some of the most interesting topics in theoretical phylogenomics.

Tuesday, May 15, 2012

Testing for common ancestry

ResearchBlogging.orgOur commentary on Douglas Theobald's test from Universal Common Ancestry (UCA) just went online. The original idea was to make a user-friendly review of his analysis described in "A formal test of the theory of universal common ancestry", but after a long e-mail exchange between Douglas and us -- actually between him and David, I didn't say much -- we decided to expand the article to include some remaining points of skepticism and spell out the basic problem with his approach.

His work

His test for UCA compares the hypothesis of one single ancestral lineage diverging into all living forms, against scenarios where more than one ancestral populations are needed to explain the current diversity of life. That is, he tries to quantify the possibility that there were two or more ancient life forms still represented today in the major domains, in comparison with the most natural possibility that only one ancestral life form prevailed (with other life forms eventually going extinct, for instance).

These scenarios can be nicely represented by phylogenies: under the UCA hypothesis all existing species can be connected by a single phylogeny (connected network or tree), while under the hypothesis of Independent Origins (IO) the species can be partitioned into disconnected groups (without a branch joining them). And then, once we have the best phylogenies under each hypothesis, we can use the arsenal of model selection methods available to chose between the hypotheses.

Using a curated set of genes highly conserved among all three domains, his test indicated that UCA is a much better explanation than IO, and that each gene is better represented by its own phylogeny than by forcing all genes to follow the same tree -- that is, horizontal gene transfer (HGT) cannot be neglected.

Our comments

Imagine that I developed a method capable of predicting a candidate's academic success based on a questionnaire. But it works only for PhDs on their mid-twenties who have at least a few high-impact single-author publications. And notice: I'm not assuming at all that such a young genius has a secured place in a university. I guess this half-baked analogy summarizes our contention with Douglas' paper.

What called our attention is the fact that the phylogenetic inference methods assume homology at each site, so it is not surprising that it favors UCA for sufficiently good alignments -- your mileage may vary on the definition of "good". These methods are delegating to the alignment the responsibility of handling homology. And he used a data set with a particularly convincing evidence of common ancestry. To transform our argument into a picture we simulated sequences under both UCA and IO scenarios (that is, using one or two phylogenies) and looked at how the resulting alignments would look like. As expected they were very different.

average sequence identity for alignments simulated under UCA and IO
(from doi:10.4081/eb.2012.e5)

The motivation for our skepticism in the UCA test is how it would perform on a blind experiment like Assemblaton or CASP: given a group of sequences of unknown homology status, can the model selection devised by Douglas Theobald tell us if they share a single ancestry? Our impression is that there would be several decisions before doing the actual test -- like optimizing the alignment, possible removal of poorly aligned regions, refusing to do the test if alignment is bad --that might undermine its applicability. So we cannot yet recommend the test for arbitrary data sets.

Frequencies of average identity per column, for alignments simulated under UCA and IO, with real data set values in gray
(from doi:10.4081/eb.2012.e5)

In our article we also wonder about the effect of HGT under the hypothesis of multiple ancestry: what if we find one, and only one gene that strongly supports independent origins? Even if all others fit nicely into a single phylogeny, wouldn't it be evidence of this otherwise extinct lineage?

The publication

We all have dreadful stories about the feedback from peer reviewers, but this manuscript was not such a case. All reviewers seemed to know very well the original work, and could make precise comments on what we were missing or mistaken. The editor, David Liberles, also joined the discussion and gave us some good advice. So we thank them for being fortiter in re, suaviter in modo.

We are preparing another manuscript that is more centered on the examples given in D. Theobald's paper. Actually we started to write this "counter-examples" manuscript prior to our present work, but we had to reorganize it since: 1) we realized that it would be harder to understand it without the current work; 2) D. Theobald recently published a reply that is relevant to our discussion, since it contains a response to a former version of our "counter-examples" manuscript. Our present paper became then necessary to minimize future misunderstandings.

The scripts necessary to reconstruct the simulations and graphics used in our study are available at our home page: http://darwin.uvigo.es/common_origin/. Please let me know if you have any trouble running the scripts, or if you want some more information.

References

Leonardo de Oliveira Martins, David Posada (2012). Proving universal common ancestry with similar sequences Trends in Evolutionary Biology, 14 (1) : 10.4081/eb.2012.e5
Theobald, D. (2010). A formal test of the theory of universal common ancestry Nature, 465 (7295), 219-222 DOI: 10.1038/nature09014

(with thanks to Jonathan Eisen for noticing the article).

Friday, April 13, 2012

The first to describe phylogenetic maximum likelihood methods, according to Joe Felsenstein

(correction: my title is incomplete, it should include for molecular sequence data :)
    Theoretical work in Berkeley has been less widely noted. It started with the Jukes­Cantor distance and the stochastic model derived from it. Shortly afterward, Berkeley's famous statistician Jerzy Neyman became involved, owing to contact with Allan Wilson. He was the first to describe maximum likelihood phylogeny methods for molecular sequence data, a development for which I often mistakenly get the credit.
Taking Variation of Evolutionary Rates Between Sites into Account in Inferring Phylogenies -- Joseph Felsenstein at the Journal of Molecular Evolution, Volume 53, Numbers 4-5

Sunday, July 24, 2011

How to summarise a collection of trees that came from a Bayesian analysis

ResearchBlogging.orgAfter running a Bayesian phylogenetic analysis we are usually left with a large collection of trees, that came from the posterior distribution of the model given our data. Then if we want to work with a single tree - that is, to have a point estimate of this posterior distribution of trees - the most usual ways are to calculate the consensus tree or to select the most frequent tree. There are other ways, but let's fix on those by now.

We might not be aware of it, but when we choose for one or another summary we are in fact deciding for the tree estimate that minimizes its distance to all other trees in the set, and in expectation this will be the closest to the true tree under this distance metric (the so called Bayes estimator). This depends on what exactly do we mean by "distance" between trees, and that's what the article "Bayes Estimators for Phylogenetic Reconstruction" (doi 10.1093/sysbio/syr021) is about.

For example, the majority-rule consensus tree is the best we can get if we assume that the Robinson-Foulds distance (RF distance) is a good way of penalizing trees far away from the true one (I won't dwell into the meaning of "truth"; for us, the True tree® is the one that originated the data). To be more explicit, the consensus tree is the one whose RF distance to all trees in the sample is the shortest possible. This will be the closest we can get to the true tree for this sample, if by "close" we mean "with a small RF distance".

Now suppose I don't like the RF metric because I can only count to two: if the trees are the same the distance is zero, but if they are different then the distance is not zero, and I don't care how different they are (think of apples and oranges). In this case the best representative of my sample is the one that appears more often, known as modal value or Maximum A Posteriori (MAP) value, since our sample comes from a posterior distribution. Is it the closest I can get to the true tree for this distribution? Yes, for this particular definition of distance: the MAP tree is the tree that maximizes the expected coincidence with the true tree.

In the article they also mention that if you want to find the tree that minimizes the expected quartet distance to the true value, then the quartet puzzling method will find this tree for you. But the quartet puzzling tree is not as easy to calculate as the consensus or MAP tree, and there is no straightforward way to find the tree that minimizes other distances in general (e.g. the dSPR, the geodesic distance or the Gene Tree Parsimony). Therefore the authors offer the well-known hill-climbing heuristics for finding the best tree, and use the squared path difference as an example of distance.

Below you can find the presentation I gave to my group last week about this paper, it contains basically some background information and a summary of their method. One thing that is absent from the slides are the results, which I briefly summarize below:
  • their method (called "Bayes" in the figures or "BE") always used the path difference as distance measure; this is the overall distance they were trying to minimize.
  • they simulated many data sets with several levels of sequence divergence, and reconstructed the phylogeny using Maximum Likelihood, Neighbor-Joining, and Bayesian analysis. From the Bayesian posterior distribution they elected as point estimates the consensus tree, MAP tree, and used their method to find the BE under the path difference.
  • Figures 3 and 5 show the distance between the inferred and the true trees, where on figure 3 this distance is the path difference and in figure 5 it is the RF distance. As expected,  the Bayes estimator is better than any other measure at minimizing the path difference distance to the true tree, while the consensus tree wins if we want the closest in terms of RF distance.
  • this result is rephrased in figures 8 and 9, which now look specifically at the distances between BE or MAP trees and the true tree. What they plot is distance(BE, true) - distance(MAP, true) for a different definition of distance(,) in each case. The MAP tree is correlated to the consensus tree (if the MAP frequency is larger than 50% they are equal, for instance). Therefore it should come as no surprise that if we define closeness to the true tree in terms of RF distance, the MAP tree will be closer than the BE as shown in figure 9. Because BE assumes that closeness to true is calculated in terms of the path difference, which is reinforced in figure 8.
  • The authors wisely avoid offering the "best" Bayes estimator, since it depends on your judgment of how to penalize trees different from the true one.




OBS: This was my first time using beamer for Latex (after all these years, I know), so the slides are not prime time material. This is also my first submission to slideshare, and I like the idea of an embedded presentation within the blog post. I use latex a lot, and I think it would be easier for me to prepare a post with figures, equations and text within a presentation, and then simply embed it here with a minimum of extra text.

Maybe I'll try this next time, a presentation but with much more text than the recommended - in real life presentations the slides should support and complement but not replace the lecturer. Then you tell me if you would like to read on such a format or if you prefer a more traditional article-ish post.

Reference:
Huggins, P., Li, W., Haws, D., Friedrich, T., Liu, J., & Yoshida, R. (2011). Bayes Estimators for Phylogenetic Reconstruction Systematic Biology, 60 (4), 528-540 DOI: 10.1093/sysbio/syr021

Saturday, July 23, 2011

Distribution of recombination distances between trees - poster at SMBE2010

(Post exported from my old blog, now defunct -- published originally on July 12, 2010)

I just came back from SMBE2010, where I presented a poster about our recombination detection software and had the chance to see awesome research other people are doing. The poster can be downloaded here (1.MB in pdf format) and I'm distributing it under the Creative Commons License. Given the great feedback I got from other participants of the meeting, I thought it might be a good opportunity to comment about the work, guided by the poster figures. Please refer to the poster to follow the figures and the explanation, I'll try to reproduce here my presentation taking into account the commentaries I received.

The motivation for the development of the model was to be able to say, for a given a mosaic structure, if the breakpoints can be explained by one recombination event or several. The recombination mosaic structure is usually inferred assuming the parental sequences (those not recombining) are known beforehand - in the figure they are the subtypes B and F - and recombination is then inferred when there is a change in the parental closest to the query sequence. Another problem is that it is common to analyze each one of the query sequences independently against the parentals - if all one wants is the "coloring", then this might be enough. For the above figure I analyzed each query sequence against one reference sequence from the subtypes B, F and C (thus comprising a quartet for each analysis). And we know that these mosaics don't tell the whole story.

If we know the topologies for both segments separated by the recombination breakpoint, then we can say, at least in theory, the minimum number of recombination events necessary to explain the difference (the real number can be much larger since we only detect those that lead to a change in the topology...). This minimum number is the Subtree Prune-and-Regraft distance, and is related to the problem of detection of horizontal gene transfers. In our case we devised an approximation to this distance based on the disagreement between all pairs of bipartitions belonging to the topologies: at each iteration we remove the smallest number of "leaves" such that the topologies will become more similar, and our approximate "uSPR distance" will be how many times we iterate this removal.

It is just an approximation, but it is better (closer to the SPR distance) than the Robinson-Foulds or the (complementary) Maximum Agreement Subtree distances, which compete in speed with our algorithm. For larger topologies it apparently works better than for smaller, but this is an artifact of the simulation - one realized SPR "neutralizes" previous ones, and it happens more often for small trees.

Our Bayesian model works with a partitioning of the alignment, where recombination can only occur between segments and never within. This doesn't pose a problem in practice since it will "shift" the recombinations to the border - the idea is that several neighboring breakpoints are equivalent to one breakpoint with a larger distance. This segments could be composed of one site each, but for computational reasons we usually set it up at five or ten base pairs. The drawback is the loss of hability to detect rate heterogeneity within the segment.

Each segment will have it own toplology (represented by Tx, Ty and Tz in the figure), which will coincide for many neighboring segments since we have the distance between them as a latent variable penalizing against too many breakpoints:

This is a truncated Poisson distribution, modified so that it can handle underdispersion - the parameter w will make the Poisson sharper around the mean - and each potential breakpoint has its own lambda and w.

The posterior distribution will have K terms for the segments (topology likelihood and evolutionary model priors) and K-1 terms for the potential breakpoints (distances between segments), as well as the hyper-priors. I use the term "potential breakpoint" because if two consecutive segments happen to have the same topology (distance equals zero) then we don't have an actual breakpoint. Again, considering only the recombinations that change the topology. This posterior distribution is calculated through MCMC sampling in the program called biomc2.

To test the algorithm, we did simulations with eight and twelve taxa datasets, simulating one (for the eight taxa datasets) or two recombinations per breakpoint. We present the output of the program biomc2.summarise, which interprets the posterior samples for one replicate: based on the posterior distribution of distances for each potential breakpoint, we neglect the actual distances and focus on whether it is larger than or equal to zero (second figure of the panel). Based on this multimodal distribution of breakpoints we infer the regions where no recombination was detected (that we call "cold spots"), credible intervals around each mode (red bars on top) or based on all values (red dots at bottom, together with the cold spots).

We also looked at the average distance for each potential breakpoint per replicate dataset, and show that indeed the software can correctly infer the location and amount of recombination for most replicates. It is worth remembering that we were generous in our simulations, in that there is still phylogenetic signal preserved on alignments with many mutations and a few recombinations. If recombination is much more frequent, then any two sites might represent distinct evolutionary histories and our method will fail.

We then analyzed HIV genomic sequences with a very similar mosaic structure, as inferred by cBrother (an implementation of DualBrothers). Here it is important to say that we used cBrother only to estimate the mosaic structure of the recombinants, doing an independent analysis for each sequence against three reference parental sequences. Therefore the figure is not a direct comparison of the programs, contrary to what its unfortunate caption might induce us to think. The distinction is between analyzing all sequences at once or independently, in quartets of sequences. If we superpose the panels it might become clearer to compare them:


click on the figure for a larger version (this one is not on the poster)

The curve in blue shows the positions where there is a change in closest parental for the query sequence, if each query sequence is analyzed neglecting the others. In red we have our algorithm estimating recombinations between all eleven sequences (eight query and three parental sequences). We can see that:

  1. all breakpoints detected by the independent analysis were also detected by our method;

  2. many recombinations were detected only when all sequences were analyzed at once, indicating that they do not involve the parental sequences - de novo recombination;

  3. if we look at the underlying topologies estimated by our method (figure S2 of the PLoS ONE paper), we see that those also detected by the independent analysis in fact involve the parentals while the others don't;

  4. biomc2 not only infers the location of recombination, but also its "strength" - given by the distance between topologies;



Finally, we show two further developments of the software: a point estimate for the recombination mosaic, and the relevance of the chosen prior over distances. The point estimate came from the need of a more easily interpretable summary of the distribution of breakpoints: instead of looking at the whole multimodal distribution, we may want to pay attention only to the peaks, or some other similar measure. This is a common problem in bioinformatics: to represent a collection of trees by a single one or to find a protein structure that represents best an ensemble of structures. In our case we have a collection of recombination mosaics (one per sample of the posterior distribution), and we elect the one with the smallest distance from all other mosaics - we had to devise a distance for this as well...

To show the importance of the prior distribution of distances, we compared it with simplified versions, like setting the penalty parameter w fixed at a low or high value. The overall behavior for all scenarios is lower resolution around breakpoints, and for weaker penalties we reconstruct the topologies better than for stronger ones, at the cost of inferring spurious breakpoints more often. We also compared the original model with a simplification where the topological distance is neglected and the prior considers only if the topologies are equal or not. This is similar to what cBrother and other programs do, and by looking at the top panel we observe that the results were also equivalent (blue lines labeled "cBrother" and "m=0"). In the same panel we plot the performance using our original ("unrestricted") model as a gray area.

I also submitted the poster to the F1000 Poster Bank.

References:
de Oliveira Martins, L., Leal, É., & Kishino, H. (2008). Phylogenetic Detection of Recombination with a Bayesian Prior on the Distance between Trees PLoS ONE, 3 (7) DOI: 10.1371/journal.pone.0002651
de Oliveira Martins, L., & Kishino, H. (2009). Distribution of distances between topologies and its effect on detection of phylogenetic recombination Annals of the Institute of Statistical Mathematics, 62 (1), 145-159 DOI: 10.1007/s10463-009-0259-8

Fault-tolerant conversion between sequence alignments

(Post exported from my old blog, now defunct -- published originally on May 17, 2010)

Despite I'm very charitable when testing my own programs, I'm not so nice when asked to scrutinize other people's work. That's why I was happy to see the announcement about the ALTER web server being published at Nucleic Acids Research (open access!).

I am not involved in the project, but I was in the very comfortable position of being one of the beta testers: all I needed to do is to find the largest and most obscure datasets I had and try them; then complain to the authors about the minimal details. I tried some big datasets (I think it was influenza H3N2 HA and HIV-1 complete genomes from South America, around 2 and 4 Mbytes each), and my simulated alignments created "by hand" from PAML. And ALTER could handle them in the end: they even sent me a report explaining how each one of my commentaries was used to improve the software, and asking me to try again until I feel satisfied.



The ALTER web server is a converter between multiple sequence alignment (MSA) formats, for DNA or protein, focused not only on the format itself (like FASTA or NEXUS) but more on the softwares that generated the alignment and the software where the alignment is going to be used in (e.g. clustal or MrBayes). They mention that this program-oriented format conversion is necessary since all useful softwares eventually violate the (outdated) format specification. In their own words
[D]uring the last years MSA's formats have `evolved' very much like the sequences they contain, with mutational events consisting of long names, extra spaces, additional carriage returns, etc.

The web service can automatically recognize the input format, and generate an output for several programs, in several formats. I found it very easy to use, as you proceed it automatically shows you the possible next steps in the same page. Another very nice feature is the possibility of collapsing duplicate (identical) sequences, working then only with the haplotypes (unique sequences). If later you need the information about the collapsed duplicates check out the "info" panel on the bottom of the screen (inside the "log" window).

The obvious case when this elimination of duplicates is useful is when doing phylogenetic reconstruction (in many cases you can safely remove identical sequences), but another option offered by ALTER is to remove very similar sequences, where you can define the threshold of similarity. Sometimes when I'm doing a preliminary analysis on a dataset, I want to discard sequences too similar in order to get an overall picture of the data, and some other times I must remove closely-related sequences since my recombination-detection program has a limitation on the number of taxa...

Besides the user-friendly web service, they also offer a geek-friendly API - if you want your program to communicate directly with the service - and the source code, licensed under the LGPL.

Reference:
Glez-Pena, D., Gomez-Blanco, D., Reboiro-Jato, M., Fdez-Riverola, F., & Posada, D. (2010). ALTER: program-oriented conversion of DNA and protein alignments Nucleic Acids Research DOI: 10.1093/nar/gkq321

The specialization of novel genes

(Post exported from my old blog, now defunct -- published originally on April 27, 2010)

Recently a paper about the software MANTiS called my attention, and I've been trying to write about it for a while. This announcement at the EvolDir list seemed like the perfect opportunity. I must warn you though that I've never used the software and I don't have any intimacy with the underlying databases, but the article is easy to follow.

The main result of the paper, published in Genome Biology and Evolution, is that there is a correlation between the mean number of anatomical systems (human tissues or cell types) where the gene is expressed and the time when the gene appeared on the phylogeny of the species. In other words, recent gene families are expressed in fewer anatomical systems (are more specific) than ancient ones. An anatomical system is a hierarchical classification of human tissues (e.g. the first level of the hierarchy: nervous, dermal, embryo, etc) available from gene expression data. So the age of appearance of a gene is an indicator of its specificity. Since the genes are subject to duplication we may have more than one member of the gene family in the same species, and the authors show that this correlation is maintained if we consider the appearance of the gene itself (as a result of duplication) or the appearance of the whole gene family to which the gene belongs.

They worked with gene families identified by MANTiS, which is a pipeline that 1) downloads data from metazoan genomes at ENSEMBL, 2) infers the gene tree based on the protein alignment of the gene family and 3) detects duplications through a reconciliation with a given species tree. Each gene tree is produced by EnsemblCompara which, as I understand, employs an extension of "reciprocal best hits" (that allow for many-to-many relations) to find the members of the family, and then maximum likelihood to find the tree itself. I will talk more about the gene tree/species tree reconciliation in the future, but it is enough to say that it's the minimal list of nodes on the gene tree that represent duplications. We have an example of such a reconciled gene tree below, where the duplications are represented by the red boxes:

MANTiS creates a new character (the brown polygons, that I think of as an orthologous group) for each duplication event, and the phylogenetic profile generated by these characters is then used to calculate the branch lengths of the species tree through a least squares approach. The phylogenetic profiles are represented by 0's and 1's in the inlet figure above, from which a distance matrix must be calculated in order to have the branch lengths.

In the study two datasets were created for the presence/absence of genes: one called "families only" composed of one character for each single gene and for each protein family, and another called "with duplications" where a new character is created for each duplication event. Both analyses were necessary since gene gain through duplications is important in explaining genome size increase.

MANTiS creates a database relating each gene to its biological function and anatomical system: the biological processes and molecular functions (ontology terms) of protein families are given by the PANTHER database for human, mouse, rat and D. melanogaster, while the gene expression data (related to the anatomical systems) comes from eGenetics, GNF and HMDEG. When comparing the time of appearance of the gene (as explained above) and the expression data for the genes we have a figure like the following:

We must notice that in this graph the X axis is inverted (that is, left is older with the present day at the right) giving the impression of a negative correlation. So older gene families - or duplications - are expressed in more cell types in humans.  Similar results were obtained using rat expression data - since the expression datasets had information for both - or using the other expression datasets.

The authors say that a possible explanation for this behaviour is the increase in the number of distinct cell types (blue line, notice the inverted axis again :D), where new genes are likely to be more specific to a cell type (which may have appeared recently itself). Associated with this explanation is the subfunctionalization of duplicated genes, and the tendency to subfunctionalize ("specialize") can explain the decreased extent of expression. The subfunctionalization process itself might be related to the generation of a new cell phenotype.

One shortcoming of the analysis is that the gene family inference might fail to detect distantly related genes, and therefore what appears to be a gene gain (the "birth" of a new gene family) might be in fact a duplication of a more ancient single gene family. For example if after the duplication number 3 on the first figure the sequences diverged too much, we might wrongly classify them as two gene families. But to be free from this problem is a tall order. The authors also call our attention to the problem of low coverage of some genomes and taxonomic bias.

References:
Milinkovitch, M., Helaers, R., & Tzika, A. (2009). Historical Constraints on Vertebrate Genome Evolution Genome Biology and Evolution, 2010, 13-18 DOI: 10.1093/gbe/evp052
Tzika, A., Helaers, R., Van de Peer, Y., & Milinkovitch, M. (2007). MANTIS: a phylogenetic framework for multi-species genome comparisons Bioinformatics, 24 (2), 151-157 DOI: 10.1093/bioinformatics/btm567