Skip to content

GitLab

Commit 2c07f592 authored by cmaffeo2's avatar cmaffeo2
Browse files

Made it so a bead can be located at multiple distant locations within a... 

Made it so a bead can be located at multiple distant locations within a segment (e.g. for looped constructs)
parent 15b4880d
Hide whitespace changes
Inline Side-by-side
mrdna/segmentmodel.py
View file @ 2c07f592
......@@ -215,14 +215,14 @@ class SegmentParticle(PointParticle):
""" Returns bead directly above self """
# assert( len(self.intrahelical_neighbors) <= 2 )
for b in self.intrahelical_neighbors:
if b.get_contour_position(self.parent) > self.contour_position:
if b.get_contour_position(self.parent, self.contour_position) > self.contour_position:
return b
def get_intrahelical_below(self):
""" Returns bead directly below self """
# assert( len(self.intrahelical_neighbors) <= 2 )
for b in self.intrahelical_neighbors:
if b.get_contour_position(self.parent) < self.contour_position:
if b.get_contour_position(self.parent, self.contour_position) < self.contour_position:
return b
def _neighbor_should_be_added(self,b):
......@@ -230,14 +230,14 @@ class SegmentParticle(PointParticle):
return True
c1 = self.contour_position
c2 = b.get_contour_position(self.parent)
c2 = b.get_contour_position(self.parent,c1)
if c2 < c1:
b0 = self.get_intrahelical_below()
else:
b0 = self.get_intrahelical_above()
if b0 is not None:
c0 = b0.get_contour_position(self.parent)
c0 = b0.get_contour_position(self.parent,c1)
if np.abs(c2-c1) < np.abs(c0-c1):
## remove b0
self.intrahelical_neighbors.remove(b0)
......@@ -256,17 +256,42 @@ class SegmentParticle(PointParticle):
self.intrahelical_neighbors.append(b)
b.intrahelical_neighbors.append(self)
def get_nt_position(self, seg):
def conceptual_get_position(self, context):
"""
context: object
Q: does this function do too much?
Danger of mixing return values
Q: does context describe the system or present an argument?
"""
## Validate Inputs
...
## Algorithm
"""
context specifies:
- kind of output: real space, nt within segment, fraction of segment
- absolute or relative
- constraints: e.g. if passing through
"""
"""
given context, provide the position
input
"""
def get_nt_position(self, seg, near_address=None):
""" Returns the "address" of the nucleotide relative to seg in
nucleotides, taking the shortest (intrahelical) contour length route to seg
"""
if seg == self.parent:
return seg.contour_to_nt_pos(self.contour_position)
pos = self.contour_position
else:
pos = self.get_contour_position(seg)
return seg.contour_to_nt_pos(pos)
pos = self.get_contour_position(seg,near_address)
return seg.contour_to_nt_pos(pos)
def get_contour_position(self,seg):
def get_contour_position(self,seg, address = None):
""" TODO: fix paradigm where a bead maps to exactly one location in a polymer
- One way: modify get_contour_position to take an optional argument that indicates where in the polymer you are looking from
"""
......@@ -315,8 +340,10 @@ class SegmentParticle(PointParticle):
results,history = descend_search_tree(self.parent, self.contour_position)
if results is None or len(results) == 0:
raise Exception("Could not find location in segment") # TODO better error
return sorted(results,key=lambda x:x[0])[0][1]
if address is not None:
return sorted(results,key=lambda x:(x[0],(x[1]-address)**2))[0][1]
else:
return sorted(results,key=lambda x:x[0])[0][1]
# nt_pos = self.get_nt_position(seg)
# return seg.nt_pos_to_contour(nt_pos)
......@@ -829,9 +856,7 @@ class Segment(ConnectableElement, Group):
ret.append( tmp )
return ret
def _add_bead(self,b,set_contour=False):
if set_contour:
b.contour_position = b.get_contour_position(self)
def _add_bead(self,b):
# assert(b.parent is None)
if b.parent is not None:
......@@ -892,27 +917,29 @@ class Segment(ConnectableElement, Group):
# if self.name == "S001":
# pdb.set_trace()
existing_beads0 = {l.particle for l in self.locations if l.particle is not None}
existing_beads = sorted( list(existing_beads0), key=lambda b: b.get_contour_position(self) )
# pdb.set_trace()
existing_beads0 = { (l.particle, l.particle.get_contour_position(self,l.address))
for l in self.locations if l.particle is not None }
existing_beads = sorted( list(existing_beads0), key=lambda bc: bc[1] )
# if self.num_nt == 1 and all([l.particle is not None for l in self.locations]):
# pdb.set_trace()
# return
for b in existing_beads:
for b,c in existing_beads:
assert(b.parent is not None)
## Add ends if they don't exist yet
## TODOTODO: test 1 nt segments?
if len(existing_beads) == 0 or existing_beads[0].get_nt_position(self) >= 0.5:
if len(existing_beads) == 0 or existing_beads[0][0].get_nt_position(self,0) >= 0.5:
# if len(existing_beads) > 0:
# assert(existing_beads[0].get_nt_position(self) >= 0.5)
b = self._generate_one_bead( self.nt_pos_to_contour(0), 0)
existing_beads = [b] + existing_beads
existing_beads = [(b,0)] + existing_beads
if existing_beads[-1].get_nt_position(self)-(self.num_nt-1) < -0.5 or len(existing_beads)==1:
if existing_beads[-1][0].get_nt_position(self,1)-(self.num_nt-1) < -0.5 or len(existing_beads)==1:
b = self._generate_one_bead( self.nt_pos_to_contour(self.num_nt-1), 0)
existing_beads.append(b)
existing_beads.append( (b,1) )
assert(len(existing_beads) > 1)
## Walk through existing_beads, add beads between
......@@ -920,15 +947,16 @@ class Segment(ConnectableElement, Group):
last = None
for I in range(len(existing_beads)-1):
eb1,eb2 = [existing_beads[i] for i in (I,I+1)]
assert( eb1 is not eb2 )
eb1,eb2 = [existing_beads[i][0] for i in (I,I+1)]
ec1,ec2 = [existing_beads[i][1] for i in (I,I+1)]
assert( (eb1,ec1) is not (eb2,ec2) )
# if np.isclose(eb1.position[2], eb2.position[2]):
# import pdb
# pdb.set_trace()
# print(" %s working on %d to %d" % (self.name, eb1.position[2], eb2.position[2]))
e_ds = eb2.get_contour_position(self) - eb1.get_contour_position(self)
e_ds = ec2-ec1
num_beads = self._get_num_beads( e_ds, max_basepairs_per_bead, max_nucleotides_per_bead )
## Ensure there is a ssDNA bead between dsDNA beads
......@@ -945,7 +973,7 @@ class Segment(ConnectableElement, Group):
if eb1.parent == self:
tmp_children.append(eb1)
s0 = eb1.get_contour_position(self)
s0 = ec1
if last is not None:
last.make_intrahelical_neighbor(eb1)
last = eb1
......@@ -1668,63 +1696,75 @@ class SegmentModel(ArbdModel):
print( "WARNING: none type found in connection, skipping" )
cabs = [e for e in cabs if e[2].particle is not None]
beads = set(s.beads + [A.particle for c,A,B in cabs])
## Add nearby beads
for c,A,B in cabs:
## TODOTODO test?
filter_fn = lambda x: x is not None and x not in beads
bs = list( filter( filter_fn, B.particle.intrahelical_neighbors ) )
beads.update(bs)
for i in range(3):
bs = list( filter( filter_fn, [n for b in bs for n in b.intrahelical_neighbors] ) )
def get_beads_and_contour_positions(s):
ret_list = []
def zip_bead_contour(beads,address=None):
if isinstance(address,list):
assert(False)
for b,a in zip(beads,address):
if b is None: continue
try:
ret_list.append((b, b.get_contour_position(s,a)))
except:
...
else:
for b in beads:
if b is None: continue
try:
ret_list.append((b, b.get_contour_position(s,address)))
except:
...
return ret_list
## Add beads from segment s
beads_contours = zip_bead_contour(s.beads)
beads_contours.extend( zip_bead_contour([A.particle for c,A,B in cabs]) )
beads = set([b for b,c in beads_contours])
## Add nearby beads
for c,A,B in cabs:
## TODOTODO test?
filter_fn = lambda x: x is not None and x not in beads
bs = list( filter( filter_fn, B.particle.intrahelical_neighbors ) )
beads_contours.extend( zip_bead_contour( bs, A.address ) )
beads.update(bs)
for i in range(3):
bs = list( filter( filter_fn, [n for b in bs for n in b.intrahelical_neighbors] ) )
beads_contours.extend( zip_bead_contour( bs, A.address ) )
beads.update(bs)
beads = list(beads)
## Skip beads that are None (some locations were not assigned a particle to avoid double-counting)
# beads = [b for b in beads if b is not None]
beads = list(filter(lambda x: x is not None, beads))
if isinstance(s, DoubleStrandedSegment):
beads = list(filter(lambda x: x.type_.name[0] == "D", beads))
contours = []
remove = []
for b in beads:
try:
# if s.name == "61-1" and b.parent.name == "60-1":
# pdb.set_trace()
## might fail if s is an entire segment away
contours.append(b.get_contour_position(s))
except:
## ignore failed attempts
remove.append(b)
beads = list(filter(lambda x: x not in remove, beads))
beads_contours = list(set(beads_contours))
beads = list(beads)
## Skip beads that are None (some locations were not assigned a particle to avoid double-counting)
# beads = [b for b in beads if b is not None]
assert( np.any([b is None for b,c in beads_contours]) == False )
# beads = list(filter(lambda x: x[0] is not None, beads))
if isinstance(s, DoubleStrandedSegment):
beads_contours = list(filter(lambda x: x[0].type_.name[0] == "D", beads_contours))
return beads_contours
beads_contours = get_beads_and_contour_positions(s)
contours = [c for b,c in beads_contours]
contours = np.array(contours, dtype=np.float16) # deliberately use low precision
contours,ids = np.unique(contours, return_index=True)
if np.any( (contours[:-1] - contours[1:])**2 < 1e-8 ):
pdb.set_trace()
contours,ids1 = np.unique(contours, return_index=True)
beads_contours = [beads_contours[i] for i in ids1]
assert( np.any( (contours[:-1] - contours[1:])**2 >= 1e-8 ) )
## TODO: keep closest beads beyond +-1.5 if there are fewer than 2 beads
tmp = []
dist = 1
while len(tmp) < 5 and dist < 3:
tmp = [ci for ci in zip(contours,ids) if np.abs(ci[0]-0.5) < dist]
tmp = list(filter(lambda bc: np.abs(bc[1]-0.5) < dist, beads_contours))
dist += 0.1
if len(tmp) <= 1:
raise Exception("Failed to fit spline into segment {}".format(s))
contours = [ci[0] for ci in tmp]
ids = [ci[1] for ci in tmp]
# cb = sorted( zip(contours,beads), key=lambda a:a[0] )
# beads = [b for c,b in cb]
# contours = [c for c,b in cb]
beads = [beads[i] for i in ids]
beads = [b for b,c in tmp]
contours = [c for b,c in tmp]
ids = [b.idx for b in beads]
if len(beads) <= 1:
......@@ -1929,11 +1969,13 @@ class SegmentModel(ArbdModel):
assert( np.isclose(a,0) or np.isclose(a,1) )
## TODO improve this for combinations of ssDNA and dsDNA (maybe a1/a2 should be calculated differently)
""" Search to see whether bead at location is already found """
b = None
if isinstance(s1,DoubleStrandedSegment):
b = s1.get_nearest_bead(a1)
b = s1.get_nearest_bead(a1)
if b is not None:
assert( b.parent is s1 )
""" if above assertion is true, no problem here """
if np.abs(b.get_nt_position(s1) - s1.contour_to_nt_pos(a1)) > 0.5:
b = None
if b is None and isinstance(s2,DoubleStrandedSegment):
......@@ -2014,6 +2056,7 @@ class SegmentModel(ArbdModel):
# pdb.set_trace()
""" Add intrahelical neighbors at connections """
for c,A,B in self.get_connections("intrahelical"):
b1,b2 = [l.particle for l in (A,B)]
if b1 is b2:
......@@ -2402,7 +2445,7 @@ class SegmentModel(ArbdModel):
c1,A1,B1,dir1 = crossovers[i]
c2,A2,B2,dir2 = crossovers[i+1]
s1,s2 = [l.container for l in (A1,A2)]
sep = A1.particle.get_nt_position(s1) - A2.particle.get_nt_position(s2)
sep = A1.particle.get_nt_position(s1,near_address=A1.address) - A2.particle.get_nt_position(s2,near_address=A2.address)
sep = np.abs(sep)
assert(sep >= 0)
......
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