/
arpeggio.py
3028 lines (2314 loc) · 112 KB
/
arpeggio.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/python
##################################################################################################################################
# #
# CREDO-BASED DEFINITIONS: #
# #
# https://bitbucket.org/harryjubb/credovi/src/ac69222542134bb86d26a24561e4f740466e2b0e/credovi/config/credo.json?at=default #
# https://bitbucket.org/harryjubb/credovi/src/ac69222542134bb86d26a24561e4f740466e2b0e/credovi/structbio/structure.py?at=default #
# #
##################################################################################################################################
###########
# IMPORTS #
###########
import argparse
import collections
import logging
#import math
import operator
try:
import resource
except ImportError:
logging.info('Resource module not available, resource usage info won\'t be logged.')
import sys
from collections import OrderedDict
import numpy as np
#from Bio.PDB import PDBIO
from Bio.PDB.PDBParser import PDBParser
from Bio.PDB import NeighborSearch
from Bio.PDB.Atom import Atom
from Bio.PDB.Residue import Residue
from Bio.PDB.Polypeptide import PPBuilder
import openbabel as ob
#############
# CONSTANTS #
#############
from config import ATOM_TYPES, CONTACT_TYPES, VDW_RADII, METALS, \
HALOGENS, CONTACT_TYPES_DIST_MAX, FEATURE_SIFT, VALENCE, \
MAINCHAIN_ATOMS, THETA_REQUIRED, STD_RES, PROT_ATOM_TYPES, \
AMIDE_SMARTS, COMMON_SOLVENTS, STANDARD_NUCLEOTIDES
###########
# CLASSES #
###########
class HydrogenError(Exception):
def __init__(self):
logging.error('Please remove all hydrogens from the structure then re-run.')
class OBBioMatchError(Exception):
def __init__(self, serial=''):
if not serial:
logging.error('An OpenBabel atom could not be matched to a BioPython counterpart.')
else:
logging.error('OpenBabel OBAtom with PDB serial number {} could not be matched to a BioPython counterpart.'.format(serial))
class AtomSerialError(Exception):
def __init__(self):
logging.error('One or more atom serial numbers are duplicated.')
class SiftMatchError(Exception):
def __init__(self):
logging.error('Seeing is not believing.')
class SelectionError(Exception):
def __init__(self, selection):
logging.error('Invalid selector: {}'.format(selection))
#############
# FUNCTIONS #
#############
def int2(x):
'''
Return integer from base 2 number.
Can accept a list/tuple of 0s and 1s.
'''
if isinstance(x, collections.Iterable):
x = ''.join([str(k) for k in x])
else:
x = str(x)
return int(x, 2)
def int3(x):
'''
Return integer from base 3 number.
Can accept a list/tuple of 0s, 1s and 2s.
'''
if isinstance(x, collections.Iterable):
x = ''.join([str(k) for k in x])
else:
x = str(x)
return int(x, 3)
def selection_parser(selection_list, atom_list):
'''
Selection syntax:
/<chain_id>/<res_num>[<ins_code>]/<atom_name>
HET:<het_id>
Other formats will be rejected, for now.
You can omit fields as long as the number of `/` is correct.
Selections are additive, so if you select chain A with /A//,
adding /A/91/C23 won't make any difference.
'''
final_atom_list = set([])
for selection in selection_list:
#selection_dict = {
# 'chain': None,
# 'residue_number': None,
# 'atom_name': None
#}
chain = None
residue_number = None
insertion_code = ' '
#residue_range = None # TODO
atom_name = None
current_atom_list = atom_list[:]
original_selection = selection
selection = selection.strip()
if selection.startswith('RESNAME:'):
selection = selection.replace('RESNAME:', '').strip()
# RESNAMES ARE MAX LENGTH 3
if len(selection) > 3:
raise SelectionError(original_selection)
current_atom_list = [x for x in current_atom_list if x.get_parent().resname.strip() == selection]
for selected_atom in current_atom_list:
final_atom_list.add(selected_atom)
# SELECT ALL ORGANIC SMALL-MOLECULE LIGANDS
elif selection.startswith('LIGANDS'):
current_atom_list = [x for x in current_atom_list if
x.get_parent().is_polypeptide == False and # MUST NOT BE POLYPEPTIDE
len(x.get_parent().child_list) >= 5 and # MIN NUMBER OF ATOMS
len(x.get_parent().child_list) <= 100 and # MAX NUMBER OF ATOMS
'C' in set([y.element for y in x.get_parent().child_list]) and # MUST CONTAIN CARBON
x.get_parent().resname.strip().upper() not in COMMON_SOLVENTS and # MUST NOT BE COMMON SOLVENT
x.get_parent().resname.strip().upper() not in STANDARD_NUCLEOTIDES and # MUST NOT BE NUCLEOTIDE
not x.get_parent().resname.startswith('+') # MUST NOT BE MODIFIED NUCLEOTIDE
]
for selected_atom in current_atom_list:
final_atom_list.add(selected_atom)
elif selection.startswith('/'):
selection = selection.lstrip('/').split('/')
#print selection
if len(selection) != 3:
raise SelectionError(original_selection)
# CHAIN
if selection[0]:
chain = selection[0]
# RESIDUE AND INS CODE
if selection[1]:
if selection[1].isdigit():
# JUST THE RESNUM
residue_number = int(selection[1])
elif selection[1].isalnum():
# CHECK FOR VALID RESNUM+INSCODE
if selection[1][-1].isalpha() and selection[1][:-1].isdigit():
residue_number = int(selection[1][:-1])
insertion_code = selection[1][-1]
else:
raise SelectionError(original_selection)
else:
raise SelectionError(original_selection)
# ATOM NAME
if selection[2]:
if not selection[2].isalnum() and "'" not in selection[2]:
raise SelectionError(original_selection)
else:
atom_name = selection[2]
# NOW MAKE THE SELECTION WITH BIOPYTHON
if chain:
current_atom_list = [x for x in current_atom_list if x.get_parent().get_parent().id == chain]
if residue_number:
current_atom_list = [x for x in current_atom_list if x.get_parent().id[1] == residue_number and x.get_parent().id[2] == insertion_code]
if atom_name:
current_atom_list = [x for x in current_atom_list if x.name == atom_name]
for selected_atom in current_atom_list:
final_atom_list.add(selected_atom)
else:
raise SelectionError(original_selection)
final_atom_list = list(final_atom_list)
if len(final_atom_list) == 0:
logging.error('Selection was empty.')
sys.exit(1)
return list(final_atom_list)
def make_pymol_string(entity):
'''
Feed me a BioPython atom or BioPython residue.
See `http://pymol.sourceforge.net/newman/user/S0220commands.html`.
chain-identifier/resi-identifier/name-identifier
chain-identifier/resi-identifier/
'''
if isinstance(entity, Atom):
chain = entity.get_parent().get_parent()
residue = entity.get_parent()
atom_name = entity.name
elif isinstance(entity, Residue):
chain = entity.get_parent()
residue = entity
atom_name = ''
else:
raise TypeError('Cannot make a PyMOL string from a non-Atom or Residue object.')
res_num = residue.id[1]
# ADD INSERTION CODE IF NEED BE
if residue.id[2] != ' ':
res_num = str(res_num) + residue.id[2]
macro = '{}/{}/{}'.format(chain.id,
res_num,
atom_name)
return macro
def get_single_bond_neighbour(ob_atom):
'''
'''
for bond in ob.OBAtomBondIter(ob_atom):
if not bond.IsSingle():
continue
current_neighbour = bond.GetNbrAtom(ob_atom)
if current_neighbour.IsHydrogen():
continue
return current_neighbour
return None
def max_mem_usage():
'''
Returns maximum memory usage of the program thus far, in megabytes, as a string.
'''
try:
return str(resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1000.0) + ' MB'
except Exception as err:
logging.warn('Resource usage information not available ().'.format(str(err)))
def get_angle(point_a, point_b, point_c):
'''
Get the angle between three points in 3D space.
Points should be supplied in Numpy array format.
http://stackoverflow.com/questions/19729831/angle-between-3-points-in-3d-space
'''
#In pseudo-code, the vector BA (call it v1) is:
#v1 = {A.x - B.x, A.y - B.y, A.z - B.z}
v1 = point_a - point_b
#Similarly the vector BC (call it v2) is:
#v2 = {C.x - B.x, C.y - B.y, C.z - B.z}
v2 = point_c - point_b
#The dot product of v1 and v2 is a function of the cosine of the angle between them
#(it's scaled by the product of their magnitudes). So first normalize v1 and v2:
#v1mag = sqrt(v1.x * v1.x + v1.y * v1.y + v1.z * v1.z)
v1_mag = np.sqrt(v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2])
#v1norm = {v1.x / v1mag, v1.y / v1mag, v1.z / v1mag}
v1_norm = np.array([v1[0] / v1_mag, v1[1] / v1_mag, v1[2] / v1_mag])
#v2mag = sqrt(v2.x * v2.x + v2.y * v2.y + v2.z * v2.z)
v2_mag = np.sqrt(v2[0] * v2[0] + v2[1] * v2[1] + v2[2] * v2[2])
#v2norm = {v2.x / v2mag, v2.y / v2mag, v2.z / v2mag}
v2_norm = np.array([v2[0] / v2_mag, v2[1] / v2_mag, v2[2] / v2_mag])
#Then calculate the dot product:
#res = v1norm.x * v2norm.x + v1norm.y * v2norm.y + v1norm.z * v2norm.z
res = v1_norm[0] * v2_norm[0] + v1_norm[1] * v2_norm[1] + v1_norm[2] * v2_norm[2]
#And finally, recover the angle:
angle = np.arccos(res)
if np.isnan(angle):
logging.warn('Angle for <{}, {}, {}> was NaN, setting to Pi.'.format(point_a, point_b, point_c))
angle = np.pi
return angle
def group_angle(group, point_coords, degrees=False, signed=False):
'''
Adapted from CREDO: `https://bitbucket.org/blundell/credovi/src/bc337b9191518e10009002e3e6cb44819149980a/credovi/structbio/aromaticring.py?at=default`
`group` should be a dict with Numpy array 'center' and 'normal' attributes.
`point_coords` should be a Numpy array.
'''
cosangle = np.dot(group['normal'], point_coords) / (np.linalg.norm(group['normal']) * np.linalg.norm(point_coords))
# GET THE ANGLE AS RADIANS
rad = np.arccos(cosangle)
if not degrees: return rad
# CONVERT RADIANS INTO DEGREES
else:
# CONVERT INTO A SIGNED ANGLE
if signed: rad = rad -np.pi if rad > np.pi / 2 else rad
# RETURN DEGREES
return rad * 180 / np.pi
def group_group_angle(group, group2, degrees=False, signed=False):
'''
Adapted from CREDO: `https://bitbucket.org/blundell/credovi/src/bc337b9191518e10009002e3e6cb44819149980a/credovi/structbio/aromaticring.py?at=default`
`group` and `group2` should be a dict with Numpy array 'center' and 'normal' attributes.
'''
cosangle = np.dot(group['normal'], group2['normal']) / (np.linalg.norm(group['normal']) * np.linalg.norm(group2['normal']))
# GET THE ANGLE AS RADIANS
rad = np.arccos(cosangle)
if not degrees: return rad
# CONVERT RADIANS INTO DEGREES
else:
# CONVERT INTO A SIGNED ANGLE
if signed: rad = rad -np.pi if rad > np.pi / 2 else rad
# RETURN DEGREES
return rad * 180 / np.pi
## GOLDEN SECTION SPIRAL
## THANK YOU BOSCOH!
## http://boscoh.com/protein/calculating-the-solvent-accessible-surface-area-asa.html
## http://boscoh.com/protein/asapy.html
#def points_on_sphere(n):
# pts = np.empty((int(n), 3))
# n = float(n)
#
# inc = math.pi * (3 - math.sqrt(5))
# off = 2 / n
#
# for k in range(0, int(n)):
#
# y = k * off - 1 + (off / 2)
# r = math.sqrt(1 - y*y)
# phi = k * inc
# pts[k][0] = math.cos(phi) * r
# pts[k][1] = y
# pts[k][2] = math.sin(phi) * r
#
# return pts
# CONTACT FUNCTIONS
def is_hbond(donor, acceptor):
'''
Feed me BioPython atoms.
'''
for hydrogen_coord in donor.h_coords:
h_dist = np.linalg.norm(hydrogen_coord - acceptor.coord)
if h_dist <= VDW_RADII['H'] + acceptor.vdw_radius + VDW_COMP_FACTOR:
if get_angle(donor.coord, hydrogen_coord, acceptor.coord) >= CONTACT_TYPES['hbond']['angle rad']:
return 1
return 0
def is_weak_hbond(donor, acceptor):
'''
Feed me BioPython atoms.
'''
for hydrogen_coord in donor.h_coords:
h_dist = np.linalg.norm(hydrogen_coord - acceptor.coord)
if h_dist <= VDW_RADII['H'] + acceptor.vdw_radius + VDW_COMP_FACTOR:
if get_angle(donor.coord, hydrogen_coord, acceptor.coord) >= CONTACT_TYPES['weak hbond']['angle rad']:
return 1
return 0
def is_halogen_weak_hbond(donor, halogen, ob_mol):
'''
Feed me BioPython atoms and the OpenBabel molecule.
'''
# `nbr` WILL BE A BIOPYTHON ATOM
# ... HOPEFULLY
nbr = ob_to_bio[get_single_bond_neighbour(ob_mol.GetAtomById(bio_to_ob[halogen])).GetId()]
for hydrogen_coord in donor.h_coords:
h_dist = np.linalg.norm(halogen.coord - hydrogen_coord)
if h_dist <= VDW_RADII['H'] + halogen.vdw_radius + VDW_COMP_FACTOR:
if CONTACT_TYPES['weak hbond']['cx angle min rad'] <= get_angle(nbr.coord, halogen.coord, hydrogen_coord) <= CONTACT_TYPES['weak hbond']['cx angle max rad']:
return 1
return 0
def is_xbond(donor, acceptor, ob_mol):
'''
Feed me BioPython atoms and the OpenBabel molecule.
'''
# `nbr` WILL BE A BIOPYTHON ATOM
# ... HOPEFULLY
nbr = ob_to_bio[get_single_bond_neighbour(ob_mol.GetAtomById(bio_to_ob[donor])).GetId()]
theta = get_angle(nbr.coord, donor.coord, acceptor.coord)
if (theta >= CONTACT_TYPES['xbond']['angle theta 1 rad']):
return 1
return 0
def update_atom_sift(atom, addition, contact_type='INTER'):
'''
'''
atom.sift = [x or y for x, y in zip(atom.sift, addition)]
if contact_type == 'INTER':
atom.sift_inter_only = [x or y for x, y in zip(atom.sift_inter_only, addition)]
if 'INTRA' in contact_type:
atom.sift_intra_only = [x or y for x, y in zip(atom.sift_intra_only, addition)]
if 'WATER' in contact_type:
atom.sift_water_only = [x or y for x, y in zip(atom.sift_water_only, addition)]
def update_atom_fsift(atom, addition, contact_type='INTER'):
'''
'''
atom.actual_fsift = [x or y for x, y in zip(atom.actual_fsift, addition)]
if contact_type == 'INTER':
atom.actual_fsift_inter_only = [x or y for x, y in zip(atom.actual_fsift_inter_only, addition)]
if 'INTRA' in contact_type:
atom.actual_fsift_intra_only = [x or y for x, y in zip(atom.actual_fsift_intra_only, addition)]
if 'WATER' in contact_type:
atom.actual_fsift_water_only = [x or y for x, y in zip(atom.actual_fsift_water_only, addition)]
def update_atom_integer_sift(atom, addition, contact_type='INTER'):
'''
'''
atom.integer_sift = [x + y for x, y in zip(atom.sift, addition)]
if contact_type == 'INTER':
atom.integer_sift_inter_only = [x + y for x, y in zip(atom.sift_inter_only, addition)]
if 'INTRA' in contact_type:
atom.integer_sift_intra_only = [x + y for x, y in zip(atom.sift_intra_only, addition)]
if 'WATER' in contact_type:
atom.integer_sift_water_only = [x + y for x, y in zip(atom.sift_water_only, addition)]
def sift_xnor(sift1, sift2):
'''
'''
out = []
for x, y in zip(sift1, sift2):
if x and not y: # TF
out.append(0)
elif not x and not y: # FF
out.append(1)
elif x and y: # TT
out.append(1)
elif not x and y: # FT
out.append(0)
else:
raise ValueError('Invalid SIFts for matching: {} and {}'.format(sift1, sift2))
return out
def sift_match_base3(sift1, sift2):
'''
0 = UNMATCHED
1 = MATCHED
2 = MATCH NOT POSSIBLE
Assuming that sift1 is the potential SIFt, and sift2 is the actual SIFt.
'''
out = []
for x, y in zip(sift1, sift2):
if x and not y: # TF
out.append(0) # UNMATCHED
elif not x and not y: # FF
out.append(2) # MATCH NOT POSSIBLE
elif x and y: # TT
out.append(1) # MATCHED
elif not x and y: # FT
raise SiftMatchError
else:
raise ValueError('Invalid SIFts for matching: {} and {}'.format(sift1, sift2))
return out
def human_sift_match(sift_match, feature_sift=FEATURE_SIFT):
'''
Takes a base-3 SIFt indicating contact matched-ness and converts it to a human readable form.
'''
terms = []
for e, fp in enumerate(sift_match):
if fp == 2: # MATCH NOT POSSIBLE
continue
elif fp == 1:
terms.append('Matched {}'.format(feature_sift[e]))
elif fp == 0:
terms.append('Unmatched {}'.format(feature_sift[e]))
else:
raise ValueError
terms.sort()
return ':'.join(terms)
########
# MAIN #
########
if __name__ == '__main__':
# ARGUMENT PARSING
parser = argparse.ArgumentParser(description='''
############
# ARPEGGIO #
############
A program for calculating interactions,
using only Open Source dependencies.
Dependencies:
- Python (v2.7)
- Numpy
- BioPython (>= v1.60)
- OpenBabel (with Python bindings)
''', formatter_class=argparse.RawTextHelpFormatter)
parser.add_argument('pdb', type=str, help='Path to the PDB file to be analysed.')
selection_group = parser.add_mutually_exclusive_group(required=False)
selection_group.add_argument('-s', '--selection', type=str, nargs='+', help='Select the "ligand" for interactions, using selection syntax: /<chain_id>/<res_num>[<ins_code>]/<atom_name> or RESNAME:<het_id>. Fields can be omitted.')
selection_group.add_argument('-sf', '--selection-file', type=str, help='Selections as above, but listed in a file.')
parser.add_argument('-wh', '--write-hydrogenated', action='store_true', help='Write a PDB file including the added hydrogen coordinates.')
parser.add_argument('-mh', '--minimise-hydrogens', action='store_true', help='Energy minimise OpenBabel added hydrogens.')
parser.add_argument('-ms', '--minimisation-steps', type=int, default=50, help='Number of hydrogen minimisation steps to perform.')
parser.add_argument('-mf', '--minimisation-forcefield', type=str, choices=('MMFF94', 'UFF', 'Ghemical'), default='MMFF94', help='Choose the forcefield to minimise hydrogens with. Ghemical is not recommended.')
parser.add_argument('-mm', '--minimisation-method', type=str, choices=('DistanceGeometry', 'SteepestDescent', 'ConjugateGradients'), default='ConjugateGradients', help='Choose the method to minimise hydrogens with. ConjugateGradients is recommended.')
parser.add_argument('-co', '--vdw-comp', type=float, default=0.1, help='Compensation factor for VdW radii dependent interaction types.')
parser.add_argument('-i', '--interacting', type=float, default=5.0, help='Distance cutoff for grid points to be \'interacting\' with the entity.')
parser.add_argument('-ph', type=float, default=7.4, help='pH for hydrogen addition.')
parser.add_argument('-sa', '--include-sequence-adjacent', action='store_true', help='For intra-polypeptide interactions, include non-bonding interactions between residues that are next to each other in sequence; this is not done by default.')
parser.add_argument('-a', '--use-ambiguities', action='store_true', help='Turn on abiguous definitions for ambiguous contacts.')
parser.add_argument('-he', '--headers', action='store_true', help='Write out column headers in output files.')
#parser.add_argument('-sr', '--solvent_radius', type=float, default=1.4, help='Radius of solvent probe for accessibility calculations.')
#parser.add_argument('-ssp', '--solvent-sphere-points', type=int, default=960, help='Number of points to use for solvent shell spheres for accessibility calculations.')
#parser.add_argument('-st', '--sasa-threshold', type=float, default=1.0, help='Floating point solvent accessible surface area threshold (squared Angstroms) for considering an atom as \'accessible\' or not.')
#parser.add_argument('-ca', '--consider-all', action='store_true', help='Consider all entity/selection atoms, not just solvent accessible ones. If this is set, SASAs won\'t be calculated.')
#parser.add_argument('-spdb', '--sasa-pdb', action='store_true', help='Store a PDB with atom b-factors set based on boolean solvent accessibility.')
parser.add_argument('-op', '--output-postfix', type=str, help='Custom text to append to output filename (but before .extension).')
parser.add_argument('-v', '--verbose', action='store_true', help='Be chatty.')
args = parser.parse_args()
pdb_filename = args.pdb
VDW_COMP_FACTOR = args.vdw_comp
INTERACTING_THRESHOLD = args.interacting
#SOLVENT_RADIUS = args.solvent_radius
#NUM_SOLVENT_SPHERE_POINTS = args.solvent_sphere_points
#SASA_THRESHOLD = args.sasa_threshold
# LOGGING
if args.verbose:
logging.basicConfig(level=logging.INFO, format='%(levelname)s//%(asctime)s.%(msecs).03d//%(message)s', datefmt='%H:%M:%S')
else:
logging.basicConfig(level=logging.WARN, format='%(levelname)s//%(asctime)s.%(msecs).03d//%(message)s', datefmt='%H:%M:%S')
logging.info('Program begin.')
# ADDRESS AMBIGUITIES
if not args.use_ambiguities:
# REMOVE IF NOT USING THE AMBIGUITIES (DEFAULT)
# REMOVE FROM SMARTS DEFINITIONS
ATOM_TYPES['hbond acceptor'].pop('NH2 terminal amide', None)
ATOM_TYPES['hbond donor'].pop('oxygen amide term', None)
ATOM_TYPES['xbond acceptor'].pop('NH2 terminal amide', None)
ATOM_TYPES['weak hbond acceptor'].pop('NH2 terminal amide', None)
# REMOVE FROM PROTEIN ATOM DEFINITIONS
PROT_ATOM_TYPES['hbond acceptor'] = [x for x in PROT_ATOM_TYPES['hbond acceptor'] if x not in ('ASNND2', 'GLNNE2', 'HISCE1', 'HISCD2')]
PROT_ATOM_TYPES['hbond donor'] = [x for x in PROT_ATOM_TYPES['hbond donor'] if x not in ('ASNOD1', 'GLNOE1', 'HISCE1', 'HISCD2')]
PROT_ATOM_TYPES['xbond acceptor'] = [x for x in PROT_ATOM_TYPES['xbond acceptor'] if x not in ('ASNND2', 'GLNNE2', 'HISCE1', 'HISCD2')]
PROT_ATOM_TYPES['weak hbond acceptor'] = [x for x in PROT_ATOM_TYPES['weak hbond acceptor'] if x not in ('ASNND2', 'GLNNE2', 'HISCE1', 'HISCD2')]
# LOAD STRUCTURE (BIOPYTHON)
pdb_parser = PDBParser()
s = pdb_parser.get_structure('structure', pdb_filename)
s_atoms = list(s.get_atoms())
logging.info('Loaded PDB structure (BioPython)')
# CHECK FOR HYDROGENS IN THE INPUT STRUCTURE
input_has_hydrogens = False
hydrogens = [x for x in s_atoms if x.element == 'H']
if hydrogens:
logging.info('Detected that the input structure contains hydrogens. Hydrogen addition will be skipped.')
input_has_hydrogens = True
# LOAD STRUCTURE (OPENBABEL)
ob_conv = ob.OBConversion()
ob_conv.SetInFormat('pdb')
mol = ob.OBMol()
ob_conv.ReadFile(mol, pdb_filename)
logging.info('Loaded PDB structure (OpenBabel)')
# RENAME FOR OUTPUTS IF REQUESTED
if args.output_postfix:
pdb_filename = pdb_filename.replace('.pdb', args.output_postfix + '.pdb')
# CHECK THAT EACH ATOM HAS A UNIQUE SERIAL NUMBER
all_serials = [x.serial_number for x in s_atoms]
if len(all_serials) > len(set(all_serials)):
raise AtomSerialError
# MAPPING OB ATOMS TO BIOPYTHON ATOMS AND VICE VERSA
# FIRST MAP PDB SERIAL NUMBERS TO BIOPYTHON ATOMS FOR SPEED LATER
# THIS AVOIDS LOOPING THROUGH `s_atoms` MANY TIMES
serial_to_bio = {x.serial_number: x for x in s_atoms}
# DICTIONARIES FOR CONVERSIONS
ob_to_bio = {}
bio_to_ob = {}
for ob_atom in ob.OBMolAtomIter(mol):
serial = ob_atom.GetResidue().GetSerialNum(ob_atom)
# MATCH TO THE BIOPYTHON ATOM BY SERIAL NUMBER
try:
biopython_atom = serial_to_bio[serial]
except KeyError:
# ERRORWORTHY IF WE CAN'T MATCH AN OB ATOM TO A BIOPYTHON ONE
raise OBBioMatchError(serial)
# `Id` IS A UNIQUE AND STABLE ID IN OPENBABEL
# CAN RECOVER THE ATOM WITH `mol.GetAtomById(id)`
ob_to_bio[ob_atom.GetId()] = biopython_atom
bio_to_ob[biopython_atom] = ob_atom.GetId()
logging.info('Mapped OB to BioPython atoms and vice-versa.')
# ADD EMPTY DATA STRUCTURES FOR TAGGED ATOM DATA
# IN A SINGLE ITERATION
for atom in s_atoms:
# FOR ATOM TYPING VIA OPENBABEL
atom.atom_types = set([])
# LIST FOR EACH ATOM TO STORE EXPLICIT HYDROGEN COORDINATES
atom.h_coords = []
# DETECT METALS
if atom.element.upper() in METALS:
atom.is_metal = True
else:
atom.is_metal = False
# DETECT HALOGENS
if atom.element.upper() in HALOGENS:
atom.is_halogen = True
else:
atom.is_halogen = False
# ADD EXPLICIT HYDROGEN COORDS FOR H-BONDING INTERACTIONS
# ADDING HYDROGENS DOESN'T SEEM TO INTERFERE WITH ATOM SERIALS (THEY GET ADDED AS 0)
# SO WE CAN STILL GET BACK TO THE PERSISTENT BIOPYTHON ATOMS THIS WAY.
if not input_has_hydrogens:
mol.AddHydrogens(False, True, args.ph) # polaronly, correctForPH, pH
logging.info('Added hydrogens.')
# ATOM TYPING VIA OPENBABEL
# ITERATE OVER ATOM TYPE SMARTS DEFINITIONS
for atom_type, smartsdict in ATOM_TYPES.items():
#logging.info('Typing: {}'.format(atom_type))
# FOR EACH ATOM TYPE SMARTS STRING
for smarts in smartsdict.values():
#logging.info('Smarts: {}'.format(smarts))
# GET OPENBABEL ATOM MATCHES TO THE SMARTS PATTERN
ob_smart = ob.OBSmartsPattern()
ob_smart.Init(str(smarts))
#logging.info('Initialised for: {}'.format(smarts))
ob_smart.Match(mol)
#logging.info('Matched for: {}'.format(smarts))
matches = [x for x in ob_smart.GetMapList()]
#logging.info('List comp matches: {}'.format(smarts))
if matches:
# REDUCE TO A SINGLE LIST
matches = set(reduce(operator.add, matches))
#logging.info('Set reduce matches: {}'.format(smarts))
for match in matches:
atom = mol.GetAtom(match)
ob_to_bio[atom.GetId()].atom_types.add(atom_type)
#logging.info('Assigned types: {}'.format(smarts))
# ALL WATER MOLECULES ARE HYDROGEN BOND DONORS AND ACCEPTORS
for atom in (x for x in s_atoms if x.get_full_id()[3][0] == 'W'):
atom.atom_types.add('hbond acceptor')
atom.atom_types.add('hbond donor')
# OVERRIDE PROTEIN ATOM TYPING FROM DICTIONARY
for residue in s.get_residues():
if residue.resname in STD_RES:
for atom in residue.child_list:
# REMOVE TYPES IF ALREADY ASSIGNED FROM SMARTS
for atom_type in PROT_ATOM_TYPES.keys():
atom.atom_types.discard(atom_type)
# ADD ATOM TYPES FROM DICTIONARY
for atom_type, atom_ids in PROT_ATOM_TYPES.iteritems():
atom_id = residue.resname.strip() + atom.name.strip()
if atom_id in atom_ids:
atom.atom_types.add(atom_type)
with open(pdb_filename.replace('.pdb', '.atomtypes'), 'wb') as fo:
if args.headers:
fo.write('{}\n'.format('\t'.join(
['atom', 'atom_types']
)))
for atom in s_atoms:
fo.write('{}\n'.format('\t'.join([str(x) for x in [make_pymol_string(atom), sorted(tuple(atom.atom_types))]])))
logging.info('Typed atoms.')
# DETERMINE ATOM VALENCES AND EXPLICIT HYDROGEN COUNTS
for ob_atom in ob.OBMolAtomIter(mol):
if not input_has_hydrogens:
if ob_atom.IsHydrogen():
continue
# `http://openbabel.org/api/2.3/classOpenBabel_1_1OBAtom.shtml`
# CURRENT NUMBER OF EXPLICIT CONNECTIONS
valence = ob_atom.GetValence()
# MAXIMUM NUMBER OF CONNECTIONS EXPECTED
implicit_valence = ob_atom.GetImplicitValence()
# BOND ORDER
bond_order = ob_atom.BOSum()
# NUMBER OF BOUND HYDROGENS
num_hydrogens = ob_atom.ExplicitHydrogenCount()
# ELEMENT NUMBER
atomic_number = ob_atom.GetAtomicNum()
# FORMAL CHARGE
formal_charge = ob_atom.GetFormalCharge()
bio_atom = ob_to_bio[ob_atom.GetId()]
bio_atom.valence = valence
bio_atom.implicit_valence = implicit_valence
bio_atom.num_hydrogens = num_hydrogens
bio_atom.bond_order = bond_order
bio_atom.atomic_number = atomic_number
bio_atom.formal_charge = formal_charge
logging.info('Determined atom explicit and implicit valences, bond orders, atomic numbers, formal charge and number of bound hydrogens.')
# INITIALISE ATOM FULL SIFT
# INITIALISE ATOM FEATURE SIFT AND
# DETERMINE ATOM POTENTIAL FEATURE SIFT
# 5: 0: HBOND
# 6: 1: WEAK_HBOND
# 7: 2: HALOGEN_BOND
# 8: 3: IONIC
# 9: 4: METAL_COMPLEX
#10: 5: AROMATIC
#11: 6: HYDROPHOBIC
#12: 7: CARBONYL
# 8: POLAR - H-BONDS WITHOUT ANGLES
# 9: WEAK POLAR - WEAK H-BONDS WITHOUT ANGLES
# ALSO INITIALISE ATOM NUMBER OF POTENTIAL HBONDS/POLAR
# AND NUMBER OF ACTUAL HBONDS/POLAR
for atom in s_atoms:
# INTEGER SIFTS
atom.integer_sift = [0] * 15
atom.integer_sift_inter_only = [0] * 15
atom.integer_sift_intra_only = [0] * 15
atom.integer_sift_water_only = [0] * 15
# BINARY SIFTS
atom.sift = [0] * 15
atom.sift_inter_only = [0] * 15
atom.sift_intra_only = [0] * 15
atom.sift_water_only = [0] * 15
atom.potential_fsift = [0] * 10
atom.actual_fsift = [0] * 10
atom.actual_fsift_inter_only = [0] * 10
atom.actual_fsift_intra_only = [0] * 10
atom.actual_fsift_water_only = [0] * 10
atom.potential_hbonds = 0
atom.potential_polars = 0
atom.actual_hbonds = 0
atom.actual_polars = 0
atom.actual_hbonds_inter_only = 0
atom.actual_polars_inter_only = 0
atom.actual_hbonds_intra_only = 0
atom.actual_polars_intra_only = 0
atom.actual_hbonds_water_only = 0
atom.actual_polars_water_only = 0
# ATOM POTENTIAL FEATURE SIFT
# 0: HBOND
if 'hbond acceptor' in atom.atom_types or 'hbond donor' in atom.atom_types:
atom.potential_fsift[0] = 1
if 'hbond acceptor' in atom.atom_types:
# NUMBER OF LONE PAIRS
lone_pairs = VALENCE[atom.atomic_number] - atom.bond_order - atom.formal_charge
if lone_pairs != 0:
lone_pairs = lone_pairs / 2
atom.potential_hbonds = atom.potential_hbonds + lone_pairs
atom.potential_polars = atom.potential_polars + lone_pairs
if 'hbond donor' in atom.atom_types:
atom.potential_hbonds = atom.potential_hbonds + atom.num_hydrogens
atom.potential_polars = atom.potential_polars + atom.num_hydrogens
# 1: WEAK HBOND
if 'weak hbond acceptor' in atom.atom_types or 'weak hbond donor' in atom.atom_types or 'hbond donor' in atom.atom_types or 'hbond acceptor' in atom.atom_types or atom.is_halogen:
atom.potential_fsift[1] = 1
# 2: HALOGEN BOND
if 'xbond acceptor' in atom.atom_types or 'xbond donor' in atom.atom_types:
atom.potential_fsift[2] = 1
# 3: IONIC
if 'pos ionisable' in atom.atom_types or 'neg ionisable' in atom.atom_types:
atom.potential_fsift[3] = 1
# 4: METAL COMPLEX