5
$\begingroup$

I am trying to assign character state changes from a presence-absence matrix to a phylogeny. Thus I want to identify the most parsimonious hypothesis for which node and branch a given "mutation" or "change" has occurred.

I have tried assigning each character to its leaf node, and then if the leaf node's sister has the same character I reassign the character to the parent node (and work back until all nodes are assigned). I am using a dummy dataset to try to achieve this:

Matrix
>Dme_001
1110000000000111
>Dme_002
1110000000000011
>Cfa_001
0110000000000011
>Mms_001
0110000000000011
>Hsa_001
0110000000000010
>Ptr_002
0110000000000011
>Mmu_002
0110000000000011
>Hsa_002
0110000000000011
>Ptr_001
0110000000000011
>Mmu_001
0110000000000011

Phylogeny
((Dme_001,Dme_002),(((Cfa_001,Mms_001),((Hsa_001,Ptr_001),Mmu_001)),(Ptr_002,(Hsa_002,Mmu_002))));

I assign internal nodes using ete3, so my output should be:

BranchID    CharacterState    Change
Node_1:    0    0->1
Hsa_001:    15    1->0 

As my code assigns character states based on their sisters if a loss is encountered it messes up the output so that:

BranchID    CharacterState    Change
Node_1:   0     0->1
Node_3    15    0->1
Node_5    15    0->1
Node_8    15    0->1

Could someone please help me with this? I'm coding in python 2.7 and developing tunnel vision. Thanks in advance

My code:

from ete3 import PhyloTree
from collections import Counter
import itertools

PAM = open('PAM','r')

gene_tree = '((Dme_001,Dme_002),(((Cfa_001,Mms_001),((Hsa_001,Ptr_001),Mmu_001)),(Ptr_002,(Hsa_002,Mmu_002))));'

NodeIDs = []

tree = PhyloTree(gene_tree)
edge = 0
for node in tree.traverse():
    if not node.is_leaf():
        node.name = "Node_%d" %edge
        edge +=1
        NodeIDs.append(node.name)
    if node.is_leaf():
        NodeIDs.append(node.name)

f = open('PAM','r')

taxa = []
pap = []

for line in f:
    term = line.strip().split('\t')
    taxa.append(term[0])
    p = [p for p in term[1]]
    pap.append(p)

statesD = dict(zip(taxa, pap))

def PlotCharacterStates():

    Plots = []

    events = []

    for key, value in statesD.iteritems():
        count = -1
        for s in value: 
            count+=1
            if s == CharacterState:
                a = key, count
                events.append(a)

    Round3_events = []
    while len(events) > 0:
        for rel in Relationships:
            node_store = []
            sis_store = []
            for event in events:
                if rel[0] == event[0]:
                    node_store.append(event[1])
                if rel[1] == event[0]:
                    sis_store.append(event[1])
            if (len(node_store) > 0) and (len(sis_store) > 0):
                place = rel, node_store, sis_store
                Round3_events.append(place)

        moved = []
        for placement in Round3_events:
            intercept = (set(placement[1]) & set(placement[2]))
            node_plot = (set(placement[1]) - set(placement[2]))
            sis_plot = (set(placement[2]) - set(placement[1]))
            if len(node_plot) > 0:
                for x in node_plot:
                    y = placement[0][0], x
                    Plots.append(y)
                    moved.append(y)
            if len(sis_plot) > 0:
                for x in sis_plot:
                    y = placement[0][1], x
                    Plots.append(y)
                    moved.append(y)
            if len(intercept) > 0:
                for x in intercept:
                    y = placement[0][2], x
                    y1 = placement[0][0], x
                    y2 = placement[0][1], x
                    moved.append(y1)
                    moved.append(y2)
                    events.append(y)

        for event in events:
            if event[0] == "Node_0":
                Plots.append(event)
                moved.append(event)

        events2 = (set(events) - set(moved))
        events = []
        for event in events2:
            events.append(event)


    pl = set(Plots)
    Plots = []
    for p in pl:
        Plots.append(p)

    print CharacterState, Plots


'''
assign sisters to leaves, internals
'''

e = []
round1b_e = []
round2a_e = []
placements = []
Relationships = []
Rounds = []
for node in tree.traverse():
    sisters = node.get_sisters()
    parent = node.up
    cycle1 = []
    if node.is_leaf():
        for sister in sisters:
            if sister.is_leaf():
                round1a = ["Round1a", node.name, sister.name, parent.name]
                node_names = node.name, sister.name
                Rounds.append(round1a)
                e.append(node_names)
                x = node.name, sister.name, parent.name, "leaf-leaf"
                Relationships.append(x)
            if not sister.is_leaf():
                round1b =  ["Round1b", node.name, sister.name, parent.name]
                node_names = node.name, sister.name
                Rounds.append(round1b)
                round1b_e.append(node_names)
                x = node.name, sister.name, parent.name, "node-leaf"
                Relationships.append(x)
    elif not node.is_leaf():
        if not node.is_root():
            for sister in sisters:
                if not sister.is_leaf():
                    node_names = node.name, sister.name
                    round2a_e.append(node_names)
                    x = node.name, sister.name, parent.name, "node-node"
                    Relationships.append(x)

x = []
CharacterStates = []                
for key, value in statesD.iteritems():
    for value in value:
        x.append(value)

y = sorted(set(x))
for x in y:
    CharacterStates.append(x)

for CharacterState in CharacterStates:
    PlotCharacterStates()
$\endgroup$
  • $\begingroup$ My code is added, thank you for the suggestion. I am using python 2.7, and I am not yet using biopython but happy to use it if there is a solution. I am looking to replicate the -apo function from TnT (Goloboff 2008) if possible (if you are familiar with this) $\endgroup$ – Gloom Oct 8 '18 at 13:57
  • $\begingroup$ Please edit the question to include this information, it might be helpful for others (I am not familiar with these area, sorry). $\endgroup$ – llrs Oct 8 '18 at 14:13
1
$\begingroup$

To give some formal words for what you are doing, I would say this sounds like you are trying to reconstruct ancestral states. There has been some work done on this, and a similar algorithm to yours was devised by Fitch in the 1970s: Fitch W: Toward defining the course of evolution: minimum change for a specific tree topology. Syst Zool 1971, 20: 406–416. 10.2307/2412116

The missing part, from your approach, is that once Fitch gets to the root of the tree, he than retraverses back to the leaves, resolving ambiguities along the way (if two sister nodes have different states, how do you decide what the parental state was? Fitch's idea is that once you get to the root, you can come back down and use the parental state to set the child state in ambiguous cases).

There is a nice implementation of Fitch's algorithm in the R phangorn package, which makes the reconstruction a one-liner (given a rooted tree object, and a matrix of leaf node states):

anc.mpr = ancestral.pars(rooted_tree, tree_samples_data_matrix.phyDat, "MPR")

If you don't want to use the R package itself, you could perhaps use it as pseudocode to port into python.

$\endgroup$
0
$\begingroup$

The OP could use the Rpy2 library to connect Python to R, but it is often better to use R directly.

The alternative is to simply use MacClade 4.0, although you will need an old Mac. This is a GUI interface that will allow you to identfy the parsimonious character state for each node in the phylogeny. The API allows you to interactively explore the phylogeny by point and click moving the branches around the tree. The tree will then be assessed in terms of "number character steps" and automatically calculate this for each topology loaded. The output is very colourful so its easy to understand a large data set.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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