I have a large set of protein sequences and would like to create a "similarity matrix" out of them (for phylo trees later). So I have to run n(n-1) sequence alignments. Wondering what library I should use for that and also what sophisticated lib would provide this all against all without me writing a for loop and parsing results myself...

Currently using Bio.pairwise2.align.globalxx from BioPython as a development placeholder but I do not think the results look really useful? It returns a score > 1 and no similarity in [0,1[ or evalue.

I have been out of the Bioinfo game for a while, so any suggestions for libraries (python preferred) much appreciated! :)


1 Answer 1


I wouldn't use Python here. I've looked at using Bio.pairwise2.align. and its has major shortcomings. In this context the problem is establishing homologous alignment positions ... if you want to make a tree.

Whats the problems? This is you start introducing the gap character '-' as an additional state in constructing your tree and the tree will bias towards not the maximum number of substitutions but the maximum gap (indel) difference between the sequences.

The second problem is changes in homologous sites between different pairwise alignments. This gets complicated fast.

Is that bad? The problem is we have little understanding of the mutation model behind generating indels and there isn't a model to account for varying descriptions of homologous alignment positions. However, we have a very clear understanding of amino acid substitutions in a single global alignment, i.e. a singular description of homologous amino acid positions. A lot of work has gone into modelling amino acid mutation, with a tree a practical output of that work.

Thus, I would simply align using muscle5 under the -super5 option. This will produce a global alignment very quickly and for vast amounts of data. This alignment will produce the best model of homologous amino acid sites.

If you just want a distance matrix then import the entire alignment into MEGAXI and request the JTT model of amino acid substitution with the gamma distribution. This will output a distance matrix nicely and will then allow you to produce a tree.

BTW estimating the gamma distribution without maximum likelihood can be difficult in distance methods (MEGA might have resolved this). This tree will be much better than the one you are proposing (you can estimate this parameter via e.g. IQTREE or ModelTest).

Why? Again alot of work has been done on modelling amino acid substations and JTT is one of the better models (assuming this is nuclear sequence not mtDNA). The LG matrix is much better BTW, but I don't think thats in MEGA.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.