CanonicalNucleotideAtoms.py 5.7 KB
Newer Older
1
import json
cmaffeo2's avatar
cmaffeo2 committed
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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():
20
21
    # DefaultOrientation = rotationAboutAxis([0,0,1], 68)
    DefaultOrientation = rotationAboutAxis([0,0,1], 68+180)
22
23
    def __init__(self, prefix, seq):
        self.sequence = seq
cmaffeo2's avatar
cmaffeo2 committed
24
25
26
27
28
29
30
31
32
33
34
35
        
        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]]

36
        self.props['beta'] = np.ones(len(d['f1']))
cmaffeo2's avatar
cmaffeo2 committed
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
        # 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

55
    def transformed(self, orientation, translation, scale=0.5, singleStranded=False):
cmaffeo2's avatar
cmaffeo2 committed
56
57
58
59
60
61
        # rotationAboutZ -= 68
        c = copy.copy(self)
        c.props = copy.deepcopy(self.props)

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

69

cmaffeo2's avatar
cmaffeo2 committed
70
71
        ## Apply transformation
        for i in range(len(c)):
72
            rvec = np.array([c.props[k][i] for k in ('x','y','z')]) - center
73
            rvec = rvec.dot( R )
cmaffeo2's avatar
cmaffeo2 committed
74
75
76
            for k,r,o in zip(('x','y','z'), rvec, translation):
                c.props[k][i] = r + o

77
        if scale != 1:
78
79
80
81
82
83
84
85
            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
86
                    c.props["beta"][i] = 0
87
88
89
90
91
92
93
94
95
96
97
98
            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
99
                        c.props["beta"][i] = 0
100
        
cmaffeo2's avatar
cmaffeo2 committed
101
102
        return c
    
103
104
105
106
    # def __iter__(self):
    #     for a in zip(self.props['name'], self.coords):
    #         yield a

cmaffeo2's avatar
cmaffeo2 committed
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
    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():
125
126
127
    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
128
        canonicalNtFwd[key] = CanonicalNucleotide( prefix, seq )
129
130
131
    # canonicalNtFwd[seq] = CanonicalNucleotide( prefix )
    # canonicalNtFwd3[seq] = CanonicalNucleotide( prefix+"3prime" )
    # canonicalNtFwd5[seq] = CanonicalNucleotide( prefix+"5prime" )
cmaffeo2's avatar
cmaffeo2 committed
132
133
134
135

canonicalNtRev = dict()
direction = 'rev'
for seq in seqComplement.keys():
136
137
138
    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
139
        canonicalNtRev[key] = CanonicalNucleotide( prefix, seq )
140
141
142
143
144
145
146
147
148
149
150
151
152

# 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)