From 7ffdf594e9eebe82ca44cbb3c1ed412571fe20fe Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Tue, 28 Jul 2020 18:41:35 -0500
Subject: [PATCH] Fixed model_from_lists for 3'-to-5' structures; fixed issue
with single-nucleotide 'crossovers'; added debug flag for tracing strands
---
mrdna/readers/segmentmodel_from_lists.py | 30 ++++++---------
mrdna/segmentmodel.py | 48 +++++++++++++++++-------
2 files changed, 46 insertions(+), 32 deletions(-)
diff --git a/mrdna/readers/segmentmodel_from_lists.py b/mrdna/readers/segmentmodel_from_lists.py
index 96f610a..c1ce49b 100644
--- a/mrdna/readers/segmentmodel_from_lists.py
+++ b/mrdna/readers/segmentmodel_from_lists.py
@@ -307,24 +307,17 @@ def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
assert( np.all(bps[ss_residues] == -1) )
for s,c in zip(strands, strand_is_circular):
- ss_residues = s[np.where(hmap[s] < 0)[0]]
- if len(ss_residues) == 0: continue
- resid_diff = np.diff(ss_residues)
-
- tmp = np.where( resid_diff != 1 )[0]
- first_res = ss_residues[0]
- for i in tmp:
- ids = np.arange(first_res, ss_residues[i]+1)
- assert( np.all(hmap[ids] == -1) )
- hmap[ids] = hid
- hrank[ids] = ids-first_res
- first_res = ss_residues[i+1]
- hid += 1
- ids = np.arange(first_res, ss_residues[-1]+1)
- assert( np.all(hmap[ids] == -1) )
- hmap[ids] = hid
- hrank[ids] = ids-first_res
- hid+=1
+ strand_segment_ends = [i for i in np.where( np.diff(hmap[s]) != 0 )[0]] + [len(s)-1]
+ seg_start = 0
+ for i in strand_segment_ends:
+ if hmap[s[i]] < 0:
+ ## Found single-stranded segment
+ ids = s[seg_start:i+1]
+ assert( np.all(hmap[ids] == -1) )
+ hmap[ids] = hid
+ hrank[ids] = np.arange(i+1-seg_start)
+ hid+=1
+ seg_start = i+1
single_stranded_helices = np.arange(double_stranded_helices[-1]+1,hid)
@@ -419,6 +412,7 @@ def model_from_basepair_stack_3prime(coordinate, basepair, stack, three_prime,
s = allSegments[hid]
s.sequence = [sequence[r] for r in sorted(resids,key=lambda x: hrank[x])]
+
## Build model
model = SegmentModel( allSegments,
max_basepairs_per_bead = max_basepairs_per_bead,
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index b7603b7..3e1c069 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -23,7 +23,6 @@ from .model.spring_from_lp import k_twist as twist_spring_from_lp
import csv
-
# import pdb
"""
TODO:
@@ -48,6 +47,9 @@ TODO:
- refactor parts of Segment into an abstract_polymer class
- make each call generate_bead_model, generate_atomic_model, generate_oxdna_model return an object with only have a reference to original object
"""
+
+_DEBUG_TRACE = False
+
class CircularDnaError(Exception):
pass
@@ -947,7 +949,8 @@ class Segment(ConnectableElement, Group):
## Stop if we found the 3prime end
if l.on_fwd_strand == is_fwd and l.type_ == "3prime" and l.connection is None:
- # print(" found end at",l)
+ if _DEBUG_TRACE:
+ print(" found end at",l)
return pos, None, None, None, None
## Check location connections
@@ -957,15 +960,17 @@ class Segment(ConnectableElement, Group):
## Found a location on the same strand?
if l.on_fwd_strand == is_fwd:
- # print(" passing through",l)
- # print("from {}, connection {} to {}".format(nt_pos,l,B))
+ if _DEBUG_TRACE:
+ print(" passing through",l)
+ print("from {}, connection {} to {}".format(nt_pos,l,B))
Bpos = B.get_nt_pos()
return pos, B.container, Bpos, B.on_fwd_strand, 0.5
## Stop at other strand crossovers so basepairs line up
elif c.type_ == "crossover":
if nt_pos == pos: continue
- # print(" pausing at",l)
+ if _DEBUG_TRACE:
+ print(" pausing at",l)
return pos, l.container, pos+(2*is_fwd-1), is_fwd, 0
raise Exception("Shouldn't be here")
@@ -1350,14 +1355,26 @@ class SingleStrandedSegment(Segment):
assert(nt < self.num_nt)
assert(other_nt < other.num_nt)
- if nt in (0,1,self.num_nt-1) and other_nt in (0,1,other.num_nt-1):
- if nt_on_5prime == True:
- other_end = other.start5 if strands_fwd[1] else other.end5
- self.connect_end3( other_end )
- else:
- other_end = other.end3 if strands_fwd[1] else other.start3
- self.connect_start5( other_end )
- return
+ if nt_on_5prime:
+ if nt == self.num_nt-1 and other_nt in (0,other.num_nt-1):
+ other_end = None
+ if strands_fwd[1] and other_nt == 0:
+ other_end = other.start5
+ elif (not strands_fwd[1]) and other_nt == other.num_nt-1:
+ other_end = other.end5
+ if other_end is not None:
+ self.connect_end3( other_end )
+ return
+ else:
+ if nt == 0 and other_nt in (0,other.num_nt-1):
+ other_end = None
+ if strands_fwd[1] and other_nt == other.num_nt-1:
+ other_end = other.end3
+ elif (not strands_fwd[1]) and other_nt == 0:
+ other_end = other.start3
+ if other_end is not None:
+ self.connect_start5( other_end )
+ return
# c1 = self.nt_pos_to_contour(nt)
# # TODOTODO
@@ -3194,7 +3211,8 @@ class SegmentModel(ArbdModel):
if locs is None: continue
# for pos, is_fwd in locs:
for l in locs:
- # print("Tracing",l)
+ if _DEBUG_TRACE:
+ print("Tracing 5prime",l)
# TODOTODO
pos = seg.contour_to_nt_pos(l.address, round_nt=True)
is_fwd = l.on_fwd_strand
@@ -3240,6 +3258,8 @@ class SegmentModel(ArbdModel):
history = []
if not strands_cover_segment(seg, is_fwd):
pos = nt = find_nt_not_in_strand(seg, is_fwd)
+ if _DEBUG_TRACE:
+ print("Tracing circular strand", pos, is_fwd)
s = Strand(segname="S{:03d}".format(len(strands)),
is_circular = True)
history = _recursively_build_strand(s, seg, pos, is_fwd)
--
GitLab