7987035156

Committed 21 Feb 2024 09:47AM UTC coverage: 46.646% (+0.4%) from 46.211%

Build # 7987035156

Build Type

push

github

Committed by

js51

Commit Message

add taxa ordering and an Alignment class for keeping track of it

Run Details

21 of 45 new or added lines in 7 files covered. (46.67%)

6 existing lines in 3 files now uncovered.

445 of 954 relevant lines covered (46.65%)

1.4 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

13.98

/splitp/simulation.py

import itertools
import numpy as np
from warnings import warn
from random import choices
from networkx import dfs_tree
from splitp import constants, alignment


def evolve_pattern(tree, model=None):
    def __evolve_on_subtree(subtree, state):
        root_node = [n for n, d in subtree.in_degree() if d == 0][0]
        children = list(subtree.successors(root_node))
        if model == None: # Get the transition matrix from the phylogeny
                M = (tree.networkx_graph.nodes[root_node])["transition_matrix"]
        else:
            branch_length = (tree.networkx_graph.nodes[root_node])["branch_length"]
            M = model.transition_matrix(branch_length)
        probs = list(M[:, constants.DNA_state_space.index(state)])
        new_state = choices(constants.DNA_state_space, probs)[0]
        if len(children) == 0:
            return f"{str(root_node)}:{new_state}"
        else:
            subtrees = [dfs_tree(tree.networkx_graph, child) for child in children]
            return ",".join(
                __evolve_on_subtree(subtree, new_state) for subtree in subtrees
            )

    root_state = choices(constants.DNA_state_space, (1/4, 1/4, 1/4, 1/4))[0]
    root_node = [n for n, d in tree.networkx_graph.in_degree() if d == 0][0]
    children = list(tree.networkx_graph.successors(root_node))
    subtrees = [dfs_tree(tree.networkx_graph, child) for child in children]
    result_string = ",".join(
        __evolve_on_subtree(subtree, root_state) for subtree in subtrees
    )
    result = {
        pair.split(":")[0]: pair.split(":")[1] for pair in result_string.split(",")
    }
    taxa = tree.taxa
    return "".join(result[k] for k in sorted(result.keys(), key=taxa.index))


def generate_alignment(tree, model, sequence_length):
    counts = {}
    for i in range(sequence_length):
        pattern = evolve_pattern(tree, model)
        if pattern not in counts:
            counts[pattern] = float(1)
        else:
            counts[pattern] += 1
    probs = {}
    for k in sorted(
        counts.keys(), key=lambda p: [constants.DNA_state_space.index(c) for c in p]
    ):
        probs[k] = counts[k] / float(sequence_length)
    
    return probs


def draw_from_multinomial(pattern_probabilities, n):
    """Use a given table of probabilities from get_pattern_probabilities() and draw from its distribution"""
    patterns, probs = zip(*pattern_probabilities.items())
    data = np.random.multinomial(n, probs)
    data = [d / n for d in data]
    results = {}
    for p, pattern in enumerate(patterns):
        if data[p] != 0:
            results[pattern] = data[p]
    return results


def get_pattern_probabilities(tree, model=None):
    """Returns a full table of site-pattern probabilities (binary character set)"""
    combinations = list(
        itertools.product(
            "".join(s for s in constants.DNA_state_space), repeat=tree.get_num_taxa()
        )
    )
    combinations = ["".join(c) for c in combinations]
    emptyArray = {
        combination: __likelihood_start(tree, combination, model) for combination in combinations
    }
    pattern_probs = alignment.Alignment(emptyArray, taxa=tree.taxa)
    return emptyArray

pattern_probabilities = get_pattern_probabilities


def __likelihood(tree, n, likelihood_table, model):
    """Recursive part of the likelihood algorithm"""
    for b in range(4):
        L = 1
        for d in tree.get_descendants(n, return_iter=True):
            d_index = tree.node_index(d)
            if model == None: # Get the transition matrix from the phylogeny
                M = (tree.networkx_graph.nodes[d])["transition_matrix"]
            else:
                branch_length = (tree.networkx_graph.nodes[d])["branch_length"]
                M = model.transition_matrix(branch_length)
            s = 0
            for b2 in range(4):
                if likelihood_table[d_index, b2] == None:
                    __likelihood(tree, d, likelihood_table, model)
                s += M[b2, b] * likelihood_table[d_index, b2]
            L *= s
        likelihood_table[tree.node_index(n), b] = L
    if not tree.is_root(n):
        __likelihood(tree, tree.get_parent(n), likelihood_table, model)


