Skip to content
GitLab
Commit 9063e58f authored by cmaffeo2's avatar cmaffeo2
Browse files

Fixed numerous bugs relating to strand routing

parent 5cf9276c
Hide whitespace changes
Inline Side-by-side
cadnano_segments.py
View file @ 9063e58f
......@@ -257,9 +257,10 @@ class cadnano_part(SegmentModel):
## Build helices
for hid in range(numHID):
# print("Working on helix",hid)
print("Working on helix",hid)
helices[hid] = []
helix_ranges[hid] = []
self.strand_occupancies[hid] = []
helixStrands = strand_list[hid]
if helixStrands is None: continue
......@@ -283,8 +284,8 @@ class cadnano_part(SegmentModel):
ends1,ends2 = [ [(e[i],e[i+1]) for i in range(0,len(e),2)] for e in (ends1,ends2) ]
regions = combineRegionLists(ends1,ends2)
if hid == 43:
import pdb
# if hid == 43:
# import pdb
for zid1,zid2 in regions:
zMid = int(0.5*(zid1+zid2))
assert( zMid in strandOccupancies[0] or zMid in strandOccupancies[1] )
......@@ -434,16 +435,10 @@ 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)
# if seg1.name == "1-1" or seg2.name == "1-1":
# # if seg1.name == "19-3" or seg2.name == "19-3":
# import pdb
# pdb.set_trace()
## TODO: validate
if zmid1 in occ[0] and zmid2 in occ[0]:
......@@ -457,7 +452,7 @@ class cadnano_part(SegmentModel):
if zmid2 in occ[0]:
seg2.connect_start3(end)
else:
seg2.connect_end3(seg1.start5)
seg2.connect_end3(end)
def _add_crossovers(self):
for h1,f1,z1,h2,f2,z2 in self.xover_list:
......@@ -473,11 +468,12 @@ class cadnano_part(SegmentModel):
def _add_prime_ends(self):
for h,fwd,z in self._5prime_list:
seg, nt = self._get_segment_nucleotide(h,z)
print(seg.name,nt,fwd)
print("adding 5prime",seg.name,nt,fwd)
seg.add_5prime(nt,fwd)
for h,fwd,z in self._3prime_list:
seg, nt = self._get_segment_nucleotide(h,z)
print("adding 3prime",seg.name,nt,fwd)
seg.add_3prime(nt,fwd)
def get_bead(self, hid, zid):
......
run.py
View file @ 9063e58f
......@@ -18,4 +18,7 @@ if __name__ == '__main__':
except:
print("WARNING: skipping unreadable json file {}".format(f))
continue
run_simulation_protocol( out, "job-id", data, gpu=0 )
try:
run_simulation_protocol( out, "job-"+out+"-", data, gpu=0 )
except:
print("ERROR: failed to simulate {}".format(f))
segmentmodel.py
View file @ 9063e58f
......@@ -68,6 +68,18 @@ class Location():
self.connection = connection # TODO weakref?
self.is_3prime_side_of_connection = is_3prime_side_of_connection
def get_nt_pos(self):
try:
pos = self.container.contour_to_nt_pos(self.address, round_nt=True)
except:
if self.address == 0:
pos = 0
elif self.address == 1:
pos = self.container.num_nts-1
else:
raise
return pos
def __repr__(self):
if self.on_fwd_strand:
on_fwd = "on_fwd_strand"
......@@ -460,7 +472,7 @@ class Segment(ConnectableElement, Group):
for a in atoms:
a.position = scale*a.position
a.beta = 0
atoms.position = pos - atoms.atoms_by_name["C1'"].collapsed_position()
atoms.position = pos - atoms.atoms_by_name["C1'"].collapsedPosition()
else:
if scale is not None and scale != 1:
if atoms.sequence in ("A","G"):
......@@ -485,9 +497,13 @@ class Segment(ConnectableElement, Group):
## TODO? Replace with abstract strand-based model?
def add_5prime(self, nt, on_fwd_strand=True):
if isinstance(self,SingleStrandedSegment):
on_fwd_strand = True
self.add_location(nt,"5prime",on_fwd_strand)
def add_3prime(self, nt, on_fwd_strand=True):
if isinstance(self,SingleStrandedSegment):
on_fwd_strand = True
self.add_location(nt,"3prime",on_fwd_strand)
def get_3prime_locations(self):
......@@ -516,7 +532,12 @@ class Segment(ConnectableElement, Group):
## 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))
def loc_rank(l):
nt = l.get_nt_pos()
## optionally add logic about type of connection
return (nt, not l.on_fwd_strand)
# locations = sorted(self.locations, key=lambda l:(l.address,not l.on_fwd_strand), reverse=(not is_fwd))
locations = sorted(self.locations, key=loc_rank, reverse=(not is_fwd))
# print(locations)
for l in locations:
......@@ -530,8 +551,6 @@ class Segment(ConnectableElement, Group):
## DEBUG
# import pdb
# pdb.set_trace()
## Skip locations encountered before our strand
# tol = 0.1
......@@ -557,15 +576,7 @@ class Segment(ConnectableElement, Group):
if l.on_fwd_strand == is_fwd:
print(" passing through",l)
print("from {}, connection {} to {}".format(nt_pos,l,B))
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
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
......@@ -935,8 +946,8 @@ class SingleStrandedSegment(Segment):
## Ensure connections occur at opposing ends
assert(np.isclose(other_nt,0) or np.isclose(other_nt,self.num_nts-1))
other_loc = other.get_location_at( c2, True )
if ("22-2" in (self.name, other.name)):
pdb.set_trace()
# if ("22-2" in (self.name, other.name)):
# pdb.set_trace()
if nt_on_5prime:
self.connect_end3( other_loc )
else:
......@@ -1017,16 +1028,17 @@ class Strand(Group):
# 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
pdb.set_trace()
print( "WARNING: segment constructed with a very small number of nts ({})".format(segment.contour_to_nt_pos(np.abs(start-end))) )
# import pdb
# pdb.set_trace()
for s in self.strand_segments:
if s.segment == segment and s.is_fwd == is_fwd:
# assert( s.start not in (start,end) )
# assert( s.end not in (start,end) )
if s.start in (start,end) or s.end in (start,end):
print(" CIRCULAR DNA")
import pdb
pdb.set_trace()
print(" CIRCULAR DNA")
s = StrandInSegment( segment, start, end, is_fwd )
self.add( s )
......@@ -1281,11 +1293,23 @@ class SegmentModel(ArbdModel):
## 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]
contours = []
for b in beads:
try:
## might fail if s is an entire segment away
contours.append(b.get_contour_position(s))
except:
## ignore failed attempts
beads.remove(b)
contours = [b.get_contour_position(s) for b in beads]
cb = sorted( zip(contours,beads), key=lambda a:a[0] )
beads = [b for c,b in cb]
contours = [c for c,b in cb]
contours,ids = np.unique(contours, return_index=True)
# 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]
ids = [b.idx for b in beads]
......@@ -1300,7 +1324,7 @@ class SegmentModel(ArbdModel):
s.position_spline_params = tck
""" Get twist """
cb = [e for e in cb if 'orientation_bead' in e[1].__dict__]
cb = [e for e in zip(contours,beads) if 'orientation_bead' in e[1].__dict__]
beads = [b for c,b in cb]
contours = [c for c,b in cb]
ids = [b.idx for b in beads]
......@@ -1380,7 +1404,7 @@ class SegmentModel(ArbdModel):
## Loop through all connections, generating beads at appropriate locations
for c,A,B in self.get_connections("intrahelical"):
s1,s2 = [l.container for l in (A,B)]
print("Working on {}".format(c))
## TODO be more elegant!
# 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
......@@ -1421,6 +1445,8 @@ class SegmentModel(ArbdModel):
A.particle = B.particle = b
b.locations.extend([A,B])
# pdb.set_trace()
""" Generate beads at other junctions """
for c,A,B in self.get_connections(exclude="intrahelical"):
s1,s2 = [l.container for l in (A,B)]
......@@ -1444,7 +1470,7 @@ class SegmentModel(ArbdModel):
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)) < 19:
if b is not None and s2.contour_to_nt_pos(np.abs(b.contour_position-a2)) < 1:
## combine beads
b.contour_position = 0.5*(b.contour_position + a2) # avg position
else:
......@@ -1482,7 +1508,8 @@ class SegmentModel(ArbdModel):
beadtype_s = dict()
for segment in segments:
for b in segment:
b.num_nts = np.around( b.num_nts, decimals=1 )
# b.num_nts = np.around( b.num_nts, decimals=1 )
b.num_nts = 3*np.around( float(b.num_nts)/3, decimals=1 )
key = (b.type_.name[0].upper(), b.num_nts)
if key in beadtype_s:
b.type_ = beadtype_s[key]
......@@ -1628,6 +1655,14 @@ class SegmentModel(ArbdModel):
pot = self.get_bond_potential(4,18.5)
self.add_bond(b1,b2, pot)
""" Add connection potentials """
for c,A,B in self.get_connections("sscrossover"):
b1,b2 = [loc.particle for loc in (c.A,c.B)]
## TODO: improve parameters
pot = self.get_bond_potential(4,6)
self.add_bond(b1,b2, pot)
crossover_bead_pots = set()
for c,A,B in self.get_connections("crossover"):
......@@ -1881,15 +1916,15 @@ class SegmentModel(ArbdModel):
""" Build strands from connectivity of helices """
def _recursively_build_strand(strand, segment, pos, is_fwd, mycounter=0, move_at_least=0.5):
mycounter+=1
if mycounter > 1000:
if mycounter > 10000:
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()
# 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)
......
Supports Markdown
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