15 #include "JmolWriter.h"
21 #include <boost/program_options.hpp>
22 namespace po = boost::program_options;
24 int main(
int argc,
char **argv)
27 for (
int i = 0; i < argc; i++) std::cerr << argv[i] <<
" ";
28 std::cerr << std::endl;
30 int profile_size = 500;
32 float background_adjustment_q = 0.0;
33 float excluded_volume_c1 = 0.0;
34 bool use_offset =
false;
36 float MAX_C2 = 4.0;
float MIN_C2 = -MAX_C2/2.0;
37 float water_layer_c2 = MAX_C2;
38 bool heavy_atoms_only =
true;
39 bool residue_level =
false;
40 bool score_log =
false;
41 bool interval_chi =
false;
42 int multi_model_pdb = 1;
43 po::options_description desc(
"Options");
45 (
"help",
"Any number of input PDBs and profiles is supported. \
46 Each PDB will be fitted against each profile.")
47 (
"version",
"FoXS (IMP applications)\nCopyright 2007-2013 IMP Inventors.\n\
48 All rights reserved. \nLicense: GNU LGPL version 2.1 or later\n\
49 <http://gnu.org/licenses/lgpl.html>.\n\
50 Written by Dina Schneidman.")
51 (
"max_q,q", po::value<float>(&max_q)->default_value(0.5),
52 "maximal q value (default = 0.5)")
53 (
"profile_size,s", po::value<int>(&profile_size)->default_value(500),
54 "number of points in the profile (default = 500)")
56 po::value<float>(&water_layer_c2)->default_value(MAX_C2),
57 "set hydration layer density. \
58 Valid range: -2.0 < c2 < 4.0 (default = 0.0)")
59 (
"hydrogens,h",
"explicitly consider hydrogens in PDB files \
61 (
"residues,r",
"perform fast coarse grained profile calculation using \
62 CA atoms only (default = false)")
64 po::value<float>(&excluded_volume_c1)->default_value(0.0),
65 "excluded volume parameter, enumerated by default. \
66 Valid range: 0.95 < c1 < 1.05")
68 po::value<float>(&background_adjustment_q)->default_value(0.0),
69 "background adjustment, not used by default. if enabled, \
70 recommended q value is 0.2")
71 (
"offset,o",
"use offset in fitting (default = false)")
72 (
"multi-model-pdb,m", po::value<int>(&multi_model_pdb)->default_value(1),
73 "1 - read the first MODEL only (default), \
74 2 - read each MODEL into a separate structure, \
75 3 - read all models into a single structure")
76 (
"score_log,l",
"use log(intensity) in fitting and scoring \
81 std::string form_factor_table_file;
82 bool ab_initio =
false;
84 bool javascript =
false;
85 po::options_description hidden(
"Hidden options");
87 (
"input-files", po::value< std::vector<std::string> >(),
88 "input PDB and profile files")
89 (
"form_factor_table,f", po::value<std::string>(&form_factor_table_file),
91 (
"ab_initio,a",
"compute profile for a bead model with \
92 constant form factor (default = false)")
93 (
"vacuum,v",
"compute profile in vacuum (default = false)")
95 "output javascript for browser viewing of the results (default = false)")
96 (
"interval_chi,i",
"compute chi for intervals (default = false)")
99 po::options_description cmdline_options;
100 cmdline_options.add(desc).add(hidden);
102 po::options_description visible(
"Usage: <pdb_file1> <pdb_file2> \
103 ... <profile_file1> <profile_file2> ... ");
106 po::positional_options_description p;
107 p.add(
"input-files", -1);
108 po::variables_map vm;
109 po::store(po::command_line_parser(argc,
110 argv).options(cmdline_options).positional(p).run(), vm);
113 std::vector<std::string> files, pdb_files, dat_files;
114 if(vm.count(
"input-files")) {
115 files = vm[
"input-files"].as< std::vector<std::string> >();
117 if(vm.count(
"help") || files.size() == 0) {
118 std::cout << visible <<
"\n";
121 if(vm.count(
"hydrogens")) heavy_atoms_only=
false;
122 if(vm.count(
"residues")) residue_level=
true;
123 if(vm.count(
"offset")) use_offset=
true;
124 if(vm.count(
"score_log")) score_log=
true;
126 if(vm.count(
"ab_initio")) { ab_initio=
true; water_layer_c2 = 0.0;
127 fit =
false; excluded_volume_c1 = 1.0; }
128 if(vm.count(
"vacuum")) { vacuum =
true; }
129 if(vm.count(
"javascript")) { javascript =
true; }
130 if(vm.count(
"interval_chi")) { interval_chi =
true; }
132 if(multi_model_pdb != 1 && multi_model_pdb != 2 && multi_model_pdb != 3) {
133 std::cerr <<
"Incorrect option for multi_model_pdb "
134 << multi_model_pdb << std::endl;
135 std::cerr <<
"Use 1 to read first MODEL only\n"
136 <<
" 2 to read each MODEL into a separate structure,\n"
137 <<
" 3 to read all models into a single structure\n";
138 std::cerr <<
"Default value of 1 is used\n";
141 float delta_q = max_q / profile_size;
143 bool reciprocal =
false;
145 if(form_factor_table_file.length() > 0) {
148 0.0, max_q, delta_q);
154 std::vector<IMP::Particles> particles_vec;
155 std::vector<IMP::saxs::Profile *> exp_profiles;
156 for(
unsigned int i=0; i<files.size(); i++) {
158 std::ifstream in_file(files[i].c_str());
160 std::cerr <<
"Can't open file " << files[i] << std::endl;
return 1;
164 IMP::atom::Hierarchies mhds;
174 if(multi_model_pdb == 2) {
177 if(multi_model_pdb == 3) {
189 for(
unsigned int h_index=0; h_index<mhds.size(); h_index++) {
190 IMP::ParticlesTemp particles =
192 if(particles.size() > 0) {
193 std::string pdb_id = files[i];
194 if(mhds.size() > 1) {
195 pdb_id = trim_extension(files[i]) +
"_m" +
196 std::string(boost::lexical_cast<std::string>(h_index+1)) +
".pdb";
198 pdb_files.push_back(pdb_id);
199 particles_vec.push_back(IMP::get_as<IMP::Particles>(particles));
200 std::cout << particles.size() <<
" atoms were read from PDB file "
202 if(mhds.size() > 1) std::cout <<
" MODEL " << h_index+1;
203 std::cout << std::endl;
205 if(water_layer_c2 != 0.0) {
208 if(residue_level) ff_type = IMP::saxs::CA_ATOMS;
209 if(!heavy_atoms_only) ff_type = IMP::saxs::ALL_ATOMS;
210 for(
unsigned int p_index=0; p_index<particles.size(); p_index++) {
211 float radius = ft->
get_radius(particles[p_index], ff_type);
217 }
catch(IMP::ValueException e) {
220 if(profile->
size() == 0) {
221 std::cerr <<
"can't parse input file " << files[i] << std::endl;
224 dat_files.push_back(files[i]);
225 if(background_adjustment_q > 0.0) {
228 exp_profiles.push_back(profile);
229 std::cout <<
"Profile read from file " << files[i] <<
" size = "
230 << profile->
size() << std::endl;
237 std::vector<IMP::saxs::Profile *> profiles;
238 std::vector<IMP::saxs::FitParameters> fps;
239 for(
unsigned int i=0; i<particles_vec.size(); i++) {
243 if(water_layer_c2 != 0.0) {
249 std::cerr <<
"Computing profile for " << pdb_files[i]
250 <<
" " << particles_vec[i].size() <<
" atoms "<< std::endl;
254 if(excluded_volume_c1 == 1.0 && water_layer_c2 == 0.0) fit =
false;
256 if(!heavy_atoms_only) ff_type = IMP::saxs::ALL_ATOMS;
257 if(residue_level) ff_type = IMP::saxs::CA_ATOMS;
260 float average_radius = 0.0;
262 for(
unsigned int k=0; k<particles_vec[i].size(); k++) {
263 average_radius += ft->
get_radius(particles_vec[i][k], ff_type);
264 volume += ft->
get_volume(particles_vec[i][k], ff_type);
266 average_radius /= particles_vec[i].size();
267 partial_profile->set_average_radius(average_radius);
269 if(dat_files.size() == 0 || !fit) {
272 calculate_profile_constant_form_factor(particles_vec[i]);
276 surface_area, ff_type);
280 ff_type, reciprocal);
285 partial_profile->calculate_profile_reciprocal_partial(particles_vec[i],
286 surface_area, ff_type);
289 surface_area, ff_type);
292 profiles.push_back(partial_profile);
294 std::string profile_file_name = std::string(pdb_files[i]) +
".dat";
297 Gnuplot::print_profile_script(pdb_files[i]);
300 for(
unsigned int j=0; j<dat_files.size(); j++) {
303 std::string fit_file_name2 = trim_extension(pdb_files[i]) +
"_" +
304 trim_extension(basename(const_cast<char *>(dat_files[j].c_str())))
307 float min_c1=0.95;
float max_c1=1.05;
308 if(excluded_volume_c1 > 0.0) { min_c1 = max_c1 = excluded_volume_c1; }
309 if(std::fabs(water_layer_c2 - MAX_C2) < 0.00000000001) {
310 }
else { MIN_C2 = MAX_C2 = water_layer_c2; }
312 IMP::saxs::FitParameters fp;
314 IMP::Pointer<IMP::saxs::ProfileFitter<IMP::saxs::ChiScoreLog> > pf =
315 new IMP::saxs::ProfileFitter<IMP::saxs::ChiScoreLog>(*exp_saxs_profile);
316 fp = pf->fit_profile(*partial_profile,
317 min_c1, max_c1, MIN_C2, MAX_C2,
318 use_offset, fit_file_name2);
320 IMP::Pointer<IMP::saxs::ProfileFitter<IMP::saxs::ChiScore> > pf =
321 new IMP::saxs::ProfileFitter<IMP::saxs::ChiScore>(*exp_saxs_profile);
322 fp = pf->fit_profile(*partial_profile,
323 min_c1, max_c1, MIN_C2, MAX_C2,
324 use_offset, fit_file_name2);
326 std::cout <<
"interval_chi " <<pdb_files[i] <<
" "
327 << pf->compute_score(*partial_profile, 0.0, 0.05) <<
" "
328 << pf->compute_score(*partial_profile, 0.0, 0.1) <<
" "
329 << pf->compute_score(*partial_profile, 0.0, 0.15) <<
" "
330 << pf->compute_score(*partial_profile, 0.0, 0.2) <<
" "
331 << pf->compute_score(*partial_profile) << std::endl;
334 std::cout << pdb_files[i] <<
" " << dat_files[j]
335 <<
" Chi = " << fp.get_chi()
336 <<
" c1 = " << fp.get_c1()
337 <<
" c2 = " << fp.get_c2()
338 <<
" default chi = " << fp.get_default_chi() << std::endl;
339 fp.set_pdb_file_name(pdb_files[i]);
340 fp.set_profile_file_name(dat_files[j]);
342 Gnuplot::print_fit_script(fp);
347 std::sort(fps.begin(), fps.end(),
348 IMP::saxs::FitParameters::compare_fit_parameters());
350 if(pdb_files.size() > 1) {
351 Gnuplot::print_profile_script(pdb_files);
352 if(dat_files.size() > 0) Gnuplot::print_fit_script(fps);
355 if(dat_files.size() > 0) {
356 Gnuplot::print_canvas_script(fps, JmolWriter::MAX_DISPLAY_NUM_);
357 JmolWriter::prepare_jmol_script(fps, particles_vec,
"jmoltable");
359 Gnuplot::print_canvas_script(pdb_files, JmolWriter::MAX_DISPLAY_NUM_);
360 JmolWriter::prepare_jmol_script(pdb_files, particles_vec,
"jmoltable");