Skip to content
GitLab
Commit 5cf9276c authored by cmaffeo2's avatar cmaffeo2
Browse files

Sorta got something to work

parent 6194b684
Hide whitespace changes
Inline Side-by-side
cadnano_segments.py
View file @ 5cf9276c
......@@ -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):
......
segmentmodel.py
View file @ 5cf9276c
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)