From ee1e32647cae9fc96d57ebcc8a9df388f8e860f0 Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Mon, 24 Feb 2020 12:14:40 -0600
Subject: [PATCH] Various changes to make it easy to run circular structures
---
mrdna/segmentmodel.py | 67 +++++++++++++++++++++++++++----------------
mrdna/simulate.py | 2 ++
2 files changed, 45 insertions(+), 24 deletions(-)
diff --git a/mrdna/segmentmodel.py b/mrdna/segmentmodel.py
index 6a8c6a6..33db573 100644
--- a/mrdna/segmentmodel.py
+++ b/mrdna/segmentmodel.py
@@ -199,8 +199,7 @@ class ConnectableElement():
self.connections.append(connection)
if other is not self:
other.connections.append(connection)
- else:
- raise NotImplementedError("Segments cannot yet be connected to themselves; if you are attempting to make a circular object, try breaking the object into multiple segments")
+
l = A.container.locations
if A not in l: l.append(A)
l = B.container.locations
@@ -977,7 +976,6 @@ class Segment(ConnectableElement, Group):
def _generate_beads(self, bead_model, max_basepairs_per_bead, max_nucleotides_per_bead):
-
""" Generate beads (positions, types, etc) and bonds, angles, dihedrals, exclusions """
## TODO: decide whether to remove bead_model argument
## (currently unused)
@@ -1003,7 +1001,7 @@ class Segment(ConnectableElement, Group):
## 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:
+ 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)
@@ -1011,7 +1009,7 @@ class Segment(ConnectableElement, Group):
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 or len(existing_beads)==1:
+ 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)
@@ -2141,6 +2139,7 @@ class SegmentModel(ArbdModel):
""" Generate beads at intrahelical junctions """
if self.DEBUG: print( "Adding intrahelical beads at junctions" )
+
## 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)]
@@ -2155,21 +2154,22 @@ class SegmentModel(ArbdModel):
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)
- """ Search to see whether bead at location is already found """
- b = None
- if isinstance(s1,DoubleStrandedSegment):
- b = s1.get_nearest_bead(a1)
+
+ def _get_nearest(seg,address):
+ b = seg.get_nearest_bead(address)
if b is not None:
- assert( b.parent is s1 )
+ assert( b.parent is seg )
""" if above assertion is true, no problem here """
- if np.abs(b.get_nt_position(s1) - s1.contour_to_nt_pos(a1)) > 0.5:
- b = None
- if b is None and 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:
+ 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
@@ -2183,6 +2183,8 @@ class SegmentModel(ArbdModel):
A.particle = B.particle = b
b.is_intrahelical = True
b.locations.extend([A,B])
+ if s1 is s2:
+ b.is_circular_bead = True
# pdb.set_trace()
@@ -2242,12 +2244,27 @@ class SegmentModel(ArbdModel):
# pdb.set_trace()
""" Add intrahelical neighbors at connections """
-
+ special_seps = dict()
for c,A,B in self.get_connections("intrahelical"):
b1,b2 = [l.particle for l in (A,B)]
if b1 is b2:
- ## already handled by Segment._generate_beads
- pass
+ ## already handled by Segment._generate_beads, unless A,B are on same segment
+ if A.container == B.container:
+ sorted_sites = sorted([(abs(L.address-b1.contour_position),L)
+ for L in (A,B)])
+ close_site, far_site = [s[1] for s in sorted_sites]
+
+ b3 = far_site.container.get_nearest_bead(far_site.address)
+ b1.intrahelical_neighbors.append(b3)
+ b3.intrahelical_neighbors.append(b1)
+
+ if far_site.address > 0.5:
+ sep = b1.contour_position + (1-b3.contour_position)
+ else:
+ sep = b3.contour_position + (1-b1.contour_position)
+ sep = A.container.num_nt * sep
+ special_seps[(b1,b3)] = special_seps[(b3,b1)] = sep
+
else:
for b in (b1,b2): assert( b is not None )
b1.make_intrahelical_neighbor(b2)
@@ -2358,10 +2375,12 @@ class SegmentModel(ArbdModel):
parent = self._getParent(b1,b2)
- seg = b2.parent
- c0 = b2.contour_position
- sep = np.abs(b1.get_nt_position(seg,c0)-b2.get_nt_position(seg))
-
+ if (b1,b2) in special_seps:
+ sep = special_seps[(b1,b2)]
+ else:
+ seg = b2.parent
+ c0 = b2.contour_position
+ sep = np.abs(b1.get_nt_position(seg,c0)-b2.get_nt_position(seg))
is_dsdna = b1.type_.name[0] == "D" and b2.type_.name[0] == "D"
if is_dsdna:
@@ -2426,7 +2445,7 @@ class SegmentModel(ArbdModel):
parent = self._getParent(b1,b2,b3)
seg = b2.parent
c0 = b2.contour_position
- sep = np.abs(b1.get_nt_position(seg,c0)-b2.get_nt_position(seg)) + np.abs(b3.get_nt_position(seg,c0)-b2.get_nt_position(seg))
+ sep = dists[b2][b1] + dists[b2][b3]
if b1.type_.name[0] == "D" and b2.type_.name[0] == "D" and b3.type_.name[0] == "D":
diff --git a/mrdna/simulate.py b/mrdna/simulate.py
index b685098..1d040ac 100644
--- a/mrdna/simulate.py
+++ b/mrdna/simulate.py
@@ -10,6 +10,8 @@ import numpy as np
def _remove_bonds_beyond(model, cutoff):
bonds = model.get_bonds()
+ for seg in model.segments:
+ bonds.extend(seg.get_bonds())
new_bonds = []
for bond in bonds:
n1,n2,b,ex = bond
--
GitLab