Skip to content

Commit 7248729

Browse files
committed
hairpin resolution bugfix
1 parent a7396a5 commit 7248729

File tree

1 file changed

+7
-1
lines changed

1 file changed

+7
-1
lines changed

src/scripts/resolve_triplets_kmerify.py

+7-1
Original file line numberDiff line numberDiff line change
@@ -243,9 +243,14 @@ def get_valid_triplets(node, edges, paths_crossing, min_edge_support, min_covera
243243

244244
def resolve_hairpins(nodelength, nodes, paths_crossing, node_seqs, node_lens, edges, maybe_resolvable, min_edge_support, min_coverage, removable_nodes):
245245
hairpins = set()
246+
unresolvable_hairpins = set()
246247
for node in iterate_deterministic(nodes):
247248
if ">" + node not in edges: continue
248249
if "<" + node not in edges: continue
250+
if len(edges[">" + node]) >= 2 and "<" + node in edges[">" + node]:
251+
unresolvable_hairpins.add(node)
252+
if len(edges["<" + node]) >= 2 and ">" + node in edges["<" + node]:
253+
unresolvable_hairpins.add(node)
249254
if len(edges[">" + node]) == 1 and getone(edges[">" + node]) == "<" + node:
250255
hairpins.add(">" + node)
251256
if len(edges["<" + node]) == 1 and getone(edges["<" + node]) == ">" + node:
@@ -257,6 +262,7 @@ def resolve_hairpins(nodelength, nodes, paths_crossing, node_seqs, node_lens, ed
257262
for node in hairpins:
258263
if revnode(node) in hairpins: continue # double hairpin resolution hard to implement, so just ignore it
259264
if len(edges[revnode(node)]) == 0: continue # don't resolve disconnected hairpins, they're probably spurious anyway
265+
if node[1:] in unresolvable_hairpins: continue # weird unresolvable structures
260266
resolutions = {}
261267
for path in iterate_paths(paths_crossing, node[1:]):
262268
if len(path) < 4: continue
@@ -310,7 +316,7 @@ def resolve_hairpins(nodelength, nodes, paths_crossing, node_seqs, node_lens, ed
310316
remove_paths.append(path)
311317
if len(path) < 4: continue
312318
add_this = []
313-
add_this.append(path[0])
319+
if path[0][1:] != node[1:]: add_this.append(path[0])
314320
for i in range(1, len(path)-2):
315321
if path[i] == node and path[i+1] == revnode(node):
316322
key = (path[i-1], path[i+2])

0 commit comments

Comments
 (0)