Commit f87c6dfe authored by cmaffeo2's avatar cmaffeo2
Browse files

Updated arbdmodel writePdb and nucleotide generation

parent 35e5c535
......@@ -20,6 +20,68 @@ def stringToIntTuples(string, tupleLen, offset):
return ret
class CanonicalNucleotideFactory(Group):
DefaultOrientation = rotationAboutAxis([0,0,1], 80)
def __init__(self, prefix, seq):
self.sequence = seq # TODO: used?
self.resname = resnames[seq]
self.children = []
## idx, type,
d = np.loadtxt("%s.dat" % prefix, dtype='i4,S4,S4,f4,f4,f4,f4,f4')
ids = d['f0']
assert( np.all(ids == sorted(ids)) )
names = d['f1']
tns = d['f2']
xs = d['f5']
ys = d['f6']
zs = d['f7']
type_list = []
for i in range(len(tns)):
tn = tns[i]
if tn not in atomic_types:
atomic_types[tn] = ParticleType(tn.decode('UTF-8'),
charge = d['f3'][i],
mass = d['f4'][i])
p = PointParticle( type_ = atomic_types[tn],
position = [xs[i],ys[i],zs[i]],
name = names[i].decode('UTF-8') )
self.children.append( p )
atoms = self.children
# minIdx = np.min(self.indices) - 1
self.bonds = np.loadtxt("%s.bonds.txt" % prefix, dtype='i4') # - minIdx
self.angles = np.loadtxt("%s.angles.txt" % prefix, dtype='i4') # - minIdx
self.dihedrals = np.loadtxt("%s.dihedrals.txt" % prefix, dtype='i4')# - minIdx
self.impropers = np.loadtxt("%s.impropers.txt" % prefix, dtype='i4')# - minIdx
def generate(self):
nt = Group(orientation = CanonicalNucleotide.DefaultOrientation)
nt.__dict__['sequence'] = self.sequence
nt.__dict__['resname'] = self.resname
atoms = nt.children = [copy.copy(c) for c in self.children]
for a in atoms:
a.parent = nt
nt.atoms_by_name = {a.name:a for a in nt.children}
b = self.bonds
a = self.angles
d = self.dihedrals
i = self.impropers
for i,j in zip(b[::2],b[1::2]):
nt.add_bond( atoms[i], atoms[j], None )
for i,j,k in zip(a[::3],a[1::3],a[2::3]):
nt.add_angle( atoms[i], atoms[j], atoms[k], None )
for i,j,k,l in zip(d[::4],d[1::4],d[2::4],d[3::4]):
nt.add_dihedral( atoms[i], atoms[j], atoms[k], atoms[l], None )
for i,j,k,l in zip(d[::4],d[1::4],d[2::4],d[3::4]):
nt.add_improper( atoms[i], atoms[j], atoms[k], atoms[l], None )
return nt
class CanonicalNucleotide(Group):
DefaultOrientation = rotationAboutAxis([0,0,1], 80)
def __init__(self, prefix, seq):
......@@ -128,10 +190,7 @@ for seq in seqComplement.keys():
for keystring,suff in zip(("5%s","%s","%s3","5%s3"),("-5prime","","-3prime","-singlet")):
prefix = 'resources/%s-%s%s' % (seq,direction,suff)
key = keystring % seq
canonicalNtFwd[key] = CanonicalNucleotide( prefix, seq )
# canonicalNtFwd[seq] = CanonicalNucleotide( prefix )
# canonicalNtFwd3[seq] = CanonicalNucleotide( prefix+"3prime" )
# canonicalNtFwd5[seq] = CanonicalNucleotide( prefix+"5prime" )
canonicalNtFwd[key] = CanonicalNucleotideFactory( prefix, seq )
canonicalNtRev = dict()
direction = 'rev'
......@@ -139,13 +198,7 @@ for seq in seqComplement.keys():
for keystring,suff in zip(("5%s","%s","%s3","5%s3"),("-5prime","","-3prime","-singlet")):
prefix = 'resources/%s-%s%s' % (seq,direction,suff)
key = keystring % seq
canonicalNtRev[key] = CanonicalNucleotide( prefix, seq )
# for seq in seqComplement.keys():
# prefix = 'resources/%s_%s' % (seq,direction)
# canonicalNtRev[seq] = CanonicalNucleotide( prefix )
# canonicalNtRev3[seq] = CanonicalNucleotide( prefix+"3prime" )
# canonicalNtRev5[seq] = CanonicalNucleotide( prefix+"5prime" )
canonicalNtRev[key] = CanonicalNucleotideFactory( prefix, seq )
with open("resources/enm-template-honeycomb.json") as ch:
enmTemplateHC = json.load(ch)
......
......@@ -382,7 +382,7 @@ class Group(Transformable, Parent, Child):
# orientation = self.orientation)
# g.children = deepcopy self.children.deepcopy() # lists are the same object
## TODO override deepcopy so parent can be excluded from copying?
# def __getstate__(self):
# return (self.children, self.parent, self.position, self.orientation)
......@@ -411,7 +411,8 @@ class PdbModel(Transformable, Parent):
fh.write("CRYST1 {:>5f} {:>5f} {:>5f} 90.00 90.00 90.00 P 1 1\n".format( *self.dimensions ))
## Write coordinates
formatString = "ATOM {:>5d} {:^4.4s}{:1.1s}{:3.3s} {:1.1s}{:>5.5s} {:8.3f}{:8.3f}{:8.3f}{:6.2f}{:6.2f}{:2.2s}{:2f}\n"
# formatString = "ATOM {:>5d} {:^4.4s}{:1.1s}{:3.3s} {:1.1s}{:>5.5s} {:8.3f}{:8.3f}{:8.3f}{:6.2f}{:6.2f}{:2.2s}{:2f}\n"
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}{:2f}\n"
for p in self.particles:
## http://www.wwpdb.org/documentation/file-format-content/format33/sect9.html#ATOM
idx = p.idx+1
......@@ -424,6 +425,12 @@ class PdbModel(Transformable, Parent):
except:
resid = idx
## TODOTODO test
if np.log10(idx) >= 5:
idx = "******"
else:
idx = "{:>6d}".format(idx)
# name = str(p.name)[:4]
name = str(p.name)[:4]
resname = resname
......@@ -431,8 +438,14 @@ class PdbModel(Transformable, Parent):
charge = 0
occ = 0
beta = 0
x,y,z = [x for x in p.collapsedPosition()]
# x,y,z = [x for x in p.collapsedPosition()]
pos = p.collapsedPosition()
dig = [max(int(np.log10(np.abs(x))//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)]
assert(resid < 1e5)
resid = "{:<4d}".format(resid)
......
......@@ -13,7 +13,7 @@ from scipy import interpolate
from CanonicalNucleotideAtoms import canonicalNtFwd, canonicalNtRev, seqComplement
from CanonicalNucleotideAtoms import enmTemplateHC, enmTemplateSQ, enmCorrectionsHC
import pdb
# import pdb
"""
TODO:
+ fix handling of crossovers for atomic representation
......@@ -313,11 +313,18 @@ class Segment(ConnectableElement, Group):
def contour_to_orientation(self,s):
assert( isinstance(s,float) or isinstance(s,int) or len(s) == 1 ) # TODO make vectorized version
orientation = None
if self.start_orientation is not None:
# axis = self.start_orientation.dot( np.array((0,0,1)) )
if self.quaternion_spline_params is None:
axis = self.contour_to_tangent(s)
orientation = rotationAboutAxis( axis, self.twist_per_nt*self.contour_to_nt_pos(s), normalizeAxis=True )
else:
q = interpolate.splev( s, self.quaternion_spline_params )
if len(q) > 1: q = np.array(q).T # TODO: is this needed?
orientation = quaternion_to_matrix(q)
return orientation
def get_contour_sorted_connections_and_locations(self,type_):
sort_fn = lambda c: c[1].address
......@@ -377,7 +384,7 @@ class Segment(ConnectableElement, Group):
nt_dict = canonicalNtFwd
else:
nt_dict = canonicalNtRev
atoms = nt_dict[ key ].duplicate() # TODO: clone
atoms = nt_dict[ key ].generate() # TODO: clone?
atoms.orientation = orientation.dot(atoms.orientation) # this one should be correct
atoms.position = pos
......@@ -693,7 +700,7 @@ class DoubleStrandedSegment(Segment):
self.twist_per_nt = float(360 * num_turns) / num_nts
if start_orientation is None:
start_orientation = np.array(((1,0,0),(0,1,0),(0,0,1)))
start_orientation = np.eye(3) # np.array(((1,0,0),(0,1,0),(0,0,1)))
self.start_orientation = start_orientation
self.nicks = []
......@@ -1123,17 +1130,18 @@ class SegmentModel(ArbdModel):
if np.any( [B.particle is None for c,A,B in cabs] ):
print( "WARNING: none type found in connection, skipping" )
cabs = [e for e in cabs if e[2].particle is not None]
beads = list(set(s.beads + [A.particle for c,A,B in cabs]))
beads = set(s.beads + [A.particle for c,A,B in cabs])
## Add nearby beads
for c,A,B in cabs:
## TODOTODO test?
filter_fn = lambda x: x is not None and x not in beads
bs = list( filter( filter_fn, B.particle.intrahelical_neighbors ) )
beads.extend(bs)
beads.update(bs)
bs = list( filter( filter_fn, [n for b in bs for n in b.intrahelical_neighbors] ) )
beads.extend(bs)
beads.update(bs)
beads = list(beads)
## Skip beads that are None (some locations were not assigned a particle to avoid double-counting)
beads = [b for b in beads if b is not None]
......@@ -1167,7 +1175,7 @@ class SegmentModel(ArbdModel):
lastq = None
for b,t in zip(beads,tangents):
o = b.orientation_bead
positions
# positions
# angleVec = o.position - b.position
angleVec = bead_coordinates[o.idx,:] - bead_coordinates[b.idx,:]
angleVec = angleVec - angleVec.dot(t)*t
......
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