CanonicalNucleotideAtoms.py 3.88 KB
Newer Older
1
import json
cmaffeo2's avatar
cmaffeo2 committed
2
3
4
import numpy as np
import copy
from  coords import rotationAboutAxis
cmaffeo2's avatar
cmaffeo2 committed
5
from arbdmodel import Group, PointParticle, ParticleType
cmaffeo2's avatar
cmaffeo2 committed
6
7

seqComplement = dict(A='T',G='C')
8
resnames = dict(A="ADE", G="GUA", C="CYT", T="THY")
cmaffeo2's avatar
cmaffeo2 committed
9
10
11
12
13
14
15
16
17
18
19
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

20
class CanonicalNucleotideFactory(Group):
cmaffeo2's avatar
cmaffeo2 committed
21
    DefaultOrientation = rotationAboutAxis([0,0,1], 180)
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
    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]
42
43
44
            type_ = ParticleType(tn.decode('UTF-8'),
                                 charge = d['f3'][i],
                                 mass   = d['f4'][i])
45
            p = PointParticle( type_ = type_,
46
                               position = [xs[i],ys[i],zs[i]],
47
48
                               name = names[i].decode('UTF-8'),
                               beta = 1)
49
50
51
52
53
54
55
56
57
58
59
            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):
60
        nt = Group(orientation = CanonicalNucleotideFactory.DefaultOrientation)
61
62
63
64
65
66
67
68
69
70
        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
cmaffeo2's avatar
cmaffeo2 committed
71
        imp = self.impropers
72
73
74
75
76
77
        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 )
cmaffeo2's avatar
cmaffeo2 committed
78
        for i,j,k,l in zip(imp[::4],imp[1::4],imp[2::4],imp[3::4]):
79
80
81
            nt.add_improper( atoms[i], atoms[j], atoms[k], atoms[l], None )
        return nt

cmaffeo2's avatar
cmaffeo2 committed
82
83
84
85
    
canonicalNtFwd = dict()
direction = 'fwd'
for seq in seqComplement.keys():
86
87
88
    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
89
        canonicalNtFwd[key] = CanonicalNucleotideFactory( prefix, seq )
cmaffeo2's avatar
cmaffeo2 committed
90
91
92
93

canonicalNtRev = dict()
direction = 'rev'
for seq in seqComplement.keys():
94
95
96
    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
97
        canonicalNtRev[key] = CanonicalNucleotideFactory( prefix, seq )
98
99
100
101
102
103
104

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)