From 76463bb3422d8fb6ae1cb9e909d4c8f745ada448 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Mon, 29 Jul 2019 13:28:46 -0500
Subject: [PATCH] Added support for circular structures in
segmentmodel_from_lists
---
mrdna/readers/segmentmodel_from_lists.py | 40 ++++++++++++++++++++++++
1 file changed, 40 insertions(+)
diff --git a/mrdna/readers/segmentmodel_from_lists.py b/mrdna/readers/segmentmodel_from_lists.py
index b13fb9a..aa12c75 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
--
GitLab