hdiff output

  0:   0:    1:   1:    2:   2:    3:   3: 

r32605/CG2FA.py 2017-05-24 13:30:21.374598387 +0100 r32604/CG2FA.py 2017-05-24 13:30:24.418638135 +0100
  1: #! /usr/bin/env python  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/OPEP/FA2CG/CG2FA.py' in revision 32604
  2: # -*- coding: UTF8 -*- 
  3:  
  4: import gzip 
  5: import sys 
  6: from string import split 
  7: from math import * 
  8:  
  9: # --------------------- Functions ------------------------# 
 10:  
 11:  
 12: def rotateB(x, y, z, sxy, a, b): 
 13:     r = sqrt(x*x + y*y) 
 14:     R = sqrt(x*x + y*y + z*z) 
 15:     st = y/r 
 16:     ct = x/r 
 17:     sf = z/R 
 18:     cf = r/R 
 19:     for i in range(a, b): 
 20:         x = (sxy[i][6]) 
 21:         y = (sxy[i][7]) 
 22:         z = (sxy[i][8]) 
 23:         X = cf*ct*x - sf*ct*z - st*y 
 24:         Y = st*cf*x - st*sf*z + ct*y 
 25:         Z = sf*x + cf*z 
 26:         sxy[i][6] = X 
 27:         sxy[i][7] = Y 
 28:         sxy[i][8] = Z 
 29:  
 30:  
 31: def rotateF(x, y, z, sxy, a, b): 
 32:     r = sqrt(x*x + y*y) 
 33:     R = sqrt(x*x + y*y + z*z) 
 34:     st = y/r 
 35:     ct = x/r 
 36:     sf = z/R 
 37:     cf = r/R 
 38:     for i in range(a, b): 
 39:         x = sxy[i][6] 
 40:         y = sxy[i][7] 
 41:         z = sxy[i][8] 
 42:         X = cf*ct*x + cf*st*y + sf*z 
 43:         Y = -st*x + ct*y 
 44:         Z = -sf*ct*x - st*sf*y + cf*z 
 45:         sxy[i][6] = X 
 46:         sxy[i][7] = Y 
 47:         sxy[i][8] = Z 
 48:  
 49:  
 50: def translate(x, y, z, sxy, a, b): 
 51:     for i in range(a, b): 
 52:         sxy[i][6] = sxy[i][6]-x 
 53:         sxy[i][7] = sxy[i][7]-y 
 54:         sxy[i][8] = sxy[i][8]-z 
 55:  
 56:  
 57: def prstr(sxy): 
 58:     print "*" 
 59:     for i in range(L): 
 60:         print '%3.3f\t %3.3f\t %3.3f' % (sxy[i][6], sxy[i][7], sxy[i][8]) 
 61:  
 62:  
 63: def toPDB(s): 
 64:     for i in range(len(s)): 
 65:         print '%-6s%5d %-4s%1s%3s %1s%4d%1s   %8.3f%8.3f%8.3f%6.2f%6.2f %2s%2s' % (s[i][0], int(s[i][1]), s[i][2], ' ', s[i][3], s[i][4], int(s[i][5]), ' ', float(s[i][6]), float(s[i][7]), float(s[i][8]), float(s[i][9]), float(s[i][10]), s[i][11], ' ') 
 66:  
 67:  
 68: def GG_sugar(sxy, qxy, PL, i, Na): 
 69:     O5x = sxy[3][6] 
 70:     O5y = sxy[3][7] 
 71:     O5z = sxy[3][8] 
 72:     C5x = sxy[4][6] 
 73:     C5y = sxy[4][7] 
 74:     C5z = sxy[4][8] 
 75:     CAx = sxy[5][6] 
 76:     CAy = sxy[5][7] 
 77:     CAz = sxy[5][8] 
 78:     CYx = sxy[11][6] 
 79:     CYy = sxy[11][7] 
 80:     CYz = sxy[11][8] 
 81:  
 82:     O5gx = qxy[1][5] 
 83:     O5gy = qxy[1][6] 
 84:     O5gz = qxy[1][7] 
 85:     C5gx = qxy[2][5] 
 86:     C5gy = qxy[2][6] 
 87:     C5gz = qxy[2][7] 
 88:     CAgx = qxy[3][5] 
 89:     CAgy = qxy[3][6] 
 90:     CAgz = qxy[3][7] 
 91:     CYgx = qxy[4][5] 
 92:     CYgy = qxy[4][6] 
 93:     CYgz = qxy[4][7] 
 94:  
 95:     DO5 = sqrt((O5x-O5gx)*(O5x-O5gx)+(O5y-O5gy)*(O5y-O5gy)+(O5z-O5gz)*(O5z-O5gz)) 
 96:     DC5 = sqrt((C5x-C5gx)*(C5x-C5gx)+(C5y-C5gy)*(C5y-C5gy)+(C5z-C5gz)*(C5z-C5gz)) 
 97:     DCA = sqrt((CAx-CAgx)*(CAx-CAgx)+(CAy-CAgy)*(CAy-CAgy)+(CAz-CAgz)*(CAz-CAgz)) 
 98:     DCY = sqrt((CYx-CYgx)*(CYx-CYgx)+(CYy-CYgy)*(CYy-CYgy)+(CYz-CYgz)*(CYz-CYgz)) 
 99:  
