IMP  2.0.0
The Integrative Modeling Platform
rg.cpp
1 /**
2  This is the program for computation of radius of gyration from SAXS profiles.
3  */
4 #include <IMP/Model.h>
5 #include <IMP/atom/pdb.h>
6 
7 #include <IMP/saxs/Profile.h>
8 #include <IMP/saxs/utility.h>
9 
10 #include <fstream>
11 #include <vector>
12 #include <string>
13 
14 #include <boost/program_options.hpp>
15 namespace po = boost::program_options;
16 
17 int main(int argc, char **argv)
18 {
19  // output arguments
20  for (int i = 0; i < argc; i++) std::cerr << argv[i] << " ";
21  std::cerr << std::endl;
22 
23  float end_q_rg = 1.3;
24  po::options_description desc("Usage: <pdb_file1> <pdb_file2> \
25 ... <profile_file1> <profile_file2> ...");
26  desc.add_options()
27  ("help", "Any number of input PDBs and profiles is supported. \
28 Each PDB will be fitted against each profile.")
29  ("input-files", po::value< std::vector<std::string> >(),
30  "input PDB and profile files")
31  ("end_q_rg,q*rg", po::value<float>(&end_q_rg)->default_value(1.3),
32  "end q*rg value for rg computation, q*rg<end_q_rg (default = 1.3), \
33 use 0.8 for elongated proteins")
34  ;
35  po::positional_options_description p;
36  p.add("input-files", -1);
37  po::variables_map vm;
38  po::store(
39  po::command_line_parser(argc,argv).options(desc).positional(p).run(), vm);
40  po::notify(vm);
41 
42  std::vector<std::string> files, pdb_files, dat_files;
43  if (vm.count("input-files")) {
44  files = vm["input-files"].as< std::vector<std::string> >();
45  }
46  if(vm.count("help") || files.size() == 0) {
47  std::cout << desc << "\n";
48  return 0;
49  }
50 
51  // 1. read pdbs and profiles, prepare particles
52  IMP::Model *model = new IMP::Model();
53  std::vector<IMP::Particles> particles_vec;
54  std::vector<IMP::saxs::Profile *> exp_profiles;
55  for(unsigned int i=0; i<files.size(); i++) {
56  // check if file exists
57  std::ifstream in_file(files[i].c_str());
58  if(!in_file) {
59  std::cerr << "Can't open file " << files[i] << std::endl;
60  exit(1);
61  }
62  // A. try as pdb
63  try {
65  IMP::atom::read_pdb(files[i], model,
67  // don't add radii
68  true, true);
69  IMP::Particles particles
70  = IMP::get_as<IMP::Particles>(get_by_type(mhd, IMP::atom::ATOM_TYPE));
71  if(particles.size() > 0) { // pdb file
72  pdb_files.push_back(files[i]);
73  particles_vec.push_back(particles);
74  std::cout << particles.size() << " atoms were read from PDB file "
75  << files[i] << std::endl;
76  }
77  } catch(IMP::ValueException e) { // not a pdb file
78  // B. try as dat file
79  IMP::saxs::Profile *profile = new IMP::saxs::Profile(files[i]);
80  if(profile->size() == 0) {
81  std::cerr << "can't parse input file " << files[i] << std::endl;
82  return 1;
83  } else {
84  dat_files.push_back(files[i]);
85  exp_profiles.push_back(profile);
86  std::cout << "Profile read from file " << files[i] << " size = "
87  << profile->size() << std::endl;
88  }
89  }
90  }
91 
92  // 2. compute rg for input profiles
93  for(unsigned int i=0; i<dat_files.size(); i++) {
94  IMP::saxs::Profile* exp_saxs_profile = exp_profiles[i];
95  double rg = exp_saxs_profile->radius_of_gyration(end_q_rg);
96  std::cout << dat_files[i] << " Rg= " << rg << std::endl;
97  }
98 
99  // 3. compute rg for input pdbs
100  for(unsigned int i=0; i<pdb_files.size(); i++) {
101  double rg = IMP::saxs::radius_of_gyration(particles_vec[i]);
102  std::cout << pdb_files[i] << " Rg= " << rg << std::endl;
103  }
104  return 0;
105 }