How to plot character state changes from a presence-absence matrix on to a phylogeny

Question

I wish 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

I'm coding in python 2.7 and welcome assistance.

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()

Answer 1

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.

Answer 2

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.