forked from mpopara/SAXS_profile
-
Notifications
You must be signed in to change notification settings - Fork 0
/
EnsAvg_SAXS_profile.py
168 lines (106 loc) · 4.23 KB
/
EnsAvg_SAXS_profile.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
# -*- coding: utf-8 -*-
"""
Created on Thu Jul 28 11:35:05 2022
@author: popara
"""
import sys
sys.path.insert(1,'C:\\Program Files\\IMP-2.17.0\\python\\')
import IMP
import IMP.atom
import IMP.core
import IMP.saxs
import numpy as np
import pathlib
import os
import tqdm
import matplotlib as mpl
import matplotlib.pyplot as plt
def get_model_profile(
pdb_fn: pathlib.Path,
model_delta_q = 0.5 / 500,
model_min_q = 0.0,
model_max_q = 0.5
) -> IMP.saxs.Profile:
m = IMP.Model()
pdb_fn = str(pdb_fn)
saxs_fn = pdb_fn + '.saxs.dat'
# calculate SAXS profile
model_profile = IMP.saxs.Profile(model_min_q, model_max_q, model_delta_q)
if os.path.exists(saxs_fn):
model_profile.read_SAXS_file(saxs_fn)
else:
mp = IMP.atom.read_pdb(pdb_fn, m, IMP.atom.NonWaterNonHydrogenPDBSelector(), True, True)
# select particles from the model
particles = IMP.atom.get_by_type(mp, IMP.atom.ATOM_TYPE)
# add radius for water layer computation
ft = IMP.saxs.get_default_form_factor_table()
for i in range(0, len(particles)):
radius = ft.get_radius(particles[i])
IMP.core.XYZR.setup_particle(particles[i], radius)
# compute surface accessibility
s = IMP.saxs.SolventAccessibleSurface()
surface_area = s.get_solvent_accessibility(IMP.core.XYZRs(particles))
model_profile.calculate_profile_partial(particles, surface_area)
# Write SAXS curve
model_profile.write_SAXS_file(saxs_fn)
return model_profile
######## input parameters and data paths ##########################
pdb_path = pathlib.Path('C:/user/SAXS_profile/example_data/PDBs/')
weights_path = pathlib.Path('C:/user/SAXS_profile/example_data/')
outfile_path = pathlib.Path('C:/user/SAXS_profile/')
# q vector range for calculation of theoretical scattering profiles
model_delta_q = 0.6 / 500
model_min_q = 0.0
model_max_q = 0.6
################### Compute SAXS profiles for PDBs and ensemble average ##########
weights = np.loadtxt(weights_path/'conformer_weights.dat', delimiter=' ', usecols=[1])
pdb_fns = sorted(list(pdb_path.glob('*.pdb')))
assert len(weights)==len(pdb_fns), "Length of weights file does not match number of conformers"
model_profiles = list()
for fn in tqdm.tqdm(pdb_fns):
model_profiles.append(
get_model_profile(fn, model_delta_q, model_min_q, model_max_q)
)
avg_saxs_fn = str(outfile_path / 'ensemble_averaged_saxs_profile.dat')
EnsAvg_profile = IMP.saxs.Profile(model_min_q, model_max_q, model_delta_q)
for p, w in zip(model_profiles, weights):
EnsAvg_profile.add(p, w)
EnsAvg_profile.write_SAXS_file(avg_saxs_fn)
################ Visualization #################################
mpl.rcParams['font.sans-serif']="Arial"
plt.rcParams['figure.figsize'] = 3.25, 4.5
params = {'legend.fontsize': 6,
'legend.handlelength': 2}
plt.rcParams['legend.frameon'] = True
plt.rcParams['legend.fancybox'] = True
mpl.rcParams['axes.linewidth'] = 0.3
mpl.rcParams['axes.labelsize'] = 8
mpl.rcParams['xtick.labelsize'] = 6
mpl.rcParams['ytick.labelsize'] = 6
mpl.rcParams['xtick.major.width'] = 0.3
mpl.rcParams['ytick.major.width'] = 0.3
mpl.rcParams['xtick.minor.width'] = 0.3
mpl.rcParams['ytick.minor.width'] = 0.3
mpl.rcParams['xtick.major.size'] = 2
mpl.rcParams['ytick.major.size'] = 2
mpl.rcParams['xtick.minor.size'] = 1 # half of the major ticks length
mpl.rcParams['ytick.minor.size'] = 1
plt.rcParams.update(params)
fig = plt.figure()
gs = fig.add_gridspec(2, 1, hspace=0, wspace=0)
ax_t = fig.add_subplot(gs[0,0]) # top panel I versus q
ax_b = fig.add_subplot(gs[1,0]) # bottom; Kratky plot I*q^2 versus q
saxs_data = np.loadtxt(avg_saxs_fn)
ax_t.plot(saxs_data[:,0],saxs_data[:,1],'-', lw=0.8)
ax_b.plot(saxs_data[:,0],saxs_data[:,1]*saxs_data[:,0]**2,'-', lw=0.8)
ax_t.set_yscale('log')
ax_t.set_ylabel(r'$I(q)$')
ax_t.set_xlim(model_min_q, None)
ax_b.set_xlim(model_min_q, None)
ax_t.axes.xaxis.set_ticklabels([])
ax_b.set_xlabel(r'$q$ $[Å^{-1}$]')
ax_b.set_ylabel(r'$I(q) \cdot q^{2}$')
fig.tight_layout()
fig.savefig(outfile_path/'EnsAvg_saxs_profile.png',dpi=300, transparent=True)
fig.savefig(outfile_path/'EnsAvg_saxs_profile.svg', transparent=True)
plt.close(fig)