import json
import numpy as np
import copy
from ..arbdmodel.coords import rotationAboutAxis
from ..arbdmodel import Group, PointParticle, ParticleType
from .. import get_resource_path
seqComplement = dict(A='T',G='C')
resnames = dict(A="ADE", G="GUA", C="CYT", T="THY")
for x,y in list(seqComplement.items()): seqComplement[y] = x
[docs]
def stringToIntTuples(string, tupleLen, offset):
l = string.split()
# for i in range(0,len(l),tupleLen):
# ret.append( (int(x)-offset for x in l[i:i+tupleLen]) )
ret = [ (int(x)-offset for x in l[i:i+tupleLen])
for i in range(0,len(l),tupleLen) ]
return ret
[docs]
class Nucleotide(Group):
[docs]
def __init__(self, children, orientation, sequence, resname, factory):
Group.__init__(self, children=children, orientation = orientation)
self.atoms_by_name = {a.name:a for a in self.children}
self._factory = factory
self.sequence = sequence
self.resname = resname
self._first_atomic_index = 0
def _get_atomic_index(self, name):
return self._factory._name_to_idx[name] + self._first_atomic_index
[docs]
def get_intrahelical_above(self):
""" Returns next nucleotide along 5'-to-3' in same strand-side of helix """
return self.get_intrahelical_nucleotide(1)
[docs]
def get_intrahelical_below(self):
""" Returns last nucleotide along 5'-to-3' in same strand-side of helix """
return self.get_intrahelical_nucleotide(-1)
[docs]
def get_intrahelical_nucleotide(self,offset=1):
""" Returns nucleotide in same strand-side of helix with some offset along 5'-to-3' """
from ..segmentmodel import SingleStrandedSegment
strand_piece = self.parent
is_fwd = strand_piece.is_fwd
lo,hi = sorted([strand_piece.start,strand_piece.end])
if is_fwd:
nt_idx = strand_piece.children.index(self)+lo + offset
else:
nt_idx = hi-strand_piece.children.index(self) - offset
seg = strand_piece.segment
if nt_idx < 0 or nt_idx >= seg.num_nt:
""" Cross intrahelical connections if possible """
ret = seg._ntpos_to_seg_and_ntpos(nt_idx, is_fwd)
if ret is None:
ret = seg._ntpos_to_seg_and_ntpos(nt_idx, is_fwd)
return None
else:
seg, nt_idx, is_fwd = ret
if isinstance(seg, SingleStrandedSegment):
return None
else:
return seg._get_atomic_nucleotide(nt_idx, is_fwd)
else:
return seg._get_atomic_nucleotide(nt_idx, is_fwd)
[docs]
class CanonicalNucleotideFactory(Group):
DefaultOrientation = rotationAboutAxis([0,0,1], 180).dot(
rotationAboutAxis([1,0,0], 180))
[docs]
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]
type_ = ParticleType(tn.decode('UTF-8'),
charge = d['f3'][i],
mass = d['f4'][i])
p = PointParticle( type_ = type_,
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
self._name_to_idx = { atoms[i].name:i for i in range(len(atoms)) }
[docs]
def generate(self):
atoms = [copy.copy(c) for c in self.children]
nt = Nucleotide(
children = atoms,
orientation = CanonicalNucleotideFactory.DefaultOrientation,
sequence = self.sequence,
resname = self.resname,
factory = self)
b = self.bonds
a = self.angles
d = self.dihedrals
imp = 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(imp[::4],imp[1::4],imp[2::4],imp[3::4]):
nt.add_improper( atoms[i], atoms[j], atoms[k], atoms[l], None )
return nt
canonicalNtFwd = dict()
direction = 'fwd'
for seq in seqComplement.keys():
for keystring,suff in zip(("5%s","%s","%s3","5%s3"),("-5prime","","-3prime","-singlet")):
prefix = get_resource_path('%s-%s%s' % (seq,direction,suff))
key = keystring % seq
canonicalNtFwd[key] = CanonicalNucleotideFactory( prefix, seq )
canonicalNtRev = dict()
direction = 'rev'
for seq in seqComplement.keys():
for keystring,suff in zip(("5%s","%s","%s3","5%s3"),("-5prime","","-3prime","-singlet")):
prefix = get_resource_path('%s-%s%s' % (seq,direction,suff))
key = keystring % seq
canonicalNtRev[key] = CanonicalNucleotideFactory( prefix, seq )
with open(get_resource_path("enm-template-honeycomb.json")) as ch:
enmTemplateHC = json.load(ch)
with open(get_resource_path("enm-template-square.json")) as ch:
enmTemplateSQ = json.load(ch)
with open(get_resource_path("enm-corrections-honecomb.json")) as ch:
enmCorrectionsHC = json.load(ch)