diff --git a/mrdna/readers/segmentmodel_from_lists.py b/mrdna/readers/segmentmodel_from_lists.py index b13fb9a24e6231d29334f611310dd387fcd8f405..aa12c75f3315aeaf5e29b4970667457ef5b5364b 100644 --- a/mrdna/readers/segmentmodel_from_lists.py +++ b/mrdna/readers/segmentmodel_from_lists.py @@ -67,6 +67,46 @@ def basepairs_and_stacks_to_helixmap(basepairs,stacks_above): is_fwd[bp] = 0 rank +=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