def __likelihood_start(tree, pattern, model):
    """Starts the likelihood algorithm.

    Starts the likelihood algorithm to determine the probability of seeing a given site pattern.

    Args:
        pattern: The site pattern to obtain the probability for.

    Returns:
        A probability (range 0-1) of observing the given site pattern given the tree.
    """

    def _to_int(p):
        return constants.DNA_state_space.index(p)

    pattern = [
        _to_int(p) for p in pattern
    ]  # A list of indices which correspond to taxa.
    taxa = tree.taxa # The list of taxa.
    # Likelihood table for dynamic prog alg ~ lTable[node_index, character]
    likelihood_table = np.array(
        [[None for _ in range(4)] for _ in range(tree.get_num_nodes())]
    )

    # Encoding pattern into likelihood table
    for i, p in enumerate(pattern):
        likelihood_table[tree.node_index(taxa[i]), :] = [
            0 for _ in range(4)
        ]
        likelihood_table[tree.node_index(taxa[i]), p] = 1

    # Starting with node which has all descendants as leaves
    initNode = None
    for n in tree.networkx_graph.nodes:
        desc = tree.get_descendants(n)
        if desc and all(
            tree.is_leaf(d) for d in desc
        ):  # If the node has descendants and they are all leaves
            initNode = n

    __likelihood(tree, initNode, likelihood_table, model)  # Should fill in the entire table

    root_index = tree.node_index(tree.root(return_index=False))
    L = 0
    for i in range(4):
        L += (1/4) * likelihood_table[root_index, i]

    return L

