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
|
#!/usr/bin/python
# Copyright (C) 2002, Thomas Hamelryck (thamelry@vub.ac.be)
# This code is part of the Biopython distribution and governed by its
# license. Please see the LICENSE file that should have been included
# as part of this package.
"""Calculates solvent exposure for a PDB file using one of 5 different methods.
-DSSP (DSSP needs to be installed)
-Residue depth (MSMS needs to be installed)
-Coordination number (ie. number of CA atoms within a sphere)
-HSEalpha half sphere exposure
-HSEbeta half sphere exposure
A PDB file can be written out with the exposure in the B factor field.
See --help for all args.
"""
import argparse
import sys
from Bio.PDB import (
DSSP,
ExposureCN,
HSExposureCA,
HSExposureCB,
PDBParser,
PDBIO,
ResidueDepth,
Selection,
)
ap = argparse.ArgumentParser(description=__doc__)
ap.add_argument("pdbfile", help="Input structure in PDB format.")
ap.add_argument(
"-t",
"--type",
dest="exp",
choices=["HSEAU", "HSEAD", "HSEBU", "HSEBD", "CN", "DSSPr", "DSSPa", "RD", "RDa"],
help="Exposure Type",
default="HSEb",
)
ap.add_argument(
"-o", "--out", dest="outfile", help="output to PDB file (B factor=exposure)",
)
ap.add_argument(
"-r",
"--radius",
dest="radius",
type=float,
help="sphere radius (default 13.0 A)",
default=13.0,
)
ap.add_argument(
"-m",
"--model",
dest="model",
type=int,
help="PDB model number (default 0)",
default=0,
)
ap.add_argument("--dssp", help="Path to the DSSP executable", default="dssp")
ap.add_argument("--msms", help="Path to the MSMS executable", default=None)
args = ap.parse_args()
# Get the structure
p = PDBParser()
s = p.get_structure("X", args.pdbfile)
# First model by default
m = s[args.model]
RADIUS = args.radius
# d=dictionary of exposures
# k=position in ntuple containing the desired exposure
format = "%4i"
args.exp = args.exp.upper()
if args.exp[0] == "H" and args.exp[3] == "A":
hse = HSExposureCA(m, RADIUS)
if args.exp[-1] == "D":
k = "EXP_HSE_A_D"
else:
k = "EXP_HSE_A_U"
elif args.exp[0] == "H" and args.exp[3] == "B":
hse = HSExposureCB(m, RADIUS)
# hse.write_pymol_script()
if args.exp[-1] == "D":
k = "EXP_HSE_B_U"
else:
k = "EXP_HSE_B_D"
elif args.exp == "CN":
hse = ExposureCN(m, RADIUS)
k = "EXP_CN"
elif args.exp == "ANGLE":
hse = HSExposureCA(m, RADIUS)
k = "EXP_CB_PCB_ANGLE"
format = "%4.1f"
elif args.exp == "DSSPR":
d = DSSP(m, args.pdbfile, dssp=args.dssp)
k = "EXP_DSSP_RASA"
format = "%.4f"
elif args.exp == "DSSPA":
d = DSSP(m, args.pdbfile, dssp=args.dssp)
k = "EXP_DSSP_ASA"
elif args.exp == "RD":
d = ResidueDepth(m, args.pdbfile, msms_exec=args.msms)
k = "EXP_RD"
format = "%4.1f"
elif args.exp == "RDA":
d = ResidueDepth(m, args.pdbfile, msms_exec=args.msms)
k = "EXP_RD_CA"
format = "%4.1f"
else:
print("ERROR: Unknown option.")
sys.exit()
residue_list = Selection.unfold_entities(m, "R")
for r in residue_list:
if k in r.xtra:
exposure = r.xtra[k]
if args.exp == "DSSPR":
# to 0=exposed, 1=buried
exposure = 1 - exposure
# Print info
hetflag, resseq, icode = r.get_id()
if icode == " ":
icode = "_"
resname = r.get_resname()
print(("%s %4i %c\t" + format) % (resname, resseq, icode, exposure))
else:
exposure = 0.0
for atom in r.get_iterator():
atom.set_bfactor(exposure)
if args.outfile:
io = PDBIO()
io.set_structure(s)
io.save(args.outfile)
|
