Source code for mrdna.model.CanonicalNucleotideAtoms

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)