CanonicalNucleotideAtoms.py 5.44 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
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
class Nucleotide(Group):
    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

    def get_intrahelical_above(self):
        """ Returns next nucleotide along 5'-to-3' in same strand-side of helix """
        return self.get_intrahelical_nucleotide(1)

    def get_intrahelical_below(self):
        """ Returns last nucleotide along 5'-to-3' in same strand-side of helix """
        return self.get_intrahelical_nucleotide(-1)

    def get_intrahelical_nucleotide(self,offset=1):
        """ Returns nucleotide in same strand-side of helix with some offset along 5'-to-3'

        TODO: Does not yet cross intrahelical connections
        """

        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:
            return None
        else:
            return seg._get_atomic_nucleotide(nt_idx, is_fwd)

61
class CanonicalNucleotideFactory(Group):
62
63
    # DefaultOrientation = rotationAboutAxis([0,0,1], 180)
    DefaultOrientation = rotationAboutAxis([1,0,0], 180)
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
    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]
84
85
86
            type_ = ParticleType(tn.decode('UTF-8'),
                                 charge = d['f3'][i],
                                 mass   = d['f4'][i])
87
            p = PointParticle( type_ = type_,
88
                               position = [xs[i],ys[i],zs[i]],
89
90
                               name = names[i].decode('UTF-8'),
                               beta = 1)
91
92
93
94
            self.children.append( p )
        
        atoms = self.children

95

96
97
98
99
100
101
        # 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

102
103
        self._name_to_idx = { atoms[i].name:i for i in range(len(atoms)) }

104
    def generate(self):
105
106
107
108
109
110
111
        atoms = [copy.copy(c) for c in self.children]
        nt = Nucleotide(
            children = atoms,
            orientation = CanonicalNucleotideFactory.DefaultOrientation,
            sequence = self.sequence,
            resname = self.resname,
            factory = self)
112
113
114
115
        
        b = self.bonds
        a = self.angles
        d = self.dihedrals
cmaffeo2's avatar
cmaffeo2 committed
116
        imp = self.impropers
117
118
119
120
121
122
        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
123
        for i,j,k,l in zip(imp[::4],imp[1::4],imp[2::4],imp[3::4]):
124
125
126
            nt.add_improper( atoms[i], atoms[j], atoms[k], atoms[l], None )
        return nt

cmaffeo2's avatar
cmaffeo2 committed
127
128
129
130
    
canonicalNtFwd = dict()
direction = 'fwd'
for seq in seqComplement.keys():
131
132
133
    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
134
        canonicalNtFwd[key] = CanonicalNucleotideFactory( prefix, seq )
cmaffeo2's avatar
cmaffeo2 committed
135
136
137
138

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

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)