LOOS  v2.3.2
eigenflucc.cpp
1 /*
2  eigenflucc
3 
4  Predict isotropic B-factors from a set of eigenpairs...
5 */
6 
7 
8 /*
9  This file is part of LOOS.
10 
11  LOOS (Lightweight Object-Oriented Structure library)
12  Copyright (c) 2010 Tod D. Romo
13  Department of Biochemistry and Biophysics
14  School of Medicine & Dentistry, University of Rochester
15 
16  This package (LOOS) is free software: you can redistribute it and/or modify
17  it under the terms of the GNU General Public License as published by
18  the Free Software Foundation under version 3 of the License.
19 
20  This package is distributed in the hope that it will be useful,
21  but WITHOUT ANY WARRANTY; without even the implied warranty of
22  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  GNU General Public License for more details.
24 
25  You should have received a copy of the GNU General Public License
26  along with this program. If not, see <http://www.gnu.org/licenses/>.
27 
28 */
29 
30 
31 
32 #include <loos.hpp>
33 #include <boost/format.hpp>
34 #include <boost/program_options.hpp>
35 
36 
37 using namespace std;
38 using namespace loos;
39 
40 namespace po = boost::program_options;
41 
42 
43 // --------------- GLOBALS
44 
45 bool verbose;
46 bool pca_input;
47 vector<uint> modes;
48 string eigvals_name, eigvecs_name, pdb_name;
49 string out_name;
50 double scale;
51 string selection;
52 
53 const double kB = 6.950356e-9; // \AA^{-1} K
54 
55 void fullHelp() {
56  //string msg =
57  cout << "\n"
58  "SYNOPSIS\n"
59  "\n"
60  "Predict isotropic B-factors from a set of eigenpairs\n"
61  "\n"
62  "DESCRIPTION\n"
63  "\n"
64  "Given the results of a network model or a simulation PCA\n"
65  "this tool will calculate the isotropic B-factors from the\n"
66  "eigenpairs. \n"
67  "\n"
68  "There are two modes of output:\n"
69  "\t- A list of B-factors numbered sequentially\n"
70  "\t- An updated PDB file containing the B-factors\n"
71  "\t NOTE: Make sure the same selection string used to\n"
72  "\t compute the ENM is used to ensure the correct\n"
73  "\t mapping of the B-factors.\n"
74  "\n"
75  //
76  "EXAMPLES\n"
77  "\n"
78  "eigenflucc anm_s.asc anm_U.asc > b_factors\n"
79  "\tCompute the B-factors of 'anm' and stream them to the\n"
80  "\tfile 'b_factors'. This outputs a sequential list of\n"
81  "\tof the values (may be more convenient for plotting).\n"
82  "\n"
83  "eigenflucc -p1 model.pdb -s 'name==\"CA\"' anm_s.asc anm_U.asc > b_factors\n"
84  "\tSame as above, but in addition we make a new pdb\n"
85  "\twhere the B-factors are modified based on our result.\n"
86  "\tThe original model.pdb is unaltered, but model-ef.pdb\n"
87  "\twill contain our results. In this case, the selection\n"
88  "\tstring includes all CA's, so they will be updated in\n"
89  "\tthe file output. \n"
90  "\n"
91  "eigenflucc -p1 model.pdb -o1 model_new_b-factors.pdb -S2 -s 'name==\"CA\"' \\\n"
92  " anm_s.asc anm_U.asc > b_factors\n"
93  "\tSame as previous, except the output pdb file is named\n"
94  "\tby the string \"model_new_b-factors.pdb\" and the\n"
95  "\tresults are scale by a factor of 2.\n"
96  "\n"
97  "eigenflucc -m1:3 -P1 pca_s.asc pca_U.asc > b_factors\n"
98  "\tComputes the B-factors from a PCA result. In\n"
99  "\taddtion, only the first 3 modes (or principal\n"
100  "\tcomponents) are used for the calculation. Only\n"
101  "\the sequential list is output (However a new PDB\n"
102  "\tfile can be written if desired).\n"
103  "\n"
104  "\n";
105 }
106 
107 
108 void parseArgs(int argc, char *argv[]) {
109 
110  try {
111  po::options_description generic("Allowed options");
112  generic.add_options()
113  ("help", "Produce this help message")
114  ("fullhelp", "Get extended help")
115  ("verbose,v", po::value<bool>(&verbose)->default_value(false), "Verbose output")
116  ("selection,s", po::value<string>(&selection)->default_value("name == 'CA'"), "Selection used to make the ENM (only when altering a PDB)")
117  ("pdb,p", po::value<string>(&pdb_name), "Alter the B-factors in a PDB")
118  ("outpdb,o", po::value<string>(&out_name), "Filename to output PDB to")
119  ("modes,m", po::value< vector<string> >(), "Modes to use (default is all)")
120  ("scale,S", po::value<double>(&scale)->default_value(1.0), "Scaling factor to apply to eigenvalues")
121  ("pca,P", po::value<bool>(&pca_input)->default_value(false), "Eigenpairs come from PCA, not ENM");
122 
123  po::options_description hidden("Hidden options");
124  hidden.add_options()
125  ("eigvals", po::value<string>(&eigvals_name), "Eigenvalues filename")
126  ("eigvecs", po::value<string>(&eigvecs_name), "Eigenvectors filename");
127 
128 
129  po::options_description command_line;
130  command_line.add(generic).add(hidden);
131 
132  po::positional_options_description p;
133  p.add("eigvals", 1);
134  p.add("eigvecs", 1);
135 
136  po::variables_map vm;
137  po::store(po::command_line_parser(argc, argv).
138  options(command_line).positional(p).run(), vm);
139  po::notify(vm);
140 
141  if (vm.count("help") || vm.count("fullhelp") || !(vm.count("eigvals") && vm.count("eigvecs"))) {
142  cout << "Usage- " << argv[0] << " [options] eigenvalues eigenvectors\n";
143  cout << generic;
144  if (vm.count("fullhelp"))
145  fullHelp();
146  exit(0);
147  }
148 
149  if (vm.count("modes")) {
150  vector<string> mode_list = vm["modes"].as< vector<string> >();
151  modes = parseRangeList<uint>(mode_list);
152  }
153 
154  if (!pdb_name.empty())
155  if (out_name.empty())
156  out_name = findBaseName(pdb_name) + "-ef.pdb";
157 
158 
159  }
160  catch(exception& e) {
161  cerr << "Error - " << e.what() << endl;
162  exit(-1);
163  }
164 
165 }
166 
167 
168 
169 int main(int argc, char *argv[]) {
170 
171  string hdr = invocationHeader(argc, argv);
172  cout << "# " << hdr << endl;
173  parseArgs(argc, argv);
174 
175  DoubleMatrix eigvals;
176  readAsciiMatrix(eigvals_name, eigvals);
177 
178  DoubleMatrix eigvecs;
179  readAsciiMatrix(eigvecs_name, eigvecs);
180 
181  if (modes.empty())
182  for (uint i=0; i<eigvals.rows(); ++i)
183  modes.push_back(i);
184 
185  uint n = modes.size();
186  uint m = eigvecs.rows();
187 
188  // V = m x n matrix of eigenvectors
189  DoubleMatrix V(m, n);
190  for (uint i=0; i<n; ++i)
191  for (uint j=0; j<m; ++j)
192  V(j, i) = eigvecs(j, modes[i]);
193 
194  // Make a copy and scale by the appropriate eigenvalue,
195  // remembering that eigenvalues are inverted for ENM,
196  // or squaring and not inverting in the case of PCA
197  DoubleMatrix VS = V.copy();
198  for (uint i=0; i<n; ++i) {
199  double e = eigvals[modes[i]];
200  double s = (pca_input) ? (scale * e * e): (scale / e);
201  for (uint j=0; j<m; ++j)
202  VS(j, i) *= s;
203  }
204 
205  // U = Covariance matrix
206  DoubleMatrix U = loos::Math::MMMultiply(VS,V,false,true);
207 
208  // B-factors come from trace of diagonal 3x3 superblocks
209  vector<double> B;
210  double prefactor = 8.0 * M_PI * M_PI / 3.0;
211  for (uint i=0; i<m; i += 3) {
212  double b = prefactor * (U(i,i) + U(i+1, i+1) + U(i+2, i+2));
213  B.push_back(b);
214  cout << boost::format("%-8d %g\n") % (i/3) % b;
215  }
216 
217  if (!pdb_name.empty()) {
218  AtomicGroup model = createSystem(pdb_name);
219  AtomicGroup subset = selectAtoms(model, selection);
220 
221  if ((unsigned int)subset.size() != B.size()) {
222  cerr << boost::format("Error- selection has %d atoms, but %d were expected.\n") % subset.size() % B.size();
223  exit(-10);
224  }
225 
226  for (uint i=0; i<B.size(); ++i)
227  subset[i]->bfactor(B[i]);
228 
229  PDB pdb = PDB::fromAtomicGroup(model);
230  pdb.remarks().add(hdr);
231  ofstream ofs(out_name.c_str());
232  ofs << pdb;
233  }
234 
235 }
236 
Remarks & remarks(void)
Accessor for the remarks object...
Definition: pdb.cpp:547
void add(const std::string s)
Add a remark.
Definition: pdb_remarks.cpp:52
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
PDB reading/writing class.
Definition: pdb.hpp:69
AtomicGroup selectAtoms(const AtomicGroup &source, const std::string selection)
Applies a string-based selection to an atomic group...
Definition: utils.cpp:195
Math::Matrix< T, P, S > readAsciiMatrix(std::istream &is)
Read in a matrix from a stream returning a newly created matrix.
Definition: MatrixRead.hpp:72
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
std::string findBaseName(const std::string &s)
Pull off the file name extension (if present)
Definition: utils.cpp:68
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
Matrix< T, OrderPolicy, StoragePolicy > copy(void) const
Deep copy...
Definition: MatrixImpl.hpp:170