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
|
from __future__ import absolute_import
from __future__ import print_function
# #START_LICENSE###########################################################
#
#
# This file is part of the Environment for Tree Exploration program
# (ETE). http://etetoolkit.org
#
# ETE is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# ETE is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
# License for more details.
#
# You should have received a copy of the GNU General Public License
# along with ETE. If not, see <http://www.gnu.org/licenses/>.
#
#
# ABOUT THE ETE PACKAGE
# =====================
#
# ETE is distributed under the GPL copyleft license (2008-2015).
#
# If you make use of ETE in published work, please cite:
#
# Jaime Huerta-Cepas, Joaquin Dopazo and Toni Gabaldon.
# ETE: a python Environment for Tree Exploration. Jaime BMC
# Bioinformatics 2010,:24doi:10.1186/1471-2105-11-24
#
# Note that extra references to the specific methods implemented in
# the toolkit may be available in the documentation.
#
# More info at http://etetoolkit.org. Contact: huerta@embl.de
#
#
# #END_LICENSE#############################################################
import logging
import os
import six
log = logging.getLogger("main")
from ..master_task import TreeMergeTask
from ..master_job import Job
from ..utils import (load_node_size, PhyloTree, SeqGroup, generate_id,
NPR_TREE_STYLE, NodeStyle, DEBUG,
faces, pjoin, GLOBALS)
from ..errors import ConfigError, TaskError
from .. import db
__all__ = ["TreeMerger"]
class TreeMerger(TreeMergeTask):
def __init__(self, nodeid, seqtype, task_tree, conf, confname):
# Initialize task
self.confname = confname
self.conf = conf
self.task_tree_file = task_tree
TreeMergeTask.__init__(self, nodeid, "treemerger", "TreeMerger",
None, self.conf[self.confname])
self.main_tree = None
self.task_tree = None
self.seqtype = seqtype
self.rf = None, None # Robinson foulds to orig partition
self.outgroup_match_dist = 0.0
self.outgroup_match = ""
self.pre_iter_support = None # support of the node pre-iteration
self.init()
def finish(self):
def euc_dist(x, y):
return len(x.symmetric_difference(y)) / float((len(x) + len(y)))
dataid = db.get_dataid(*self.task_tree_file.split("."))
ttree = PhyloTree(db.get_data(dataid))
mtree = self.main_tree
ttree.dist = 0
cladeid, target_seqs, out_seqs = db.get_node_info(self.threadid, self.nodeid)
self.out_seqs = out_seqs
self.target_seqs = target_seqs
ttree_content = ttree.get_cached_content()
if mtree and not out_seqs:
mtree_content = mtree.get_cached_content()
log.log(24, "Finding best scoring outgroup from previous iteration.")
for _n in mtree_content:
if _n.cladeid == cladeid:
orig_target = _n
target_left = set([_n.name for _n in mtree_content[orig_target.children[0]]])
target_right = set([_n.name for _n in mtree_content[orig_target.children[1]]])
partition_pairs = []
everything = set([_n.name for _n in ttree_content[ttree]])
for n, content in six.iteritems(ttree_content):
if n is ttree:
continue
left = set([_n.name for _n in content])
right = everything - left
d1 = euc_dist(left, target_left)
d2 = euc_dist(left, target_right)
best_match = min(d1, d2)
partition_pairs.append([best_match, left, right, n])
partition_pairs.sort()
self.outgroup_match_dist = partition_pairs[0][0]
#self.outgroup_match = '#'.join( ['|'.join(partition_pairs[0][1]),
# '|'.join(partition_pairs[0][2])] )
outgroup = partition_pairs[0][3]
ttree.set_outgroup(outgroup)
ttree.dist = orig_target.dist
ttree.support = orig_target.support
# Merge task and main trees
parent = orig_target.up
orig_target.detach()
parent.add_child(ttree)
elif mtree and out_seqs:
log.log(26, "Rooting tree using %d custom seqs" %
len(out_seqs))
self.outgroup_match = '|'.join(out_seqs)
#log.log(22, "Out seqs: %s", len(out_seqs))
#log.log(22, "Target seqs: %s", target_seqs)
if len(out_seqs) > 1:
#first root to a single seqs outside the outgroup
#(should never fail and avoids random outgroup split
#problems in unrooted trees)
ttree.set_outgroup(ttree & list(target_seqs)[0])
# Now tries to get the outgroup node as a monophyletic clade
outgroup = ttree.get_common_ancestor(out_seqs)
if set(outgroup.get_leaf_names()) ^ out_seqs:
msg = "Monophyly of the selected outgroup could not be granted! Probably constrain tree failed."
#dump_tree_debug(msg, self.taskdir, mtree, ttree, target_seqs, out_seqs)
raise TaskError(self, msg)
else:
outgroup = ttree & list(out_seqs)[0]
ttree.set_outgroup(outgroup)
orig_target = self.main_tree.get_common_ancestor(target_seqs)
found_target = outgroup.get_sisters()[0]
ttree = ttree.get_common_ancestor(target_seqs)
outgroup.detach()
self.pre_iter_support = orig_target.support
# Use previous dist and support
ttree.dist = orig_target.dist
ttree.support = orig_target.support
parent = orig_target.up
orig_target.detach()
parent.add_child(ttree)
else:
# ROOTS FIRST ITERATION
log.log(24, "Getting outgroup for first NPR split")
# if early split is provided in the command line, it
# overrides config file
mainout = GLOBALS.get("first_split_outgroup", "midpoint")
if mainout.lower() == "midpoint":
log.log(26, "Rooting to midpoint.")
best_outgroup = ttree.get_midpoint_outgroup()
if best_outgroup:
ttree.set_outgroup(best_outgroup)
else:
log.warning("Midpoint outgroup could not be set!")
ttree.set_outgroup(next(ttree.iter_leaves()))
else:
if mainout.startswith("~"):
# Lazy defined outgroup. Will trust in the common
# ancestor of two or more OTUs
strict_common_ancestor = False
outs = set(mainout[1:].split())
if len(outs) < 2:
raise TaskError(self, "First split outgroup error: common "
"ancestor calculation requires at least two OTU names")
else:
strict_common_ancestor = True
outs = set(mainout.split())
if outs - target_seqs:
raise TaskError(self, "Unknown seqs cannot be used to set first split rooting:%s" %(outs - target_seqs))
if len(outs) > 1:
anchor = list(set(target_seqs) - outs)[0]
ttree.set_outgroup(ttree & anchor)
common = ttree.get_common_ancestor(outs)
out_seqs = common.get_leaf_names()
if common is ttree:
msg = "First split outgroup could not be granted:%s" %out_seqs
#dump_tree_debug(msg, self.taskdir, mtree, ttree, target_seqs, outs)
raise TaskError(self, msg)
if strict_common_ancestor and set(out_seqs) ^ outs:
msg = "Monophyly of first split outgroup could not be granted:%s" %out_seqs
#dump_tree_debug(msg, self.taskdir, mtree, ttree, target_seqs, outs)
raise TaskError(self, msg)
log.log(26, "@@8:First split rooting to %d seqs@@1:: %s" %(len(out_seqs),out_seqs))
ttree.set_outgroup(common)
else:
single_out = outs.pop()
common = ttree.set_outgroup(single_out)
log.log(26, "@@8:First split rooting to 1 seq@@1:: %s" %(single_out))
self.main_tree = ttree
orig_target = ttree
tn = orig_target.copy()
self.pre_iter_task_tree = tn
self.rf = orig_target.robinson_foulds(ttree)
self.pre_iter_support = orig_target.support
# Reloads node2content of the rooted tree and generate cladeids
ttree_content = self.main_tree.get_cached_content()
for n, content in six.iteritems(ttree_content):
cid = generate_id([_n.name for _n in content])
n.add_feature("cladeid", cid)
#ttree.write(outfile=self.pruned_tree)
self.task_tree = ttree
def dump_tree_debug(msg, taskdir, mtree, ttree, target_seqs, out_seqs):
try:
if out_seqs is None: out_seqs = set()
if target_seqs is None: target_seqs = set()
if ttree:
for n in ttree.get_leaves():
if n.name in out_seqs:
n.name = n.name + " *__OUTGROUP__*"
if mtree:
for n in mtree.get_leaves():
if n.name in out_seqs:
n.name = n.name + " *__OUTGROUP__*"
if n.name in target_seqs:
n.name = n.name + " [ TARGET ]"
OUT = open(pjoin(taskdir, "__debug__"), "w")
print(msg, file=OUT)
print("MainTree:", mtree, file=OUT)
print("TaskTree:", ttree, file=OUT)
print("Expected outgroups:", out_seqs, file=OUT)
OUT.close()
except Exception as e:
print(e)
|
