How can I apply proportional (p) distances (Nucleotide) using bioPython

Question

I am working on phylogenetic tree generation. I used BioPython and ClustalW 2.x.x for this purpose. I have generate the tree using BioPython but when I try to generate tree using "MegaSoftware GUI" tool my tree does not matched with the output of MegaSoftware.

When I use "No. of difference" model then results are matched with biopython tree Like below image shown.

But whenever I choose "p-distance" model my results are not matched.

Output MegaSoftware using p-distance

Output BioPython

Here is my code for generating phylogenetic tree:

def conversion_fasta_dnd(clustalw, filepath, type, clustering, output, pwmatrix):
    """
        ClustalW tool to conversion the file format.
        :param clustalw: path of the executable programm.
        :param filepath: path of `fasta` extension file.
        :param type: dna, protein
        :param clustering: upgma, nj
        :param output: clustal, phylip, fasta
        :param pwmatrix: blosum, pam, id
        :return: dnd, aln (files)
        """
    # esblish the full path to the clustalW program
    clustalw_exe = clustalw

    # check to make sure the program is there
    assert os.path.isfile(clustalw_exe), "Clustal W executable missing"

    # create the appropriate command line
    clustalw_cline = ClustalwCommandline(clustalw_exe,
                                         infile=filepath,
                                         type=type,
                                         clustering=clustering,
                                         output=output,
                                         pwmatrix=pwmatrix)

    # print the command line
    print(clustalw_cline)

    # execute the command
    stdout, stderr = clustalw_cline()

# This function draw the tree on image file
def simple_tree(filepath):
    """
    Generate the simple tree image.
    :param filepath: path of dnd file format.
    """
    def get_label(leaf):
        return leaf.name

    file = filepath
    tree = Phylo.read(file, 'newick')
    Phylo.draw(tree, label_func=get_label, do_show=False, branch_labels=lambda c: c.branch_length)

    # axis management
    pylab.axis('on')
    pylab.savefig(OUTPUT_IMAGE, format='png', bbox_inches='tight', dpi=100, bbox_extra_artists=[])

In above code I have used two things first BioPhython lib and the other is Clustal W for align the sequences.

Could you please clarify what p-distances are?
I only have these reference for p-distance:

• p-distance (Nucleotide).
• Proportional (p) distances

I need to add p-distance model in my code, someone please help me out How can I do that.

Answer 1

Based on the sites you mention, it seems that p-distances are pairwise distances between aligned sequences obtained by counting the number of differences and dividing by the length of the alignment. This is the proportion of differences between two sequences.

I'm not sure exactly what you mean by "add p-distance model in my code", but here is how you could compute those distances from an existing alignment (here I use as an example an alignment I have in fasta format):

from Bio import AlignIO

def p_dist(seq1, seq2):
    ali_len = len(seq1)
    assert ali_len == len(seq2)
    return sum(0 if nuc1 == nuc2 else 1 for (nuc1, nuc2) in zip(seq1, seq2)) / ali_len

ali = AlignIO.read("most_targeted_histones/transcript_alignments/H4_top4_loss_spliced_NT.fa", "fasta")
nb_seq = len(ali)
distances = {}
for i in range(nb_seq - 1):
    rec1 = ali[i]
    j = i + 1
    while j < nb_seq:
        rec2 = ali[j]
        distances[(rec1.id, rec2.id)] = p_dist(rec1.seq, rec2.seq)
        j += 1

print(*distances.items(), sep="\n")

Output:

(('ZK131.4_spliced', 'ZK131.8_spliced'), 0.0)
(('ZK131.4_spliced', 'T10C6.14_spliced'), 0.10897435897435898)
(('ZK131.4_spliced', 'ZK131.1_spliced'), 0.0)
(('ZK131.8_spliced', 'T10C6.14_spliced'), 0.10897435897435898)
(('ZK131.8_spliced', 'ZK131.1_spliced'), 0.0)
(('T10C6.14_spliced', 'ZK131.1_spliced'), 0.10897435897435898)

If you want to build a tree from those distances, I'm not sure you can do this with BioPython. It seems however that there are some possibilities with DendroPy.

You would need to write the distances in a suitable format, and to choose the method between UPGMA and neighbour-joining.