Commit b631e30f authored by cmaffeo2's avatar cmaffeo2
Browse files

Added support for circular strands

parent 7bea2ccd
......@@ -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:
......
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