Commit 76463bb3 authored by cmaffeo2's avatar cmaffeo2
Browse files

Added support for circular structures in segmentmodel_from_lists

parent 93b5c9da
...@@ -67,6 +67,46 @@ def basepairs_and_stacks_to_helixmap(basepairs,stacks_above): ...@@ -67,6 +67,46 @@ def basepairs_and_stacks_to_helixmap(basepairs,stacks_above):
is_fwd[bp] = 0 is_fwd[bp] = 0
rank +=1 rank +=1
hid += 1 hid += 1
## Create "helix" for each circular segment
processed = set()
unclaimed_bases = np.where( (basepairs >= 0)*(helixmap == -1) )[0]
for nt0 in unclaimed_bases:
if nt0 in processed: continue
nt = nt0
all_nts = [nt]
counter = 0
while stacks_above[nt] >= 0:
nt = stacks_above[nt]
bp = basepairs[nt]
assert( bp >= 0 )
assert(helixmap[nt] == -1)
assert(helixmap[bp] == -1)
all_nts.append(nt)
if nt == nt0:
break
if counter > 1e6:
raise Exception("DNA is apparently too long; probably there is something wrong with the structure")
counter += 1
## Split circular helix into two groups (TODO: fix segmentmodel so that circular helices work and this is not needed)
for group in (all_nts[:len(all_nts)//2], all_nts[len(all_nts)//2:]):
rank = 0
for nt in group:
bp = basepairs[nt]
is_fwd[bp] = 0
rank +=1
helixmap[nt] = helixmap[bp] = hid
helixrank[nt] = helixrank[bp] = rank
processed.add(nt)
processed.add(bp)
rank +=1
hid += 1
return helixmap, helixrank, is_fwd return helixmap, helixrank, is_fwd
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment