CanonicalNucleotideAtoms.py

import json
import numpy as np
import copy
from  coords import rotationAboutAxis

seqComplement = dict(A='T',G='C')
for x,y in list(seqComplement.items()): seqComplement[y] = x


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

class CanonicalNucleotide():
    # DefaultOrientation = rotationAboutAxis([0,0,1], 68)
    DefaultOrientation = rotationAboutAxis([0,0,1], 68+180)
    def __init__(self, prefix, seq):
        self.sequence = seq
        
        d = np.loadtxt("%s.dat" % prefix, dtype='i4,S4,S4,f4,f4,f4,f4,f4')

        self.indices = d['f0']
        self.props = dict()
        for k1,k2 in zip(('charge','mass','x','y','z'),
                         ('f3','f4','f5','f6','f7')):
            self.props[k1] = d[k2]

        for k1,k2 in zip(('name','type'),('f1','f2')):
            self.props[k1] = [str(x,'utf-8') for x in d[k2]]

        self.props['beta'] = np.ones(len(d['f1']))
        # self.props['name'] = [str(x) for x in d['f1']]
        # self.props['type'] = [str(x) for x in d['f2']]

        minIdx = np.min(self.indices) - 1
        b = np.loadtxt("%s.bonds.txt" % prefix, dtype='i4') - minIdx
        self.bonds = np.array( [x for x in zip(b[::2],b[1::2])] )

        a = np.loadtxt("%s.angles.txt" % prefix, dtype='i4') - minIdx
        self.angles = np.array( [x for x in zip(a[::3],a[1::3],a[2::3])] )

        d = np.loadtxt("%s.dihedrals.txt" % prefix, dtype='i4') - minIdx
        self.dihedrals = np.array( [x for x in zip(d[::4],d[1::4],d[2::4],d[3::4])] )

        i = np.loadtxt("%s.impropers.txt" % prefix, dtype='i4') - minIdx
        self.impropers = np.array( [x for x in zip(i[::4],i[1::4],i[2::4],i[3::4])] )
        
        self.indices = self.indices - minIdx

    def transformed(self, orientation, translation, scale=0.5, singleStranded=False):
        # rotationAboutZ -= 68
        c = copy.copy(self)
        c.props = copy.deepcopy(self.props)

        ## Find transformation
        R = CanonicalNucleotide.DefaultOrientation.dot( orientation )
        
        ## Find offset of nt (shift for ssDNA)
        center = np.array((0,0,0))
        if singleStranded:
            i0 = self.index("C3'")-1
            center = np.array([c.props[k][i0] for k in ('x','y','z')])


        ## Apply transformation
        for i in range(len(c)):
            rvec = np.array([c.props[k][i] for k in ('x','y','z')]) - center
            rvec = rvec.dot( R )
            for k,r,o in zip(('x','y','z'), rvec, translation):
                c.props[k][i] = r + o

        if scale != 1:
            if singleStranded:
                i0 = self.index("C3'")-1
                r0 = np.array([c.props[k][i0] for k in ('x','y','z')])
                for i in range(len(c)):
                    rvec = np.array([c.props[k][i] for k in ('x','y','z')])
                    rvec = scale*(rvec-r0) + r0
                    for k,r in zip(('x','y','z'), rvec):
                        c.props[k][i] = r
                    c.props["beta"][i] = 0
            else:
                if self.sequence in ("A","G"):
                    i0 = self.index("N9")-1
                else:
                    i0 = self.index("N1")-1
                r0 = np.array([c.props[k][i0] for k in ('x','y','z')])
                for i in range(len(c)):
                    if c.props['name'][i][-1] in ("'","P","T"):
                        rvec = np.array([c.props[k][i] for k in ('x','y','z')])
                        rvec = scale*(rvec-r0) + r0
                        for k,r in zip(('x','y','z'), rvec):
                            c.props[k][i] = r
                        c.props["beta"][i] = 0
        
        return c
    
    # def __iter__(self):
    #     for a in zip(self.props['name'], self.coords):
    #         yield a

    def __len__(self):
        return len(self.indices)

    def index(self, name):
        for i,n in zip(self.indices, self.props['name']):
            if n == name: return i
        raise Exception("No atomic index found for atom '%s'!" % name)

    def atomicProps(self):
        for i in range(len(self)):
            tmp = dict()
            for k,a in self.props.items():
                tmp[k] = a[i]
            yield tmp
    
canonicalNtFwd = dict()
direction = 'fwd'
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" )

canonicalNtRev = dict()
direction = 'rev'
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" )

with open("resources/enm-template-honeycomb.json") as ch:
    enmTemplateHC = json.load(ch)
with open("resources/enm-template-square.json") as ch:
    enmTemplateSQ = json.load(ch)
with open("resources/enm-corrections-honecomb.json") as ch:
    enmCorrectionsHC = json.load(ch)