From 2c07f592f83b1237f77775481faf6b88d6d94e30 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Sun, 20 Jan 2019 00:55:46 -0600
Subject: [PATCH] Made it so a bead can be located at multiple distant
locations within a segment (e.g. for looped constructs)
---
mrdna/segmentmodel.py | 185 ++++++++++++++++++++++++++----------------
1 file changed, 114 insertions(+), 71 deletions(-)
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index 00cb905..5e698c4 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -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)
--
GitLab