100:     dist = sqrt(DO5*DO5+DC5*DC5+DCA*DCA+DCY*DCY) 
101:     return dist 
102:  
103:  
104: def GG_A(sxy, qxy, PL, i, Na): 
105:     B1x = (sxy[12][6]+sxy[13][6]+sxy[14][6]+sxy[15][6]+sxy[21][6])/5 
106:     B1y = (sxy[12][7]+sxy[13][7]+sxy[14][7]+sxy[15][7]+sxy[21][7])/5 
107:     B1z = (sxy[12][8]+sxy[13][8]+sxy[14][8]+sxy[15][8]+sxy[21][8])/5 
108:     B2x = (sxy[15][6]+sxy[16][6]+sxy[18][6]+sxy[19][6]+sxy[20][6]+sxy[21][6])/6 
109:     B2y = (sxy[15][7]+sxy[16][7]+sxy[18][7]+sxy[19][7]+sxy[20][7]+sxy[21][7])/6 
110:     B2z = (sxy[15][8]+sxy[16][8]+sxy[18][8]+sxy[19][8]+sxy[20][8]+sxy[21][8])/6 
111:     G1x = qxy[5][5] 
112:     G1y = qxy[5][6] 
113:     G1z = qxy[5][7] 
114:     G2x = qxy[6][5] 
115:     G2y = qxy[6][6] 
116:     G2z = qxy[6][7] 
117:  
118:     DB1 = sqrt((B1x-G1x)**2+(B1y-G1y)**2+(B1z-G1z)**2) 
119:     DB2 = sqrt((B2x-G2x)**2+(B2y-G2y)**2+(B2z-G2z)**2) 
120:  
121:     dist = sqrt(DB1**2+DB2**2) 
122:     return dist 
123:  
124:  
125: def GG_G(sxy, qxy, PL, i, Na): 
126:     B1x = (sxy[12][6]+sxy[13][6]+sxy[14][6]+sxy[15][6]+sxy[22][6])/5 
127:     B1y = (sxy[12][7]+sxy[13][7]+sxy[14][7]+sxy[15][7]+sxy[22][7])/5 
128:     B1z = (sxy[12][8]+sxy[13][8]+sxy[14][8]+sxy[15][8]+sxy[22][8])/5 
129:     B2x = (sxy[15][6]+sxy[16][6]+sxy[18][6]+sxy[19][6]+sxy[21][6]+sxy[22][6])/6 
130:     B2y = (sxy[15][7]+sxy[16][7]+sxy[18][7]+sxy[19][7]+sxy[21][7]+sxy[22][7])/6 
131:     B2z = (sxy[15][8]+sxy[16][8]+sxy[18][8]+sxy[19][8]+sxy[21][8]+sxy[22][8])/6 
132:  
133:     G1x = qxy[5][5] 
134:     G1y = qxy[5][6] 
135:     G1z = qxy[5][7] 
136:     G2x = qxy[6][5] 
137:     G2y = qxy[6][6] 
138:     G2z = qxy[6][7] 
139:  
140:     DB1 = sqrt((B1x-G1x)*(B1x-G1x)+(B1y-G1y)*(B1y-G1y)+(B1z-G1z)*(B1z-G1z)) 
141:     DB2 = sqrt((B2x-G2x)*(B2x-G2x)+(B2y-G2y)*(B2y-G2y)+(B2z-G2z)*(B2z-G2z)) 
142:  
143:     dist = sqrt(DB1*DB1+DB2*DB2) 
144:     return dist 
145:  
146:  
147: def GG_CU(sxy, qxy, PL, i, Na): 
148:     B1x = (sxy[12][6]+sxy[13][6]+sxy[15][6]+sxy[16][6]+sxy[18][6]+sxy[19][6])/6 
149:     B1y = (sxy[12][7]+sxy[13][7]+sxy[15][7]+sxy[16][7]+sxy[18][7]+sxy[19][7])/6 
150:     B1z = (sxy[12][8]+sxy[13][8]+sxy[15][8]+sxy[16][8]+sxy[18][8]+sxy[19][8])/6 
151:     G1x = qxy[5][5] 
152:     G1y = qxy[5][6] 
153:     G1z = qxy[5][7] 
154:     DB1 = sqrt((B1x-G1x)*(B1x-G1x)+(B1y-G1y)*(B1y-G1y)+(B1z-G1z)*(B1z-G1z)) 
155:  
156:     dist = sqrt(DB1*DB1) 
157:     return dist 
158:  
159:  
160: def rebuildA(qxy, Sxy, PL, i, Na): 
161:     L = len(Sxy) 
162:     LG = len(qxy) 
163:     wxy = [] 
164:     N = 33 
165:     for i in range(0, L): 
166:         sxy = [] 
167:         if(Sxy[i][2] == 'P'): 
168:             for j in range(N): 
169:                 sxy.append(Sxy[i+j]) 
170:  
171:             x = sxy[0][6] 
172:             y = sxy[0][7] 
173:             z = sxy[0][8] 
174:             translate(x, y, z, sxy, 0, N) 
175:  
176:             # Align CA 
177:             x = sxy[5][6] 
178:             y = sxy[5][7] 
179:             z = sxy[5][8] 
180:             rotateF(x, y, z, sxy, 0, N) 
181:             x = qxy[3][5]-qxy[0][5] 
182:             y = qxy[3][6]-qxy[0][6] 
183:             z = qxy[3][7]-qxy[0][7] 
184:             rotateB(x, y, z, sxy, 0, N) 
185:  
186:             # Center Cy and B2 
187:             x = sxy[11][6] 
188:             y = sxy[11][7] 
189:             z = sxy[11][8] 
190:             translate(x, y, z, sxy, 0, N) 
191:             x = (sxy[15][6]+sxy[16][6]+sxy[18][6]+sxy[19][6]+sxy[20][6]+sxy[21][6])/6 
192:             y = (sxy[15][7]+sxy[16][7]+sxy[18][7]+sxy[19][7]+sxy[20][7]+sxy[21][7])/6 
193:             z = (sxy[15][8]+sxy[16][8]+sxy[18][8]+sxy[19][8]+sxy[20][8]+sxy[21][8])/6 
194:             # rotateF(x, y, z, sxy, 12, 22) 
195:             # rotateF(x, y, z, sxy, 29, N) 
196:  
197:             # Align B2 
198:             x = qxy[6][5] - qxy[4][5] 
199:             y = qxy[6][6] - qxy[4][6] 
200:             z = qxy[6][7] - qxy[4][7] 
201:             # rotateB(x, y, z, sxy, 12, 22) 
202:             # rotateB(x, y, z, sxy, 29, N) 
203:  
204:             # translation ref base to base coordinates: P atoms coincide 
205:             x = sxy[0][6]-qxy[0][5] 
206:             y = sxy[0][7]-qxy[0][6] 
207:             z = sxy[0][8]-qxy[0][7] 
208:             translate(x, y, z, sxy, 0, N) 
209:  
210:             if(i == 0): 
211:                 wxy = sxy 
212:                 dist = 100 
213:  
214:             ndist = GG_sugar(sxy, qxy, PL, i, Na) 
215:             if ndist < dist: 
216:                 wxy = sxy 
217:                 dist = ndist 
218:  
219:     for i in range(0, L): 
220:         # print 'A', i 
221:         hxy = [] 
222:         if Sxy[i][2] == 'P': 
223:             for j in range(N): 
224:                 hxy.append(Sxy[i+j]) 
225:  
226:             # Center Cy and B2 
227:             x = hxy[11][6] 
228:             y = hxy[11][7] 
229:             z = hxy[11][8] 
230:             translate(x, y, z, hxy, 0, N) 
231:             x = (hxy[15][6]+hxy[16][6]+hxy[18][6]+hxy[19][6]+hxy[20][6]+hxy[21][6])/6 
232:             y = (hxy[15][7]+hxy[16][7]+hxy[18][7]+hxy[19][7]+hxy[20][7]+hxy[21][7])/6 
233:             z = (hxy[15][8]+hxy[16][8]+hxy[18][8]+hxy[19][8]+hxy[20][8]+hxy[21][8])/6 
234:             rotateF(x, y, z, hxy, 12, 22) 
235:             rotateF(x, y, z, hxy, 29, N) 
236:  
237:             # Align B2 
238:             x = qxy[6][5] - qxy[4][5] 
239:             y = qxy[6][6] - qxy[4][6] 
240:             z = qxy[6][7] - qxy[4][7] 
241:             rotateB(x, y, z, hxy, 12, 22) 
242:             rotateB(x, y, z, hxy, 29, N) 
243:  
244:             # print 'hxy-B' 
245:             # for k in hxy: 
246:             #     print k 
247:  
248:             # translation ref base to base coordinates: CY atoms coincide 
249:             x = hxy[11][6]-qxy[4][5] 
250:             y = hxy[11][7]-qxy[4][6] 
251:             z = hxy[11][8]-qxy[4][7] 
252:             translate(x, y, z, hxy, 0, N) 
253:  
254:             if(i == 0): 
255:                 dist = 100 
256:  
257:             # print 'hxy' 
258:             # for k in hxy: 
259:             #     print k 
260:  
261:             ndist = GG_A(hxy, qxy, PL, i, Na) 
262:             # print dist 
263:             if ndist < dist: 
264:                 for i in range(12, 22): 
265:                     wxy[i] = hxy[i] 
266:                 for i in range(29, N): 
267:                     wxy[i] = hxy[i] 
268:                 dist = ndist 
269:  
270:             # print 'wxy' 
271:             # for k in wxy: 
272:             #     print k 
273:  
274:     return wxy 
275:  
276:  
277: def rebuildG(qxy, Sxy, PL, i, Na): 
278:     L = len(Sxy) 
279:     LG = len(qxy) 
280:     wxy = [] 
281:     N = 34 
282:     for i in range(0, L): 
283:         sxy = [] 
284:         if(Sxy[i][2] == 'P'): 
285:             for j in range(N): 
286:                 sxy.append(Sxy[i+j]) 
287:  
288:             # print i 
289:             x = sxy[0][6] 
290:             y = sxy[0][7] 
291:             z = sxy[0][8] 
292:             translate(x, y, z, sxy, 0, N) 
293:  
294:             # Align CA 
295:             x = sxy[5][6] 
296:             y = sxy[5][7] 
297:             z = sxy[5][8] 
298:  
299:             rotateF(x, y, z, sxy, 0, N) 
300:             x = qxy[3][5]-qxy[0][5] 
301:             y = qxy[3][6]-qxy[0][6] 
302:             z = qxy[3][7]-qxy[0][7] 
303:             rotateB(x, y, z, sxy, 0, N) 
304:  
305:             # Center Cy and B2 
306:             x = sxy[11][6] 
307:             y = sxy[11][7] 
308:             z = sxy[11][8] 
309:             translate(x, y, z, sxy, 0, N) 
310:             x = (sxy[14][6]+sxy[16][6]+sxy[18][6]+sxy[19][6]+sxy[21][6]+sxy[22][6])/6 
311:             y = (sxy[14][7]+sxy[16][7]+sxy[18][7]+sxy[19][7]+sxy[21][7]+sxy[22][7])/6 
312:             z = (sxy[14][8]+sxy[16][8]+sxy[18][8]+sxy[19][8]+sxy[21][8]+sxy[22][8])/6 
313:             # rotateF(x, y, z, sxy, 12, 23) 
314:             # rotateF(x, y, z, sxy, 29, N) 
315:  
316:             # Align B2 
317:             x = qxy[6][5] - qxy[4][5] 
318:             y = qxy[6][6] - qxy[4][6] 
319:             z = qxy[6][7] - qxy[4][7] 
320:             # rotateB(x, y, z, sxy, 12, 23) 
321:             # rotateB(x, y, z, sxy, 29, N) 
322:  
323:             # translation ref base to base coordinates: P atoms coincide 
324:             x = sxy[0][6]-qxy[0][5] 
325:             y = sxy[0][7]-qxy[0][6] 
326:             z = sxy[0][8]-qxy[0][7] 
327:             translate(x, y, z, sxy, 0, N) 
328:  
329:             if(i == 0): 
330:                 wxy = sxy 
331:                 dist = 100 
332:  
333:             ndist = GG_sugar(sxy, qxy, PL, i, Na) 
334:             if ndist < dist: 
335:                 wxy = sxy 
336:                 dist = ndist 
337:  
338:             # print 'wxy - sugar', ndist 
339:             # for k in wxy: 
340:             #     print k 
341:  
342:     for i in range(0, L): 
343:         # print 'G', i 
344:         hxy = [] 
345:         if Sxy[i][2] == 'P': 
346:             for j in range(N): 
347:                 hxy.append(Sxy[i+j]) 
348:  
349:             # Center Cy and B2 
350:             x = hxy[11][6] 
351:             y = hxy[11][7] 
352:             z = hxy[11][8] 
353:             translate(x, y, z, hxy, 0, N) 
354:             x = (hxy[14][6]+hxy[16][6]+hxy[18][6]+hxy[19][6]+hxy[21][6]+hxy[22][6])/6 
355:             y = (hxy[14][7]+hxy[16][7]+hxy[18][7]+hxy[19][7]+hxy[21][7]+hxy[22][7])/6 
356:             z = (hxy[14][8]+hxy[16][8]+hxy[18][8]+hxy[19][8]+hxy[21][8]+hxy[22][8])/6 
357:             rotateF(x, y, z, hxy, 12, 23) 
358:             rotateF(x, y, z, hxy, 29, N) 
359:  
360:             # Align B2 
361:             x = qxy[6][5] - qxy[4][5] 
362:             y = qxy[6][6] - qxy[4][6] 
363:             z = qxy[6][7] - qxy[4][7] 
364:             rotateB(x, y, z, hxy, 12, 23) 
365:             rotateB(x, y, z, hxy, 29, N) 
366:  
367:             # print 'hxy-prima' 
368:             # for k in hxy: 
369:             #     print k 
370:  
371:             # translation ref base to base coordinates: P atoms coincide 
372:             x = hxy[11][6]-qxy[4][5] 
373:             y = hxy[11][7]-qxy[4][6] 
374:             z = hxy[11][8]-qxy[4][7] 
375:             translate(x, y, z, hxy, 0, N) 
376:  
377:             # print 'hxy' 
378:             # for k in hxy: 
379:             #     print k 
380:  
381:             if i == 0: 
382:                 dist = 100 
383:  
384:             ndist = GG_G(hxy, qxy, PL, i, Na) 
385:             # print ndist 
386:             if ndist < dist: 
387:                 for i in range(12, 23): 
388:                     wxy[i] = hxy[i] 
389:                 for i in range(29, N): 
390:                     wxy[i] = hxy[i] 
391:                 dist = ndist 
392:  
393:             # print 'wxy - final', ndist 
394:             # for k in wxy: 
395:             #     print k 
396:  
397:     return wxy 
398:  
399:  
400: def rebuildC(qxy, Sxy, PL, i, Na): 
401:     L = len(Sxy) 
402:     LG = len(qxy) 
403:     wxy = [] 
404:     N = 31 
405:     for i in range(0, L): 
406:         sxy = [] 
407:         if Sxy[i][2] == 'P': 
408:             for j in range(N): 
409:                 sxy.append(Sxy[i+j]) 
410:  
411:             # print i 
412:             x = sxy[0][6] 
413:             y = sxy[0][7] 
414:             z = sxy[0][8] 
415:             translate(x, y, z, sxy, 0, N) 
416:  
417:             # Align CA 
418:             x = sxy[5][6] 
419:             y = sxy[5][7] 
420:             z = sxy[5][8] 
421:  
422:             rotateF(x, y, z, sxy, 0, N) 
423:             x = qxy[3][5]-qxy[0][5] 
424:             y = qxy[3][6]-qxy[0][6] 
425:             z = qxy[3][7]-qxy[0][7] 
426:             rotateB(x, y, z, sxy, 0, N) 
427:  
428:             # Center Cy and B2 
429:             x = sxy[11][6] 
430:             y = sxy[11][7] 
431:             z = sxy[11][8] 
432:             translate(x, y, z, sxy, 0, N) 
433:             x = (sxy[12][6]+sxy[13][6]+sxy[15][6]+sxy[16][6]+sxy[18][6]+sxy[19][6])/6 
434:             y = (sxy[12][7]+sxy[13][7]+sxy[15][7]+sxy[16][7]+sxy[18][7]+sxy[19][7])/6 
435:             z = (sxy[12][8]+sxy[13][8]+sxy[15][8]+sxy[16][8]+sxy[18][8]+sxy[19][8])/6 
436:             # rotateF(x, y, z, sxy, 12, 20) 
437:             # rotateF(x, y, z, sxy, 27, N) 
438:  
439:             # Align B2 
440:             x = qxy[5][5] - qxy[4][5] 
441:             y = qxy[5][6] - qxy[4][6] 
442:             z = qxy[5][7] - qxy[4][7] 
443:             # rotateB(x, y, z, sxy, 12, 20) 
444:             # rotateB(x, y, z, sxy, 27, N) 
445:  
446:             # translation ref base to base coordinates: P atoms coincide 
447:             x = sxy[0][6]-qxy[0][5] 
448:             y = sxy[0][7]-qxy[0][6] 
449:             z = sxy[0][8]-qxy[0][7] 
450:             translate(x, y, z, sxy, 0, N) 
451:  
452:             if(i == 0): 
453:                 wxy = sxy 
454:                 dist = 100 
455:  
456:             ndist = GG_sugar(sxy, qxy, PL, i, Na) 
457:             if ndist < dist: 
458:                 wxy = sxy 
459:                 dist = ndist 
460:  
461:     for i in range(0, L): 
462:         # print 'C', i 
463:         hxy = [] 
464:         if Sxy[i][2] == 'P': 
465:             for j in range(N): 
466:                 hxy.append(Sxy[i+j]) 
467:  
468:             # Center Cy and B2 
469:             x = hxy[11][6] 
470:             y = hxy[11][7] 
471:             z = hxy[11][8] 
472:             translate(x, y, z, hxy, 0, N) 
473:             x = (hxy[12][6]+hxy[13][6]+hxy[15][6]+hxy[16][6]+hxy[18][6]+hxy[19][6])/6 
474:             y = (hxy[12][7]+hxy[13][7]+hxy[15][7]+hxy[16][7]+hxy[18][7]+hxy[19][7])/6 
475:             z = (hxy[12][8]+hxy[13][8]+hxy[15][8]+hxy[16][8]+hxy[18][8]+hxy[19][8])/6 
476:             rotateF(x, y, z, hxy, 12, 20) 
477:             rotateF(x, y, z, hxy, 27, N) 
478:  
479:             # Align B2 
480:             x = qxy[5][5] - qxy[4][5] 
481:             y = qxy[5][6] - qxy[4][6] 
482:             z = qxy[5][7] - qxy[4][7] 
483:             rotateB(x, y, z, hxy, 12, 20) 
484:             rotateB(x, y, z, hxy, 27, N) 
485:  
486:             # translation ref base to base coordinates: P atoms coincide 
487:             x = hxy[11][6]-qxy[4][5] 
488:             y = hxy[11][7]-qxy[4][6] 
489:             z = hxy[11][8]-qxy[4][7] 
490:             translate(x, y, z, hxy, 0, N) 
491:  
492:             if(i == 0): 
493:                 dist = 100 
494:  
495:             ndist = GG_CU(hxy, qxy, PL, i, Na) 
496:             if ndist < dist: 
497:                 for i in range(12, 20): 
498:                     wxy[i] = hxy[i] 
499:                 for i in range(27, N): 
500:                     wxy[i] = hxy[i] 
501:                 dist = ndist 
502:  
503:     return wxy 
504:  
505:  
506: def rebuildU(qxy, Sxy, PL, i, Na): 
507:     L = len(Sxy) 
508:     LG = len(qxy) 
509:     wxy = [] 
510:     N = 30 
511:     for i in range(0, L): 
512:         sxy = [] 
513:         if(Sxy[i][2] == 'P'): 
514:             for j in range(N): 
515:                 sxy.append(Sxy[i+j]) 
516:  
517:             x = sxy[0][6] 
518:             y = sxy[0][7] 
519:             z = sxy[0][8] 
520:             translate(x, y, z, sxy, 0, N) 
521:  
522:             # Align CA 
523:             x = sxy[5][6] 
524:             y = sxy[5][7] 
525:             z = sxy[5][8] 
526:  
527:             rotateF(x, y, z, sxy, 0, N) 
528:             x = qxy[3][5]-qxy[0][5] 
529:             y = qxy[3][6]-qxy[0][6] 
530:             z = qxy[3][7]-qxy[0][7] 
531:             rotateB(x, y, z, sxy, 0, N) 
532:  
533:             # Center Cy and B2 
534:             x = sxy[11][6] 
535:             y = sxy[11][7] 
536:             z = sxy[11][8] 
537:             translate(x, y, z, sxy, 0, N) 
538:             x = (sxy[12][6]+sxy[13][6]+sxy[15][6]+sxy[16][6]+sxy[18][6]+sxy[19][6])/6 
539:             y = (sxy[12][7]+sxy[13][7]+sxy[15][7]+sxy[16][7]+sxy[18][7]+sxy[19][7])/6 
540:             z = (sxy[12][8]+sxy[13][8]+sxy[15][8]+sxy[16][8]+sxy[18][8]+sxy[19][8])/6 
541:             # rotateF(x, y, z, sxy, 12, 20) 
542:             # rotateF(x, y, z, sxy, 27, N) 
543:  
544:             # Align B2 
545:             x = qxy[5][5] - qxy[4][5] 
546:             y = qxy[5][6] - qxy[4][6] 
547:             z = qxy[5][7] - qxy[4][7] 
548:             # rotateB(x, y, z, sxy, 12, 20) 
549:             # rotateB(x, y, z, sxy, 27, N) 
550:  
551:             # translation ref base to base coordinates: P atoms coincide 
552:             x = sxy[0][6]-qxy[0][5] 
553:             y = sxy[0][7]-qxy[0][6] 
554:             z = sxy[0][8]-qxy[0][7] 
555:             translate(x, y, z, sxy, 0, N) 
556:  
557:             if(i == 0): 
558:                 wxy = sxy 
559:                 dist = 100 
560:  
561:             ndist = GG_sugar(sxy, qxy, PL, i, Na) 
562:             if ndist < dist: 
563:                 wxy = sxy 
564:                 dist = ndist 
565:  
566:     for i in range(0, L): 
567:         # print 'U', i 
568:         hxy = [] 
569:         if(Sxy[i][2] == 'P'): 
570:             for j in range(N): 
571:                 hxy.append(Sxy[i+j]) 
572:  
573:             # Center Cy and B2 
574:             x = hxy[11][6] 
575:             y = hxy[11][7] 
576:             z = hxy[11][8] 
577:             translate(x, y, z, hxy, 0, N) 
578:             x = (hxy[12][6]+hxy[13][6]+hxy[15][6]+hxy[16][6]+hxy[18][6]+hxy[19][6])/6 
579:             y = (hxy[12][7]+hxy[13][7]+hxy[15][7]+hxy[16][7]+hxy[18][7]+hxy[19][7])/6 
580:             z = (hxy[12][8]+hxy[13][8]+hxy[15][8]+hxy[16][8]+hxy[18][8]+hxy[19][8])/6 
581:             rotateF(x, y, z, hxy, 12, 20) 
582:             rotateF(x, y, z, hxy, 27, N) 
583:  
584:             # Align B2 
585:             x = qxy[5][5] - qxy[4][5] 
586:             y = qxy[5][6] - qxy[4][6] 
587:             z = qxy[5][7] - qxy[4][7] 
588:             rotateB(x, y, z, hxy, 12, 20) 
589:             rotateB(x, y, z, hxy, 27, N) 
590:  
591:             # translation ref base to base coordinates: P atoms coincide 
592:             x = hxy[11][6]-qxy[4][5] 
593:             y = hxy[11][7]-qxy[4][6] 
594:             z = hxy[11][8]-qxy[4][7] 
595:             translate(x, y, z, hxy, 0, N) 
596:  
597:             if(i == 0): 
598:                 dist = 100 
599:  
600:             ndist = GG_CU(hxy, qxy, PL, i, Na) 
601:             if ndist < dist: 
602:                 for i in range(12, 20): 
603:                     wxy[i] = hxy[i] 
604:                 for i in range(27, N): 
605:                     wxy[i] = hxy[i] 
606:                 dist = ndist 
607:  
608:     return wxy 
609:  
610:  
611: def main(filename): 
612:     # --------------------- Main Routine ------------------------ # 
613:  
614:     Na = 0 
615:  
616:     Gxy = [] 
617:     q = [] 
618:     for l in open(filename, "r"): 
619:         if l[0:4] == "ATOM": 
620:             q.append(l) 
621:             Gxy.append(split(l)) 
622:     LG = len(Gxy) 
623:  
624:     PL = [] 
625:     for i in range(len(Gxy)): 
626:         # print i 
627:         Gxy[i][5] = float(q[i][30:38]) 
628:         Gxy[i][6] = float(q[i][38:46]) 
629:         Gxy[i][7] = float(q[i][46:54]) 
630:         Na = max(int(q[i][22:26]), Na) 
631:         if(Gxy[i][2] == 'P'): 
632:             PL.append([Gxy[i][5], Gxy[i][6], Gxy[i][7]]) 
633:  
634:     SxyA = [] 
635:     SxyG = [] 
636:     SxyC = [] 
637:     SxyU = [] 
638:     fname_list = ["baseA.lib.gz", "baseG.lib.gz", "baseC.lib.gz", "baseU.lib.gz"] 
639:     Sxy_list = [SxyA, SxyG, SxyC, SxyU] 
640:     for Sxyi, fi in zip(Sxy_list, fname_list): 
641:         for l in gzip.open(fi, 'r'): 
642:             Sxyi.append(split(l)) 
643:             Sxyi[-1][6] = float(Sxyi[-1][6]) 
644:             Sxyi[-1][7] = float(Sxyi[-1][7]) 
645:             Sxyi[-1][8] = float(Sxyi[-1][8]) 
646:  
647:     Fxy = [] 
648:     ATOM = 1 
649:     for i in range(2, Na+1): 
650:         qxy = [] 
651:         for j in range(LG): 
652:             if int(q[j][22:26]) == i: 
653:             &nError running /usr/bin/diff -r -N -U 10 r32605 r32604 bsp;   qxy.append(Gxy[j]) 
654:  
655:         if 'G' in qxy[0][3]: 
656:             wxy = rebuildG(qxy, SxyG, PL, i, Na) 
657:         elif 'A' in qxy[0][3]: 
658:             wxy = rebuildA(qxy, SxyA, PL, i, Na) 
659:         elif 'C' in qxy[0][3]: 
660:             wxy = rebuildC(qxy, SxyC, PL, i, Na) 
661:         elif 'U' in qxy[0][3]: 
662:             wxy = rebuildU(qxy, SxyU, PL, i, Na) 
663:  
664:     #    for k in range(len(wxy)): 
665:     #        wxy[k][5] = i 
666:     #        Fxy.append(wxy[k]) 
667:  
668:         for k in range(len(wxy)): 
669:             wxy[k][1] = ATOM+k 
670:             wxy[k][5] = i 
671:  
672:     #    for k in wxy: 
673:     #        print k 
674:  
675:         ATOM = ATOM+len(wxy) 
676:         toPDB(wxy) 
677:  
678:     # for i in range(len(Fxy)): 
679:     #      Fxy[i][1] = i 
680:  
681:     # toPDB(Fxy) 
682:  
683: if __name__ == "__main__": 
684:     if len(sys.argv) < 2: 
685:         fname = "6TNA_fit.pdb" 
686:     else: 
687:         fname = sys.argv[1] 
688:  
689:     main(fname)