4
$\begingroup$

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.

$\endgroup$

2 Answers 2

2
$\begingroup$

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')
$\endgroup$
2
  • 1
    $\begingroup$ Thank you very much i think this was the issue. I will report it to their github if you haven't already raised an issue. $\endgroup$
    – Mirko
    Dec 5, 2023 at 18:34
  • $\begingroup$ I haven't yet, so, go for it! Should probably be a simple fix for them. $\endgroup$
    – Jesse
    Dec 5, 2023 at 18:47
0
$\begingroup$

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! :)

$\endgroup$
4
  • $\begingroup$ I don't think this is the case. This is a Minimal working example to show the problem. The real data is hundreds of species with longer sequences. When I run multiple times manually the clades in the tree change, but when I use biopython consensus tree it does not happen, so it must definitely be a bug or an error in my code. $\endgroup$
    – Mirko
    Dec 5, 2023 at 18:29
  • $\begingroup$ When you mean you ran the samples "multiple times manually" as opposed to "when I use the biopython consensus tree"? This sentence is confusing and misleading potentially. "Consensus trees are used to summarize a set of trees defined on the same set of taxa. A consensus algorithm takes the trees as input, so that the method of producing the input trees is not part of the consensus algorithm" as cited on: ncbi.nlm.nih.gov/pmc/articles/PMC2909780/…. $\endgroup$ Dec 8, 2023 at 14:00
  • $\begingroup$ A consensus tree uses the consensus algorithm to generate a tree. That makes me wonder if the algorithm you are using is appropriate for the amount of samples you have and the variation of species from different axons that you are using. What taxons are you using? Are you looking at family, class, order, genus, etc...? Family would defiantly be a broader range of values than genus. If you can maybe narrow down to a more exclusive taxon like order or genus, maybe your values may be more consistent? $\endgroup$ Dec 8, 2023 at 14:06
  • $\begingroup$ In the example i gave (that is completely sufficient to understand the problem), there is no mention of taxa. There is no need to invoke more complex processes when it is evident that when you run a consensus tree either on a very random sequence set (MWE) or in a very big dataset (my actual case) you do not expect to see 100% on all the nodes. It's like saying that you made the best possible alignment, also all the assumptions of your models are met and there is no measure error in the sequences and so on... $\endgroup$
    – Mirko
    Dec 8, 2023 at 17:49

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.