1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
|
#! /usr/bin/env python
##############################################################################
## DendroPy Phylogenetic Computing Library.
##
## Copyright 2010-2015 Jeet Sukumaran and Mark T. Holder.
## All rights reserved.
##
## See "LICENSE.rst" for terms and conditions of usage.
##
## If you use this work or any portion thereof in published work,
## please cite it as:
##
## Sukumaran, J. and M. T. Holder. 2010. DendroPy: a Python library
## for phylogenetic computing. Bioinformatics 26: 1569-1571.
##
##############################################################################
import math
import unittest
import collections
import dendropy
from dendropy.calculate import treecompare
from dendropy.calculate.treecompare import TreeShapeKernel
from dendropy.calculate.treecompare import AssemblageInducedTreeManager
from dendropy.calculate.treecompare import AssemblageInducedTreeShapeKernel
class TreeShapeKernelBasicCalculationTest(unittest.TestCase):
def test_small_single(self):
## Test based on:
## KAMPHIR
## By Art F.Y. Poon and Rosemary McCloskey
## https://github.com/ArtPoon/kamphir.git
T1_str = "[&R] ( ( A:0.5, B:0.25 )E:0.5, ( C:0.25, D:0.25 )F:0.5 )G;"
T2_str = "[&R] ( ( ( A:0.25, B:0.25 )E:0.5, C:0.25 )F:0.5, D:0.25 )G;"
taxon_namespace = dendropy.TaxonNamespace()
T1 = dendropy.Tree.get_from_string(T1_str, "newick", taxon_namespace=taxon_namespace)
T2 = dendropy.Tree.get_from_string(T2_str, "newick", taxon_namespace=taxon_namespace)
tree_shape_kernel = TreeShapeKernel(decay_factor=0.5, gauss_factor=1)
assert tree_shape_kernel(T1, T2) == 1.125 * (1+math.exp(-0.0625))
def test_different_size_trees(self):
# sigma = 1
# gaussFactor = 1
# decayFactor = 0.1
trees_str = """\
[&R] (T1:12.5150033896,((T2:9.56663932282,(((T3:1.07963880769,T4:1.07963880769):4.33934073323,(T5:0.649236454849,T6:0.649236454849):4.76974308607):1.9951772698,T7:7.41415681072):2.1524825121):0.294448382741,T8:9.86108770556):2.65391568409):1.88016636219;
[&R] ((((T1:2.63300806671,(T2:1.36260395686,T3:1.36260395686):1.27040410984):13.0805656885,(((T4:6.49642565476,(T5:4.76391492099,(T6:2.31456992456,(T7:1.43347670305,(T8:0.845099591454,(T9:0.749944158281,T10:0.749944158281):0.0951554331727):0.588377111599):0.881093221511):2.44934499642):1.73251073377):1.78625805849,(T11:1.99279678602,T12:1.99279678602):6.28988692723):5.93122188014,T13:14.2139055934):1.49966816178):9.5804936484,(((((((T14:0.205658120519,T15:0.205658120519):7.75482186971,T16:7.96047999023):0.460360751119,(T17:2.05148090483,T18:2.05148090483):6.36935983652):7.37034613696,T19:15.7911868783):0.64452345349,T20:16.4357103318):3.12259410597,(T21:13.3918100022,T22:13.3918100022):6.16649443557):3.77060549196,(((((T23:0.117188712769,T24:0.117188712769):0.350994673936,T25:0.468183386705):0.434444802184,T26:0.902628188889):15.8023737799,T27:16.7050019688):5.11792644629,(T28:18.0216678715,T29:18.0216678715):3.80126054359):1.50598151462):1.96515747383):2.27498485797,T30:27.5690522615):2.7054878136;
[&R] ((T1:18.4726052336,(T2:16.7820588123,T3:16.7820588123):1.69054642133):21.0361145203,(((T4:22.1454645535,T5:22.1454645535):1.47756990124,((((((T6:3.72910479068,T7:3.72910479068):0.635795000427,(T8:3.12137642206,T9:3.12137642206):1.24352336905):4.02980332076,T10:8.39470311187):0.544576038938,T11:8.93927915081):4.15398513837,((T12:1.29491896299,T13:1.29491896299):3.95016532419,(T14:3.55950898913,T15:3.55950898913):1.68557529805):7.848180002):6.79081967153,(T16:13.3961088469,((T17:7.5123548288,T18:7.5123548288):4.47081545633,(T19:5.42073528802,T20:5.42073528802):6.56243499711):1.41293856181):6.48797511377):3.73895049406):0.434178411336,(T21:6.05244405939,(T22:2.47879429586,T23:2.47879429586):3.57364976352):18.0047688067):15.4515068878):33.9171079393;
[&R] ((((T1:12.2171017101,(T2:5.80709236615,T3:5.80709236615):6.4100093439):0.163771078375,(T4:1.48920874144,T5:1.48920874144):10.891664047):6.14111109308,((T6:18.0816707702,T7:18.0816707702):0.315080693117,(T8:9.85665216391,((T9:3.58056555459,T10:3.58056555459):0.704080460611,T11:4.2846460152):5.57200614872):8.54009929938):0.125232418213):0.445776242893,(((T12:10.3677145836,(T13:1.58010326492,T14:1.58010326492):8.78761131866):2.28511142688,T15:12.6528260105):1.47666702282,(T16:12.4604750725,(((T17:1.66562876613,T18:1.66562876613):2.60275247994,T19:4.26838124607):0.527073064007,(T20:0.771151073421,T21:0.771151073421):4.02430323665):7.66502076239):1.66901796081):4.83826709113):19.4421800373;
[&R] (((T1:10.2887625813,((T2:3.40673675066,T3:3.40673675066):4.35878457692,T4:7.76552132758):2.52324125374):8.94208354363,((((((T5:1.25729509129,T6:1.25729509129):1.98815007493,T7:3.24544516622):9.92816453819,T8:13.1736097044):0.131783344535,((T9:5.12529436775,T10:5.12529436775):2.8254015761,T11:7.95069594385):5.3546971051):0.31354548999,((T12:6.23707496437,T13:6.23707496437):1.40621592147,T14:7.64329088584):5.97564765309):2.09487953589,(((T15:4.48342349888,T16:4.48342349888):1.33113258373,T17:5.81455608261):0.722465080174,(T18:3.26170665225,T19:3.26170665225):3.27531451053):9.17679691204):3.51702805012):12.2978992269,(T20:26.4312399558,((((((T21:0.115067790962,T22:0.115067790962):0.614386167869,T23:0.729453958831):10.0469766167,(((T24:5.54290468454,T25:5.54290468454):0.24279629634,((T26:0.0173560007622,T27:0.0173560007622):1.03837416324,T28:1.055730164):4.72997081688):0.935147415819,T29:6.7208483967):4.05558217886):5.34149070025,T30:16.1179212758):0.249537503423,((T31:7.68676243321,T32:7.68676243321):0.224359722764,T33:7.91112215597):8.45633662326):4.68102050585,(((T34:6.33634296608,T35:6.33634296608):11.7427817803,T36:18.0791247464):1.88140825859,(T37:3.77606556293,T38:3.77606556293):16.1844674421):1.08794628008):5.38276067077):5.097505396):12.8756560997;
[&R] (((T1:21.6283980853,(T2:10.2523638386,T3:10.2523638386):11.3760342467):5.43593515808,(((T4:14.5815467291,((T5:10.2802760256,(T6:2.13975995776,T7:2.13975995776):8.14051606783):1.62385561861,T8:11.9041316442):2.67741508492):5.60705586175,(T9:5.28631833362,(T10:5.19018203656,((T11:0.654535739431,(T12:0.451288383503,T13:0.451288383503):0.203247355928):0.3240577623,T14:0.978593501732):4.21158853482):0.0961362970605):14.9022842572):6.30642464811,(T15:6.96742024513,T16:6.96742024513):19.5276069938):0.56930600442):9.39112113938,(((T17:9.86297817832,(T18:1.09548421955,(T19:0.163658446308,T20:0.163658446308):0.931825773239):8.76749395877):13.4127381553,(((((T21:0.162641367112,T22:0.162641367112):7.60612335657,(T23:1.19220280515,T24:1.19220280515):6.57656191854):4.64147681058,T25:12.4102415343):2.18305782204,(T26:13.4944558875,((T27:4.0278122205,(T28:3.83892525427,T29:3.83892525427):0.188886966234):6.19053229632,(T30:0.0720558467085,T31:0.0720558467085):10.1462886701):3.27611137066):1.09884346882):4.83967942511,(T32:13.108711714,(T33:1.21934185331,T34:1.21934185331):11.8893698607):6.32426706741):3.84273755221):12.2937144592,((T35:3.69731694953,(T36:1.08556862997,T37:1.08556862997):2.61174831956):2.01430134794,T38:5.71161829747):29.8578124954):0.886023589904):13.2943061023;
[&R] (((((T1:2.69139513211,T2:2.69139513211):8.551485817,T3:11.2428809491):6.39974203014,((T4:17.0587824705,(((((T5:1.0038538602,(T6:0.65693286312,(T7:0.372810170449,T8:0.372810170449):0.284122692671):0.346920997081):0.992649847145,T9:1.99650370735):8.59545004131,T10:10.5919537487):3.93618536755,(T11:11.4771452793,(T12:11.3117342188,(((T13:1.44496178243,T14:1.44496178243):1.39593941106,T15:2.84090119349):7.36176225558,T16:10.2026634491):1.10907076976):0.16541106051):3.05099383686):0.446923840522,((T17:2.07495063443,(T18:1.03526596889,T19:1.03526596889):1.03968466554):3.02573649405,(T20:3.09572655531,T21:3.09572655531):2.00496057316):9.87437582825):2.08371951376):0.535710348712,T22:17.5944928192):0.0481301600476):0.485204130118,(T23:2.10622863693,T24:2.10622863693):16.0215984724):8.64518545213,((T25:12.6080734092,((T26:8.70251718564,T27:8.70251718564):1.05101612319,((T28:6.04918606932,((T29:3.85573365482,T30:3.85573365482):0.0757038480577,T31:3.93143750288):2.11774856645):2.49018043616,(T32:4.24643709587,T33:4.24643709587):4.29292940962):1.21416680334):2.85454010035):3.79696048295,((((T34:2.04492296218,T35:2.04492296218):1.53035052582,T36:3.57527348801):0.0643892487425,T37:3.63966273675):12.6279014452,T38:16.267564182):0.137469710183):10.3679786694):5.11271917475;
[&R] ((((T1:7.07157972517,(T2:0.0218808463127,T3:0.0218808463127):7.04969887886):0.8655582016,(T4:3.02712438178,T5:3.02712438178):4.91001354499):1.21415733877,(((T6:0.0583489618229,T7:0.0583489618229):1.01650984984,T8:1.07485881166):4.07712736061,T9:5.15198617228):3.99930909326):7.26950003749,(((T10:0.424456263841,T11:0.424456263841):8.04936250082,T12:8.47381876466):3.5117568466,((T13:6.98368939693,(T14:0.777662333494,T15:0.777662333494):6.20602706344):2.06669325023,T16:9.05038264716):2.9351929641):4.43521969176):10.6934883223;
[&R] (((T1:14.4454814308,((T2:3.56688601323,T3:3.56688601323):5.84751651584,T4:9.41440252907):5.03107890172):20.1795615891,(((((T5:4.15054680055,T6:4.15054680055):11.295603376,(T7:6.43442739414,T8:6.43442739414):9.01172278238):0.10084264253,(T9:8.4004080966,T10:8.4004080966):7.14658472244):2.73330732663,(((T11:6.10577944378,(T12:4.04456805007,(T13:2.80049127557,T14:2.80049127557):1.2440767745):2.06121139371):1.90192017764,(T15:7.79084487337,T16:7.79084487337):0.216854748047):7.40306903527,T17:15.4107686567):2.86953148899):11.2997061766,(T18:1.00480866062,T19:1.00480866062):28.5751976616):5.04503669769):0.0473375886983,((((T20:5.77874649173,(T21:5.42355087634,(T22:3.52279709572,(T23:2.90310692644,(T24:2.09879096434,T25:2.09879096434):0.804315962107):0.619690169282):1.90075378062):0.355195615385):4.90915591989,(T26:4.01715333576,T27:4.01715333576):6.67074907586):3.94107350362,(T28:12.4346731779,(T29:3.0690908858,T30:3.0690908858):9.36558229213):2.1943027373):12.9196714424,(T31:14.8264557361,(T32:0.0945338917738,T33:0.0945338917738):14.7319218443):12.7221916216):7.123733251):12.8703596545;
"""
expected = [[ [] for i in range(10) ] for i in range(10)]
expected[0][0] = 1.02943649816
expected[0][1] = 0.629696880372
expected[0][2] = 0.226436166785
expected[0][3] = 0.504437998991
expected[0][4] = 0.532944026681
expected[0][5] = 1.06675009986
expected[0][6] = 0.63997331036
expected[0][7] = 0.533570048403
expected[0][8] = 0.316006296153
expected[0][9] = 0.763121914566
expected[1][0] = 0.629696880372
expected[1][1] = 5.73381840453
expected[1][2] = 0.609230566852
expected[1][3] = 1.36716715734
expected[1][4] = 1.66142414246
expected[1][5] = 3.3424585382
expected[1][6] = 2.81647328252
expected[1][7] = 1.6426813156
expected[1][8] = 1.71534697686
expected[1][9] = 2.11186097095
expected[2][0] = 0.226436166785
expected[2][1] = 0.609230566852
expected[2][2] = 3.41707388515
expected[2][3] = 0.965818487686
expected[2][4] = 1.8148659483
expected[2][5] = 1.09814847947
expected[2][6] = 1.70646748651
expected[2][7] = 0.543992787817
expected[2][8] = 1.89068146553
expected[2][9] = 1.2301930343
expected[3][0] = 0.504437998991
expected[3][1] = 1.36716715734
expected[3][2] = 0.965818487686
expected[3][3] = 3.52734267013
expected[3][4] = 1.6121784729
expected[3][5] = 1.9713400052
expected[3][6] = 1.968844455
expected[3][7] = 0.727319654436
expected[3][8] = 1.07080235987
expected[3][9] = 2.06317190944
expected[4][0] = 0.532944026681
expected[4][1] = 1.66142414246
expected[4][2] = 1.8148659483
expected[4][3] = 1.6121784729
expected[4][4] = 7.13717804604
expected[4][5] = 2.47204082769
expected[4][6] = 2.17747252512
expected[4][7] = 1.58536532878
expected[4][8] = 2.31434377631
expected[4][9] = 2.23786059877
expected[5][0] = 1.06675009986
expected[5][1] = 3.3424585382
expected[5][2] = 1.09814847947
expected[5][3] = 1.9713400052
expected[5][4] = 2.47204082769
expected[5][5] = 7.77913926686
expected[5][6] = 2.87296840413
expected[5][7] = 2.30072202059
expected[5][8] = 2.12831286543
expected[5][9] = 3.00420928696
expected[6][0] = 0.63997331036
expected[6][1] = 2.81647328252
expected[6][2] = 1.70646748651
expected[6][3] = 1.968844455
expected[6][4] = 2.17747252512
expected[6][5] = 2.87296840413
expected[6][6] = 7.32395549398
expected[6][7] = 1.3099686298
expected[6][8] = 2.72462663407
expected[6][9] = 2.60804726918
expected[7][0] = 0.533570048403
expected[7][1] = 1.6426813156
expected[7][2] = 0.543992787817
expected[7][3] = 0.727319654436
expected[7][4] = 1.58536532878
expected[7][5] = 2.30072202059
expected[7][6] = 1.3099686298
expected[7][7] = 3.00348381586
expected[7][8] = 1.05294859724
expected[7][9] = 1.30614199254
expected[8][0] = 0.316006296153
expected[8][1] = 1.71534697686
expected[8][2] = 1.89068146553
expected[8][3] = 1.07080235987
expected[8][4] = 2.31434377631
expected[8][5] = 2.12831286543
expected[8][6] = 2.72462663407
expected[8][7] = 1.05294859724
expected[8][8] = 6.04463164133
expected[8][9] = 2.07119628549
expected[9][0] = 0.763121914566
expected[9][1] = 2.11186097095
expected[9][2] = 1.2301930343
expected[9][3] = 2.06317190944
expected[9][4] = 2.23786059877
expected[9][5] = 3.00420928696
expected[9][6] = 2.60804726918
expected[9][7] = 1.30614199254
expected[9][8] = 2.07119628549
expected[9][9] = 4.7667655343
trees = dendropy.TreeList.get(data=trees_str, schema="newick")
tree_shape_kernel = TreeShapeKernel(
sigma=1,
gauss_factor=1,
decay_factor=0.1,
)
for idx1, t1 in enumerate(trees):
for idx2, t2 in enumerate(trees):
self.assertAlmostEqual(tree_shape_kernel(t1, t2), expected[idx1][idx2])
# print("{}, {} = {}".format(idx1+1, idx2+1, tree_shape_kernel(t1, t2)))
class AssemblageInducedTreeManagerTestBase(unittest.TestCase):
GROUP_IDS = ("a", "b", "c", "d", "e")
@staticmethod
def get_random_tree():
from dendropy.model import birthdeath
tns = dendropy.TaxonNamespace()
for group_id in AssemblageInducedTreeManagerTests.GROUP_IDS:
for group_member in range(10):
t = tns.require_taxon(label="{}{}".format(
group_id,
group_member))
tree = dendropy.simulate.birth_death_tree(
birth_rate=0.1,
death_rate=0.0,
taxon_namespace=tns)
tree.assemblage_leaf_sets = []
tree.assemblage_classification_regime_subtrees = []
for group_id in AssemblageInducedTreeManagerTests.GROUP_IDS:
node_filter_fn=lambda nd: nd.taxon is None or nd.taxon.label.startswith(group_id)
subtree1 = tree.extract_tree(
node_filter_fn=node_filter_fn)
assemblage_leaf_set = set()
for leaf_nd in tree.leaf_node_iter():
if node_filter_fn(leaf_nd):
assert leaf_nd.taxon.label.startswith(group_id)
assemblage_leaf_set.add(leaf_nd)
tree.assemblage_leaf_sets.append(assemblage_leaf_set)
tree.assemblage_classification_regime_subtrees.append(subtree1)
assert len(tree.assemblage_classification_regime_subtrees) == len(AssemblageInducedTreeManagerTests.GROUP_IDS)
return tree
def validate_managed_trees(self, test_target, trees):
self.assertEqual(test_target._num_assemblage_classifications, len(AssemblageInducedTreeManagerTests.GROUP_IDS))
self.assertEqual(len(test_target._tree_assemblage_induced_trees_map), len(trees))
for tree in trees:
self.assertIn(tree, test_target._tree_assemblage_induced_trees_map)
self.assertEqual(len(test_target._tree_assemblage_induced_trees_map[tree]), len(AssemblageInducedTreeManagerTests.GROUP_IDS))
self.assertEqual(len(test_target._tree_assemblage_induced_trees_map[tree]), len(tree.assemblage_leaf_sets))
induced_trees = test_target._tree_assemblage_induced_trees_map[tree]
for ( induced_tree, group_id, original_leafset_nodes) in zip( induced_trees, AssemblageInducedTreeManagerTests.GROUP_IDS, tree.assemblage_leaf_sets):
original_leafset = set(original_leafset_nodes)
for leaf_nd in induced_tree.leaf_node_iter():
self.assertTrue(leaf_nd.taxon.label.startswith(group_id), leaf_nd.taxon.label)
original_node = leaf_nd.extraction_source
self.assertIn(original_node, original_leafset)
original_leafset.remove(original_node)
self.assertEqual(len(original_leafset), 0)
labels=[x.taxon.label for x in original_leafset_nodes]
t2 = tree.extract_tree_with_taxa_labels(labels=labels)
self.assertEqual(treecompare.weighted_robinson_foulds_distance(t2, induced_tree), 0.0)
t3 = dendropy.Tree(tree)
t3.retain_taxa_with_labels(labels=labels)
# print(t3.as_string("newick"))
# print(induced_tree.as_string("newick"))
self.assertAlmostEqual(treecompare.weighted_robinson_foulds_distance(t3, induced_tree), 0.0)
class AssemblageInducedTreeManagerTests(AssemblageInducedTreeManagerTestBase):
def test_basic_subtree_generation_and_caching(self):
test_target = AssemblageInducedTreeManager()
trees = []
for idx in range(10):
tree = self.get_random_tree()
tree.induced_trees = test_target.generate_induced_trees(
tree=tree,
assemblage_leaf_sets=tree.assemblage_leaf_sets)
self.assertEqual(len(tree.induced_trees), len(AssemblageInducedTreeManagerTests.GROUP_IDS))
self.assertEqual(len(tree.induced_trees), len(tree.assemblage_leaf_sets))
for induced_tree in tree.induced_trees:
self.assertIn(induced_tree, test_target._tree_assemblage_induced_trees_map[tree])
trees.append(tree)
self.validate_managed_trees(
test_target=test_target,
trees=trees)
class AssemblageInducedTreeShapeKernelTests(AssemblageInducedTreeManagerTestBase):
def validate_managed_and_cached_trees(self, test_target, trees):
self.validate_managed_trees(
test_target=test_target,
trees=trees)
for tree in trees:
self.assertIn(tree, test_target._tree_cache)
self.assertEqual(len(test_target._tree_assemblage_induced_trees_map[tree]), len(AssemblageInducedTreeManagerTests.GROUP_IDS))
# for induced_tree in test_target._tree_assemblage_induced_trees_map[tree]:
# self.assertIn(induced_tree, test_target._tree_cache)
def test_comparison_add_trees_directly(self):
test_target = AssemblageInducedTreeShapeKernel()
trees = [self.get_random_tree() for idx in range(10)]
for tree in trees:
test_target.update_assemblage_induced_tree_cache(
tree=tree,
assemblage_leaf_sets=tree.assemblage_leaf_sets)
## to check for inadvertent double adding
test_target.update_assemblage_induced_tree_cache(
tree=tree,
assemblage_leaf_sets=tree.assemblage_leaf_sets)
self.validate_managed_and_cached_trees(
test_target=test_target,
trees=trees)
def test_comparison_add_trees_by_comparison(self):
test_target = AssemblageInducedTreeShapeKernel()
trees = [self.get_random_tree() for idx in range(10)]
tree1 = trees[0]
for tree2 in trees[1:]:
x = test_target(
tree1=tree1,
tree2=tree2,
tree1_assemblage_leaf_sets=tree1.assemblage_leaf_sets,
tree2_assemblage_leaf_sets=tree2.assemblage_leaf_sets,
)
self.validate_managed_and_cached_trees(
test_target=test_target,
trees=trees)
if __name__ == "__main__":
unittest.main()
|
