BioPython bootstrap is not reliable?

Question

Here i will show you a minimal working example of code and as you can see the support values for the tree is always 100.

I am using synthetic sequences of 100bp for 6 elements. The sequences have been generated at random choosing from ATCG for each position with equal probabilities.

The input fasta file is

>S1
CTCAGGCACTAGGAGCTTCTCCAGGGCAAAGTTGTCTACGAATATGCCGACTCAGAAGGTTATCAATACGGTTTACTTATCTGCACGCCATTTTCCTATG
>S2
TAACAAATGTTTCTCGCTCGAACCCCGTGGTGCGAGGGGTCACGATAAAAGGTCCCTCCTTCCACGGCATAATTGCTCCCTTTCTTTTCCGTGGGCAGGA
>S3
TTCCTAAAGGACGACTGGAAGCCGGGTACACGTCAACGGAGCTTCTAGCCCGAGGTCTATAGACCGGTTAATGAACTAGACTTAATGTGGAGCTCGTAGA
>S4
ACTTCGTGCCCAACCACTCTATGAGCAGAGTGCTCGAATGAACCTCAAAAGGATATCCGCATTTACTCTATATAACAAACGGCCGTCCCCCCATCTGTCA
>S5
TGACACTGGGTCATTTTACCCGCACTGATCGCTGGGCAGGTCGGCAATTTCGTCAGAATGCCGTGGCCGCACTGAAAATATCATTACCGGACTGAGTATG
>S6
ATAAGCCGAGGGGTAGCCCTCTATTTCGCACCGTATAGACGAAGTGATAAACTTTCTAACCACGTGTCGCCTATCTCTACCTAGCCACATTTGAGGTGCG

The following python code makes the multiple alignment using MUSCLE algorithm

### alignment
from Bio.Align.Applications import MuscleCommandline
from Bio import AlignIO, SeqIO
import os
muscle_exe = "C:/Scripts/muscle3.exe"
path = os.getcwd()
with open('tmp_in.fa', 'w') as fp:
    fp.write(fasta)
