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

Improved scripts for building cadnano structures; seems it will be neccesary...

Improved scripts for building cadnano structures; seems it will be neccesary to make contour go from nt[0]-0.5 to nt[-1]+0.5
parent dc623af9
This diff is collapsed.
# -*- coding: utf-8 -*-
import sys
from glob import glob
import re
from cadnano_segments import read_json_file, run_simulation_protocol
if __name__ == '__main__':
for f in glob('json/*'):
if len(sys.argv) > 1:
json_files = sys.argv[1:]
else:
print(len(sys.argv))
json_files = glob('json/*')
for f in json_files:
print("Working on {}".format(f))
out = re.match('json/(.*).json',f).group(1)
data = read_json_file(f)
try:
data = read_json_file(f)
except:
print("WARNING: skipping unreadable json file {}".format(f))
continue
run_simulation_protocol( out, "job-id", data, gpu=0 )
......@@ -47,6 +47,7 @@ class Location():
self.type_ = type_
self.particle = None
self.connection = None
self.is_3prime_side_of_connection = None
self.prev_in_strand = None
self.next_in_strand = None
......@@ -59,8 +60,9 @@ class Location():
else:
return self.connection.other(self)
def set_connection(self,connection):
def set_connection(self, connection, is_3prime_side_of_connection):
self.connection = connection # TODO weakref?
self.is_3prime_side_of_connection = is_3prime_side_of_connection
def __repr__(self):
if self.on_fwd_strand:
......@@ -106,7 +108,27 @@ class ConnectableElement():
else:
counter[l] = 1
assert( np.all( [counter[l] == 1 for l in locs] ) )
return locs
return locs
def get_location_at(self, address, on_fwd_strand=True, new_type="crossover"):
loc = None
if (self.num_nts == 1):
# import pdb
# pdb.set_trace()
## Assumes that intrahelical connections have been made before crossovers
for l in self.locations:
if l.on_fwd_strand == on_fwd_strand and l.connection is None:
assert(loc is None)
loc = l
assert( loc is not None )
else:
for l in self.locations:
if l.address == address and l.on_fwd_strand == on_fwd_strand:
assert(loc is None)
loc = l
if loc is None:
loc = Location( self, address=address, type_=new_type, on_fwd_strand=on_fwd_strand )
return loc
def get_connections_and_locations(self, connection_type=None, exclude=[]):
""" Returns a list with each entry of the form:
......@@ -125,9 +147,13 @@ class ConnectableElement():
raise Exception("Object contains connection that fails to refer to object")
return ret
def _connect(self, other, connection):
def _connect(self, other, connection, in_3prime_direction=None):
## TODO fix circular references
A,B = [connection.A, connection.B]
if in_3prime_direction is not None:
A.is_3prime_side_of_connection = not in_3prime_direction
B.is_3prime_side_of_connection = in_3prime_direction
A.connection = B.connection = connection
self.connections.append(connection)
other.connections.append(connection)
......@@ -296,10 +322,19 @@ class Segment(ConnectableElement, Group):
for c,loc,other in self.get_connections_and_locations():
loc.particle = None
def contour_to_nt_pos(self,contour_pos):
return contour_pos*(self.num_nts-1)
def contour_to_nt_pos(self, contour_pos, round_nt=False):
nt = contour_pos*(self.num_nts-1)
if round_nt:
assert( (np.around(nt) - nt)**2 < 1e-3 )
nt = np.around(nt)
return nt
def nt_pos_to_contour(self,nt_pos):
return nt_pos/(self.num_nts-1)
if self.num_nts == 1:
assert(nt_pos == 0)
return 0
else:
return nt_pos/(self.num_nts-1)
def contour_to_position(self,s):
p = interpolate.splev( s, self.position_spline_params )
......@@ -439,25 +474,36 @@ class Segment(ConnectableElement, Group):
yield c
## TODO rename
def get_strand_segment(self, nt_pos, is_fwd):
def get_strand_segment(self, nt_pos, is_fwd, move_at_least=0.5):
""" Walks through locations, checking for crossovers """
# if self.name in ("6-1","1-1"):
# import pdb
# pdb.set_trace()
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))
# print(locations)
for l in locations:
pos = self.contour_to_nt_pos(l.address)
pos = self.contour_to_nt_pos(l.address, round_nt=True)
## DEBUG
## Skip locations encountered before our strand
tol = 0.1
if is_fwd:
if pos-nt_pos <= tol: continue
elif pos-nt_pos >= -tol: continue
# tol = 0.1
# if is_fwd:
# if pos-nt_pos <= tol: continue
# elif nt_pos-pos <= tol: continue
if (pos-nt_pos)*(2*is_fwd-1) < move_at_least: continue
## TODO: remove move_at_least
if np.isclose(pos,nt_pos):
if l.is_3prime_side_of_connection: continue
## Stop if we found the 3prime end
if l.on_fwd_strand == is_fwd and l.type_ == "3prime":
return pos, None, None, None
print(" found end at",l)
return pos, None, None, None, None
## Check location connections
c = l.connection
......@@ -466,16 +512,19 @@ class Segment(ConnectableElement, Group):
## 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))
Bpos = B.container.contour_to_nt_pos(B.address)
return pos, B.container, Bpos, B.on_fwd_strand
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)
return pos, B.container, Bpos, B.on_fwd_strand, 0.5
## Stop at other strand crossovers so basepairs line up
elif c.type_ == "crossover":
# print(" pausing at",l)
return pos, l.container, pos+(2*is_fwd-1), is_fwd
if nt_pos == pos: continue
print(" pausing at",l)
return pos, l.container, pos+(2*is_fwd-1), is_fwd, 0
import pdb
pdb.set_trace()
raise Exception("Shouldn't be here")
# print("Shouldn't be here")
## Made it to the end of the segment without finding a connection
......@@ -584,10 +633,13 @@ class Segment(ConnectableElement, Group):
if True:
print("WARNING: DEBUG")
## Remove duplicates, preserving order
tmp = []
for c in new_children:
if c not in tmp:
tmp.append(c)
else:
print(" duplicate particle found!")
new_children = tmp
for b in new_children:
......@@ -606,7 +658,7 @@ class Segment(ConnectableElement, Group):
def _generate_beads(self, bead_model, max_basepairs_per_bead, max_nucleotides_per_bead):
""" Generate beads (positions, types, etcl) and bonds, angles, dihedrals, exclusions """
""" Generate beads (positions, types, etc) and bonds, angles, dihedrals, exclusions """
## TODO: decide whether to remove bead_model argument
## (currently unused)
......@@ -616,7 +668,6 @@ class Segment(ConnectableElement, Group):
# existing_beads = [l.particle for l in locs if l.particle is not None]
existing_beads = {l.particle for l in self.locations if l.particle is not None}
existing_beads = sorted( list(existing_beads), key=lambda b: b.get_contour_position(self) )
if len(existing_beads) != len(set(existing_beads)):
pdb.set_trace()
......@@ -627,13 +678,13 @@ class Segment(ConnectableElement, Group):
## 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) > 1.5)
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) > 1.5)
assert(self.num_nts-1-existing_beads[0].get_nt_position(self) > 1.5)
# 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)
b = self._generate_one_bead(1, 0)
existing_beads.append(b)
assert(len(existing_beads) > 1)
......@@ -764,26 +815,23 @@ class DoubleStrandedSegment(Segment):
## Validate other, nt, other_nt
## TODO
## Create locations, connections and add to segments
c = self.nt_pos_to_contour(nt)
assert(c >= 0 and c <= 1)
def get_loc(seg, address, on_fwd_strand):
loc = None
for l in seg.locations:
if l.address == address and l.on_fwd_strand == on_fwd_strand:
assert(loc is None)
loc = l
if loc is None:
loc = Location( seg, address=address, type_="crossover", on_fwd_strand=on_fwd_strand )
return loc
if isinstance(other,SingleStrandedSegment):
other.add_crossover(other_nt, self, nt, strands_fwd[::-1])
else:
loc = get_loc(self, c, strands_fwd[0])
## Create locations, connections and add to segments
c = self.nt_pos_to_contour(nt)
assert(c >= 0 and c <= 1)
c = other.nt_pos_to_contour(other_nt)
assert(c >= 0 and c <= 1)
other_loc = get_loc( other, c, strands_fwd[1] )
self._connect(other, Connection( loc, other_loc, type_="crossover" ))
loc = self.get_location_at(c, strands_fwd[0])
c = other.nt_pos_to_contour(other_nt)
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]
## Real work
def _connect_ends(self, end1, end2, type_, force_connection):
......@@ -794,8 +842,10 @@ class DoubleStrandedSegment(Segment):
assert( end.type_ in ("end3","end5") )
assert( end1.type_ != end2.type_ )
## Create and add connection
end1.container._connect( end2.container, Connection( end1, end2, type_=type_ ) )
if end2.type_ == "end3":
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 )
def _get_num_beads(self, contour, max_basepairs_per_bead, max_nucleotides_per_bead):
return int(contour*self.num_nts // max_basepairs_per_bead)
......@@ -820,7 +870,6 @@ class DoubleStrandedSegment(Segment):
contour_position=contour_position )
self._add_bead(bead)
return bead
class SingleStrandedSegment(Segment):
......@@ -842,10 +891,10 @@ class SingleStrandedSegment(Segment):
for l in (self.start5,self.end3):
self.locations.append(l)
def connect_3end(self, end5, force_connection=False):
def connect_end3(self, end5, force_connection=False):
self._connect_end( end5, _5_to_3 = False, force_connection = force_connection )
def connect_5end(self, end3, force_connection=False):
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 )
def _connect_end(self, other, _5_to_3, force_connection):
......@@ -853,11 +902,47 @@ class SingleStrandedSegment(Segment):
if _5_to_3 == True:
my_end = self.end5
assert( other.type_ == "end3" )
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" )
conn = Connection( other, my_end, type_="intrahelical" )
other.container._connect( self, conn, in_3prime_direction=True )
self._connect( other.container, Connection( my_end, other, type_="intrahelical" ) )
def add_crossover(self, nt, other, other_nt, strands_fwd=[True,False]):
""" 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))
loc = self.get_location_at(c1, True)
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))
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":
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":
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 )
def _get_num_beads(self, contour, max_basepairs_per_bead, max_nucleotides_per_bead):
return int(contour*self.num_nts // max_nucleotides_per_bead)
......@@ -925,11 +1010,17 @@ class Strand(Group):
def add_dna(self, segment, start, end, is_fwd):
""" 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()
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) )
# 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):
import pdb
pdb.set_trace()
print(" CIRCULAR DNA")
s = StrandInSegment( segment, start, end, is_fwd )
self.add( s )
self.num_nts += s.num_nts
......@@ -1307,27 +1398,29 @@ class SegmentModel(ArbdModel):
s1,s2 = [l.container for l in (A,B)]
if A.particle is not None and B.particle is not None:
continue
assert( A.particle is None )
assert( B.particle is None )
# assert( A.particle is None )
# assert( B.particle is None )
## TODO: offload the work here to s1/s2 (?)
a1,a2 = [l.address for l in (A,B)]
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:
## 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 = s2.get_nearest_bead(a2)
if b is not None and s2.contour_to_nt_pos(np.abs(b.contour_position-a2)) < 1.9:
## 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)
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:
## 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)
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:
## 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)
""" Some tests """
for c,A,B in self.get_connections("intrahelical"):
......@@ -1758,27 +1851,27 @@ class SegmentModel(ArbdModel):
# if c[0]
""" Build strands from connectivity of helices """
def _recursively_build_strand(strand, segment, pos, is_fwd, mycounter=0):
def _recursively_build_strand(strand, segment, pos, is_fwd, mycounter=0, move_at_least=0.5):
mycounter+=1
if mycounter > 1000:
import pdb
pdb.set_trace()
s,seg = [strand, segment]
end_pos, next_seg, next_pos, next_dir = seg.get_strand_segment(pos, is_fwd)
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)
if next_seg is not None:
# print(" next_dir: {}".format(next_dir))
_recursively_build_strand(s, next_seg, next_pos, next_dir, mycounter)
_recursively_build_strand(s, next_seg, next_pos, next_dir, mycounter, move_at_least)
for seg in self.segments:
locs = seg.get_5prime_locations()
if locs is None: continue
# for pos, is_fwd in locs:
for l in locs:
# print("Tracing",l)
pos = seg.contour_to_nt_pos(l.address)
print("Tracing",l)
pos = seg.contour_to_nt_pos(l.address, round_nt=True)
is_fwd = l.on_fwd_strand
s = Strand()
_recursively_build_strand(s, seg, 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