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.
