-
Notifications
You must be signed in to change notification settings - Fork 45
/
buildPhosOxy.py
162 lines (127 loc) · 6.37 KB
/
buildPhosOxy.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
#!/usr/bin/env python
"""Functions for building non-bridging phosphoryl oxygens"""
# Copyright 2010 Kevin Keating
#
# Licensed under the Educational Community License, Version 2.0 (the
# "License"); you may not use this file except in compliance with the
# License. You may obtain a copy of the License at
#
# http://www.osedu.org/licenses/ECL-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an "AS IS"
# BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
# or implied. See the License for the specific language governing
# permissions and limitations under the License.
from math import radians, cos, sin
from strucCalc import minus, plus, dotProd, crossProd, scalarProd, magnitude, normalize
#constants
PHOSBONDLENGTH = 1.485
PHOSBONDANGLE = 119.6
INITPHOSANGLE = 108 #the angle between the O5' and the phosphoryl oxygens in the first phosphate of a segment
def buildPhosOxy(curResAtoms, prevResAtoms):
"""Calculate non-bridging phosphoryl oxygen coordinates using the coordinates of the current and previous nucleotides
ARGUMENTS:
curResAtoms - a dictionary of the current nucleotide coordinates (i.e. the nucleotide to build the phosphoryl oxygens on) in the format atomName: [x, y, z]
prevResAtoms - a dictionary of the previous nucleotide coordinates in the format atomName: [x, y, z]
RETURNS:
phosOxyCoords - a dictionary of the phosphoryl oxygen coordinates in the format atomName: [x, y, z]
or None if the residue is missing the P, O5', or O3' atoms
"""
try:
P = curResAtoms ["P"]
O5 = curResAtoms ["O5'"]
O3 = prevResAtoms["O3'"]
except KeyError:
return None
#calculate a line from O5' to O3'
norm = minus(O5, O3)
#calculate the intersection of a plane (with normal $norm and point $P) and a line (from O5' to O3')
#using formula from http://local.wasp.uwa.edu.au/~pbourke/geometry/planeline/
# norm dot (P - O3)
# i = ------------------------
# norm dot (O5 - O3)
#intersecting point = O3 + i(O5 - O3)
i = dotProd(norm, minus(P, O3)) / dotProd(norm, minus(O5, O3))
interPoint = plus(O3, scalarProd(i, minus(O5, O3)))
#move $interPoint so that the distance from $P to $interPoint is 1.485 (the length of the P-OP1 bond)
#we also reflect the point about P
PIline = minus(P, interPoint) #here's is where the reflection occurs, because we do $P-$interPoint instead of $interPoint-$P
#scaledPoint = scalarProd(1/magnitude(PIline) * PHOSBONDLENGTH, PIline)
scaledPoint = normalize(PIline, PHOSBONDLENGTH)
#to get the new point location, we would do P + scaledPoint
#but we need to rotate the point first before translating it back
#rotate this new point by 59.8 and -59.8 degrees to determine the phosphoryl oxygen locations
#we rotate about the axis defined by $norm
angle = radians(PHOSBONDANGLE / 2)
(x, y, z) = scaledPoint
#unitnorm = scalarProd( 1/magnitude(norm), norm)
unitnorm = normalize(norm)
(u, v, w) = unitnorm
a = u*x + v*y + w*z
phosOxyCoords = {}
for (atomName, theta) in (("OP1", angle), ("OP2", -angle)):
cosTheta = cos(theta)
sinTheta = sin(theta)
#perform the rotation, and then add $P to the coordinates
newX = a*u + (x-a*u)*cosTheta + (v*z-w*y)*sinTheta + P[0]
newY = a*v + (y-a*v)*cosTheta + (w*x-u*z)*sinTheta + P[1]
newZ = a*w + (z-a*w)*cosTheta + (u*y-v*x)*sinTheta + P[2]
phosOxyCoords[atomName] = [newX, newY, newZ]
return phosOxyCoords
def buildInitOrTerminalPhosOxy(curResAtoms, prevResAtoms = None):
"""build phosphoryl oxygens using only the O3' or O5'atom (intended for the first or last nucleotide of a chain/segment)
ARGUMENTS:
curResAtoms - a dictionary of the current nucleotide coordinates (i.e. the nucleotide to build the phosphoryl oxygens on) in the format atomName: [x, y, z]
this dictionary must contain at least the phosphate
OPTIONAL ARGUMENTS:
prevResAtoms - a dictionary of the previous nucleotide coordinates in the format atomName: [x, y, z]
if provided, the curResAtoms["P"] and prevResAtoms["O3'"] will be used to place the phosphoryl oxygens
if not provided, the curResAtoms["P"] and curResAtoms["O5'"]
RETURNS:
phosOxyCoords - a dictionary of the phosphoryl oxygen coordinates in the format atomName: [x, y, z]
or None if the residue is missing the necessary atoms
"""
try:
P = curResAtoms ["P"]
if prevResAtoms is not None:
O = prevResAtoms["O3'"]
C = prevResAtoms["C3'"]
else:
O = curResAtoms["O5'"]
C = curResAtoms["C5'"]
except KeyError:
return None
#place atom along the P-O5' bond that is the appropriate distance from P
phosOxy = minus(O, P)
phosOxy = normalize(phosOxy, PHOSBONDLENGTH)
#define plane with C5'-O5'-P
norm = crossProd( minus(C, O), minus(P, O))
norm = normalize(norm)
#rotate dummy atom in plane about P by INITPHOSANGLE (which is the appropriate O5'-OP1 angle)
(x, y, z) = phosOxy
(u, v, w) = norm
theta = -radians(INITPHOSANGLE) #use the negative angle so that we rotate away from the C5'
cosTheta = cos(theta)
sinTheta = sin(theta)
#perform the rotation
a = u*x + v*y + w*z
phosOxy[0] = a*u + (x-a*u)*cosTheta + (v*z-w*y)*sinTheta
phosOxy[1] = a*v + (y-a*v)*cosTheta + (w*x-u*z)*sinTheta
phosOxy[2] = a*w + (z-a*w)*cosTheta + (u*y-v*x)*sinTheta
#rotate dummy atom about O5'-P axis by 59.8 and -59.8 degrees
norm = minus(O, P)
norm = normalize(norm)
(x, y, z) = phosOxy
(u, v, w) = norm
angle = radians(PHOSBONDANGLE / 2)
phosOxyCoords = {}
for (atomName, theta) in (("OP1", angle), ("OP2", -angle)):
cosTheta = cos(theta)
sinTheta = sin(theta)
#perform the rotation, and then add $P to the coordinates
newX = a*u + (x-a*u)*cosTheta + (v*z-w*y)*sinTheta + P[0]
newY = a*v + (y-a*v)*cosTheta + (w*x-u*z)*sinTheta + P[1]
newZ = a*w + (z-a*w)*cosTheta + (u*y-v*x)*sinTheta + P[2]
phosOxyCoords[atomName] = [newX, newY, newZ]
return phosOxyCoords