import numpy as np
from arbdmodel import PointParticle, ParticleType, Group, ArbdModel
from nonbonded import *
from copy import copy, deepcopy
from nbPot import nbDnaScheme
"""
# TODO:
"""
class Location():
""" Site for connection within an object """
def __init__(self, container, address, type_):
self.container = container
self.address = address
self.type_ = type_
self.particle = None
class Connection():
""" Abstract base class for connection between two elements """
def __init__(self, A, B, type_ = None):
assert( isinstance(A,Location) )
assert( isinstance(B,Location) )
self.A = A
self.B = B
self.type_ = type_
# class ConnectableElement(Transformable):
class ConnectableElement():
""" Abstract base class """
def __init__(self, connections=[]):
self.connections = connections
def _connect(self, other, connection):
self.connections.append(connection)
other.connections.append(connection)
def _find_connections(self, loc):
return [c for c in self.connections if c.A == loc or c.B == loc]
class Segment(ConnectableElement, Group):
""" Base class that describes a segment of DNA. When built from
cadnano models, should not span helices """
"""Define basic particle types"""
dsDNA_particle = ParticleType("D",
diffusivity = 43.5,
mass = 300,
radius = 3,
)
ssDNA_particle = ParticleType("S",
diffusivity = 43.5,
mass = 150,
radius = 3,
)
def __init__(self, name, num_nts,
start_position = np.array((0,0,0)),
end_position = None,
segment_model = None):
Group.__init__(self, name, children=[])
ConnectableElement.__init__(self, connections=[])
self._bead_model_generation = 0 # TODO: remove?
self.segment_model = segment_model # TODO: remove?
# self.end5 = Location( self, address=0, type_= "end5" )
# self.end3 = Location( self, address=-1, type_ = "end3" )
self.num_nts = num_nts
if end_position is None:
end_position = np.array((0,0,self.distance_per_nt*num_nts)) + start_position
self.start_position = start_position
self.end_position = end_position
def _generate_one_bead(self, pos, nts):
raise NotImplementedError
def _assign_particles_to_locations(self):
raise NotImplementedError
def get_all_consecutive_beads(self, number):
assert(number >= 1)
## Assume that consecutive beads in self.children are bonded
ret = []
for i in range(len(self.children)-number+1):
tmp = [self.children[i+j] for j in range(0,number)]
ret.append( tmp )
return ret
def get_beads_before_bead(self, bead, number, inclusive=True):
## Assume that consecutive beads in self.children are bonded
i = self.children.index(bead)
l = len(self.children)
if i-number < 0:
raise Exception("Not enough beads after bead")
start = 1
if inclusive: start = 0
return [self.children[i-j] for j in range(start,number)]
def get_beads_after_bead(self, bead, number, inclusive=True):
## Assume that consecutive beads in self.children are bonded
i = self.children.index(bead)
l = len(self.children)
if i+number >= l:
raise Exception("Not enough beads after bead")
start = 1
if inclusive: start = 0
return [self.children[i+i] for j in range(start,number)]
def _generate_beads(self, bead_model, max_nts_per_bead=4):
""" Generate beads (positions, types, etcl) and bonds, angles, dihedrals, exclusions """
# self._bead_model_generation += 1
# self._bead_model_max_nts_per_bead = max_nts_per_bead
direction = self.end_position - self.start_position
num_beads = (self.num_nts // max_nts_per_bead) + 1
nts_per_bead = float(self.num_nts)/num_beads
last = None
for i in range(num_beads+1):
nts = nts_per_bead
if i == 0 or i == num_beads:
nts *= 0.5
s = i*float(nts_per_bead)/(self.num_nts) # contour
pos = direction * s + self.start_position
b = self._generate_one_bead(pos,nts)
self.children.append(b)
# if last is not None:
# self.add_bond( i=last, j=b, bond="ssdna" )
# last = b
self._assign_particles_to_locations()
def _regenerate_beads(self, max_nts_per_bead=4, ):
...
def _generate_atomic(self, atomic_model):
...
class DoubleStrandedSegment(Segment):
""" Class that describes a segment of ssDNA. When built from
cadnano models, should not span helices """
def __init__(self, name, num_nts, start_position = np.array((0,0,0)),
end_position = None,
segment_model = None,
twist = None,
start_orientation = None):
self.distance_per_nt = 5
Segment.__init__(self, name, num_nts,
start_position,
end_position,
segment_model)
self.nicks = []
self.start5 = Location( self, address=0, type_= "end5" )
self.start3 = Location( self, address=0, type_ = "end3" )
self.end5 = Location( self, address=-1, type_= "end5" )
self.end3 = Location( self, address=-1, type_ = "end3" )
## Convenience methods
def connect_start5(self, end3, force_connection=False):
if isinstance(end3, SingleStrandedSegment):
end3 = end3.end3
self._connect_ends( self.start5, end3, force_connection = force_connection )
def connect_start3(self, end5, force_connection=False):
if isinstance(end5, SingleStrandedSegment):
end5 = end5.end5
self._connect_ends( self.start3, end5, force_connection = force_connection )
def connect_end3(self, end5, force_connection=False):
if isinstance(end5, SingleStrandedSegment):
end5 = end5.end5
self._connect_ends( self.end3, end5, force_connection = force_connection )
def connect_end5(self, end3, force_connection=False):
if isinstance(end3, SingleStrandedSegment):
end3 = end3.end3
self._connect_ends( self.end5, end3, force_connection = force_connection )
## Real work
def _connect_ends(self, end1, end2, force_connection):
## validate the input
for end in (end1, end2):
assert( isinstance(end, Location) )
assert( end.type_ in ("end3","end5") )
assert( end1.type_ != end2.type_ )
end1.container._connect( end2.container, Connection( end1, end2, type_="intrahelical" ) )
def _generate_one_bead(self, pos, nts):
return PointParticle( Segment.dsDNA_particle, pos, nts,
num_nts=nts, parent=self )
def _assign_particles_to_locations(self):
self.start3.particle = self.start5.particle = self.children[0]
self.end3.particle = self.end5.particle = self.children[-1]
def _generate_atomic(self, atomic_model):
...
# def add_crossover(self, locationInA, B, locationInB):
# j = Crossover( [self, B], [locationInA, locationInB] )
# self._join(B,j)
# def add_internal_crossover(self, locationInA, B, locationInB):
# j = Crossover( [self, B], [locationInA, locationInB] )
# self._join(B,j)
# def stack_end(self, myEnd):
# ## Perhaps this should not really be possible; these ends should be part of same helix
# ...
# def connect_strand(self, other):
# ...
# def break_apart(self):
# """Break into smaller pieces so that "crossovers" are only at the ends"""
# ...
class SingleStrandedSegment(Segment):
""" Class that describes a segment of ssDNA. When built from
cadnano models, should not span helices """
def __init__(self, name, num_nts, start_position = np.array((0,0,0)),
end_position = None,
segment_model = None):
self.distance_per_nt = 5
Segment.__init__(self, name, num_nts,
start_position,
end_position,
segment_model)
self.start = self.end5 = Location( self, address=0, type_= "end5" )
self.end = self.end3 = Location( self, address=-1, type_ = "end3" )
def connect_3end(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):
self._connect_end( end3, _5_to_3 = True, force_connection = force_connection )
def _connect_end(self, other, _5_to_3, force_connection):
assert( isinstance(other, Location) )
if _5_to_3 == True:
my_end = self.end5
assert( other.type_ == "end3" )
else:
my_end = self.end3
assert( other.type_ == "end5" )
self._connect( other.container, Connection( my_end, other, type_="intrahelical" ) )
def _generate_one_bead(self, pos, nts):
return PointParticle( Segment.ssDNA_particle, pos, nts,
num_nts=nts, parent=self )
def _assign_particles_to_locations(self):
self.start.particle = self.children[0]
self.end.particle = self.children[-1]
def _generate_atomic(self, atomic_model):
...
class SegmentModel(ArbdModel):
def __init__(self, segments = [],
max_basepairs_per_bead = 7,
max_nucleotides_per_bead = 4,
dimensions=(1000,1000,1000), temperature=291,
timestep=50e-6, cutoff=50,
decompPeriod=10000, pairlistDistance=None,
nonbondedResolution=0):
ArbdModel.__init__(self,segments,
dimensions, temperature, timestep, cutoff,
decompPeriod, pairlistDistance=None,
nonbondedResolution=0)
self.max_basepairs_per_bead = max_basepairs_per_bead # dsDNA
self.max_nucleotides_per_bead = max_nucleotides_per_bead # ssDNA
self.children = self.segments = segments
self._bonded_potential = dict() # cache bonded potentials
self._generate_bead_model(segments, max_nucleotides_per_bead, max_nucleotides_per_bead)
def _get_intrahelical_beads(self):
ret = []
for s in self.segments:
ret.extend( s.get_all_consecutive_beads(2) )
for s in self.segments:
for c in s.connections:
if c.type_ == "intrahelical":
if c.A.container == s: # avoid double-counting
b1,b2 = [loc.particle for loc in (c.A,c.B)]
for b in (b1,b2): assert( b is not None )
ret.append( [b1,b2] )
return ret
def _get_intrahelical_angle_beads(self):
ret = []
for s in self.segments:
ret.extend( s.get_all_consecutive_beads(3) )
for s1 in self.segments:
for c in s1.connections:
if c.A.container != s1: continue
s2 = c.B.container
if c.type_ == "intrahelical":
b1,b2 = [loc.particle for loc in (c.A,c.B)]
for b in (b1,b2): assert( b is not None )
try:
b0 = s1.get_beads_before_bead(b1,1)[0]
assert( b0 is not None )
ret.append( [b0,b1,b2] )
except:
...
try:
b0 = s1.get_beads_after_bead(b1,1)[0]
assert( b0 is not None )
ret.append( [b2,b1,b0] )
except:
...
try:
b3 = s2.get_beads_before_bead(b2,1)[0]
assert( b3 is not None )
ret.append( [b3,b2,b1] )
except:
...
try:
b3 = s2.get_beads_after_bead(b2,1)[0]
assert( b3 is not None )
ret.append( [b1,b2,b3] )
except:
...
return ret
def _get_potential(self, type_, kSpring, d):
key = (type_,kSpring,d)
if key not in self._bonded_potential:
if type_ == "bond":
self._bonded_potential[key] = HarmonicBond(kSpring,d)
elif type_ == "angle":
self._bonded_potential[key] = HarmonicAngle(kSpring,d)
elif type_ == "dihedral":
self._bonded_potential[key] = HarmonicDihedral(kSpring,d)
else:
raise Exception("Unhandled potential type '%s'" % type_)
return self._bonded_potential[key]
def get_bond_potential(self, kSpring, d):
return self._get_potential("bond", kSpring, d)
def get_angle_potential(self, kSpring, d):
return self._get_potential("angle", kSpring, d)
def get_dihedral_potential(self, kSpring, d):
return self._get_potential("dihedral", kSpring, d)
def _generate_bead_model(self, segments,
max_basepairs_per_bead = 7,
max_nucleotides_per_bead = 4):
""" Generate beads """
for s in segments:
s._generate_beads( max_nucleotides_per_bead )
""" Combine beads at junctions as needed """
for s in segments:
for c in s.connections:
if c.A.container == s:
...
""" Reassign bead types """
beadtype_s = dict()
for segment in segments:
for b in segment:
key = (b.type_.name[0].upper(), b.num_nts)
if key in beadtype_s:
b.type_ = beadtype_s[key]
else:
t = deepcopy(b.type_)
if key[0] == "D":
t.__dict__["nts"] = b.num_nts*2
elif key[0] == "S":
t.__dict__["nts"] = b.num_nts
else:
raise Exception("TODO")
print(t.nts)
t.name = t.name + "%03d" % (100*t.nts)
beadtype_s[key] = b.type_ = t
""" Add intrahelical potentials """
## First replace intrahelical bonds
for b1,b2 in self._get_intrahelical_beads():
if b1.parent == b2.parent:
sep = 0.5*(b1.num_nts+b2.num_nts)
parent = b1.parent
else:
sep = 1
parent = self
if b1.type_.name[0] == "D" and b1.type_.name[0] == "D":
k = 10.0/np.sqrt(sep) # TODO: determine from simulations
d = 3.4*sep
else:
## TODO: get correct numbers from ssDNA model
k = 1.0/np.sqrt(sep)
d = 5*sep
bond = self.get_bond_potential(k,d)
parent.add_bond( b1, b2, bond, exclude=True )
""" Add connection potentials """
## TODO
# def get_bead(self, location):
# if type(location.container) is not list:
# s = self.segments.index(location.container)
# s.get_bead(location.address)
# else:
# r
# ...
if __name__ == "__main__":
seg1 = DoubleStrandedSegment("strand", num_nts = 46)
seg2 = SingleStrandedSegment("strand",
start_position = seg1.end_position + np.array((1,0,1)),
num_nts = 12)
seg3 = SingleStrandedSegment("strand",
start_position = seg1.start_position + np.array((-1,0,-1)),
end_position = seg1.end_position + np.array((-1,0,1)),
num_nts = 128)
seg1.start3
seg1.start5
seg1.end3
seg1.end5
seg1.connect_end3(seg2)
seg1.connect_end5(seg3)
seg1.connect_start3(seg3)
model = SegmentModel( [seg1, seg2, seg3],
dimensions=(5000,5000,5000),
)
model.useNonbondedScheme( nbDnaScheme )
model.simulate( outputPrefix = 'strand-test', outputPeriod=1e4, numSteps=1e6, gpu=1 )
# seg = SingleStrandedSegment("strand", num_nts = 21)
# generate_bead_model( [seg] )
# for b in seg:
# print(b.num_nts, b.position)