Skip to content

GitLab

Commit 8fc5e3da authored by cmaffeo2's avatar cmaffeo2
Browse files

Added one_bea_per_monomer option to generate_bead_model

parent a5651d9d
Hide whitespace changes
Inline Side-by-side
mrdna/segmentmodel.py
View file @ 8fc5e3da
......@@ -1052,97 +1052,110 @@ class Segment(ConnectableElement, Group):
assert(len(old_beads) == len(self.beads))
def _generate_beads(self, bead_model, max_basepairs_per_bead, max_nucleotides_per_bead):
def _generate_beads(self, bead_model, max_basepairs_per_bead, max_nucleotides_per_bead, one_bead_per_monomer=False):
""" Generate beads (positions, types, etc) and bonds, angles, dihedrals, exclusions """
## TODO: decide whether to remove bead_model argument
## (currently unused)
## First find points between-which beads must be generated
# conn_locs = self.get_contour_sorted_connections_and_locations()
# locs = [A for c,A,B in conn_locs]
# existing_beads = [l.particle for l in locs if l.particle is not None]
# if self.name == "S001":
# pdb.set_trace()
tmp_children = []
# 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 one_bead_per_monomer:
last = None
for i in range(self.num_nt):
s = self.nt_pos_to_contour(i)
b = self._generate_one_bead(s,1)
# if self.num_nt == 1 and all([l.particle is not None for l in self.locations]):
# pdb.set_trace()
# return
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][0].get_nt_position(self,0) >= 0.5 and 'is_circular_bead' not in existing_beads[0][0].__dict__):
# if len(existing_beads) > 0:
# assert(existing_beads[0].get_nt_position(self) >= 0.5)
c = self.nt_pos_to_contour(0)
if self.num_nt == 1: c -= 0.4
b = self._generate_one_bead(c, 0)
existing_beads = [(b,0)] + existing_beads
if (existing_beads[-1][0].get_nt_position(self,1)-(self.num_nt-1) < -0.5 and 'is_circular_bead' not in existing_beads[0][0].__dict__) or len(existing_beads)==1:
c = self.nt_pos_to_contour(self.num_nt-1)
if self.num_nt == 1: c += 0.4
b = self._generate_one_bead(c, 0)
existing_beads.append( (b,1) )
assert(len(existing_beads) > 1)
## Walk through existing_beads, add beads between
tmp_children = [] # build list of children in nice order
last = None
if last is not None:
last.make_intrahelical_neighbor(b)
last = b
tmp_children.append(b)
else:
## First find points between-which beads must be generated
# conn_locs = self.get_contour_sorted_connections_and_locations()
# locs = [A for c,A,B in conn_locs]
# existing_beads = [l.particle for l in locs if l.particle is not None]
# if self.name == "S001":
# pdb.set_trace()
for I in range(len(existing_beads)-1):
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) )
# 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 np.isclose(eb1.position[2], eb2.position[2]):
# import pdb
# if self.num_nt == 1 and all([l.particle is not None for l in self.locations]):
# pdb.set_trace()
# return
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][0].get_nt_position(self,0) >= 0.5 and 'is_circular_bead' not in existing_beads[0][0].__dict__):
# if len(existing_beads) > 0:
# assert(existing_beads[0].get_nt_position(self) >= 0.5)
c = self.nt_pos_to_contour(0)
if self.num_nt == 1: c -= 0.4
b = self._generate_one_bead(c, 0)
existing_beads = [(b,0)] + existing_beads
if (existing_beads[-1][0].get_nt_position(self,1)-(self.num_nt-1) < -0.5 and 'is_circular_bead' not in existing_beads[0][0].__dict__) or len(existing_beads)==1:
c = self.nt_pos_to_contour(self.num_nt-1)
if self.num_nt == 1: c += 0.4
b = self._generate_one_bead(c, 0)
existing_beads.append( (b,1) )
assert(len(existing_beads) > 1)
## Walk through existing_beads, add beads between
last = None
for I in range(len(existing_beads)-1):
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 = 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
if num_beads == 0 and isinstance(self,SingleStrandedSegment) and isinstance(eb1.parent,DoubleStrandedSegment) and isinstance(eb2.parent,DoubleStrandedSegment):
num_beads = 1
## TODO similarly ensure there is a dsDNA bead between ssDNA beads
ds = e_ds / (num_beads+1)
nts = ds*self.num_nt
eb1.num_nt += 0.5*nts
eb2.num_nt += 0.5*nts
## Add beads
if eb1.parent == self:
tmp_children.append(eb1)
s0 = ec1
if last is not None:
last.make_intrahelical_neighbor(eb1)
last = eb1
for j in range(num_beads):
s = ds*(j+1) + s0
# if self.name in ("51-2","51-3"):
# if self.name in ("31-2",):
# print(" adding bead at {}".format(s))
b = self._generate_one_bead(s,nts)
last.make_intrahelical_neighbor(b)
last = b
tmp_children.append(b)
# print(" %s working on %d to %d" % (self.name, eb1.position[2], eb2.position[2]))
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
if num_beads == 0 and isinstance(self,SingleStrandedSegment) and isinstance(eb1.parent,DoubleStrandedSegment) and isinstance(eb2.parent,DoubleStrandedSegment):
num_beads = 1
## TODO similarly ensure there is a dsDNA bead between ssDNA beads
ds = e_ds / (num_beads+1)
nts = ds*self.num_nt
eb1.num_nt += 0.5*nts
eb2.num_nt += 0.5*nts
## Add beads
if eb1.parent == self:
tmp_children.append(eb1)
last.make_intrahelical_neighbor(eb2)
s0 = ec1
if last is not None:
last.make_intrahelical_neighbor(eb1)
last = eb1
for j in range(num_beads):
s = ds*(j+1) + s0
# if self.name in ("51-2","51-3"):
# if self.name in ("31-2",):
# print(" adding bead at {}".format(s))
b = self._generate_one_bead(s,nts)
last.make_intrahelical_neighbor(b)
last = b
tmp_children.append(b)
last.make_intrahelical_neighbor(eb2)
if eb2.parent == self:
tmp_children.append(eb2)
if eb2.parent == self:
tmp_children.append(eb2)
# if self.name in ("31-2",):
# pdb.set_trace()
self._rebuild_children(tmp_children)
......@@ -2249,10 +2262,17 @@ class SegmentModel(ArbdModel):
escapable_twist=escapable_twist)
def generate_bead_model(self,
max_basepairs_per_bead = 7,
max_nucleotides_per_bead = 4,
local_twist=False,
escapable_twist=True):
max_basepairs_per_bead = 7,
max_nucleotides_per_bead = 4,
local_twist=False,
escapable_twist=True,
sequence_dependent = False,
one_bead_per_monomer = False
):
if sequence_dependent and not local_twist:
raise ValueError("Incompatible options: if sequence_dependent==True, then local_twist must be True")
if sequence_dependent and not one_bead_per_monomer:
raise ValueError("Incompatible options: if sequence_dependent==True, then one_bead_per_monomer must be True")
self.children = self.segments # is this okay?
self.clear_beads()
......@@ -2291,97 +2311,120 @@ class SegmentModel(ArbdModel):
""" Generate beads at intrahelical junctions """
if self.DEBUG: print( "Adding intrahelical beads at junctions" )
if not one_bead_per_monomer:
## 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)]
## 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)]
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
# TODOTODO
a1,a2 = [l.address for l in (A,B)]
for a in (a1,a2):
assert( np.isclose(a,0) or np.isclose(a,1) )
## TODO: offload the work here to s1
# TODOTODO
a1,a2 = [l.address for l in (A,B)]
for a in (a1,a2):
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)
def _get_nearest(seg,address):
b = seg.get_nearest_bead(address)
if b is not None:
assert( b.parent is seg )
""" if above assertion is true, no problem here """
if np.abs(b.get_nt_position(seg) - seg.contour_to_nt_pos(address)) > 0.5:
b = None
return b
""" Search to see whether bead at location is already found """
b = None
if isinstance(s1,DoubleStrandedSegment):
b = _get_nearest(s1,a1)
if b is None and isinstance(s2,DoubleStrandedSegment):
b = _get_nearest(s2,a2)
if b is not None and b.parent not in (s1,s2):
b = None
## TODO improve this for combinations of ssDNA and dsDNA (maybe a1/a2 should be calculated differently)
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.is_intrahelical = True
b.locations.extend([A,B])
if s1 is s2:
b.is_circular_bead = True
def _get_nearest(seg,address):
b = seg.get_nearest_bead(address)
if b is not None:
assert( b.parent is seg )
""" if above assertion is true, no problem here """
if np.abs(b.get_nt_position(seg) - seg.contour_to_nt_pos(address)) > 0.5:
b = None
return b
# pdb.set_trace()
""" Search to see whether bead at location is already found """
b = None
if isinstance(s1,DoubleStrandedSegment):
b = _get_nearest(s1,a1)
if b is None and isinstance(s2,DoubleStrandedSegment):
b = _get_nearest(s2,a2)
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.is_intrahelical = True
b.locations.extend([A,B])
if s1 is s2:
b.is_circular_bead = True
""" Generate beads at other junctions """
for c,A,B in self.get_connections(exclude="intrahelical"):
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 )
# pdb.set_trace()
## TODO: offload the work here to s1/s2 (?)
a1,a2 = [l.address for l in (A,B)]
""" Generate beads at other junctions """
for c,A,B in self.get_connections(exclude="intrahelical"):
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 )
## TODO: offload the work here to s1/s2 (?)
a1,a2 = [l.address for l in (A,B)]
def maybe_add_bead(location, seg, address, ):
if location.particle is None:
b = seg.get_nearest_bead(address)
try:
distance = seg.contour_to_nt_pos(np.abs(b.contour_position-address))
max_distance = min(max_basepairs_per_bead, max_nucleotides_per_bead)*0.5
if "is_intrahelical" in b.__dict__:
max_distance = 0.5
if distance >= max_distance:
raise Exception("except")
## combine beads
b.update_position( 0.5*(b.contour_position + address) ) # avg position
except:
b = seg._generate_one_bead(address,0)
location.particle = b
b.locations.append(location)
maybe_add_bead(A,s1,a1)
maybe_add_bead(B,s2,a2)
""" Some tests """
for c,A,B in self.get_connections("intrahelical"):
for l in (A,B):
if l.particle is None: continue
assert( l.particle.parent is not None )
def maybe_add_bead(location, seg, address, ):
if location.particle is None:
b = seg.get_nearest_bead(address)
try:
distance = seg.contour_to_nt_pos(np.abs(b.contour_position-address))
max_distance = min(max_basepairs_per_bead, max_nucleotides_per_bead)*0.5
if "is_intrahelical" in b.__dict__:
max_distance = 0.5
if distance >= max_distance:
raise Exception("except")
## combine beads
b.update_position( 0.5*(b.contour_position + address) ) # avg position
except:
b = seg._generate_one_bead(address,0)
location.particle = b
b.locations.append(location)
maybe_add_bead(A,s1,a1)
maybe_add_bead(B,s2,a2)
""" Some tests """
for c,A,B in self.get_connections("intrahelical"):
for l in (A,B):
if l.particle is None: continue
assert( l.particle.parent is not None )
""" Generate beads in between """
if self.DEBUG: print("Generating beads")
for s in segments:
s._generate_beads( self, max_basepairs_per_bead, max_nucleotides_per_bead )
""" Generate beads in between """
if self.DEBUG: print("Generating beads")
for s in segments:
s._generate_beads( self, max_basepairs_per_bead, max_nucleotides_per_bead )
else:
""" Generate beads in each segment """
if self.DEBUG: print("Generating beads")
for s in segments:
s._generate_beads( self, max_basepairs_per_bead, max_nucleotides_per_bead, one_bead_per_monomer )
""" Associate beads with junctions """
# for c,A,B in self.get_connections("intrahelical"):
for c,A,B in self.get_connections():
for l in (A,B):
seg = l.container
a = l.address
b = seg.get_nearest_bead(a)
l.particle = b
b.locations.append(l)
"""
b.is_intrahelical = True
b.locations.extend([A,B])
if s1 is s2:
b.is_circular_bead = True
"""
# """ Combine beads at junctions as needed """
# for c,A,B in self.get_connections():
......
Markdown is supported
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