LOOS  v2.3.2
anm.cpp
1 /*
2  anm
3 
4  (c) 2008 Tod D. Romo, Grossfield Lab
5  Department of Biochemistry
6  University of Rochster School of Medicine and Dentistry
7 
8  Computes the anisotropic network model for a structure. It does
9  this by building a hessian for the structure, then computing the SVD
10  of it and the corresponding pseudo-inverse (ignoring the 6 lowest
11  modes).
12 
13  Usage:
14  anm [selection string] radius model-name output-prefix
15 
16  Examples:
17  anm 'resid >= 10 && resid <= 50 && name == "CA"' foo.pdb foo
18 
19  This creates the following files:
20  foo_H.asc == The hessian
21  foo_U.asc == Left singular vectors
22  foo_s.asc == Singular values
23  foo_V.asc == Right singular vectors
24  foo_Hi.asc == Pseudo-inverse of H
25 
26  Notes:
27  o The default selection (if none is specified) is to pick CA's
28  o The output is in ASCII format suitable for use with Matlab/Octave/Gnuplot
29 
30 
31 */
32 
33 
34 /*
35  This file is part of LOOS.
36 
37  LOOS (Lightweight Object-Oriented Structure library)
38  Copyright (c) 2008,2009 Tod D. Romo
39  Department of Biochemistry and Biophysics
40  School of Medicine & Dentistry, University of Rochester
41 
42  This package (LOOS) is free software: you can redistribute it and/or modify
43  it under the terms of the GNU General Public License as published by
44  the Free Software Foundation under version 3 of the License.
45 
46  This package is distributed in the hope that it will be useful,
47  but WITHOUT ANY WARRANTY; without even the implied warranty of
48  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
49  GNU General Public License for more details.
50 
51  You should have received a copy of the GNU General Public License
52  along with this program. If not, see <http://www.gnu.org/licenses/>.
53 
54 */
55 
56 
57 
58 #include <loos.hpp>
59 
60 #include <limits>
61 
62 
63 #include "hessian.hpp"
64 #include "enm-lib.hpp"
65 #include "anm-lib.hpp"
66 
67 
68 using namespace std;
69 using namespace loos;
70 using namespace ENM;
71 
72 
73 namespace opts = loos::OptionsFramework;
74 namespace po = loos::OptionsFramework::po;
75 
76 
77 // Globals... Icky poo!
78 
79 string prefix;
80 
81 int verbosity;
82 bool debug;
83 
84 string spring_desc;
85 string bound_spring_desc;
86 
87 string fullHelpMessage() {
88 
89  string s =
90  "\n"
91  "SYNOPSIS\n"
92  "\n"
93  "Compute the anisotropic network model for a structure.\n"
94  "\n"
95  "DESCRIPTION\n"
96  "\n"
97  "An anisotropic network model predicts the motions of a structure\n"
98  "using a harmonic contact (spring) potential. (See Atilgan et al. 2001)\n"
99  "It does this by building a hessian for the structure, then computing\n"
100  "the SVD of it and the corresponding pseudo-inverse (ignoring the 6\n"
101  "lowest modes).\n"
102  "\n"
103  "This creates the following files:\n"
104  "\tfoo_H.asc - The hessian\n"
105  "\tfoo_U.asc - Left singular vectors\n"
106  "\tfoo_s.asc - Singular values\n"
107  "\tfoo_V.asc - Right singular vectors\n"
108  "\tfoo_Hi.asc - Pseudo-inverse of H\n"
109  "\n"
110  "\n"
111  "* Spring Constant Control *\n"
112  "Contacts between beads in an ANM are connected by a single potential\n"
113  "which is described by a hookean spring. The stiffness of each connection\n"
114  "can be modified using various definitions of the spring constant.\n"
115  "The spring constant used is controlled by the --spring option.\n"
116  "If only the name for the spring function is given, then the default\n"
117  "parameters are used. Alternatively, the name may include a\n"
118  "comma-separated list of parameters to be passed to the spring\n"
119  "function, i.e. --spring=distance,15.0\n\n"
120  "Available spring functions:\n";
121 
122  vector<string> names = springNames();
123  for (vector<string>::const_iterator i = names.begin(); i != names.end(); ++i)
124  s = s + "\t" + *i + "\n";
125 
126  s +=
127  "\n\n"
128  "* Adding \"Connectivity\" *\n"
129  "ANM also supports construction of spring connections based on\n"
130  "pseudo-connectivity. This allows beads neighboring in sequence\n"
131  "to be connected by a separate \"bound\" spring, chosen using the\n"
132  "--bound option. In this case the other or \"non-bound\" spring is\n"
133  "chosen with the --spring option.\n"
134  "\n"
135  "\n\n"
136  "EXAMPLES\n\n"
137  "anm --selection 'resid >= 10 && resid <= 50 && name == \"CA\"' foo.pdb foo\n"
138  "\tCompute the ANM for residues #10 through #50 with a 15 Angstrom cutoff\n"
139  "\ti.e. construct contacts using only the CA's that are within 15 Angstroms\n"
140  "\n"
141  "anm -S=exponential,-1.3 foo.pdb foo\n"
142  "\tCompute an ANM using an spring function where the magnitude of\n"
143  "\tthe connection decays exponentially with distance at a rate of\n"
144  "\texp(-1.3*r) where r is the distance between contacts. Note:\n"
145  "\tin this case all beads are connected - which can eliminate\n"
146  "\tan error in the numeric eigendecomposition.\n"
147  "\n"
148  "anm --bound=constant,100 --spring=exponential,-1.3 foo.pdb foo\n"
149  "\tSimilar to the example above, but using connectivity. Here\n"
150  "\tresidues that are adjacent in sequence are connected by\n"
151  "\tsprings with a constant stiffness of \"100\" and all other\n"
152  "\tresidues are connected by springs that decay exponentially\n"
153  "\twith distance\n"
154  "\n";
155 
156  return(s);
157 }
158 
159 // @cond TOOLS_INTERNAL
160 class ToolOptions : public opts::OptionsPackage {
161 public:
162 
163  void addGeneric(po::options_description& o) {
164  o.add_options()
165  ("debug", po::value<bool>(&debug)->default_value(false), "Turn on debugging (output intermediate matrices)")
166  ("spring,S", po::value<string>(&spring_desc)->default_value("distance"),"Spring function to use")
167  ("bound", po::value<string>(&bound_spring_desc), "Bound spring");
168  }
169 
170  string print() const {
171  ostringstream oss;
172  oss << boost::format("debug=%d, spring='%s', bound='%s'") % debug % spring_desc % bound_spring_desc;
173  return(oss.str());
174  }
175 };
176 
177 // @endcond
178 
179 
180 
181 loos::Math::Matrix<int> buildConnectivity(const AtomicGroup& model) {
182  uint n = model.size();
183  loos::Math::Matrix<int> M(n, n);
184 
185  for (uint j=0; j<n-1; ++j)
186  for (uint i=j; i<n; ++i)
187  if (i == j)
188  M(j, i) = 1;
189  else {
190  M(j, i) = model[j]->isBoundTo(model[i]);
191  M(i, j) = M(j, i);
192  }
193 
194  return(M);
195 }
196 
197 
198 
199 int main(int argc, char *argv[]) {
200 
201  string header = invocationHeader(argc, argv);
202 
203  opts::BasicOptions* bopts = new opts::BasicOptions(fullHelpMessage());
204  opts::BasicSelection* sopts = new opts::BasicSelection("name == 'CA'");
206  ToolOptions* topts = new ToolOptions;
207  opts::RequiredArguments* ropts = new opts::RequiredArguments("prefix", "output-prefix");
208 
209  opts::AggregateOptions options;
210  options.add(bopts).add(sopts).add(mopts).add(topts).add(ropts);
211  if (!options.parse(argc, argv))
212  exit(-1);
213 
214  AtomicGroup model = mopts->model;
215  AtomicGroup subset = selectAtoms(model, sopts->selection);
216 
217  verbosity = bopts->verbosity;
218  prefix = ropts->value("prefix");
219 
220 
221  if (verbosity > 0)
222  cerr << boost::format("Selected %d atoms from %s\n") % subset.size() % mopts->model_name;
223 
224  // Determine which kind of scaling to apply to the Hessian...
225  vector<SpringFunction*> springs;
226  SpringFunction* spring = 0;
227  spring = springFactory(spring_desc);
228  springs.push_back(spring);
229 
230  vector<SuperBlock*> blocks;
231  SuperBlock* blocker = new SuperBlock(spring, subset);
232  blocks.push_back(blocker);
233 
234 
235  // Handle Decoration (if necessary)
236  if (!bound_spring_desc.empty()) {
237  if (! model.hasBonds()) {
238  cerr << "Error- cannot use bound springs unless the model has connectivity\n";
239  exit(-10);
240  }
241  loos::Math::Matrix<int> M = buildConnectivity(subset);
242  SpringFunction* bound_spring = springFactory(bound_spring_desc);
243  springs.push_back(bound_spring);
244 
245  BoundSuperBlock* decorator = new BoundSuperBlock(blocker, bound_spring, M);
246  blocks.push_back(decorator);
247 
248  blocker = decorator;
249  }
250 
251 
252  ANM anm(blocker);
253  anm.debugging(debug);
254  anm.prefix(prefix);
255  anm.meta(header);
256  anm.verbosity(verbosity);
257 
258  anm.solve();
259 
260 
261  // Write out the LSVs (or eigenvectors)
262  writeAsciiMatrix(prefix + "_U.asc", anm.eigenvectors(), header, false);
263  writeAsciiMatrix(prefix + "_s.asc", anm.eigenvalues(), header, false);
264 
265  writeAsciiMatrix(prefix + "_Hi.asc", anm.inverseHessian(), header, false);
266 
267  for (vector<SuperBlock*>::iterator i = blocks.begin(); i != blocks.end(); ++i)
268  delete *i;
269  for (vector<SpringFunction*>::iterator i = springs.begin(); i != springs.end(); ++i)
270  delete *i;
271 
272 }
Provides simple way to add command-line arguments (required options)
Namespace for encapsulating options processing.
This class creates superblocks in a hessian.
Definition: hessian.hpp:52
Provides a single LOOS selection (–selection)
bool hasBonds(void) const
Does any atom in the group have bond information???
STL namespace.
bool parse(int argc, char *argv[])
Parses a command line, returning true if parsing was ok.
Options common to all tools (including –fullhelp)
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
Decorator for switching spring functions based on a matrix of flags.
Definition: hessian.hpp:186
std::string value(const std::string &s) const
Retrieve the value for an argument.
Interface for ENM spring functions.
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.
SpringFunction * springFactory(const std::string &spring_desc)
Factory function for generating new SpringFunction instances based on a user string.
Namespace to encapsulate Elastic Network Model routines.
Definition: anm-lib.hpp:32
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Definition: utils.cpp:195
std::vector< std::string > springNames()
List of possible names for springFactory()
Request a model with coordinates.
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
Anisotropic network model.
Definition: anm-lib.hpp:35
Combines a set of OptionsPackages.
AggregateOptions & add(OptionsPackage *pack)
Add a pointer to an OptionsPackage that will be used for options.