Commit 7ffdf594 authored by cmaffeo2's avatar cmaffeo2
Browse files

Fixed model_from_lists for 3'-to-5' structures; fixed issue with...

Fixed model_from_lists for 3'-to-5' structures; fixed issue with single-nucleotide 'crossovers'; added debug flag for tracing strands
parent 0befffbe
......@@ -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,
......
......@@ -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)
......
Markdown is supported
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