Skip to content

Commit

Permalink
version 0.1.8
Browse files Browse the repository at this point in the history 
bugfixes to interior cycle visualization
  • Loading branch information
@jluebeck
jluebeck committed Oct 27, 2023
1 parent 954a6c2 commit 313c1d7
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 16 deletions.
14 changes: 11 additions & 3 deletions CycleViz.py
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#!/usr/bin/env python

__author__ = "Jens Luebeck (jluebeck [at] ucsd.edu)"
__version__ = "0.1.7"
__version__ = "0.1.8"

import argparse
from collections import defaultdict
Expand Down Expand Up @@ -782,7 +782,10 @@ def plot_ref_genome(ref_placements, cycle, total_length, imputed_status, label_s
# if it's very small, just put one marker in the center
if skinny and not refObj.prev_is_adjacent and not refObj.next_is_adjacent:
skinny = False
tm = (str(ts[0]) + "-" + str(te[0]), np.mean((ts[1],te[1])), 0)
if cycle[ind][1] == "+":
tm = (str(int(round(ts[0] - ts[2]))) + "-" + str(int(round(te[0] - te[2]))), np.mean((ts[1],te[1])), 0)
else:
tm = (str(int(round(te[0] - te[2]))) + "-" + str(int(round(ts[0] - ts[2]))), np.mean((ts[1],te[1])), 0)
newposns.append(tm)

else:
Expand Down Expand Up @@ -1304,6 +1307,7 @@ def construct_cycle_ref_placements(cycle, segSeqD, raw_cycle_length, prev_seg_in

n_cycles = len(interior_cycles.keys())
visualized_cycles = 0
visualized_elements = 0
for IS_cnum in sorted(interior_cycles.keys()):
interior_cycle, IS_isCircular = interior_cycles[IS_cnum], IS_circular_D[IS_cnum]

Expand All @@ -1318,6 +1322,7 @@ def construct_cycle_ref_placements(cycle, segSeqD, raw_cycle_length, prev_seg_in

if IS_rObj_placements:
visualized_cycles+=1
visualized_elements+=len(IS_rObj_placements)
plot_ref_genome(IS_rObj_placements, new_interior_cycle, total_length, [False] * len(new_interior_cycle),
False, 'none',
bar_width)
Expand All @@ -1331,15 +1336,17 @@ def construct_cycle_ref_placements(cycle, segSeqD, raw_cycle_length, prev_seg_in

else: # strict mode
print("\nTrying to match", curr_interior_cycle)
ic_pad = total_length*0.01
for IS_rObj_placements, new_interior_cycle, new_IS_links in vu.handle_IS_data(ref_placements, cycle, isCycle,
curr_interior_cycle, interior_segSeqD, IS_isCircular, IS_bh):
# print(new_interior_cycle)
# print(len(IS_rObj_placements), len(new_interior_cycle))
visualized_elements+=1
hit = True
max_olaps = 0
used_rows = set()
for currObj in IS_rObj_placements.values():
abs_start, abs_end = int(currObj.abs_start_pos)-1, int(currObj.abs_end_pos) + 1
abs_start, abs_end = int(currObj.abs_start_pos)-ic_pad, int(currObj.abs_end_pos) + ic_pad
olaps = covered_posns[abs_start:abs_end]
for o in olaps:
used_rows.add(o.data)
Expand Down Expand Up @@ -1379,6 +1386,7 @@ def construct_cycle_ref_placements(cycle, segSeqD, raw_cycle_length, prev_seg_in
print(IS_cnum, "no hit")

print(str(visualized_cycles) + " paths of the " + str(n_cycles) + " given interior paths were visualized")
print(str(visualized_elements) + " elements were placed counting all occurrences of the interior cycles")

# cytobanding
if args.annotate_structure != "genes":
Expand Down
49 changes: 36 additions & 13 deletions VizUtil.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -822,11 +822,12 @@ def cycles_match(ref_placements, cycle, isCycle, start_ref_ind, interior_cycle,

print("found initial matching seg: ", rObj.to_string())
if intdir == "+":
curr_interior_offset = rObj.ref_end - s
curr_interior_offset = s - rObj.ref_start
else:
curr_interior_offset = e - rObj.ref_start
curr_interior_offset = rObj.ref_end - e

start_offset = curr_interior_offset
print(start_offset, "start offset")

ref_ind = start_ref_ind
while left_interior_cycle_index < len(interior_cumulative_bounds) - 1:
Expand All @@ -836,49 +837,58 @@ def cycles_match(ref_placements, cycle, isCycle, start_ref_ind, interior_cycle,
ref_ind = -1

right_interior_cycle_index_hit = bisect.bisect(interior_cumulative_bounds, curr_interior_offset) - 1
print("Right index", right_interior_cycle_index_hit)
# print("Right index", right_interior_cycle_index_hit, curr_interior_offset, interior_cumulative_bounds)

ref_end_seg = ref_ind
segID, intdir = interior_cycle[left_interior_cycle_index]
c, s, e = interior_segSeqD[segID]
if rObj.direction == "+":
end_offset = e - rObj.ref_start
else:
end_offset = s - rObj.ref_start
end_offset = e - rObj.ref_start

interior_contains_ref = s <= rObj.ref_start <= rObj.ref_end <= e
overlaps = (rObj.ref_start <= s <= rObj.ref_end or rObj.ref_start <= e <= rObj.ref_end)
# print(c, rObj.chrom, intdir, rObj.direction, overlaps)
if c != rObj.chrom or intdir != rObj.direction or not any([interior_contains_ref, overlaps]):
print("different chrom, dir or bounds")
return False, ref_start_seg, ref_end_seg, start_offset, end_offset

# print(right_interior_cycle_index_hit, left_interior_cycle_index)
if right_interior_cycle_index_hit != left_interior_cycle_index:
# print(right_interior_cycle_index_hit, len(interior_cumulative_bounds) - 1)
if right_interior_cycle_index_hit == len(interior_cumulative_bounds) - 1:
for x in range(left_interior_cycle_index, right_interior_cycle_index_hit-1):
if not next_seg_index_is_adj[x]:
print(x, "was not adjacent")
return False, ref_start_seg, ref_end_seg, start_offset, end_offset

if rObj.direction == "+":
# print(s,e, rObj.ref_start, rObj.ref_end)
end_offset = e - rObj.ref_start
print(end_offset, "+ end offset")
else:
end_offset = s - rObj.ref_start
# print(s,e, rObj.ref_start, rObj.ref_end)
end_offset = e - rObj.ref_start
print(end_offset, "- end offset")

print("Reached end of interior cycle")
return True, ref_start_seg, ref_end_seg, start_offset, end_offset

else:
# print("not hitting end yet")
if rObj.direction == "+" and abs(rObj.ref_end - e) < 2:
print("both uptick positive")
left_interior_cycle_index+=1
curr_interior_offset += (rObj.ref_end - rObj.ref_start)
ref_ind += 1
rObj = ref_placements[ref_ind]
curr_interior_offset += (rObj.ref_end - rObj.ref_start)
elif rObj.direction == "-" and abs(rObj.ref_start - s) < 2:
print("both uptick negative")
left_interior_cycle_index+=1
curr_interior_offset += (rObj.ref_end - rObj.ref_start)
ref_ind += 1
rObj = ref_placements[ref_ind]
curr_interior_offset += (rObj.ref_end - rObj.ref_start)
else:
for x in range(left_interior_cycle_index, right_interior_cycle_index_hit):
if not next_seg_index_is_adj[x]:
Expand All @@ -891,7 +901,10 @@ def cycles_match(ref_placements, cycle, isCycle, start_ref_ind, interior_cycle,
else:
# it was the same interior seg
print("same interior segment")
print(c,s,e)
# print(c,s,e)
#
# print(rObj.direction, s, e, rObj.ref_start, rObj.ref_end)

if rObj.direction == "+" and abs(rObj.ref_end - e) < 2:
print("uptick both positive")
left_interior_cycle_index += 1
Expand All @@ -902,15 +915,15 @@ def cycles_match(ref_placements, cycle, isCycle, start_ref_ind, interior_cycle,
left_interior_cycle_index += 1
tb = True

curr_interior_offset += (rObj.ref_end - rObj.ref_start)
ref_ind += 1
rObj = ref_placements[ref_ind]
curr_interior_offset += (rObj.ref_end - rObj.ref_start)

# if ref_ind == start_ref_ind:
# print("went full circle")
# return False, ref_start_seg, ref_end_seg, start_offset, end_offset
# print(left_interior_cycle_index, right_interior_cycle_index_hit, len(interior_cumulative_bounds))
if right_interior_cycle_index_hit == len(interior_cumulative_bounds) - 1 and tb:
# print(right_interior_cycle_index_hit, len(interior_cumulative_bounds)-1, tb)
if right_interior_cycle_index_hit == len(interior_cumulative_bounds) - 2 and tb:
print("Finished matching")
return True, ref_start_seg, ref_end_seg, start_offset, end_offset

Expand Down Expand Up @@ -944,6 +957,7 @@ def handle_IS_data(ref_placements, cycle, isCycle, interior_cycle, interior_segS
isCycle, ind, interior_cycle, interior_segSeqD, IS_isCircular)

print(matched, ref_start_ind, ref_end_ind, ref_start_offset, ref_end_offset)
print("")
if not matched:
continue

Expand All @@ -957,7 +971,11 @@ def handle_IS_data(ref_placements, cycle, isCycle, interior_cycle, interior_segS
for matched_ref_ind in ref_inds_hit:
curr_rObj = ref_placements[matched_ref_ind]
if matched_ref_ind == ref_start_ind:
normStart = curr_rObj.abs_start_pos + ref_start_offset
if curr_rObj.direction == "+":
normStart = curr_rObj.abs_start_pos + ref_start_offset
else:
normStart = curr_rObj.abs_start_pos + ref_start_offset

if ref_start_ind == ref_end_ind:
normEnd = curr_rObj.abs_start_pos + ref_end_offset
else:
Expand All @@ -970,6 +988,11 @@ def handle_IS_data(ref_placements, cycle, isCycle, interior_cycle, interior_segS
elif matched_ref_ind == ref_end_ind:
normStart = curr_rObj.abs_start_pos
normEnd = curr_rObj.abs_start_pos + ref_end_offset
if curr_rObj.direction == "+":
normEnd = curr_rObj.abs_start_pos + ref_end_offset
else:
normStart = curr_rObj.abs_start_pos + ref_end_offset

seg_chrom, seg_s, seg_e = interior_segSeqD[interior_cycle[-1][0]]

else:
Expand Down Expand Up @@ -1013,7 +1036,7 @@ def handle_IS_data(ref_placements, cycle, isCycle, interior_cycle, interior_segS
else:
new_IS_links.append(False)

if (new_interior_cycle[-1], new_interior_cycle[0]) in valid_link_pairs or (new_interior_cycle[0], new_interior_cycle[-1]) in valid_link_pairs:
if ((new_interior_cycle[-1], new_interior_cycle[0]) in valid_link_pairs or (new_interior_cycle[0], new_interior_cycle[-1]) in valid_link_pairs) and IS_isCircular:
new_IS_links.append(True)
else:
new_IS_links.append(False)
Expand Down

0 comments on commit 313c1d7

Please sign in to comment.