LOOS  v2.3.2
water-inside.cpp
1 /*
2  water-inside.cpp
3 
4 
5  (c) 2008 Tod D. Romo,
6  Grossfield Lab,
7  University of Rochester Medical and Dental School
8 
9 
10  Applies a given set of criteria to determine whether or not a water
11  is inside a protein. A matrix is then built up where each column
12  represents a timepoint in the trajectory and each row is the
13  internal water state (i.e. 1 = water is inside, 0 = water is not
14  inside)
15 
16  Also tracks the volume of the probe region (i.e. what's defined as
17  inside, if possible) and writes out a list of atomids that describe
18  which atoms go with which rows of the matrix.
19 
20 */
21 
22 #include <loos.hpp>
23 
24 
25 
26 #include <loos.hpp>
27 #include <cmath>
28 #include <boost/format.hpp>
29 
30 #include <DensityGrid.hpp>
31 #include <DensityTools.hpp>
32 #include <DensityOptions.hpp>
33 
34 using namespace std;
35 using namespace loos;
36 using namespace loos::DensityTools;
37 
38 
39 
41 
42 
43 string fullHelpMessage(void) {
44  string msg =
45  "\n"
46  "SYNOPSIS\n"
47  "\n"
48  "\tClassify waters as inside a protein or not in a trajectory\n"
49  "\n"
50  "DESCRIPTION\n"
51  "\n"
52  "\twater-inside applies a user-specified set of criteria to determine\n"
53  "whether or not water is inside a protein. A matrix is constructed where\n"
54  "each column is a time-series for each water. A 1 means the corresponding\n"
55  "water is inside the protein, and 0 means it's not. The volume of the probe region\n"
56  "is also tracked (if possible), and written out separately. In addition, a file\n"
57  "is written that maps the water atomids to the columns of the matrix.\n"
58  "See water-hist for more information about internal-water criteria.\n"
59  "\n"
60  "\nEXAMPLES\n"
61  "\twater-inside --prefix water foo.pdb foo.dcd\n"
62  "This example will use the axis filter for water (i.e. water atoms within\n"
63  "the default radius of 10 Angstroms from the first principal axis of the protein\n"
64  "selection. The default water selection (name == 'OH2') and protein selection\n"
65  "(name == 'CA') are used. The output prefix is set to 'water', so 'water.asc',\n"
66  "'water.vol', and 'water.atoms' will be created containing the time-series matrix,\n"
67  "the internal water region volume, and the atom mapping respectively.\n\n"
68  "\twater-inside --mode radius --radius 5 --prot 'resid == 65'\\\n"
69  "\t --prefix pocket foo.pdb foo.dcd\n"
70  "This example will find water atoms (using the default selection) that are within\n"
71  "5 Angstroms of any atom in residue 65, and use the output prefix 'pocket'.\n"
72  "\n"
73  "NOTES\n"
74  "\tLOOS does not care what is called a protein or water. You can use any selection,\n"
75  "for example, to track ligands, or lipids, etc.\n"
76  "\n"
77  "SEE ALSO\n"
78  "\twater-hist\n"
79  ;
80 
81  return(msg);
82 }
83 
84 
85 
86 
87 void writeAtomIds(const string& fname, const AtomicGroup& grp, const string& hdr) {
88  ofstream ofs(fname.c_str());
89  AtomicGroup::const_iterator ci;
90 
91  uint t = 0;
92  ofs << "# " << hdr << endl;
93  ofs << "# i\tatomid(i)\tresidue(i)\n";
94  for (ci = grp.begin(); ci != grp.end(); ++ci) {
95  ofs << boost::format("%u\t%d\t%s-%d\n") % t++ % (*ci)->id() % (*ci)->name() % (*ci)->resid();
96  }
97 }
98 
99 
100 int main(int argc, char *argv[]) {
101  string hdr = invocationHeader(argc, argv);
102 
103  opts::BasicOptions* basopts = new opts::BasicOptions(fullHelpMessage());
104  opts::OutputPrefix* prefopts = new opts::OutputPrefix;
106  opts::BasicWater* watopts = new opts::BasicWater;
107 
108  opts::AggregateOptions options;
109  options.add(basopts).add(prefopts).add(tropts).add(watopts);
110  if (!options.parse(argc, argv))
111  exit(-1);
112 
113 
114 
115  AtomicGroup model = tropts->model;
116  pTraj traj = tropts->trajectory;
117  vector<uint> frames = tropts->frameList();
118 
119  AtomicGroup subset = selectAtoms(model, watopts->prot_string);
120  AtomicGroup waters = selectAtoms(model, watopts->water_string);
121 
122  uint m = waters.size();
123  uint n = traj->nframes();
124  Matrix M(m,n);
125  Math::Matrix<double> V(n, 1);
126  cerr << boost::format("Water matrix is %d x %d.\n") % m % n;
127 
128  uint i = 0;
129  cerr << "Processing - ";
130 
131  for (vector<uint>::iterator t = frames.begin(); t != frames.end(); ++t) {
132  if (i % 100 == 0)
133  cerr << ".";
134 
135  traj->readFrame(*t);
136  traj->updateGroupCoords(model);
137 
138  vector<int> mask = watopts->filter_func->filter(waters, subset);
139  if (mask.size() != m) {
140  cerr << boost::format("ERROR - returned mask has size %u but expected %u.\n") % mask.size() % m;
141  exit(-10);
142  }
143 
144  for (uint j=0; j<m; ++j)
145  M(j,i) = mask[j];
146 
147  V(i,0) = watopts->filter_func->volume();
148  ++i;
149  }
150 
151  cerr << " done\n";
152  writeAsciiMatrix(prefopts->prefix + ".asc", M, hdr);
153  writeAsciiMatrix(prefopts->prefix + ".vol", V, hdr);
154  writeAtomIds(prefopts->prefix + ".atoms", waters, hdr);
155 }
Gets a string as prefix for output files (–prefix)
Trajectory with either a –range or –skip.
Options specific to tools that work with water/internal-water.
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
STL namespace.
bool parse(int argc, char *argv[])
Parses a command line, returning true if parsing was ok.
AtomicGroup model
Model that describes the trajectory.
virtual std::vector< int > filter(const loos::AtomicGroup &, const loos::AtomicGroup &)=0
Given a molecule and a set of waters, pick which waters are inside.
std::string water_string
User-specified strings.
Options common to all tools (including –fullhelp)
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
Namespace for Density package.
Definition: DensityGrid.hpp:48
std::ostream & writeAsciiMatrix(std::ostream &os, const Math::Matrix< T, P, S > &M, const std::string &meta, const Math::Range &start, const Math::Range &end, const bool trans=false, F fmt=F())
Write a submatrix to a stream.
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Definition: utils.cpp:195
DensityTools::WaterFilterBase * filter_func
Filter for determining internal waters.
std::vector< uint > frameList() const
Returns the list of frames the user requested.
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
virtual double volume(void)=0
Calculate the volume of the region we can pick waters from...
Combines a set of OptionsPackages.
AggregateOptions & add(OptionsPackage *pack)
Add a pointer to an OptionsPackage that will be used for options.