def limb(D, j, n):
min_i = None
min_k = None
min_length = float('inf')
for i in range(n):
for k in range(n):
if i != j and k != j:
length = (D[j][k] + D[i][j] - D[i][k]) / 2
if length < min_length:
min_length = length
min_i = i
min_k = k
return (min_i, min_k, min_length)
i = 0
j = 1
k = 2
l = 3
limb(D, l, 3)
G = {
'a': [('b', 3)],
'b': [('a', 3), ('c', 2), ('d', 1)],
'c': [('b', 2)],
'd': [('b', 1)]
}
def get_neighbors(G, v):
if v in G:
return G[v]
else:
[]
get_neighbors(G, 'a')
def find_path(G, source, destination):
stack = [source]
visited = set([source])
while len(stack) > 0:
v = stack[-1]
if v == destination:
return stack
has_neighbors = False
for (w, weight) in get_neighbors(G, v):
if w not in visited:
stack.append(w)
visited.add(w)
has_neighbors = True
break
if not has_neighbors:
stack.pop()
print('Path not found')
return []
find_path(G, 'a', 'c')
def distance_between_nodes(G, node_i, node_j):
neighbors = get_neighbors(G, node_i)
for (w, weight) in neighbors:
if w == node_j:
return weight
return None
distance_between_nodes(G, 'a', 'b')
def add_node_on_path(G, source, destination, distance):
path = find_path(G, source, destination)
i = 0
j = 1
node_i = path[i]
node_j = path[j]
current_distance = distance_between_nodes(G, node_i, node_j)
while current_distance < distance:
i += 1
j += 1
node_i = path[i]
node_j = path[j]
current_distance += distance_between_nodes(G, node_i, node_j)
if current_distance == distance:
return node_j
else:
distance_between = current_distance - distance
return new_node_between_nodes(G, node_j, node_i, distance_between)
def add_neighbor(G, node, neighbor, distance):
G[node].append((neighbor, distance))
if neighbor not in G:
G[neighbor] = []
G[neighbor].append((node, distance))
def remove_neighbor(G, node, neighbor, distance):
G[node].remove((neighbor, distance))
G[neighbor].remove((node, distance))
def new_node_between_nodes(G, node_i, node_j, distance_between):
new_node = 'X{}{}'.format(node_i, node_j)
distance_ij = distance_between_nodes(G, node_i, node_j)
distance_diff = distance_ij - distance_between
remove_neighbor(G, node_i, node_j, distance_ij)
add_neighbor(G, node_i, new_node, distance_between)
add_neighbor(G, new_node, node_j, distance_diff)
return new_node
G = {
'a': [('b', 3)],
'b': [('a', 3), ('c', 2), ('d', 1)],
'c': [('b', 2)],
'd': [('b', 1)]
}
new_node = new_node_between_nodes(G, 'b', 'c', 1.5)
print(new_node)
print(G)
def additive_phylogeny(D, n):
if n == 2:
return {
0: [(1, D[0][1])],
1: [(0, D[0][1])]
}
(i, k, limb_length) = limb(D, n - 1, n)
for j in range(n - 1):
D[j][n - 1] -= limb_length
D[n - 1][j] = D[j][n - 1]
x = D[i][n - 1]
T = additive_phylogeny(D, n - 1)
v = add_node_on_path(T, i, k, x)
add_neighbor(T, v, n-1, limb_length)
return T
D = [[ 0, 13, 21, 22],
[13, 0, 12, 13],
[21, 12, 0, 13],
[22, 13, 13, 0]]
n = 4
additive_phylogeny(D, n)
class Cluster:
def __init__(self, elements):
self.elements = elements
self.age = 0
self.left = None
self.right = None
self.label = str(elements)
def add_children(self, left, right):
self.left = left
self.right = right
def update_age(self, age):
self.age = age
def __str__(self):
return 'Label: {}, Age: {}'.format(self.label, self.age)
c1 = Cluster([1,2])
c2 = Cluster([3,4])
print(c1)
print(c2)
c1.update_age(10)
print(c1)
c3 = Cluster([])
c3.add_children(c1, c2)
print('Left: ')
print(c3.left)
print('Right: ')
print(c3.right)
def distance(D, cluster_1, cluster_2):
d = 0
n = len(cluster_1.elements)
m = len(cluster_2.elements)
for i in cluster_1.elements:
for j in cluster_2.elements:
d += D[i][j]
d = d / (n * m)
return d
def two_closest(clusters, D):
min_i = None
min_j = None
min_distance = float('inf')
for cluster_i in clusters:
for cluster_j in clusters:
if cluster_i != cluster_j:
current_distance = distance(D, cluster_i, cluster_j)
if current_distance < min_distance:
min_distance = current_distance
min_i = cluster_i
min_j = cluster_j
return (min_i, min_j, min_distance)
def create_clusters(D):
n = len(D)
clusters = []
for i in range(n):
clusters.append(Cluster([i]))
return clusters
clusters = create_clusters(D)
for c in clusters:
print(c)
print()
print(D)
print()
(c1, c2, dist) = two_closest(clusters, D)
print('C1: ')
print(c1)
print('C2: ')
print(c2)
print('Distance: {}'.format(dist))
def UPGMA(D, n):
clusters = create_clusters(D)
T = {}
for c in clusters:
T[c] = []
while len(clusters) > 1:
(c_i, c_j, dist) = two_closest(clusters, D)
c_new = Cluster(c_i.elements + c_j.elements)
c_new.left = c_i
c_new.right = c_j
c_new.update_age(dist / 2)
T[c_new] = [(c_i, dist/2), (c_j, dist/2)]
new_clusters = [c_new]
for c in clusters:
if c != c_i and c != c_j:
new_clusters.append(c)
clusters = new_clusters
root = clusters[0]
return (T, root)
D = [[ 0, 13, 21, 22],
[13, 0, 12, 13],
[21, 12, 0, 13],
[22, 13, 13, 0]]
n = 4
graph, root = UPGMA(D, n)
#print(graph)
for v, neighbors in graph.items():
print(v)
print('Neighbors: ')
for w in graph[v]:
print('{} '.format(w))
def create_map(D):
D_map = {}
n = len(D)
for i in range(n):
D_map[i] = {}
for j in range(n):
D_map[i][j] = D[i][j]
return D_map
D_map = create_map(D)
D_map['a'] = {}
D_map['a']['b'] = 3
print(D_map)
def add_node(D_map, m, i, j):
nodes = [x[0] for x in list(D_map.items())]
D_map[m] = {}
for k in nodes:
if k != i and k != j:
D_map[k][m] = 1/2 * (D_map[k][i] + D_map[k][j] - D_map[i][j])
D_map[m][k] = D_map[k][m]
D_map[m][m] = 0
return D_map
def remove_node(D_map, removed_node):
del D_map[removed_node]
for node in D_map:
if removed_node in D_map[node]:
del D_map[node][removed_node]
return D_map
D_map = create_map(D)
print(D_map)
remove_node(D_map, 3)
def total_distance(D_map, node_i):
total = 0
for node_j in D_map[node_i]:
total += D_map[node_i][node_j]
return total
D_map = create_map(D)
print(D_map)
total_distance(D_map, 0)
import copy
def create_map_star(D_map):
D_map_star = copy.deepcopy(D_map)
n = len(D_map)
for i in D_map:
for j in D_map[i]:
if i != j:
D_map_star[i][j] = (n - 2) * D_map[i][j] - total_distance(D_map, i) - total_distance(D_map, j)
else:
D_map_star[i][j] = 0
return D_map_star
print(D)
D_map = create_map(D)
create_map_star(D_map)
def minimal_pair(D_star):
min_distance = float('inf')
min_i = None
min_j = None
for i in D_star:
for j in D_star[i]:
if i != j:
current_distance = D_star[i][j]
if current_distance < min_distance:
min_i = i
min_j = j
min_distance = current_distance
return (min_i, min_j)
def neighbor_joining(D, n):
if n == 2:
(i, j) = [x[0] for x in list(D.items())]
return {
i: [(j, D[i][j])],
j: [(i, D[j][i])]
}
D_star = create_map_star(D)
(i, j) = minimal_pair(D_star)
delta = (total_distance(D, i) - total_distance(D, j)) / (n - 2)
limb_length_i = 1/2 * (D[i][j] + delta)
limb_length_j = 1/2 * (D[i][j] - delta)
m = '{}+{}'.format(i, j)
D = add_node(D, m, i, j)
D = remove_node(D, i)
D = remove_node(D, j)
T = neighbor_joining(D, n - 1)
T[m].append((i, limb_length_i))
T[m].append((j, limb_length_j))
T[i] = [(m, limb_length_i)]
T[j] = [(m, limb_length_j)]
return T
D_map = create_map(D)
n = len(D_map)
T = neighbor_joining(D_map, n)
print(T)