-
Notifications
You must be signed in to change notification settings - Fork 1
/
saxs_scoring.py
128 lines (95 loc) · 3.69 KB
/
saxs_scoring.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
# -*- 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 typing
import tqdm
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
def score_weighted_profiles(
model_weights: typing.List[float],
fit_file_name: str
) -> float:
# Compute weighted average over all models
model_profile = IMP.saxs.Profile(model_min_q, model_max_q, model_delta_q)
for p, w in zip(model_profiles, model_weights):
model_profile.add(p, w)
fit_settings = {
'min_c1': 0.95, # c1 adjusts the excluded volume
'max_c1': 1.05,
'min_c2': -2.0, # c2 adjusts the density of hydration layer
'max_c2': 4.0,
'use_offset': True,
'fit_file_name': fit_file_name
}
# calculate chi-square score
saxs_score = IMP.saxs.ProfileFitterChi(exp_profile)
fit_parameters = saxs_score.fit_profile(model_profile, *fit_settings.values())
chi2 = fit_parameters.get_chi_square()
return chi2
########## input data path and input parameters
basePath = pathlib.Path('C:/user/folder')
expPath = basePath/'experimental_SAXS_profile.dat'
pdbPath = basePath /'PDBs/'
weightsPath = basePath/'conformer_weights.dat'
model_delta_q = 0.6 / 500
model_min_q = 0.0
model_max_q = 0.6
################################################################
# read experimental profile: IMP.saxs.Profile(file_name, fit_file, max_q, units)
# file_name: profile file name
# fit_file: if true, intensities are read from column 3
# max_q: read till maximal q value = max_q, or all if max_q<=0
# units: gets 1, 2, or 3 for unknown q units, 1/A, or 1/nm
exp_profile = IMP.saxs.Profile(str(expPath), False, 0.55, 3)
print('min_q = ' + str(exp_profile.get_min_q()))
print('max_q = ' + str(exp_profile.get_max_q()))
print('delta_q = ' + str(exp_profile.get_delta_q()))
# Compute model profiles for PDBs or if saxs profiles already exist in pdbPath, it will read those
pdb_fns = sorted(list(pdbPath.glob('*.pdb')))
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)
)
weights = np.loadtxt(weightsPath)[:,1]
fit_file_name = str(basePath / 'saxs.fit') # file containing fit parameters, as well as the exp and fitted theoretical saxs profile
chi2 = score_weighted_profiles(weights, fit_file_name)
print('chi2 = ' + str(chi2))