Commit 34556ed5 authored by cmaffeo2's avatar cmaffeo2
Browse files

Improved psf/pdb writer

parent 5b0294b3
......@@ -382,19 +382,63 @@ class PointParticle(Transformable, Child):
else:
raise AttributeError("'{}' object has no attribute '{}'".format(type(self).__name__, name))
def _get_psfpdb_dictionary(self):
p = self
try:
segname = p.segname
except:
segname = "A"
try:
resname = p.resname
except:
resname = p.name[:3]
try:
resid = p.resid
except:
resid = p.idx+1
try:
occ = p.occupancy
except:
occ = 0
try:
beta = p.beta
except:
beta = 0
data = dict(segname = segname,
resname = resname,
name = str(p.name)[:4],
chain = "A",
resid = int(resid),
idx = p.idx+1,
type = p.type_.name[:4],
charge = p.charge,
mass = p.mass,
occupancy = occ,
beta = beta
)
return data
class Group(Transformable, Parent, Child):
def __init__(self, name=None, children = None, parent=None,
position = np.array((0,0,0)),
orientation = np.array(((1,0,0),(0,1,0),(0,0,1))),
remove_duplicate_bonded_terms = False,
):
**kwargs):
Transformable.__init__(self, position, orientation)
Child.__init__(self, parent) # Initialize Child first
Parent.__init__(self, children, remove_duplicate_bonded_terms)
self.name = name
self.isClone = False
for key,val in kwargs.items():
self.__dict__[key] = val
def clone(self):
return Clone(self)
g = copy(self)
......@@ -445,49 +489,32 @@ class PdbModel(Transformable, Parent):
fh.write("CRYST1{:>9.3f}{:>9.3f}{:>9.3f} 90.00 90.00 90.00 P 1 1\n".format( *self.dimensions ))
## Write coordinates
formatString = "ATOM {:>6.6s} {:^4.4s}{:1.1s}{:3.3s} {:1.1s}{:>5.5s} {:8.8s}{:8.8s}{:8.8s}{:6.2f}{:6.2f}{:2.2s}{:2d}{:>6s}\n"
formatString = "ATOM {idx:>6.6s} {name:^4.4s} {resname:3.3s} {chain:1.1s}{resid:>5.5s} {x:8.8s}{y:8.8s}{z:8.8s}{occupancy:6.2f}{beta:6.2f} {charge:2d}{segname:>6s}\n"
for p in self.particles:
## http://www.wwpdb.org/documentation/file-format-content/format33/sect9.html#ATOM
idx = p.idx+1
try:
segname = p.segname
except:
segname = "A"
try:
resname = p.resname
except:
resname = p.name[:3]
try:
resid = p.resid
except:
resid = idx
## TODOTODO test
data = p._get_psfpdb_dictionary()
idx = data['idx']
if np.log10(idx) >= 5:
idx = " *****"
else:
idx = "{:>6d}".format(idx)
data['idx'] = idx
# name = str(p.name)[:4]
name = str(p.name)[:4]
resname = resname
chain = "A"
charge = 0
occ = 0
beta = 0
# x,y,z = [x for x in p.collapsedPosition()]
pos = p.collapsedPosition()
dig = [max(int(np.log10(np.abs(x)+1e-6)//1),0)+1 for x in pos]
for d in dig: assert( d <= 7 )
# assert( np.all(dig <= 7) )
fs = ["{: %d.%df}" % (8,7-d) for d in dig]
x,y,z = [f.format(x) for f,x in zip(fs,pos)]
data['x'] = x
data['y'] = y
data['z'] = z
assert(data['resid'] < 1e5)
data['charge'] = int(data['charge'])
data['resid'] = "{:<4d}".format(data['resid'])
fh.write( formatString.format(**data) )
assert(resid < 1e5)
resid = "{:<4d}".format(resid)
fh.write( formatString.format(
idx, name, "", resname, chain, resid, x, y, z, occ, beta, "", charge, segname ))
return
def writePsf(self, filename):
......@@ -506,31 +533,9 @@ class PdbModel(Transformable, Parent):
formatString = "{idx:>8d} {segname:7.7s} {resid:<10.10s} {resname:7.7s}" + \
" {name:7.7s} {type:7.7s} {charge:f} {mass:f}\n"
for p in self.particles:
idx = p.idx + 1
try:
segname = p.segname
except:
segname = "A"
try:
resname = p.resname
except:
resname = str(p.name)[:3]
try:
resid = p.resid
except:
resid = idx
data = dict(
idx = idx,
segname = segname,
resid = "%d%c%c" % (resid," "," "), # TODO: work with large indices
name = str(p.name)[:4],
resname = resname,
type = p.type_.name[:4],
charge = p.charge,
mass = p.mass
)
fh.write(formatString.format( **data ))
data = p._get_psfpdb_dictionary()
data['resid'] = "%d%c%c" % (data['resid']," "," ") # TODO: work with large indices
fh.write( formatString.format(**data) )
fh.write("\n")
## Write out bonds
......
......@@ -285,26 +285,28 @@ class cadnano_part(SegmentModel):
zMid = int(0.5*(zid1+zid2))
assert( zMid in strandOccupancies[0] or zMid in strandOccupancies[1] )
numBps = zid2-zid1+1
for ins_idx,length in insertions:
## TODO: ensure placement of insertions is correct
## (e.g. are insertions at the ends handled correctly?)
if ins_idx < zid1:
continue
if ins_idx > zid2:
break
numBps += length
bp_to_zidx = []
insertion_dict = {idx:length for idx,length in insertions}
for i in range(zid1,zid2+1):
if i in insertion_dict:
l = insertion_dict[i]
else:
l = 0
for j in range(i,i+1+l):
bp_to_zidx.append(i)
numBps = len(bp_to_zidx)
# print("Adding helix with length",numBps,zid1,zid2)
kwargs = dict(name="%d-%d" % (hid,len(helices[hid])))
name = "%d-%d" % (hid,len(helices[hid]))
# "H%03d" % hid
kwargs = dict(name=name, segname=name, occupancy=hid)
posargs1 = dict( start_position = self._get_cadnano_position(hid,zid1),
end_position = self._get_cadnano_position(hid,zid2) )
posargs2 = dict( start_position = posargs1['end_position'],
end_position = posargs1['start_position'])
## TODO get sequence from cadnano api
if zMid in strandOccupancies[0] and zMid in strandOccupancies[1]:
kwargs['num_bp'] = numBps
......@@ -317,7 +319,28 @@ class cadnano_part(SegmentModel):
seg = SingleStrandedSegment(**kwargs,**posargs2)
else:
raise Exception("Segment could not be found")
def callback(segment, bp_to_zidx=bp_to_zidx):
for b in segment.beads:
bp = int(round(b.get_nt_position(segment)))
if bp < 0: bp = 0
if bp >= segment.num_nt: bp = segment.num_nt-1
b.beta = bp_to_zidx[bp]
if 'orientation_bead' in b.__dict__:
b.orientation_bead.beta = bp_to_zidx[bp]
seg._generate_bead_callbacks.append(callback)
def atomic_callback(nucleotide, bp_to_zidx=bp_to_zidx):
nt = nucleotide
segment = nucleotide.parent.segment
bp = int(round(segment.contour_to_nt_pos( nt.contour_position )))
if bp < 0: bp = 0
if bp >= segment.num_nt: bp = segment.num_nt-1
nt.beta = bp_to_zidx[bp]
nt.parent.occupancy = segment.occupancy
seg._generate_nucleotide_callbacks.append(atomic_callback)
helices[hid].append( seg )
helix_ranges[hid].append( (zid1,zid2) )
......
......@@ -196,10 +196,10 @@ class ConnectableElement():
# return [c for c in self.connections if c.A == loc or c.B == loc]
class SegmentParticle(PointParticle):
def __init__(self, type_, position, name="A", segname="A", **kwargs):
def __init__(self, type_, position, name="A", **kwargs):
self.name = name
self.contour_position = None
PointParticle.__init__(self, type_, position, name=name, segname=segname, **kwargs)
PointParticle.__init__(self, type_, position, name=name, **kwargs)
self.intrahelical_neighbors = []
self.other_neighbors = []
self.locations = []
......@@ -355,14 +355,17 @@ class Segment(ConnectableElement, Group):
def __init__(self, name, num_nt,
start_position = None,
end_position = None,
segment_model = None):
segment_model = None,
**kwargs):
if start_position is None: start_position = np.array((0,0,0))
Group.__init__(self, name, children=[])
Group.__init__(self, name, children=[], **kwargs)
ConnectableElement.__init__(self, connection_locations=[], connections=[])
self.resname = name
if 'segname' not in kwargs:
self.segname = name
# self.resname = name
self.start_orientation = None
self.twist_per_nt = 0
......@@ -381,6 +384,10 @@ class Segment(ConnectableElement, Group):
self.start_position = start_position
self.end_position = end_position
## Used to assign cadnano names to beads
self._generate_bead_callbacks = []
self._generate_nucleotide_callbacks = []
## Set up interpolation for positions
self._set_splines_from_ends()
......@@ -477,7 +484,7 @@ class Segment(ConnectableElement, Group):
def _generate_one_bead(self, contour_position, nts):
raise NotImplementedError
def _generate_atomic_nucleotide(self, contour_position, is_fwd, seq, scale):
def _generate_atomic_nucleotide(self, contour_position, is_fwd, seq, scale, strand_segment):
""" Seq should include modifications like 5T, T3 Tsinglet; direction matters too """
# print("Generating nucleotide at {}".format(contour_position))
......@@ -507,7 +514,6 @@ class Segment(ConnectableElement, Group):
nt_dict = canonicalNtFwd if is_fwd else canonicalNtRev
atoms = nt_dict[ key ].generate() # TODO: clone?
atoms.orientation = orientation.dot(atoms.orientation)
if isinstance(self, SingleStrandedSegment):
if scale is not None and scale != 1:
......@@ -526,6 +532,12 @@ class Segment(ConnectableElement, Group):
a.position = scale*(a.position-r0) + r0
a.beta = 0
atoms.position = pos
atoms.contour_position = contour_position
strand_segment.add(atoms)
for callback in self._generate_nucleotide_callbacks:
callback(atoms)
return atoms
......@@ -820,6 +832,9 @@ class Segment(ConnectableElement, Group):
# pdb.set_trace()
self._rebuild_children(tmp_children)
for callback in self._generate_bead_callbacks:
callback(self)
def _regenerate_beads(self, max_nts_per_bead=4, ):
...
......@@ -835,14 +850,16 @@ class DoubleStrandedSegment(Segment):
local_twist = False,
num_turns = None,
start_orientation = None,
twist_persistence_length = 90 ):
twist_persistence_length = 90,
**kwargs):
self.helical_rise = 10.44
self.distance_per_nt = 3.4
Segment.__init__(self, name, num_bp,
start_position,
end_position,
segment_model)
segment_model,
**kwargs)
self.num_bp = self.num_nt
self.local_twist = local_twist
......@@ -979,14 +996,16 @@ class SingleStrandedSegment(Segment):
def __init__(self, name, num_nt, start_position = None,
end_position = None,
segment_model = None):
segment_model = None,
**kwargs):
if start_position is None: start_position = np.array((0,0,0))
self.distance_per_nt = 5
Segment.__init__(self, name, num_nt,
start_position,
end_position,
segment_model)
segment_model,
**kwargs)
self.start = self.start5 = Location( self, address=0, type_= "end5" ) # TODO change type_?
self.end = self.end3 = Location( self, address=1, type_ = "end3" )
......@@ -1255,12 +1274,9 @@ class Strand(Group):
if strand_segment_count == len(s.strand_segments) and c == 1 and not self.is_circular:
seq = seq+"3"
nt = seg._generate_atomic_nucleotide( c, s.is_fwd, seq, scale )
# if s.is_fwd:
# else:
# nt = seg._generate_atomic_nucleotide( c, s.is_fwd, "A" )
nt = seg._generate_atomic_nucleotide( c, s.is_fwd, seq, scale, s )
s.add(nt)
## Join last basepairs
if last is not None:
self.link_nucleotides(last,nt)
......
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