replacing \ with /
inp_fp = path.replace('\', '/') + '/tmp_in.fa'
out_fp = path.replace('\', '/') + '/tmp_out.afa'
muscle_cline = MuscleCommandline(muscle_exe, input=inp_fp, out=out_fp)

The following code block takes the previous alignment out_fp and creates a tree using UPGMA algorithm.
The number of bootstrap trees is set to 50 and the consensus is majority based.

from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio.Phylo.Consensus import *
read the multiple sequence alignment (msa)
with open(out_fp,'r') as afa:
    alignment = AlignIO.read(afa, 'fasta')
calculator = DistanceCalculator('identity')
constructor = DistanceTreeConstructor(method='upgma', distance_calculator=calculator)
get the list of bootstrap trees
consensus_tree = bootstrap_consensus(alignment=alignment, times=50, tree_constructor=constructor, consensus=majority_consensus)
print(consensus_tree)

The final consensus tree looks like this

Tree(rooted=True)
    Clade()
        Clade(branch_length=0.350925925925926, name='S5')
        Clade(branch_length=0.05833333333333341, confidence=100.0)
            Clade(branch_length=0.05879629629629617, confidence=100.0)
                Clade(branch_length=0.26481481481481484, name='S3')
                Clade(branch_length=0.020370370370370372, confidence=100.0)
                    Clade(branch_length=0.24444444444444446, name='S1')
                    Clade(branch_length=0.24444444444444446, name='S4')
            Clade(branch_length=0.031018518518518404, confidence=100.0)
                Clade(branch_length=0.2925925925925926, name='S6')
                Clade(branch_length=0.2925925925925926, name='S2')

My question is: why does the tree have all support values equal to 100? I am sure there must be an error either on my side or in the package module itself.

Answer 1

I think this is a bug.

It seems to work if you do this, creating an equivalent Alignment object instead of a MultipleSeqAlignment to give the bootstrap step:

from Bio.Align import Alignment
alignment2 = Alignment(list(alignment))
consensus_tree = bootstrap_consensus(alignment=alignment2, times=50, tree_constructor=constructor, consensus=majority_consensus)

---

bootstrap_consensus calls bootstrap_trees, which makes however many randomly shuffled alignments you asked for. But if the input is a MultipleSeqAlignment object like you're using it builds a tree over and over with that same input to the function as the tree function input, so it always gets the same tree over and over. The other branch of bootstrap_trees, which runs if you have an Aignment object, builds trees based on the randomly-selected columns of the input, like you'd expect. Pretty sure that's just a variable naming screw-up one of us should report to Biopython.

Self-contained version of your tree part:

#!/usr/bin/env python
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio import AlignIO
from Bio.Phylo.Consensus import *
out_fp = "tmp_out.fa"
read the multiple sequence alignment (msa)
with open(out_fp,'r') as afa:
    alignment = AlignIO.read(afa, 'fasta')
calculator = DistanceCalculator('identity')
constructor = DistanceTreeConstructor(method='upgma', distance_calculator=calculator)
get the list of bootstrap trees
consensus_tree = bootstrap_consensus(alignment=alignment, times=50, tree_constructor=constructor, consensus=majority_consensus)
print(consensus_tree)
Give bootstrap_consensus an Alignment object instead of MultipleSeqAlignment
object
from Bio.Align import Alignment
alignment2 = Alignment(list(alignment))
consensus_tree = bootstrap_consensus(alignment=alignment2, times=50, tree_constructor=constructor, consensus=majority_consensus)
print(consensus_tree)

Example outputs for the two different input alignment objects:

Tree(rooted=True)
    Clade()
        Clade(branch_length=0.012239583333333349, confidence=100.0)
            Clade(branch_length=0.31666666666666665, name='S3')
            Clade(branch_length=0.033333333333333284, confidence=100.0)
                Clade(branch_length=0.2833333333333333, name='S2')
                Clade(branch_length=0.2833333333333333, name='S4')
        Clade(branch_length=0.033072916666666716, confidence=100.0)
            Clade(branch_length=0.29583333333333334, name='S5')
            Clade(branch_length=0.03750000000000007, confidence=100.0)
                Clade(branch_length=0.2583333333333333, name='S1')
                Clade(branch_length=0.2583333333333333, name='S6')
Tree(rooted=True)
    Clade()
        Clade(branch_length=0.2583333333333333, name='S2')
        Clade(branch_length=0.03405145202020203, confidence=66.0)
            Clade(branch_length=0.275, name='S5')
            Clade(branch_length=0.049892241379310344, confidence=58.0)
                Clade(branch_length=0.21250000000000002, name='S1')
                Clade(branch_length=0.21250000000000002, name='S6')
        Clade(branch_length=0.039416666666666676, confidence=50.0)
            Clade(branch_length=0.2583333333333333, name='S4')
            Clade(branch_length=0.321875, name='S3')

Answer 2

Initially, 100 looks like a good value. However, if you lower to bootstrap value, let's say maybe to 50 as you did, you will probably get an erroneous support value. Bootstrapping (you probably already know this) is how many times a branch in a phylogenetic tree is observed. A good bootstrapping value would be anywhere between 80 and 100. 50 is not a good bootstrapping value. This means the branches are only being observed 50 times as opposed to 85, 90, etc... If I were you, try changing the bootstrapping value of 50 to at least 80. See what that gets you. You may be satisfied with it, but you may want to increase it depending on your output. By doing this, you will see a more realistic and phylogenetically accurate output of your tree. Most bioinformaticians use 100 as a bootstrap value. The only issue, I wouldn't really call it an issue, but it will take longer for an output to occur due to computing time and computing memory available. BUT you will still get an output, you just may need to wait a bit. My bioinformatics background is rooted in phylogenetic and horizontal gene transfer research involving phages and bacteria, so I have seen this issue a lot and have experienced myself a lot! I hope this information provides some insight to your project and bioinformatics journey! :)