diff --git a/cadnano_segments.py b/cadnano_segments.py
index e9d0f578bce5418077723d36087e8270ed754316..682b8b67d1a24b1058305fd6768396a441b8dccf 100644
--- a/cadnano_segments.py
+++ b/cadnano_segments.py
@@ -344,167 +344,6 @@ class cadnano_part(SegmentModel):
endList.append(lastStrand[1])
return endLists
- def combineRegionLists(self,loHi1,loHi2,intersect=False):
-
- """Combines two lists of (lo,hi) pairs specifying integer
- regions a single list of regions. """
-
- ## Validate input
- for l in (loHi1,loHi2):
- ## Assert each region in lists is sorted
- for pair in l:
- assert(len(pair) == 2)
- assert(pair[0] <= pair[1])
-
- if len(loHi1) == 0:
- if intersect:
- return []
- else:
- return loHi2
- if len(loHi2) == 0:
- if intersect:
- return []
- else:
- return loHi1
-
- ## Break input into lists of compact regions
- compactRegions1,compactRegions2 = [[],[]]
- for compactRegions,loHi in zip(
- [compactRegions1,compactRegions2],
- [loHi1,loHi2]):
- tmp = []
- lastHi = loHi[0][0]-1
- for lo,hi in loHi:
- if lo-1 != lastHi:
- compactRegions.append(tmp)
- tmp = []
- tmp.append((lo,hi))
- lastHi = hi
- if len(tmp) > 0:
- compactRegions.append(tmp)
-
- ## Build result
- result = []
- region = []
- i,j = [0,0]
- compactRegions1.append([[1e10]])
- compactRegions2.append([[1e10]])
- while i < len(compactRegions1)-1 or j < len(compactRegions2)-1:
- cr1 = compactRegions1[i]
- cr2 = compactRegions2[j]
-
- ## initialize region
- if len(region) == 0:
- if cr1[0][0] <= cr2[0][0]:
- region = cr1
- i += 1
- continue
- else:
- region = cr2
- j += 1
- continue
-
- if region[-1][-1] >= cr1[0][0]:
- region = self.combineCompactRegionLists(region, cr1, intersect=False)
- i+=1
- elif region[-1][-1] >= cr2[0][0]:
- region = self.combineCompactRegionLists(region, cr2, intersect=False)
- j+=1
- else:
- result.extend(region)
- region = []
-
- assert( len(region) > 0 )
- result.extend(region)
- result = sorted(result)
-
- # print("loHi1:",loHi1)
- # print("loHi2:",loHi2)
- # print(result,"\n")
-
- if intersect:
- lo = max( [loHi1[0][0], loHi2[0][0]] )
- hi = min( [loHi1[-1][1], loHi2[-1][1]] )
- result = [r for r in result if r[0] >= lo and r[1] <= hi]
-
- return result
-
- def combineCompactRegionLists(self,loHi1,loHi2,intersect=False):
- """Combines two lists of (lo,hi) pairs specifying regions within a
- compact integer set into a single list of regions.
-
- examples:
- loHi1 = [[0,4],[5,7]]
- loHi2 = [[2,4],[5,9]]
- out = [(0, 1), (2, 4), (5, 7), (8, 9)]
-
- loHi1 = [[0,3],[5,7]]
- loHi2 = [[2,4],[5,9]]
- out = [(0, 1), (2, 3), (4, 4), (5, 7), (8, 9)]
- """
-
- ## Validate input
- for l in (loHi1,loHi2):
- ## Assert each region in lists is sorted
- for pair in l:
- assert(len(pair) == 2)
- assert(pair[0] <= pair[1])
- ## Assert lists are compact
- for pair1,pair2 in zip(l[::2],l[1::2]):
- assert(pair1[1]+1 == pair2[0])
-
- if len(loHi1) == 0:
- if intersect:
- return []
- else:
- return loHi2
- if len(loHi2) == 0:
- if intersect:
- return []
- else:
- return loHi1
-
- ## Find the ends of the region
- lo = min( [loHi1[0][0], loHi2[0][0]] )
- hi = max( [loHi1[-1][1], loHi2[-1][1]] )
-
- ## Make a list of indices where each region will be split
- splitAfter = []
- for l,h in loHi2:
- if l != lo:
- splitAfter.append(l-1)
- if h != hi:
- splitAfter.append(h)
-
- for l,h in loHi1:
- if l != lo:
- splitAfter.append(l-1)
- if h != hi:
- splitAfter.append(h)
- splitAfter = sorted(list(set(splitAfter)))
-
- # print("splitAfter:",splitAfter)
-
- split=[]
- last = -2
- for s in splitAfter:
- split.append(s)
- last = s
-
- # print("split:",split)
- returnList = [(i+1,j) if i != j else (i,j) for i,j in zip([lo-1]+split,split+[hi])]
-
- if intersect:
- lo = max( [loHi1[0][0], loHi2[0][0]] )
- hi = min( [loHi1[-1][1], loHi2[-1][1]] )
- returnList = [r for r in returnList if r[0] >= lo and r[1] <= hi]
-
- # print("loHi1:",loHi1)
- # print("loHi2:",loHi2)
- # print(returnList,"\n")
- return returnList
-
-
def _helix_strands_to_segment_ranges(self, helix_strands):
"""Utility method to convert cadnano strand lists into list of
indices of terminal points"""
@@ -595,24 +434,40 @@ class cadnano_part(SegmentModel):
## TODO: handle nicks that are at intrahelical connections(?)
zmid1 = int(0.5*(r1[0]+r1[1]))
zmid2 = int(0.5*(r2[0]+r2[1]))
+ if seg1.name == "19-3" or seg2.name == "19-3":
+ import pdb
+ pdb.set_trace()
+
+ # if zMid in strandOccupancies[0] and zMid in strandOccupancies[1]:
+ # seg = DoubleStrandedSegment(**kwargs,**posargs1)
+ # elif zMid in strandOccupancies[0]:
+ # seg = SingleStrandedSegment(**kwargs,**posargs1)
+ # elif zMid in strandOccupancies[1]:
+ # seg = SingleStrandedSegment(**kwargs,**posargs2)
+
## TODO: validate
if zmid1 in occ[0] and zmid2 in occ[0]:
seg1.connect_end3(seg2.start5)
if zmid1 in occ[1] and zmid2 in occ[1]:
if zmid1 in occ[0]:
- seg2.connect_end3(seg1.end5)
+ end = seg1.end5
+ else:
+ end = seg1.start5
+ if zmid2 in occ[0]:
+ seg2.connect_start3(end)
else:
seg2.connect_end3(seg1.start5)
def _add_crossovers(self):
for h1,f1,z1,h2,f2,z2 in self.xover_list:
- if (h1 == 52 or h2 == 52) and z1 == 221:
- import pdb
- pdb.set_trace()
+ # if (h1 == 52 or h2 == 52) and z1 == 221:
+ # import pdb
+ # pdb.set_trace()
seg1, nt1 = self._get_segment_nucleotide(h1,z1)
seg2, nt2 = self._get_segment_nucleotide(h2,z2)
## TODO: use different types of crossovers
+ ## fwd?
seg1.add_crossover(nt1,seg2,nt2,[f1,f2])
def _add_prime_ends(self):
diff --git a/segmentmodel.py b/segmentmodel.py
index 58a2767c93ee2bf0fcb922c5be8b2708adc4fa3f..ca71dae9033434a79e4292286e7202713b127b09 100644
--- a/segmentmodel.py
+++ b/segmentmodel.py
@@ -1,3 +1,4 @@
+import pdb
import numpy as np
import random
from arbdmodel import PointParticle, ParticleType, Group, ArbdModel
@@ -36,12 +37,15 @@ TODO:
- add unit test of helices connected to themselves
"""
+class ParticleNotConnectedError(Exception):
+ pass
+
class Location():
""" Site for connection within an object """
def __init__(self, container, address, type_, on_fwd_strand = True):
## TODO: remove cyclic references(?)
self.container = container
- self.address = address # represents position along contour length in segments
+ self.address = address # represents position along contour length in segment
# assert( type_ in ("end3","end5") ) # TODO remove or make conditional
self.on_fwd_strand = on_fwd_strand
self.type_ = type_
@@ -69,8 +73,6 @@ class Location():
on_fwd = "on_fwd_strand"
else:
on_fwd = "on_rev_strand"
- # return "<Location in {} at contour {} {} with connection {}>".format( self.container.name, self.address, self.on_fwd_strand, self.connection )
- # return "<Location {} in {} at contour {} {} with connection {}>".format( self.type_, self.container.name, self.address, on_fwd, self.connection )
return "<Location {}.{}[{:.2f},{:d}]>".format( self.container.name, self.type_, self.address, self.on_fwd_strand)
class Connection():
@@ -89,6 +91,11 @@ class Connection():
return self.A
else:
raise Exception("OutOfBoundsError")
+
+ def __repr__(self):
+ return "<Connection {}--{}--{}]>".format( self.A, self.type_, self.B )
+
+
# class ConnectableElement(Transformable):
class ConnectableElement():
@@ -120,7 +127,7 @@ class ConnectableElement():
if l.on_fwd_strand == on_fwd_strand and l.connection is None:
assert(loc is None)
loc = l
- assert( loc is not None )
+ # assert( loc is not None )
else:
for l in self.locations:
if l.address == address and l.on_fwd_strand == on_fwd_strand:
@@ -173,31 +180,7 @@ class SegmentParticle(PointParticle):
PointParticle.__init__(self, type_, position, name=name, segname=segname, **kwargs)
self.intrahelical_neighbors = []
self.other_neighbors = []
-
- # def get_contour_position(self,seg):
- # assert( isinstance(seg,Segment) )
- # if seg == self.parent:
- # return self.contour_position
- # else:
- # ## TODO replace with something more elegant
- # for c,A,B in self.parent.get_connections_and_locations():
- # if A.particle is self and B.container is seg:
- # nt = np.abs( (self.contour_position - A.address)*(A.container.num_nts-1) )
- # if B.address < 0.5:
- # return B.address-nt/(seg.num_nts-1)
- # else:
- # return B.address+nt/(seg.num_nts-1)
- # ## ERROR
- # print("")
- # for c,A,B in self.parent.get_connections_and_locations():
- # print(" ",c.type_)
- # print(A,B)
- # print(A.particle,self)
- # print(B.container,seg)
- # print("")
- # import pdb
- # pdb.set_trace()
- # raise Exception("Did not find location for particle {} in Segment {}".format(self,seg))
+ self.locations = []
def get_intrahelical_above(self):
""" Returns bead directly above self """
@@ -212,45 +195,99 @@ class SegmentParticle(PointParticle):
for b in self.intrahelical_neighbors:
if b.get_contour_position(self.parent) < self.contour_position:
return b
-
- def get_nt_position(self,seg):
- if seg == self.parent:
- return seg.contour_to_nt_pos(self.contour_position)
+ def _neighbor_should_be_added(self,b):
+ c1 = self.contour_position
+ c2 = b.get_contour_position(self.parent)
+ if c2 < c1:
+ b0 = self.get_intrahelical_below()
else:
- cl = [e for e in self.parent.get_connections_and_locations() if e[2].container is seg]
- dc = [(self.contour_position - A.address)**2 for c,A,B in cl]
+ b0 = self.get_intrahelical_above()
+
+ if b0 is not None:
+ c0 = b0.get_contour_position(self.parent)
+ if np.abs(c2-c1) < np.abs(c0-c1):
+ ## remove b0
+ self.intrahelical_neighbors.remove(b0)
+ b0.intrahelical_neighbors.remove(self)
+ return True
+ else:
+ return False
+ return True
+
+ def make_intrahelical_neighbor(self,b):
+ add1 = self._neighbor_should_be_added(b)
+ add2 = b._neighbor_should_be_added(self)
+ if add1 and add2:
+ assert(len(b.intrahelical_neighbors) <= 1)
+ assert(len(self.intrahelical_neighbors) <= 1)
+ self.intrahelical_neighbors.append(b)
+ b.intrahelical_neighbors.append(self)
- if len(dc) == 0:
- pdb.set_trace()
+
+ # def get_nt_position(self,seg):
+ # if seg == self.parent:
+ # return seg.contour_to_nt_pos(self.contour_position)
+ # else:
+ # cl = [e for e in self.parent.get_connections_and_locations() if e[2].container is seg]
- i = np.argmin(dc)
- c,A,B = cl[i]
- ## TODO: generalize, removing np.abs and conditional
- delta_nt = np.abs( A.container.contour_to_nt_pos(self.contour_position - A.address) )
- B_nt_pos = seg.contour_to_nt_pos(B.address)
- if B.address < 0.5:
- return B_nt_pos-delta_nt
- else:
- return B_nt_pos+delta_nt
+ # dc = [(self.contour_position - A.address)**2 for c,A,B in cl]
+
+ # if len(dc) == 0:
+ # import pdb
+ # pdb.set_trace()
- def get_contour_position_old(self,seg):
+ # i = np.argmin(dc)
+ # c,A,B = cl[i]
+ # ## TODO: generalize, removing np.abs and conditional
+ # delta_nt = np.abs( A.container.contour_to_nt_pos(self.contour_position - A.address) )
+ # B_nt_pos = seg.contour_to_nt_pos(B.address)
+ # if B.address < 0.5:
+ # return B_nt_pos-delta_nt
+ # else:
+ # return B_nt_pos+delta_nt
+
+ def get_nt_position(self,seg):
if seg == self.parent:
- return self.contour_position
+ return seg.contour_to_nt_pos(self.contour_position)
else:
- cl = [e for e in self.parent.get_connections_and_locations() in B.container is seg]
- dc = [(self.contour_position - A.address)**2 for c,A,B in e]
- if len(dc) == 0:
- pdb.set_trace()
+ def get_nt_pos(contour1, seg1, seg2):
+ cl = [e for e in seg1.get_connections_and_locations() if e[2].container is seg2]
+ dc = [(contour1 - A.address)**2 for c,A,B in cl]
+ if len(dc) == 0: return None
- i = np.argmin(dc)
+ i = np.argmin(dc)
+ c,A,B = cl[i]
+
+ ## TODO: generalize, removing np.abs and conditional
+ delta_nt = np.abs( seg1.contour_to_nt_pos(contour1 - A.address) )
+ B_nt_pos = seg2.contour_to_nt_pos(B.address)
+ if B.address < 0.5:
+ return B_nt_pos-delta_nt
+ else:
+ return B_nt_pos+delta_nt
+
+ pos = get_nt_pos(self.contour_position, self.parent, seg)
+ if pos is None:
+ ## Particle is not directly connected
+ visited_segs = set(seg)
+ positions = []
+ for l in self.locations:
+ if l.container == self.parent: continue
+ pos0 = get_nt_pos(self.contour_position, self.parent, l.container)
+ assert(pos0 is not None)
+ pos0 = l.container.nt_pos_to_contour(pos0)
+ pos = get_nt_pos( pos0, l.container, seg )
+ if pos is not None:
+ positions.append( pos )
+ assert( len(positions) > 0 )
+ if len(positions) > 1:
+ import pdb
+ pdb.set_trace()
+ pos = positions[0]
+ return pos
- nt = np.abs( (self.contour_position - A.address)*(A.container.num_nts-1) )
- if B.address < 0.5:
- return seg.nt_pos_to_contour(B.address-nt)
- else:
- return seg.nt_pos_to_contour(B.address+nt)
def get_contour_position(self,seg):
if seg == self.parent:
@@ -323,18 +360,14 @@ class Segment(ConnectableElement, Group):
loc.particle = None
def contour_to_nt_pos(self, contour_pos, round_nt=False):
- nt = contour_pos*(self.num_nts-1)
+ nt = contour_pos*(self.num_nts) - 0.5
if round_nt:
- assert( (np.around(nt) - nt)**2 < 1e-3 )
+ assert( np.isclose(np.around(nt),nt) )
nt = np.around(nt)
return nt
def nt_pos_to_contour(self,nt_pos):
- if self.num_nts == 1:
- assert(nt_pos == 0)
- return 0
- else:
- return nt_pos/(self.num_nts-1)
+ return (nt_pos+0.5)/(self.num_nts)
def contour_to_position(self,s):
p = interpolate.splev( s, self.position_spline_params )
@@ -482,14 +515,24 @@ class Segment(ConnectableElement, Group):
move_at_least = 0
## Iterate through locations
+ # locations = sorted(self.locations, key=lambda l:(l.address,not l.on_fwd_strand), reverse=(not is_fwd))
locations = sorted(self.locations, key=lambda l:(l.address,not l.on_fwd_strand), reverse=(not is_fwd))
# print(locations)
for l in locations:
- pos = self.contour_to_nt_pos(l.address, round_nt=True)
+ # TODOTODO probably okay
+ if l.address == 0:
+ pos = 0.0
+ elif l.address == 1:
+ pos = self.num_nts-1
+ else:
+ pos = self.contour_to_nt_pos(l.address, round_nt=True)
## DEBUG
+ # import pdb
+ # pdb.set_trace()
+
## Skip locations encountered before our strand
# tol = 0.1
# if is_fwd:
@@ -514,7 +557,15 @@ class Segment(ConnectableElement, Group):
if l.on_fwd_strand == is_fwd:
print(" passing through",l)
print("from {}, connection {} to {}".format(nt_pos,l,B))
- Bpos = B.container.contour_to_nt_pos(B.address, round_nt=True)
+ try:
+ Bpos = B.container.contour_to_nt_pos(B.address, round_nt=True)
+ except:
+ if B.address == 0:
+ Bpos = 0
+ elif B.address == 1:
+ Bpos = B.container.num_nts-1
+ else:
+ raise
return pos, B.container, Bpos, B.on_fwd_strand, 0.5
## Stop at other strand crossovers so basepairs line up
@@ -530,66 +581,6 @@ class Segment(ConnectableElement, Group):
## Made it to the end of the segment without finding a connection
return 1*is_fwd, None, None, None
-
- def get_end_of_strand_old(self, contour_pos, is_fwd):
- """ Walks through locations, checking for crossovers """
-
- ## Iterate through locations
- # for l in self.locations:
- def loc_iter():
- locations = sorted(self.locations, key=lambda l:(l.address,not l.on_fwd_strand), reverse=(not is_fwd))
- # if is_fwd:
- for l in locations:
- yield l
- # else:
- # for l in locations[::-1]:
- # yield l
-
- for l in loc_iter():
- # if l.particle is None:
- # pos = l.address
- # else:
- # pos = l.particle.get_contour_position()
- pos = l.address
-
- ## DEBUG
- # if self.name == "1-0" and is_fwd == False:
- # import pdb
- # pdb.set_trace()
-
- ## Skip locations encountered before our strand
- if is_fwd:
- if pos <= contour_pos: continue
- elif pos >= contour_pos: continue
-
- # print(" ?",l)
-
- ## Stop if we found the 3prime end
- if l.on_fwd_strand == is_fwd and l.type_ == "3prime":
- return pos, None, None, None
-
- ## Check location connections
- c = l.connection
- if c is None: continue
- B = c.other(l)
-
- ## Found a location on the same strand?
- if l.on_fwd_strand == is_fwd:
- # print(" passing through",l)
- # print("from {}, connection {} to {}".format(contour_pos,l,B))
- return pos, B.container, B.address, B.on_fwd_strand
-
- ## Stop at other strand crossovers so basepairs line up
- elif c.type_ == "crossover":
- # print(" pausing at",l)
- # print("pausing at {}".format(l))
- return pos, l.container, pos, is_fwd
-
- raise Exception("Shouldn't be here")
- # print("Shouldn't be here")
- ## Made it to the end of the segment without finding a connection
- return 1*is_fwd, None, None, None
-
def get_nearest_bead(self, contour_position):
if len(self.beads) < 1: return None
cs = np.array([b.contour_position for b in self.beads]) # TODO: cache
@@ -674,17 +665,19 @@ class Segment(ConnectableElement, Group):
for b in existing_beads:
assert(b.parent is not None)
+ # if self.name == "1-1":
+ # import pdb
+ # pdb.set_trace()
+
## Add ends if they don't exist yet
## TODOTODO: test 1 nt segments?
- if len(existing_beads) == 0 or existing_beads[0].get_contour_position(self) > 0:
- if len(existing_beads) > 0:
- assert(existing_beads[0].get_nt_position(self) >= 0.5)
-
+ if len(existing_beads) == 0 or existing_beads[0].get_nt_position(self) > 0.5:
+ # if len(existing_beads) > 0:
+ # assert(existing_beads[0].get_nt_position(self) >= 0.5)
b = self._generate_one_bead(0, 0)
existing_beads = [b] + existing_beads
- if existing_beads[-1].get_contour_position(self) < 1:
- # assert((1-existing_beads[0].get_contour_position(self))*(self.num_nts-1) >= 0.5)
- assert(self.num_nts-1-existing_beads[0].get_nt_position(self) >= 0.5)
+
+ if existing_beads[-1].get_nt_position(self)-(self.num_nts-1) < -0.5:
b = self._generate_one_bead(1, 0)
existing_beads.append(b)
assert(len(existing_beads) > 1)
@@ -694,11 +687,13 @@ 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)]
- if eb1 is eb2:
- pdb.set_trace()
assert( eb1 is not eb2 )
- # print(" %s working on %d to %d" % (self.name, eb1.position[2], eb2.position[2]))
+ # 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)
num_beads = self._get_num_beads( e_ds, max_basepairs_per_bead, max_nucleotides_per_bead )
ds = e_ds / (num_beads+1)
@@ -712,26 +707,17 @@ class Segment(ConnectableElement, Group):
s0 = eb1.get_contour_position(self)
if last is not None:
- last.intrahelical_neighbors.append(eb1)
- eb1.intrahelical_neighbors.append(last)
- assert(len(last.intrahelical_neighbors) <= 2)
- assert(len(eb1.intrahelical_neighbors) <= 2)
+ last.make_intrahelical_neighbor(eb1)
last = eb1
for j in range(num_beads):
s = ds*(j+1) + s0
b = self._generate_one_bead(s,nts)
- last.intrahelical_neighbors.append(b)
- b.intrahelical_neighbors.append(last)
- assert(len(last.intrahelical_neighbors) <= 2)
- assert(len(b.intrahelical_neighbors) <= 2)
+ last.make_intrahelical_neighbor(b)
last = b
tmp_children.append(b)
- last.intrahelical_neighbors.append(eb2)
- eb2.intrahelical_neighbors.append(last)
- assert(len(last.intrahelical_neighbors) <= 2)
- assert(len(eb2.intrahelical_neighbors) <= 2)
+ last.make_intrahelical_neighbor(eb2)
if eb2.parent == self:
tmp_children.append(eb2)
@@ -778,8 +764,8 @@ class DoubleStrandedSegment(Segment):
self.end = self.end3 = Location( self, address=1, type_ = "end3" )
self.end5 = Location( self, address=1, type_= "end5", on_fwd_strand=False )
- for l in (self.start5,self.start3,self.end3,self.end5):
- self.locations.append(l)
+ # for l in (self.start5,self.start3,self.end3,self.end5):
+ # self.locations.append(l)
## Set up interpolation for azimuthal angles
a = np.array([self.start_position,self.end_position]).T
@@ -810,13 +796,13 @@ class DoubleStrandedSegment(Segment):
end3 = end3.end3
self._connect_ends( self.end5, end3, type_, force_connection = force_connection )
- def add_crossover(self, nt, other, other_nt, strands_fwd=[True,False]):
+ def add_crossover(self, nt, other, other_nt, strands_fwd=[True,False], nt_on_5prime=True):
""" Add a crossover between two helices """
## Validate other, nt, other_nt
## TODO
if isinstance(other,SingleStrandedSegment):
- other.add_crossover(other_nt, self, nt, strands_fwd[::-1])
+ other.add_crossover(other_nt, self, nt, strands_fwd[::-1], not nt_on_5prime)
else:
## Create locations, connections and add to segments
@@ -826,12 +812,16 @@ class DoubleStrandedSegment(Segment):
loc = self.get_location_at(c, strands_fwd[0])
c = other.nt_pos_to_contour(other_nt)
+ # TODOTODO: may need to subtract or add a little depending on 3prime/5prime
assert(c >= 0 and c <= 1)
other_loc = other.get_location_at(c, strands_fwd[1])
self._connect(other, Connection( loc, other_loc, type_="crossover" ))
- loc.is_3prime_side_of_connection = not strands_fwd[0]
- other_loc.is_3prime_side_of_connection = not strands_fwd[1]
-
+ if nt_on_5prime:
+ loc.is_3prime_side_of_connection = False
+ other_loc.is_3prime_side_of_connection = True
+ else:
+ loc.is_3prime_side_of_connection = True
+ other_loc.is_3prime_side_of_connection = False
## Real work
def _connect_ends(self, end1, end2, type_, force_connection):
@@ -842,7 +832,7 @@ class DoubleStrandedSegment(Segment):
assert( end.type_ in ("end3","end5") )
assert( end1.type_ != end2.type_ )
## Create and add connection
- if end2.type_ == "end3":
+ if end2.type_ == "end5":
end1.container._connect( end2.container, Connection( end1, end2, type_=type_ ), in_3prime_direction=True )
else:
end2.container._connect( end1.container, Connection( end2, end1, type_=type_ ), in_3prime_direction=True )
@@ -888,58 +878,74 @@ class SingleStrandedSegment(Segment):
self.start = self.start5 = Location( self, address=0, type_= "end5" ) # TODO change type_?
self.end = self.end3 = Location( self, address=1, type_ = "end3" )
- for l in (self.start5,self.end3):
- self.locations.append(l)
+ # for l in (self.start5,self.end3):
+ # self.locations.append(l)
def connect_end3(self, end5, force_connection=False):
- self._connect_end( end5, _5_to_3 = False, force_connection = force_connection )
+ self._connect_end( end5, _5_to_3 = True, force_connection = force_connection )
def connect_5end(self, end3, force_connection=False): # TODO: change name or possibly deprecate
- self._connect_end( end3, _5_to_3 = True, force_connection = force_connection )
+ self._connect_end( end3, _5_to_3 = False, force_connection = force_connection )
def _connect_end(self, other, _5_to_3, force_connection):
assert( isinstance(other, Location) )
if _5_to_3 == True:
- my_end = self.end5
- assert( other.type_ == "end3" )
+ my_end = self.end3
+ # assert( other.type_ == "end5" )
+ if (other.type_ is not "end5"):
+ print("Warning: code does not prevent connecting 3prime to 3prime, etc")
conn = Connection( my_end, other, type_="intrahelical" )
self._connect( other.container, conn, in_3prime_direction=True )
else:
- my_end = self.end3
- assert( other.type_ == "end5" )
+ my_end = self.end5
+ # assert( other.type_ == "end3" )
+ if (other.type_ is not "end3"):
+ print("Warning: code does not prevent connecting 3prime to 3prime, etc")
conn = Connection( other, my_end, type_="intrahelical" )
other.container._connect( self, conn, in_3prime_direction=True )
- def add_crossover(self, nt, other, other_nt, strands_fwd=[True,False]):
+ def add_crossover(self, nt, other, other_nt, strands_fwd=[True,False], nt_on_5prime=True):
""" Add a crossover between two helices """
## Validate other, nt, other_nt
## TODO
## TODO: fix direction
- c1 = self.nt_pos_to_contour(nt)
- ## Ensure connections occur at ends, otherwise the structure doesn't make sense
- assert(np.isclose(c1,0) or np.isclose(c1,1))
+ # c1 = self.nt_pos_to_contour(nt)
+ # # TODOTODO
+ # ## Ensure connections occur at ends, otherwise the structure doesn't make sense
+ # # assert(np.isclose(c1,0) or np.isclose(c1,1))
+ # assert(np.isclose(nt,0) or np.isclose(nt,self.num_nts-1))
+ if nt == 0:
+ c1 = 0
+ elif nt == self.num_nts-1:
+ c1 = 1
+ else:
+ raise Exception("Crossovers can only be at the ends of an ssDNA segment")
loc = self.get_location_at(c1, True)
- c2 = other.nt_pos_to_contour(other_nt)
+ if other_nt == 0:
+ c2 = 0
+ elif other_nt == other.num_nts-1:
+ c2 = 1
+ else:
+ c2 = other.nt_pos_to_contour(other_nt)
+
if isinstance(other,SingleStrandedSegment):
## Ensure connections occur at opposing ends
- assert(np.isclose(c2,0) or np.isclose(c2,1))
+ assert(np.isclose(other_nt,0) or np.isclose(other_nt,self.num_nts-1))
other_loc = other.get_location_at( c2, True )
- assert( loc.type_ in ("end3","end5"))
- assert( other_loc.type_ in ("end3","end5"))
- if loc.type_ == "end3":
+ if ("22-2" in (self.name, other.name)):
+ pdb.set_trace()
+ if nt_on_5prime:
self.connect_end3( other_loc )
else:
- assert( other_loc.type_ == "end3" )
other.connect_end3( self )
else:
- assert( loc.type_ in ("end3","end5"))
assert(c2 >= 0 and c2 <= 1)
other_loc = other.get_location_at( c2, strands_fwd[1] )
- if loc.type_ == "end3":
+ if nt_on_5prime:
self._connect(other, Connection( loc, other_loc, type_="sscrossover" ), in_3prime_direction=True )
else:
other._connect(self, Connection( other_loc, loc, type_="sscrossover" ), in_3prime_direction=True )
@@ -972,7 +978,7 @@ class StrandInSegment(Group):
nts = np.abs(end-start)+1
self.num_nts = int(round(nts))
- assert( np.abs(self.num_nts-nts) < 1e-5 )
+ assert( np.isclose(self.num_nts,nts) )
# print(" Creating {}-nt StrandInSegment in {} from {} to {} {}".format(self.num_nts, segment.name, start, end, is_fwd))
@@ -1008,6 +1014,7 @@ class Strand(Group):
## TODO disambiguate names of functions
def add_dna(self, segment, start, end, is_fwd):
+ # TODOTODO use nt pos ?
""" start/end should be provided expressed as contour_length, is_fwd tuples """
if not (segment.contour_to_nt_pos(np.abs(start-end)) > 0.9):
import pdb
@@ -1375,23 +1382,44 @@ class SegmentModel(ArbdModel):
s1,s2 = [l.container for l in (A,B)]
## TODO be more elegant!
- if A.on_fwd_strand == False: continue # TODO verify this avoids double-counting
+ # if isinstance(s1, DoubleStrandedSegment) and isinstance(s2, DoubleStrandedSegment) and A.on_fwd_strand == False: continue
+ if isinstance(s1, DoubleStrandedSegment) and isinstance(s2, DoubleStrandedSegment) and A.on_fwd_strand == False: continue
+ ## if A.on_fwd_strand == False: continue # TODO verify this avoids double-counting
assert( A.particle is None )
assert( B.particle is None )
## TODO: offload the work here to s1
+ # TODOTODO
a1,a2 = [l.address for l in (A,B)]
# a1,a2 = [a - s.nt_pos_to_contour(0.5) if a == 0 else a + s.nt_pos_to_contour(0.5) for a,s in zip((a1,a2),(s1,s2))]
for a in (a1,a2):
- assert( a in (0,1,0.0,1.0) )
- a1,a2 = [a - s.nt_pos_to_contour(0.5) if a == 0 else a + s.nt_pos_to_contour(0.5) for a,s in zip((a1,a2),(s1,s2))]
+ assert( np.isclose(a,0) or np.isclose(a,1) )
+ # a1,a2 = [a - s.nt_pos_to_contour(0) if a == 0 else a + s.nt_pos_to_contour(0) for a,s in zip((a1,a2),(s1,s2))]
## TODO improve this for combinations of ssDNA and dsDNA (maybe a1/a2 should be calculated differently)
- if isinstance(s2,DoubleStrandedSegment):
- b = s2._generate_one_bead(a2,0)
- else:
- b = s1._generate_one_bead(a1,0)
+
+ if isinstance(s1,DoubleStrandedSegment):
+ b = s1.get_nearest_bead(a1)
+ if b is not None:
+ if np.abs(b.get_nt_position(s1) - s1.contour_to_nt_pos(a1)) > 0.5:
+ b = None
+ elif isinstance(s2,DoubleStrandedSegment):
+ b = s2.get_nearest_bead(a2)
+ if b is not None:
+ if np.abs(b.get_nt_position(s2) - s2.contour_to_nt_pos(a2)) > 0.5:
+ b = None
+
+ if b is not None and b.parent not in (s1,s2):
+ b = None
+
+ if b is None:
+ ## need to generate a bead
+ if isinstance(s2,DoubleStrandedSegment):
+ b = s2._generate_one_bead(a2,0)
+ else:
+ b = s1._generate_one_bead(a1,0)
A.particle = B.particle = b
+ b.locations.extend([A,B])
""" Generate beads at other junctions """
for c,A,B in self.get_connections(exclude="intrahelical"):
@@ -1406,21 +1434,23 @@ class SegmentModel(ArbdModel):
if A.particle is None:
b = s1.get_nearest_bead(a1)
- if b is not None and s1.contour_to_nt_pos(np.abs(b.contour_position-a1)) < 1.9:
+ if b is not None and s1.contour_to_nt_pos(np.abs(b.contour_position-a1)) < 1:
## combine beads
b.contour_position = 0.5*(b.contour_position + a1) # avg position
- A.particle = b
else:
- A.particle = s1._generate_one_bead(a1,0)
+ b = s1._generate_one_bead(a1,0)
+ A.particle = b
+ b.locations.append(A)
if B.particle is None:
b = s2.get_nearest_bead(a2)
- if b is not None and s2.contour_to_nt_pos(np.abs(b.contour_position-a2)) < 1.9:
+ if b is not None and s2.contour_to_nt_pos(np.abs(b.contour_position-a2)) < 19:
## combine beads
b.contour_position = 0.5*(b.contour_position + a2) # avg position
- B.particle = b
else:
- B.particle = s2._generate_one_bead(a2,0)
+ b = s2._generate_one_bead(a2,0)
+ B.particle = b
+ b.locations.append(B)
""" Some tests """
for c,A,B in self.get_connections("intrahelical"):
@@ -1445,11 +1475,7 @@ class SegmentModel(ArbdModel):
continue
else:
for b in (b1,b2): assert( b is not None )
- if b2 not in b1.intrahelical_neighbors:
- b1.intrahelical_neighbors.append(b2)
- b2.intrahelical_neighbors.append(b1)
- assert(len(b1.intrahelical_neighbors) <= 2)
- assert(len(b2.intrahelical_neighbors) <= 2)
+ b1.make_intrahelical_neighbor(b2)
""" Reassign bead types """
if self.DEBUG: print("Assigning bead types")
@@ -1760,6 +1786,7 @@ class SegmentModel(ArbdModel):
for i in range(len(conn_locs)):
c,A,B = conn_locs[i]
## TODO: handle change of direction
+ # TODOTODO
address = 1*(direction==-1)
if A.address == address and segment_is_new_helix(B.container):
new_s = B.container
@@ -1809,39 +1836,40 @@ class SegmentModel(ArbdModel):
for seg in self.segments:
## TODO move into Segment calls
import pdb
- if seg.start5.connection is None:
- add_end = True
- for l in seg.get_locations("5prime"):
- if l.address == 0 and l.on_fwd_strand:
- add_end = False
- break
- if add_end:
- seg.add_5prime(0)
- if 'end5' in seg.__dict__ and seg.end5.connection is None:
- add_end = True
- for l in seg.get_locations("5prime"):
- if l.address == 1 and (l.on_fwd_strand is False):
- add_end = False
- break
- if add_end:
- seg.add_5prime(seg.num_nts-1,on_fwd_strand=False)
-
- if 'start3' in seg.__dict__ and seg.start3.connection is None:
- add_end = True
- for l in seg.get_locations("3prime"):
- if l.address == 0 and (l.on_fwd_strand is False):
- add_end = False
- break
- if add_end:
- seg.add_3prime(0,on_fwd_strand=False)
- if seg.end3.connection is None:
- add_end = True
- for l in seg.get_locations("3prime"):
- if l.address == 1 and l.on_fwd_strand:
- add_end = False
- break
- if add_end:
- seg.add_3prime(seg.num_nts-1)
+ if False: # TODO: Make this happen conditionally
+ if seg.start5.connection is None:
+ add_end = True
+ for l in seg.get_locations("5prime"):
+ if l.address == 0 and l.on_fwd_strand:
+ add_end = False
+ break
+ if add_end:
+ seg.add_5prime(0)
+ if 'end5' in seg.__dict__ and seg.end5.connection is None:
+ add_end = True
+ for l in seg.get_locations("5prime"):
+ if l.address == 1 and (l.on_fwd_strand is False):
+ add_end = False
+ break
+ if add_end:
+ seg.add_5prime(seg.num_nts-1,on_fwd_strand=False)
+
+ if 'start3' in seg.__dict__ and seg.start3.connection is None:
+ add_end = True
+ for l in seg.get_locations("3prime"):
+ if l.address == 0 and (l.on_fwd_strand is False):
+ add_end = False
+ break
+ if add_end:
+ seg.add_3prime(0,on_fwd_strand=False)
+ if seg.end3.connection is None:
+ add_end = True
+ for l in seg.get_locations("3prime"):
+ if l.address == 1 and l.on_fwd_strand:
+ add_end = False
+ break
+ if add_end:
+ seg.add_3prime(seg.num_nts-1)
# print( [(l,l.get_connected_location()) for l in seg.locations] )
# addresses = np.array([l.address for l in seg.get_locations("5prime")])
@@ -1857,7 +1885,12 @@ class SegmentModel(ArbdModel):
import pdb
pdb.set_trace()
s,seg = [strand, segment]
-
+
+ #if seg.name == "22-1" and pos > 140:
+ if seg.name == "22-2":
+ import pdb
+ pdb.set_trace()
+
end_pos, next_seg, next_pos, next_dir, move_at_least = seg.get_strand_segment(pos, is_fwd, move_at_least)
s.add_dna(seg, pos, end_pos, is_fwd)
@@ -1871,6 +1904,7 @@ class SegmentModel(ArbdModel):
# for pos, is_fwd in locs:
for l in locs:
print("Tracing",l)
+ # TODOTODO
pos = seg.contour_to_nt_pos(l.address, round_nt=True)
is_fwd = l.on_fwd_strand
s = Strand()