From b631e30fcf56229b1543a2e4db0664c5115721ed Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Tue, 4 Sep 2018 14:14:01 -0500
Subject: [PATCH] Added support for circular strands
---
dnarbd/segmentmodel.py | 116 +++++++++++++++++++++++++++++++++--------
1 file changed, 95 insertions(+), 21 deletions(-)
diff --git a/dnarbd/segmentmodel.py b/dnarbd/segmentmodel.py
index 30d01bc..6dbe417 100644
--- a/dnarbd/segmentmodel.py
+++ b/dnarbd/segmentmodel.py
@@ -36,6 +36,8 @@ TODO:
- document
- add unit test of helices connected to themselves
"""
+class CircularDnaError(Exception):
+ pass
class ParticleNotConnectedError(Exception):
pass
@@ -1117,11 +1119,12 @@ class StrandInSegment(Group):
class Strand(Group):
""" Represents an entire ssDNA strand from 5' to 3' as it routes through segments """
- def __init__(self, segname = None):
+ def __init__(self, segname = None, is_circular = False):
Group.__init__(self)
self.num_nt = 0
self.children = self.strand_segments = []
self.segname = segname
+ self.is_circular = is_circular
## TODO disambiguate names of functions
def add_dna(self, segment, start, end, is_fwd):
@@ -1135,9 +1138,7 @@ class Strand(Group):
# assert( s.start not in (start,end) )
# assert( s.end not in (start,end) )
if s.start in (start,end) or s.end in (start,end):
- print(" CIRCULAR DNA")
- import pdb
- pdb.set_trace()
+ raise CircularDnaError("Found circular DNA")
s = StrandInSegment( segment, start, end, is_fwd )
self.add( s )
@@ -1176,6 +1177,21 @@ class Strand(Group):
# assert( len(ret) == self.num_nt )
# return ret
+ def link_nucleotides(self, nt5, nt3):
+ parent = nt5.parent if nt5.parent is nt3.parent else self
+ o3,c3,c4,c2,h3 = [nt5.atoms_by_name[n]
+ for n in ("O3'","C3'","C4'","C2'","H3'")]
+ p,o5,o1,o2,c5 = [nt3.atoms_by_name[n]
+ for n in ("P","O5'","O1P","O2P","C5'")]
+ parent.add_bond( o3, p, None )
+ parent.add_angle( c3, o3, p, None )
+ for x in (o5,o1,o2):
+ parent.add_angle( o3, p, x, None )
+ parent.add_dihedral(c3, o3, p, x, None )
+ for x in (c4,c2,h3):
+ parent.add_dihedral(x, c3, o3, p, None )
+ parent.add_dihedral(o3, p, o5, c5, None)
+
def generate_atomic_model(self, scale, first_atomic_index):
last = None
resid = 1
@@ -1190,9 +1206,9 @@ class Strand(Group):
# assert(s.end != s.start)
assert( s.num_nt == 1 or (np.linalg.norm( seg.contour_to_position(contour[-1]) - seg.contour_to_position(contour[0]) ) > 0.1) )
for c,seq in zip(contour,s.get_sequence()):
- if last is None:
+ if last is None and not self.is_circular:
seq = "5"+seq
- if strand_segment_count == len(s.strand_segments) and c == 1:
+ if strand_segment_count == len(s.strand_segments) and c == 1 and not self.is_circular:
seq = seq+"3"
nt = seg._generate_atomic_nucleotide( c, s.is_fwd, seq, scale )
@@ -1203,23 +1219,17 @@ class Strand(Group):
s.add(nt)
## Join last basepairs
if last is not None:
- o3,c3,c4,c2,h3 = [last.atoms_by_name[n]
- for n in ("O3'","C3'","C4'","C2'","H3'")]
- p,o5,o1,o2,c5 = [nt.atoms_by_name[n]
- for n in ("P","O5'","O1P","O2P","C5'")]
- self.add_bond( o3, p, None )
- self.add_angle( c3, o3, p, None )
- for x in (o5,o1,o2):
- self.add_angle( o3, p, x, None )
- self.add_dihedral(c3, o3, p, x, None )
- for x in (c4,c2,h3):
- self.add_dihedral(x, c3, o3, p, None )
- self.add_dihedral(o3, p, o5, c5, None)
+ self.link_nucleotides(last,nt)
+
nt.__dict__['resid'] = resid
resid += 1
last = nt
nt._first_atomic_index = first_atomic_index
first_atomic_index += len(nt.children)
+
+ if self.is_circular:
+ self.link_nucleotides(last,self.strand_segments[0].children[0])
+
return first_atomic_index
def update_atomic_orientations(self,default_orientation):
@@ -2230,7 +2240,15 @@ class SegmentModel(ArbdModel):
# pdb.set_trace()
end_pos, next_seg, next_pos, next_dir, move_at_least = seg.get_strand_segment(pos, is_fwd, move_at_least)
- strand.add_dna(seg, pos, end_pos, is_fwd)
+
+ try:
+ strand.add_dna(seg, pos, end_pos, is_fwd)
+ except CircularDnaError:
+ ## Circular DNA was found
+ break
+ except:
+ print("Unexpected error:", sys.exc_info()[0])
+ raise
if next_seg is None:
return
@@ -2252,8 +2270,64 @@ class SegmentModel(ArbdModel):
_recursively_build_strand(s, seg, pos, is_fwd)
# print("{} {}".format(seg.name,s.num_nt))
strands.append(s)
-
- self.strands = sorted(strands, key=lambda s:s.num_nt)[::-1] # or something
+
+ ## Trace circular DNA
+ def strands_cover_segment(segment, is_fwd=True):
+ direction = 'fwd' if is_fwd else 'rev'
+ nt = 0
+ for sp in segment.strand_pieces[direction]:
+ nt += sp.num_nt
+ return nt == segment.num_nt
+
+ def find_nt_not_in_strand(segment, is_fwd=True):
+ fwd_str = 'fwd' if is_fwd else 'rev'
+
+ def check(val):
+ assert(val >= 0 and val < segment.num_nt)
+ print("find_nt_not_in_strand({},{}) returning {}".format(
+ segment, is_fwd, val))
+ return val
+
+ if is_fwd:
+ last = -1
+ for sp in segment.strand_pieces[fwd_str]:
+ if sp.start-last > 1:
+ return check(last+1)
+ last = sp.end
+ return check(last+1)
+ else:
+ last = segment.num_nt
+ for sp in segment.strand_pieces[fwd_str]:
+ if last-sp.end > 1:
+ return check(last-1)
+ last = sp.start
+ return check(last-1)
+
+ def add_strand_if_needed(seg,is_fwd):
+ if not strands_cover_segment(seg, is_fwd):
+ pos = nt = find_nt_not_in_strand(seg, is_fwd)
+ s = Strand(is_circular = True)
+ _recursively_build_strand(s, seg, pos, is_fwd)
+ strands.append(s)
+
+ for seg in self.segments:
+ add_strand_if_needed(seg,True)
+ if isinstance(seg, DoubleStrandedSegment):
+ add_strand_if_needed(seg,False)
+
+ self.strands = sorted(strands, key=lambda s:s.num_nt)[::-1]
+ def check_strands():
+ dsdna = filter(lambda s: isinstance(s,DoubleStrandedSegment), self.segments)
+ for s in dsdna:
+ nt_fwd = nt_rev = 0
+ for sp in s.strand_pieces['fwd']:
+ nt_fwd += sp.num_nt
+ for sp in s.strand_pieces['rev']:
+ nt_rev += sp.num_nt
+ assert( nt_fwd == s.num_nt and nt_rev == s.num_nt )
+ # print("{}: {},{} (fwd,rev)".format(s.name, nt_fwd/s.num_nt,nt_rev/s.num_nt))
+ check_strands()
+
## relabel segname
counter = 0
for s in self.strands:
--
GitLab