1	import itertools	3✔
2	import numpy as np	3✔
3	from warnings import warn	3✔
4	from random import choices	3✔
5	from networkx import dfs_tree	3✔
6	from splitp import constants, alignment	3✔
7
8
9	def evolve_pattern(tree, model=None):	3✔
10	def __evolve_on_subtree(subtree, state):	×
11	root_node = [n for n, d in subtree.in_degree() if d == 0][0]	×
12	children = list(subtree.successors(root_node))	×
13	if model == None: # Get the transition matrix from the phylogeny	×
14	M = (tree.networkx_graph.nodes[root_node])["transition_matrix"]	×
15	else:
16	branch_length = (tree.networkx_graph.nodes[root_node])["branch_length"]	×
17	M = model.transition_matrix(branch_length)	×
18	probs = list(M[:, constants.DNA_state_space.index(state)])	×
19	new_state = choices(constants.DNA_state_space, probs)[0]	×
20	if len(children) == 0:	×
21	return f"{str(root_node)}:{new_state}"	×
22	else:
23	subtrees = [dfs_tree(tree.networkx_graph, child) for child in children]	×
24	return ",".join(	×
25	__evolve_on_subtree(subtree, new_state) for subtree in subtrees
26	)
27
28	root_state = choices(constants.DNA_state_space, (1/4, 1/4, 1/4, 1/4))[0]	×
29	root_node = [n for n, d in tree.networkx_graph.in_degree() if d == 0][0]	×
30	children = list(tree.networkx_graph.successors(root_node))	×
31	subtrees = [dfs_tree(tree.networkx_graph, child) for child in children]	×
32	result_string = ",".join(	×
33	__evolve_on_subtree(subtree, root_state) for subtree in subtrees
34	)
35	result = {	×
36	pair.split(":")[0]: pair.split(":")[1] for pair in result_string.split(",")
37	}
NEW 38	taxa = tree.taxa	×
39	return "".join(result[k] for k in sorted(result.keys(), key=taxa.index))	×
40
41
42	def generate_alignment(tree, model, sequence_length):	3✔
43	counts = {}	×
44	for i in range(sequence_length):	×
45	pattern = evolve_pattern(tree, model)	×
46	if pattern not in counts:	×
47	counts[pattern] = float(1)	×
48	else:
49	counts[pattern] += 1	×
50	probs = {}	×
51	for k in sorted(	×
52	counts.keys(), key=lambda p: [constants.DNA_state_space.index(c) for c in p]
53	):
54	probs[k] = counts[k] / float(sequence_length)	×
55
UNCOV 56	return probs	×
57
58
59	def draw_from_multinomial(pattern_probabilities, n):	3✔
60	"""Use a given table of probabilities from get_pattern_probabilities() and draw from its distribution"""
61	patterns, probs = zip(*pattern_probabilities.items())	×
62	data = np.random.multinomial(n, probs)	×
63	data = [d / n for d in data]	×
64	results = {}	×
65	for p, pattern in enumerate(patterns):	×
66	if data[p] != 0:	×
67	results[pattern] = data[p]	×
68	return results	×
69
70
71	def get_pattern_probabilities(tree, model=None):	3✔
72	"""Returns a full table of site-pattern probabilities (binary character set)"""
UNCOV 73	combinations = list(	×
74	itertools.product(
75	"".join(s for s in constants.DNA_state_space), repeat=tree.get_num_taxa()
76	)
77	)
78	combinations = ["".join(c) for c in combinations]	×
79	emptyArray = {	×
80	combination: __likelihood_start(tree, combination, model) for combination in combinations
81	}
NEW 82	pattern_probs = alignment.Alignment(emptyArray, taxa=tree.taxa)	×
UNCOV 83	return emptyArray	×
84
85	pattern_probabilities = get_pattern_probabilities	3✔
86
87
88	def __likelihood(tree, n, likelihood_table, model):	3✔
89	"""Recursive part of the likelihood algorithm"""
90	for b in range(4):	×
91	L = 1	×
92	for d in tree.get_descendants(n, return_iter=True):	×
93	d_index = tree.node_index(d)	×
94	if model == None: # Get the transition matrix from the phylogeny	×
95	M = (tree.networkx_graph.nodes[d])["transition_matrix"]	×
96	else:
97	branch_length = (tree.networkx_graph.nodes[d])["branch_length"]	×
98	M = model.transition_matrix(branch_length)	×
99	s = 0	×
100	for b2 in range(4):	×
101	if likelihood_table[d_index, b2] == None:	×
102	__likelihood(tree, d, likelihood_table, model)	×
103	s += M[b2, b] * likelihood_table[d_index, b2]	×
104	L *= s	×
105	likelihood_table[tree.node_index(n), b] = L	×
106	if not tree.is_root(n):	×
107	__likelihood(tree, tree.get_parent(n), likelihood_table, model)	×
108
109
110	def __likelihood_start(tree, pattern, model):	3✔
111	"""Starts the likelihood algorithm.
112
113	Starts the likelihood algorithm to determine the probability of seeing a given site pattern.
114
115	Args:
116	pattern: The site pattern to obtain the probability for.
117
118	Returns:
119	A probability (range 0-1) of observing the given site pattern given the tree.
120	"""
121
122	def _to_int(p):	×
123	return constants.DNA_state_space.index(p)	×
124
125	pattern = [	×
126	_to_int(p) for p in pattern
127	] # A list of indices which correspond to taxa.
NEW 128	taxa = tree.taxa # The list of taxa.	×
129	# Likelihood table for dynamic prog alg ~ lTable[node_index, character]
130	likelihood_table = np.array(	×
131	[[None for _ in range(4)] for _ in range(tree.get_num_nodes())]
132	)
133
134	# Encoding pattern into likelihood table
135	for i, p in enumerate(pattern):	×
136	likelihood_table[tree.node_index(taxa[i]), :] = [	×
137	0 for _ in range(4)
138	]
139	likelihood_table[tree.node_index(taxa[i]), p] = 1	×
140
141	# Starting with node which has all descendants as leaves
142	initNode = None	×
143	for n in tree.networkx_graph.nodes:	×
144	desc = tree.get_descendants(n)	×
145	if desc and all(	×
146	tree.is_leaf(d) for d in desc
147	): # If the node has descendants and they are all leaves
148	initNode = n	×
149
150	__likelihood(tree, initNode, likelihood_table, model) # Should fill in the entire table	×
151
152	root_index = tree.node_index(tree.root(return_index=False))	×
153	L = 0	×
154	for i in range(4):	×
155	L += (1/4) * likelihood_table[root_index, i]	×
156
157	return L	×

js51 / SplitP / 7987035156

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous