LOOS  v2.3.2
gnm.cpp
1 /*
2  gnm
3 
4  Computes the gaussian normal mode (elastic network model)
5  decomposition for a structure. Does this by building the Kirchoff
6  matrix given a PDB and a selection, then computes the SVD of the
7  matrix and finally the pseudo-inverse.
8 
9 
10  Usage:
11 
12  gnm [selection-string] radius model-name output-prefix
13 
14  Examples:
15 
16  gnm 'resid >= 10 && resid <= 50 && name == "CA"' 7.0 foo.pdb foo
17 
18  This will create the following files:
19  foo_K.asc == Kirchoff matrix
20  foo_U.asc == Left singular vectors
21  foo_s.asc == singular values
22  foo_V.asc == Right singular vectors
23  foo_Ki.asc == Pseudo-inverse of K
24 
25  Notes:
26  o The default selection (if none is specified) is to pick CA's
27  o The output is in ASCII format suitable for use with Matlab/Octave/Gnuplot
28 
29 */
30 
31 
32 /*
33  This file is part of LOOS.
34 
35  LOOS (Lightweight Object-Oriented Structure library)
36  Copyright (c) 2008 Tod D. Romo
37  Department of Biochemistry and Biophysics
38  School of Medicine & Dentistry, University of Rochester
39 
40  This package (LOOS) is free software: you can redistribute it and/or modify
41  it under the terms of the GNU General Public License as published by
42  the Free Software Foundation under version 3 of the License.
43 
44  This package is distributed in the hope that it will be useful,
45  but WITHOUT ANY WARRANTY; without even the implied warranty of
46  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
47  GNU General Public License for more details.
48 
49  You should have received a copy of the GNU General Public License
50  along with this program. If not, see <http://www.gnu.org/licenses/>.
51 
52 */
53 
54 
55 
56 #include <loos.hpp>
57 
58 #include <boost/format.hpp>
59 #include <boost/program_options.hpp>
60 
61 using namespace std;
62 using namespace loos;
63 namespace po = boost::program_options;
64 
66 
67 
68 // Globals... Fah!
69 
70 string selection;
71 string model_name;
72 string prefix;
73 double cutoff;
74 
75 void fullHelp() {
76  //string msg =
77  cout << "\n"
78  "\n"
79  "SYNOPSIS\n"
80  "\n"
81  "Compute the normal modes of a guassian network model\n"
82  "\n"
83  "DESCRIPTION\n"
84  "Computes the gaussian normal mode analysis of an ENM\n"
85  "This is done by building the Kirchoff matrix given a PDB\n"
86  "and a selection, then computing the SVD of the matrix and\n"
87  "finally computing the pseudo-inverse.\n"
88  "See: Bahar, et al., Folding and Design 2, 173-181, (1997).\n"
89  "\n"
90  "This will create the following files:\n"
91  "\tfoo_K.asc - Kirchoff matrix\n"
92  "\tfoo_U.asc - Left singular vectors\n"
93  "\tfoo_s.asc - singular values\n"
94  "\tfoo_V.asc - Right singular vectors\n"
95  "\tfoo_Ki.asc - Pseudo-inverse of K\n"
96  "\n"
97  "Notes:\n"
98  "- The default selection (if none is specified) is to pick CA's\n"
99  "- The output is ASCII format suitable for use with Matlab/Octave/Gnuplot\n"
100  //
101  "\n"
102  "EXAMPLES\n"
103  "\n"
104  "gnm -c8.2 -s 'resid >= 10 && resid <= 50 && name == \"CA\"' model.pdb foo\n"
105  "\tCompute the GNM of model.pdb for residues #10 through #50 with\n"
106  "\tan 8.2 Angstrom cutoff i.e. construct contacts using only the CA's\n"
107  "\tthat are within 8.2 Angstroms. Write out the files to foo_X.asc\n"
108  "\t\n"
109  "SEE ALSO\n"
110  "\n"
111  "Packages/ElasticNetworks/anm - \n"
112  "The anisotropic version of this tool. Here eigenvectors predicting\n"
113  "the direction of movements are written out as well.\n"
114  "\t\n"
115  "Packages/ElasticNetworks/vsa - \n"
116  "This is an extension of the ANM method mentioned above that splits\n"
117  "the calculation into two parts - a subsystem and an environment.\n"
118  "These eigendecompositions of these two parts are performed separately\n"
119  "and the environment can then be 'subtracted' off the subsystem.\n"
120  "\n";
121  }
122 
123 void parseOptions(int argc, char *argv[]) {
124 
125  try {
126  po::options_description generic("Allowed options");
127  generic.add_options()
128  ("help", "Produce this help message")
129  ("fullhelp", "Get extended help")
130  ("selection,s", po::value<string>(&selection)->default_value("name == 'CA'"), "Which atoms to use for the network")
131  ("cutoff,c", po::value<double>(&cutoff)->default_value(7.0), "Cutoff distance for node contact");
132 
133  po::options_description hidden("Hidden options");
134  hidden.add_options()
135  ("model", po::value<string>(&model_name), "Model filename")
136  ("prefix", po::value<string>(&prefix), "Output prefix");
137 
138  po::options_description command_line;
139  command_line.add(generic).add(hidden);
140 
141  po::positional_options_description p;
142  p.add("model", 1);
143  p.add("prefix", 1);
144 
145  po::variables_map vm;
146  po::store(po::command_line_parser(argc, argv).
147  options(command_line).positional(p).run(), vm);
148  po::notify(vm);
149 
150  if (vm.count("help") || vm.count("fullhelp") || !(vm.count("model") && vm.count("prefix"))) {
151  cerr << "Usage- gnm [options] model-name output-prefix\n";
152  cerr << generic;
153  if (vm.count("fullhelp"))
154  fullHelp();
155  exit(-1);
156  }
157  }
158  catch(exception& e) {
159  cerr << "Error - " << e.what() << endl;
160  exit(-1);
161  }
162 }
163 
164 // This is the Kirchoff normalization constant (see Bahar, Atilgan,
165 // and Erman. Folding & Design 2:173)
166 double normalization = 1.0;
167 
168 
169 Matrix kirchoff(AtomicGroup& group, const double cutoff) {
170  int n = group.size();
171  Matrix M(n, n);
172  double r2 = cutoff * cutoff;
173 
174 
175  for (int j=1; j<n; j++)
176  for (int i=0; i<j; i++)
177  if (group[i]->coords().distance2(group[j]->coords()) <= r2)
178  M(i, j) = M(j, i) = -normalization;
179 
180  for (int j=0; j<n; j++) {
181  double sum = 0;
182  for (int i=0; i<n; i++) {
183  if (i == j)
184  continue;
185  sum += M(j, i);
186  }
187  M(j, j) = -sum;
188  }
189 
190  return(M);
191 }
192 
193 
194 
195 int main(int argc, char *argv[]) {
196 
197  string header = invocationHeader(argc, argv);
198  parseOptions(argc, argv);
199 
200  AtomicGroup model = createSystem(model_name);
201  AtomicGroup subset = selectAtoms(model, selection);
202 
203  cout << boost::format("Selected %d atoms from %s\n") % subset.size() % model_name;
204  Timer<WallTimer> timer;
205  cerr << "Computing Kirchoff matrix - ";
206  timer.start();
207  Matrix K = kirchoff(subset, cutoff);
208  timer.stop();
209  cerr << "done.\n" << timer << endl;
210 
211 
212  writeAsciiMatrix(prefix + "_K.asc", K, header);
213 
214  boost::tuple<DoubleMatrix, DoubleMatrix, DoubleMatrix> result = svd(K);
215  Matrix U = boost::get<0>(result);
216  Matrix S = boost::get<1>(result);
217  Matrix Vt = boost::get<2>(result);
218 
219  uint n = S.rows();
220 
221  reverseRows(S);
222  reverseColumns(U);
223  reverseRows(Vt);
224 
225  // Write out the LSV (or eigenvectors)
226  writeAsciiMatrix(prefix + "_U.asc", U, header);
227  writeAsciiMatrix(prefix + "_s.asc", S, header);
228 
229  // Now go ahead and compute the pseudo-inverse...
230 
231  // Vt = Vt * diag(1./diag(S))
232  // Remember, Vt is stored col-major but transposed, hence the
233  // somewhat funky indices...
234  //
235  // Note: We have to toss the first term (see Chennubhotla et al (Phys Biol 2(2005:S173-S180))
236  for (uint i=1; i<n; i++) {
237  double s = 1.0/S[i];
238  for (uint j=0; j<n; j++)
239  Vt(i, j) *= s;
240  }
241 
242  Matrix Ki = MMMultiply(Vt, U, true, true);
243  writeAsciiMatrix(prefix + "_Ki.asc", Ki, header);
244 }
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
STL namespace.
AtomicGroup createSystem(const std::string &filename)
Factory function for reading in structure files.
Definition: sfactories.cpp:115
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
void start(void)
Starts the timer.
Definition: timer.hpp:85
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
def svd(traj)
Returns a tuple containing the SVD result for a trajectory, along with the average structure...
Definition: ensembles.py:66
double stop(void)
Stops the timer.
Definition: timer.hpp:94
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Class for tracking time.
Definition: timer.hpp:80
Namespace for most things not already encapsulated within a class.
RealMatrix MMMultiply(const RealMatrix &A, const RealMatrix &B, const bool transa, const bool transb)
Matrix-matrix multiply (using BLAS)
Definition: MatrixOps.cpp:170