-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathpross.py
executable file
·2285 lines (1838 loc) · 72.8 KB
/
pross.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
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/env python
"""
Calulate all torsion angles and assign secondary structure and mesostate code.
This module uses much of the code from the original BIOMOL collection of
utilities written by Raj Srinivasani with enhancements by Nick Fitzkee.
The script was put together by Pat Fleming so that a user would not need
to have the BIOMOL distribution installed to run PROSS.
Note that since Raj's time the definitions of mesostates has been superceded
by the fine grained 30 deg x 30 deg grid for most purposes. Either mesostate
grid will work for PROSS. Give your choice as an argument (see USAGE below).
Date: September 2004
Author: Pat Fleming, [email protected]
"""
USAGE = """
python PROSS.py pdbfile [oldmeso | fgmeso]
Choose old mesostate definitions or finegrained mesostate definitions.
Default is finegrained.
"""
import sys
import string, re
import copy
import gzip
import math
import types
# This script does not require Numeric
HAVE_NUMPY = 0
_RAD_TO_DEG = 180.0/math.pi
_DEG_TO_RAD = math.pi/180.0
#Residue codes to match those in a PDB file.
RESIDUES = ['ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'CYX', 'GLN', 'GLU',
'GLX', 'GLY', 'HIS', 'HIP', 'ILE', 'LEU', 'LYS', 'MET',
'PHE', 'PRO', 'SER', 'THR', 'TRP', 'TYR', 'VAL', 'ACE',
'NME', 'PCA', 'FOR', 'ASX', 'AMD']
#One-Letter residue codes.
ONE_LETTER_CODES = {
'ALA': 'A', 'ARG': 'R', 'ASN': 'N', 'ASP': 'D', 'CYS': 'C',
'GLN': 'Q', 'GLU': 'E', 'GLY': 'G', 'HIS': 'H', 'ILE': 'I',
'LEU': 'L', 'MET': 'M', 'PHE': 'F', 'PRO': 'P', 'SER': 'S',
'THR': 'T', 'TRP': 'W', 'TYR': 'Y', 'VAL': 'V', 'LYS': 'K'}
_CHI1_DEFAULT = ('N', 'CA', 'CB', 'CG')
CHI1_ATOMS = {
'ARG': _CHI1_DEFAULT,
'ASN': _CHI1_DEFAULT,
'ASP': _CHI1_DEFAULT,
'CYS': ('N', 'CA', 'CB', 'SG'),
'SER': ('N', 'CA', 'CB', 'OG'),
'GLN': _CHI1_DEFAULT,
'GLU': _CHI1_DEFAULT,
'HIS': _CHI1_DEFAULT,
'ILE': ('N', 'CA', 'CB', 'CG1'),
'LEU': _CHI1_DEFAULT,
'LYS': _CHI1_DEFAULT,
'MET': _CHI1_DEFAULT,
'PHE': _CHI1_DEFAULT,
'PRO': _CHI1_DEFAULT,
'THR': ('N', 'CA', 'CB', 'OG1'),
'TRP': _CHI1_DEFAULT,
'TYR': _CHI1_DEFAULT,
'VAL': ('N', 'CA', 'CB', 'CG1'),
'GLX': _CHI1_DEFAULT,
'ASX': _CHI1_DEFAULT,
'CYX': _CHI1_DEFAULT,
'PCA': _CHI1_DEFAULT
}
CHI2_ATOMS = {
'ARG': ('CA', 'CB', 'CG', 'CD'),
'ASN': ('CA', 'CB', 'CG', 'OD1'),
'ASP': ('CA', 'CB', 'CG', 'OD1'),
'ASX': ('CA', 'CB', 'CG', 'OD1'),
'GLN': ('CA', 'CB', 'CG', 'CD'),
'GLU': ('CA', 'CB', 'CG', 'CD'),
'GLX': ('CA', 'CB', 'CG', 'CD'),
'HIS': ('CA', 'CB', 'CG', 'ND1'),
'ILE': ('CA', 'CB', 'CG1', 'CD1'),
'LEU': ('CA', 'CB', 'CG', 'CD1'),
'LYS': ('CA', 'CB', 'CG', 'CD'),
'MET': ('CA', 'CB', 'CG', 'SD'),
'PHE': ('CA', 'CB', 'CG', 'CD1'),
'PRO': ('CA', 'CB', 'CG', 'CD'),
'TRP': ('CA', 'CB', 'CG', 'CD1'),
'TYR': ('CA', 'CB', 'CG', 'CD1')
}
CHI3_ATOMS = {
'ARG': ('CB', 'CG', 'CD', 'NE'),
'GLN': ('CB', 'CG', 'CD', 'OE1'),
'GLU': ('CB', 'CG', 'CD', 'OE1'),
'GLX': ('CB', 'CG', 'CD', 'OE1'),
'LYS': ('CB', 'CG', 'CD', 'CE'),
'MET': ('CB', 'CG', 'SD', 'CE')
}
CHI4_ATOMS = {
'ARG': ('CG', 'CD', 'NE', 'CZ'),
'LYS': ('CG', 'CD', 'CE', 'NZ')
}
# Unpack PDB Line Constants - used for simple index changes
UP_SERIAL = 0
UP_NAME = 1
UP_ALTLOC = 2
UP_RESNAME = 3
UP_CHAINID = 4
UP_RESSEQ = 5
UP_ICODE = 6
UP_X = 7
UP_Y = 8
UP_Z = 9
UP_OCCUPANCY = 10
UP_TEMPFACTOR = 11
##########################
# Start Mesostate Definition Code
##########################
# This is redundant but kept for compatibility
default = 'fgmeso'
MSDEFS = {}
# Setup Fine Grain Mesostate Bins (ala Pat Fleming)
MSDEFS['fgmeso'] = {}
mydefs = MSDEFS['fgmeso']
mydefs['RC_DICT'] = {( 180,-165): 'Aa', ( 180,-135): 'Ab', ( 180,-105): 'Ac',
( 180, -75): 'Ad', ( 180, -45): 'Ae', ( 180, -15): 'Af',
( 180, 15): 'Ag', ( 180, 45): 'Ah', ( 180, 75): 'Ai',
( 180, 105): 'Aj', ( 180, 135): 'Ak', ( 180, 165): 'Al',
(-150,-165): 'Ba', (-150,-135): 'Bb', (-150,-105): 'Bc',
(-150, -75): 'Bd', (-150, -45): 'Be', (-150, -15): 'Bf',
(-150, 15): 'Bg', (-150, 45): 'Bh', (-150, 75): 'Bi',
(-150, 105): 'Bj', (-150, 135): 'Bk', (-150, 165): 'Bl',
(-120,-165): 'Ca', (-120,-135): 'Cb', (-120,-105): 'Cc',
(-120, -75): 'Cd', (-120, -45): 'Ce', (-120, -15): 'Cf',
(-120, 15): 'Cg', (-120, 45): 'Ch', (-120, 75): 'Ci',
(-120, 105): 'Cj', (-120, 135): 'Ck', (-120, 165): 'Cl',
( -90,-165): 'Da', ( -90,-135): 'Db', ( -90,-105): 'Dc',
( -90, -75): 'Dd', ( -90, -45): 'De', ( -90, -15): 'Df',
( -90, 15): 'Dg', ( -90, 45): 'Dh', ( -90, 75): 'Di',
( -90, 105): 'Dj', ( -90, 135): 'Dk', ( -90, 165): 'Dl',
( -60,-165): 'Ea', ( -60,-135): 'Eb', ( -60,-105): 'Ec',
( -60, -75): 'Ed', ( -60, -45): 'Ee', ( -60, -15): 'Ef',
( -60, 15): 'Eg', ( -60, 45): 'Eh', ( -60, 75): 'Ei',
( -60, 105): 'Ej', ( -60, 135): 'Ek', ( -60, 165): 'El',
( -30,-165): 'Fa', ( -30,-135): 'Fb', ( -30,-105): 'Fc',
( -30, -75): 'Fd', ( -30, -45): 'Fe', ( -30, -15): 'Ff',
( -30, 15): 'Fg', ( -30, 45): 'Fh', ( -30, 75): 'Fi',
( -30, 105): 'Fj', ( -30, 135): 'Fk', ( -30, 165): 'Fl',
( 0,-165): 'Ja', ( 0,-135): 'Jb', ( 0,-105): 'Jc',
( 0, -75): 'Jd', ( 0, -45): 'Je', ( 0, -15): 'Jf',
( 0, 15): 'Jg', ( 0, 45): 'Jh', ( 0, 75): 'Ji',
( 0, 105): 'Jj', ( 0, 135): 'Jk', ( 0, 165): 'Jl',
( 30,-165): 'Ha', ( 30,-135): 'Hb', ( 30,-105): 'Hc',
( 30, -75): 'Hd', ( 30, -45): 'He', ( 30, -15): 'Hf',
( 30, 15): 'Hg', ( 30, 45): 'Hh', ( 30, 75): 'Hi',
( 30, 105): 'Hj', ( 30, 135): 'Hk', ( 30, 165): 'Hl',
( 60,-165): 'Ia', ( 60,-135): 'Ib', ( 60,-105): 'Ic',
( 60, -75): 'Id', ( 60, -45): 'Ie', ( 60, -15): 'If',
( 60, 15): 'Ig', ( 60, 45): 'Ih', ( 60, 75): 'Ii',
( 60, 105): 'Ij', ( 60, 135): 'Ik', ( 60, 165): 'Il',
( 90,-165): 'Ja', ( 90,-135): 'Jb', ( 90,-105): 'Jc',
( 90, -75): 'Jd', ( 90, -45): 'Je', ( 90, -15): 'Jf',
( 90, 15): 'Jg', ( 90, 45): 'Jh', ( 90, 75): 'Ji',
( 90, 105): 'Jj', ( 90, 135): 'Jk', ( 90, 165): 'Jl',
( 120,-165): 'Ka', ( 120,-135): 'Kb', ( 120,-105): 'Kc',
( 120, -75): 'Kd', ( 120, -45): 'Ke', ( 120, -15): 'Kf',
( 120, 15): 'Kg', ( 120, 45): 'Kh', ( 120, 75): 'Ki',
( 120, 105): 'Kj', ( 120, 135): 'Kk', ( 120, 165): 'Kl',
( 150,-165): 'La', ( 150,-135): 'Lb', ( 150,-105): 'Lc',
( 150, -75): 'Ld', ( 150, -45): 'Le', ( 150, -15): 'Lf',
( 150, 15): 'Lg', ( 150, 45): 'Lh', ( 150, 75): 'Li',
( 150, 105): 'Lj', ( 150, 135): 'Lk', ( 150, 165): 'Ll'}
# What mesostate to use when the angles are invalid (e.g. 999.99)
mydefs['INVALID'] = '??'
# What mesostate to use when the omega angle is odd (e.g. < 90.0)
mydefs['OMEGA'] = '**'
# The size of mesostate codes used in this set.
mydefs['CODE_LENGTH'] = 2
# Geometric parameters: DELTA is the size of the bins, XXX_OFF is
# the offset from 0 degrees of the centers of the bins.
mydefs['DELTA'] = 30.0
mydefs['PHI_OFF'] = 0.0
mydefs['PSI_OFF'] = 15.0
# Set up turns and turn dictionary. Dictionary contains the number of
# occurences in the PDB. This number isn't used, so for new turn types
# '1' is sufficient (as is used for PII, below)
mydefs['TURNS'] = {'EfDf': 5226, 'EeEf': 4593, 'EfEf': 4061, 'EfDg': 3883,
'EeDg': 2118, 'EeEe': 1950, 'EfCg': 1932, 'EeDf': 1785,
'EkJf': 1577, 'EkIg': 1106, 'EfEe': 995, 'EkJg': 760,
'EeCg': 553, 'DfDf': 479, 'EfCf': 395, 'DgDf': 332,
'DfDg': 330, 'IhIg': 310, 'EfDe': 309, 'EkIh': 298,
'DgCg': 275, 'DfCg': 267, 'IbDg': 266, 'DfEe': 260,
'FeEf': 250, 'IbEf': 249, 'DfEf': 219, 'IhJf': 216,
'IhJg': 213, 'IgIg': 207, 'EfCh': 188, 'DgEe': 180,
'DgEf': 176, 'EeEg': 172, 'IhIh': 153, 'EeDe': 150,
'IgJg': 147, 'EkKf': 147, 'EeCh': 147, 'IbDf': 131,
'DgDg': 128, 'EgDf': 127, 'FeDg': 114, 'ElIg': 111,
'IgIh': 107, 'DfDe': 107, 'EjIg': 101, 'EeCf': 100,
'DfCh': 94, 'DgCf': 91, 'DfCf': 91, 'DeEe': 91,
'DkIh': 88, 'FeDf': 79, 'EkIf': 78, 'EeDh': 76,
'DgCh': 74, 'IgJf': 71, 'EjJg': 71, 'FeEe': 69,
'DlIh': 66, 'EgCg': 65, 'ElIh': 62, 'EjJf': 62,
'FeCg': 59, 'DlIg': 56, 'IbCg': 54, 'EfEg': 54,
'EkJe': 53, 'FkJf': 52, 'ElJg': 51, 'DgDe': 49,
'DlJg': 46, 'EgCf': 45, 'IaEf': 40, 'FkIg': 39,
'JaEf': 38, 'EjIh': 38, 'EgEf': 38, 'DkJg': 36,
'DeEf': 34, 'EeCi': 31, 'JgIh': 29, 'IcEf': 29,
'EkKe': 29, 'DkIg': 29, 'IbEe': 27, 'EgDg': 27,
'EeFe': 27, 'EjKf': 26, 'IaDf': 25, 'HhIg': 24,
'HbDg': 24, 'ElJf': 24, 'EfDh': 24, 'IcDf': 23,
'EfBh': 23, 'IcDg': 22, 'IcCg': 22, 'FkJg': 21,
'FeCh': 21, 'IgKf': 20, 'FdDg': 20, 'EkHh': 20,
'DfDh': 20, 'DgBh': 19, 'DfBh': 19, 'DeDf': 19,
'DfFe': 18, 'EfFe': 17, 'EgEe': 16, 'EgDe': 16,
'DkJf': 16, 'JgJg': 15, 'IbEg': 15, 'IbCh': 15,
'EfBg': 15, 'DgCe': 15, 'JlEf': 14, 'CgCg': 14,
'HhJf': 13, 'EeBi': 13, 'DfBi': 13, 'IhIf': 12,
'FeEg': 12, 'FdEf': 12, 'EdEf': 12, 'DlJf': 12,
'DhCg': 12, 'JgIg': 11, 'IeBg': 11, 'FjIg': 11,
'FdCh': 11, 'EdEe': 11, 'JfIh': 10, 'JaEe': 10,
'HhJg': 10, 'HbEf': 10, 'HbCh': 10, 'FkIh': 10,
'FjJf': 10, 'ElJe': 10, 'DhDf': 10, 'CgDf': 10}
# Set up the PII defitions, similar to dictionary above.
mydefs['PII'] = {'Dk':1, 'Dl':1, 'Ek':1, 'El':1}
# Set up the codes that define helix and strand. Here, rather than storing
# the dictionary like we did above, we'll store the regular expression
# matcher directly. This prevents us from recompiling it every time we
# want to find a helix or strand.
helix = ('De', 'Df', 'Ed', 'Ee', 'Ef', 'Fd', 'Fe')
strand = ('Bj', 'Bk', 'Bl', 'Cj', 'Ck', 'Cl', 'Dj', 'Dk', 'Dl')
pat_helix = "(%s){5,}" % string.join(map(lambda x: "(%s)" % x, helix), '|')
pat_strand = "(%s){3,}" % string.join(map(lambda x: "(%s)" % x, strand), '|')
mydefs['HELIX'] = re.compile(pat_helix)
mydefs['STRAND'] = re.compile(pat_strand)
###########################
# Setup Old Mesostate Bins (ala Raj Srinivasan)
##########################
MSDEFS['oldmeso'] = {}
mydefs = MSDEFS['oldmeso']
mydefs['RC_DICT'] = {( 180, 180): 'A', ( 180,-120): 'B', ( 180, -60): 'C',
( 180, 0): 'D', ( 180, 60): 'E', ( 180, 120): 'F',
(-120, 180): 'G', (-120,-120): 'H', (-120, -60): 'I',
(-120, 0): 'J', (-120, 60): 'K', (-120, 120): 'L',
( -60, 180): 'M', ( -60,-120): 'N', ( -60, -60): 'O',
( -60, 0): 'P', ( -60, 60): 'Q', ( -60, 120): 'R',
( 0, 180): 'S', ( 0,-120): 'T', ( 0, -60): 'U',
( 0, 0): 'V', ( 0, 60): 'W', ( 0, 120): 'X',
( 60, 180): 'm', ( 60,-120): 'r', ( 60, -60): 'q',
( 60, 0): 'p', ( 60, 60): 'o', ( 60, 120): 'n',
( 120, 180): 'g', ( 120,-120): 'l', ( 120, -60): 'k',
( 120, 0): 'j', ( 120, 60): 'i', ( 120, 120): 'h'}
# What mesostate to use when the angles are invalid (e.g. 999.99)
mydefs['INVALID'] = '?'
# What mesostate to use when the omega angle is odd (e.g. < 90.0)
mydefs['OMEGA'] = '*'
# The size of mesostate codes used in this set.
mydefs['CODE_LENGTH'] = 1
# Geometric parameters: DELTA is the size of the bins, XXX_OFF is
# the offset from 0 degrees of the centers of the bins.
mydefs['DELTA'] = 60.0
mydefs['PHI_OFF'] = 0.0
mydefs['PSI_OFF'] = 0.0
# Set up turns and turn dictionary. Dictionary contains the type
# of the turn (no primes)
mydefs['TURNS'] = {'OO': 1, 'OP': 1, 'OJ': 1, 'PO': 1, 'PP': 1, 'PJ': 1,
'JO': 1, 'JP': 1, 'JJ': 1,
'Mo': 2, 'Mp': 2, 'Mj': 2, 'Ro': 2, 'Rp': 2, 'Rj': 2,
'oo': 3, 'op': 3, 'oj': 3, 'po': 3, 'pp': 3, 'pj': 3,
'jo': 3, 'jp': 3, 'jj': 3,
'mO': 4, 'mP': 4, 'mJ': 4, 'rO': 4, 'rP': 4, 'rJ': 4}
# Set up the PII defitions, similar to dictionary above.
mydefs['PII'] = {'M':1, 'R':1}
# Set up the codes that define helix and strand. Here, rather than storing
# the dictionary like we did above, we'll store the regular expression
# matcher directly. This prevents us from recompiling it every time we
# want to find a helix or strand.
helix = ('O', 'P')
strand = ('L', 'G', 'F', 'A')
pat_helix = "(%s){5,}" % string.join(map(lambda x: "(%s)" % x, helix), '|')
pat_strand = "(%s){3,}" % string.join(map(lambda x: "(%s)" % x, strand), '|')
mydefs['HELIX'] = re.compile(pat_helix)
mydefs['STRAND'] = re.compile(pat_strand)
##########################
# End Mesostate Definition Code
##########################
def res_rc(r1, r2, r3=180, mcodes=None):
"""res_rc(r1, r2, r3) - get mesostate code for a residue
Given a phi (r1), psi (r2), and omega (r3) torsion angle, calculate
the mesostate code that describes that residue.
A mesostate will be returned in
all but two cases: if omega deviates from planarity by more than 90
degrees, '*' is returned. Also, if any torsions are greater than
180.0 (biomol assignes 999.0 degrees to angles that that are
indeterminate), prosst.INVALID is returned. Here, r3 defaults to
180.0.
"""
if not mcodes: mcodes = default
ms = MSDEFS[mcodes]
OMEGA = ms['OMEGA']
INVALID = ms['INVALID']
PHI_OFF = ms['PHI_OFF']
PSI_OFF = ms['PSI_OFF']
DELTA = ms['DELTA']
RC_DICT = ms['RC_DICT']
if (abs(r3) <= 90.0):
return OMEGA
elif r1>180.0 or r2>180.0 or r3>180.0:
return INVALID
ir1 = -int(PHI_OFF) + int(round((r1+PHI_OFF)/DELTA )) * int(DELTA)
ir2 = -int(PSI_OFF) + int(round((r2+PSI_OFF)/DELTA )) * int(DELTA)
while ir1 <= -180: ir1 = ir1 + 360
while ir1 > 180: ir1 = ir1 - 360
while ir2 <= -180: ir2 = ir2 + 360
while ir2 > 180: ir2 = ir2 - 360
return RC_DICT[(ir1,ir2)]
def rc_codes(chain, phi=None, psi=None, ome=None, mcodes=None):
"""rc_codes(chain, phi, psi, ome) - return rotamer codes
Given a protein chain (and optionally phi, psi, omega), this
function will return a list of mesostate codes that
applies to the chain, as determined by res_rc.
"""
if not mcodes: mcodes = default
n = range(len(chain))
if phi is None: phi = map(chain.phi, n)
if psi is None: psi = map(chain.psi, n)
if ome is None: ome = map(chain.omega, n)
return map(lambda x, y, z: res_rc(x, y, z, mcodes), phi, psi, ome)
def rc_ss(chain, phi=None, psi=None, ome=None, mcodes=None):
"""rc_ss(chain, phi, psi, ome) - calculate secondary structure
This function calculates the secondary structure using the PROSS method
with rotamer codes. Given a chain, and optionally a list of phi,
psi, and omega, it calculates the backbone secondary structure of
the chain. The return value is (phi, psi, ome, sst), where
phi, psi, and ome are calculated if not specified, and sst is the
secondary structure codes: H = helix, E = strand, P = PII, C = coil.
"""
if not mcodes: mcodes = default
ms = MSDEFS[mcodes]
PII = ms['PII']
TURNS = ms['TURNS']
HELIX = ms['HELIX']
STRAND = ms['STRAND']
if phi is None:
chain.gaps()
nres = len(chain)
phi = map(chain.phi, xrange(nres))
else:
nres = len(chain)
if psi is None: psi = map(chain.psi, xrange(nres))
if ome is None: ome = map(chain.omega, xrange(nres))
codes = rc_codes(chain, phi, psi, ome, mcodes)
chain.gaps()
sst = ['C']*nres
is_PII = PII.has_key
for i in xrange(nres-1):
code = codes[i]
if is_PII(code):
sst[i] = 'P'
is_turn = TURNS.has_key
for i in xrange(nres-1):
code = codes[i] + codes[i+1]
if is_turn(code):
sst[i] = sst[i+1] = 'T'
helices = _rc_find(codes, HELIX, mcodes)
strands = _rc_find(codes, STRAND, mcodes)
for helix in helices:
i, j = helix
for k in range(i, j):
sst[k] = 'H'
print strands
for strand in strands:
i, j = strand
for k in range(i, j):
if sst[k] in ('C', 'P'): sst[k] = 'E'
# if sst[k] == 'C': sst[k] = 'E'
return phi, psi, ome, sst
def _rc_find(codes, pattern, mcodes=None):
"""_rc_find(codes, pat_obj) - find a endpoints of a regexp
Given a list of mesostate codes, this function identifies a endpoints
of a match <pattern>. <pat_obj> is a compiled regular expression
pattern whose matches will be returned as pairs indicated start,
end in <codes>
"""
if not mcodes: mcodes = default
CODE_LENGTH = MSDEFS[mcodes]['CODE_LENGTH']
if not type(codes) == type(''):
codes = string.join(codes, '')
matches = []
it = pattern.finditer(codes)
try:
while 1:
mat = it.next()
matches.append((mat.start()/CODE_LENGTH, mat.end()/CODE_LENGTH))
except StopIteration:
pass
return matches
##############################
# Protein objects and methods
##############################
def is_atom(o):
"""Determine if an object in a TypedList is an atom."""
return hasattr(o, '_is_an_atom3d')
def is_residue(o):
"""Determine if an object in a TypedList is a residue."""
return hasattr(o, '_is_a_residue')
def is_chain(o):
"""Determine if an object in a TypedList is a chain."""
return hasattr(o, '_is_a_chain')
def is_mol(o):
"""Determine if an object in a TypedList is a molecule."""
return hasattr(o, '_is_a_mol')
####################
HAVE_POP = hasattr([], 'pop')
HAVE_EXTEND = hasattr([], 'extend')
class TypedList:
"""A Python list restricted to having objects of the same type.
An instance of a TypedList is created as follows:
mylist = TypedList(function, [elements])
where function is a python function which takes an argument
and returns 1 or 0 indicating whether the object represented
by the argument is of the correct type, and elements is an optional
list of elements to be added into the instance. Here is a
full blown example.
def is_int(o):
return type(o) == type(0)
mylist = TypedList(is_int, [0, 1, 2, 3])
New elements are added to the list as follows:
mylist.append(25)
Instances of TypedList support all operations available for
Python Lists (as of Python version 1.5.2a2)
"""
_is_a_typed_list = 1
def __init__(self, function, elements=None):
self._func = function
self.elements = []
if not elements is None:
if self._func(elements):
self.elements.append(elements)
else:
for el in elements:
self.append(el)
def append(self, el):
"""Add an element to the TypedList."""
if self._func(el):
self.elements.append(el)
else:
raise TypeError, 'Element to be added to list has incorrect type.'
def __len__(self):
"""Returns the length of the TypedList."""
return len(self.elements)
def __repr__(self):
"""Returns a printable representation of the TypedList."""
return "%s(%s, %s)" % (self.__class__, self._func.__name__,
`self.elements`)
def __str__(self):
"""Returns all elements of the TypedList."""
return `self.elements`
def __getitem__(self, i):
"""Returns an item from the TypedList located at index i."""
return self.elements[i]
def __setitem__(self, i, v):
"""Sets item at index i to the value of v."""
if self._func(v):
self.elements[i] = v
else:
raise TypeError, 'Item not of correct type in __setitem__'
def __delitem__(self, i):
"""Removes an item from the TypedList."""
del self.elements[i]
def __getslice__(self, i, j):
"""Returns a slice of the typed list between indexes i and j."""
new = self.__class__(self._func)
new.elements = self.elements[i:j]
return new
def __setslice__(self, i, j, v):
"""Sets the value of a slice of items between i and j to v."""
if self._alltrue(v):
self.elements[i:j] = v
def __delslice__(self, i, j):
"""Deleted a slice of items located between indices i and j."""
del self.elements[i:j]
def __add__(self, other):
"""Adds a list to a TypedList, as long as it's of the same type."""
if not hasattr(other, '_is_a_typed_list'):
raise TypeError,'List to be concatenated not instance of %s' %\
self.__class__
if self._func <> other._func:
raise TypeError, 'Lists to be added not of same type'
new = self.__class__(self._func)
new.elements = self.elements + other.elements
return new
def __mul__(self, other):
"""Multiply set of elements and another, and store in new.elements."""
if type(other) == type(0):
new = self.__class__(self._func)
new.elements = self.elements * other
return new
else:
raise TypeError, "can't multiply list with non-int"
__rmul__ = __mul__
def __copy__(self):
"""Returns a copy of a TypedList."""
new = self.__class__(self._func)
for el in self.elements:
new.elements.append(el.__copy__())
have = new.__dict__.has_key
for key in self.__dict__.keys():
if not have(key):
new.__dict__[key] = copy.deepcopy(self.__dict__[key])
return new
__deepcopy__ = clone = __copy__
def _alltrue(self, els):
"""Determine if all TypedList functions are set."""
return len(els) == len(filter(None, map(self._func, els)))
def sort(self):
"""Sort the elements in a TypedList."""
self.elements.sort()
def reverse(self):
"""Reverse the order of elements in a TypedList."""
self.elements.reverse()
def count(self, el):
"""Returns the number of elements in a TypedList."""
return self.elements.count(el)
def extend(self, els):
"""Add an attribute to an element in a TypedList."""
if self._alltrue(els):
if HAVE_EXTEND:
self.elements.extend(els)
else:
for el in els:
self.elements.append(el)
else:
raise TypeError, 'One or more elements of list not of correct type'
def pop(self):
"""Remove the last element in a TypedList."""
if HAVE_POP:
return self.elements.pop()
else:
el = self.elements[-1]
del self.elements[-1]
return el
def index(self, el):
"""Return an index if the items in a TypedList."""
return self.elements.index(el)
def remove(self, el):
"""Remove an element from a TypedList."""
self.elements.remove(el)
def insert(self, pos, el):
"""Insert an element in a TypedList at position pos."""
if self._func(el):
self.elements.insert(pos, el)
else:
raise TypeError, 'Item not of correct type in insert'
def indices(self, len=len):
"""Returns the length of a range of indices."""
return xrange(len(self.elements))
def reverse_indices(self, len=len):
"""Reverses indices of elements in a TypedList."""
return xrange(len(self.elements)-1, -1, -1)
####################
class molResidue(TypedList):
"""A TypedList containing only objects of the residue
type, set with `is_a_residue`
"""
_is_a_residue = 1
def __init__(self, name='', atoms=None, **kw):
TypedList.__init__(self, is_atom, atoms)
self.name = name
for key, value in kw.items():
setattr(self, key, value)
def num_atoms(self):
"""returns the number of atoms in residue"""
return len(self.elements)
def has_atom(self, name):
"""returns true if residue has an atom named 'name'"""
for atom in self.elements:
if atom.name == name:
return 1
return 0
def atoms(self):
"""Returns a list of atoms within the molResidue list."""
return self.elements
def atoms_with_name(self, *names):
"""returns atoms in residue with specified names"""
ret = []
Append = ret.append
for name in names:
for atom in self.elements:
if atom.name == name:
Append(atom)
break
return ret
def delete_atoms_with_name(self, *names):
"""delete atoms in residue with specified names"""
els = self.elements
for i in self.reverse_indices():
atom = els[i]
if atom.name in names:
del els[i]
def atoms_not_with_name(self, *names):
"""returns atoms in residue excluding specified names"""
ret = []
for atom in self.elements:
if not atom.name in names:
ret.append(atom)
return ret
def atom_coordinates(self, *names):
"""returns coordinates of named atoms.
If names is omitted all atom coordinates are returned."""
if len(names)==0:
atoms = self.elements
else:
atoms = apply(self.atoms_with_name, names)
na = len(atoms)
if na == 0: return
if HAVE_NUMPY:
a = Numeric.zeros((na, 3), 'd')
else:
a = [None]*na
pos = 0
for atom in atoms:
a[pos] = atom.coords()
pos = pos + 1
return a
def assign_radii(self):
raise AttributeError, 'Should be defined in specialized class'
def type(self):
return 'residue'
class molChain(TypedList):
"""A TypedList containing only items of the chain
type, set with `is_a_chain`
"""
_is_a_chain = 1
def __init__(self, name='', residues=None, **kw):
self.name = name
TypedList.__init__(self, is_residue, residues)
for key, value in kw.items():
setattr(self, key, value)
def num_residues(self):
"""Returns the number of residues within a chain."""
return len(self)
def num_atoms(self):
"""Returns the number of atoms in each residue within a chain."""
na = 0
for res in self.elements:
na = na + len(res.elements)
return na
def atoms(self):
"""Returns a list of atoms within each residue in a chain."""
ret = []
Append = ret.append
for res in self.elements:
for atom in res.elements:
Append(atom)
return ret
def residues_with_name(self, *names):
"""returns specified residues as a python list"""
if len(names) == 0:
return
l = []
for res in self.elements:
if res.name in names:
l.append(res)
return l
def delete_residues_with_name(self, *names):
"""delete specified residues from Chain"""
if len(names) == 0:
return
els = self.elements
for i in self.reverse_indices():
if els[i].name in names:
del els[i]
def residues_not_with_name(self, *names):
"""returns residues excluding specified names as a python list"""
ret = []
for res in self.elements:
if not res.name in names:
ret.append(res)
return ret
def atoms_with_name(self, *names):
"""Returns a list of atoms that match the specified names."""
ret = []
Append = ret.append
if len(names) > 0:
for res in self.elements:
for name in names:
for atom in res.elements:
if atom.name == name:
Append(atom)
break
else:
for res in self.elements:
for atom in res.elements:
Append(atom)
return ret
def delete_atoms_with_name(self, *names):
"""delete atoms in residue with specified names"""
for res in self.elements:
apply(res.delete_atoms_with_name, names)
def atoms_not_with_name(self, *names):
"""returns atoms in residue excluding specified names"""
ret = []
for res in self.elements:
ret[len(ret):] = apply(res.atoms_not_with_name, names)
return ret
def atom_coordinates(self, *names):
"""returns coordinates of named atoms. if names is None
all atom coordinates are returned."""
coords = []
if len(names) > 0:
for res in self.elements:
for atom in res.elements:
if atom.name in names:
coords.append(atom.coords())
else:
atoms = apply(self.atoms_with_name, names)
coords = map(lambda a:a.coords(), atoms)
if HAVE_NUMPY:
return Numeric.array(coords)
else:
return coords
def atoms(self):
"""returns a python list of atoms in a residue."""
atoms = []
for res in self.elements:
for atom in res: atoms.append(atom)
return atoms
def delete_alt_locs(self):
"""delete_alt_locs - remove secondary conformations in the chain
In a chain with multiple occupancy and alternate location identifiers,
it is often desirable to eliminate the secondary conformations for
use in simulation, etc. This function (abitrarily) finds and selects
the first given conformation and deletes all other conformations.
"""
AtomCount = self.present
chain = self.elements
delete = []
for i in xrange(len(chain)):
residue = chain[i]
rid = (residue.idx, residue.icode)
for j in xrange(len(residue)):
atom = residue[j]
anam = atom.name
try:
acnt = AtomCount[rid][anam]
except KeyError:
print "Unable to locate %s %s %s in present dictionary."%\
(rid[0], rid[1], anam)
return
if acnt == 1:
continue
if acnt < 1:
AtomCount[rid][anam] = acnt + 1
delete.append((i, j))
continue
atom.alt = ' '
AtomCount[rid][anam] = -acnt + 2
delete.reverse()
for r, a in delete:
del chain[r][a]
def assign_radii(self):
"""Assigns the radius of residues in a molChain"""
for res in self:
res.assign_radii()
def preserve_chain_hetero(self):
"""prevent hetero residues from being deleted as hetero atoms
Normally, delete_hetero will delete all hetero atoms from a
molecule. This includes waters and heterogroups (hemes, etc.),
but it also includes hetero residues--nonstandard residues
that have backbone connectivity but perhaps extra atoms (e.g.
S-hydroxy-cysteine). Deleting these residues may disrupt an
otherwise continuous chain and may be undesirable.
Given that a chain has a valid SEQRES entry, this function will
'unset' the hetero flag for heterogroups that are involved in
the sequence itself. When delete_hetero is run, these atoms
will be preserved.
"""
for i in xrange(self.num_residues()):
if hasattr(self[i], 'chain_het') and hasattr(self[i], 'het'):
delattr(self[i], 'het')
def delete_hetero(self):
"""delete all hetero atoms from a protein
This function removes all HETATM records from the molecular
structure. This include waters, heme groups, as well as residues
with nonstandard structures
"""
for i in xrange(self.num_residues()-1, -1, -1):
if hasattr(self[i], 'het'):
del self[i]
def delete_waters(self):
"""Delete water molecules from residues in the molChain."""
for i in xrange(self.num_residues()-1, -1, -1):
if self[i].name == 'HOH':
del self[i]
def translate(self, dx, dy=None, dz=None):
for res in self.elements:
res.translate(dx, dy, dz)
def rotatex(self, theta):
"""Rotates residues in the molChain by the X-axis."""
for res in self.elements:
res.rotatex(theta)
def rotatey(self, theta):
"""Rotates residues in the molChain by the Y-axis."""
res.rotatey(theta)
def rotatez(self, theta):
"""Rotates residues in the molChain by the Z-axis."""
for res in self.elements:
res.rotatez(theta)
def type(self):
"""Returns the type of the molChain. Should always be `Chain`"""
return 'Chain'
class molMol(TypedList):
"""A typed list containing only molecule objects,
as set by `is_a_mol
"""
_is_a_mol = 1
def __init__(self, name='', chains=None):
self.name = name