LOOS  v2.3.2
fidpick.cpp
1 /*
2  Picks fiducial structures for a structural histogram, a la Lyman &
3  Zuckerman, Biophys J (2006) 91:164-172
4 
5 
6  Usage- fidpick model trajectory selection output-name cutoff [seed]
7 
8 
9  Note: This is the older method of partitioning based on a distance
10  cutoff only, rather than only taking the closest N structures (used
11  in the decorr-time and Neff methods
12 */
13 
14 
15 /*
16 
17  This file is part of LOOS.
18 
19  LOOS (Lightweight Object-Oriented Structure library)
20  Copyright (c) 2010, Tod D. Romo
21  Department of Biochemistry and Biophysics
22  School of Medicine & Dentistry, University of Rochester
23 
24  This package (LOOS) is free software: you can redistribute it and/or modify
25  it under the terms of the GNU General Public License as published by
26  the Free Software Foundation under version 3 of the License.
27 
28  This package is distributed in the hope that it will be useful,
29  but WITHOUT ANY WARRANTY; without even the implied warranty of
30  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31  GNU General Public License for more details.
32 
33  You should have received a copy of the GNU General Public License
34  along with this program. If not, see <http://www.gnu.org/licenses/>.
35 */
36 
37 
38 
39 
40 #include <loos.hpp>
41 
42 using namespace std;
43 using namespace loos;
44 
45 typedef vector<AtomicGroup> vGroup;
46 
47 
48 
49 
50 string fullHelpMessage(void) {
51  string msg =
52  "\n"
53  "SYNOPSIS\n"
54  "\tPick fiducial structures for a structural histogram using a distance cutoff\n"
55  "\n"
56  "DESCRIPTION\n"
57  "\n"
58  "\tThis tool implements the older method of constructing a structural histogram\n"
59  "where bins are defined in terms of the closest N structures to a randomly picked\n"
60  "fiducial. See Lyman and Zuckerman, Biophys J (2006) 91:164-72 for more information.\n"
61  "\n"
62  "EXAMPLES\n"
63  "\n"
64  "\tfidpick model.pdb simulation.dcd all 'name == \"CA\" fiducials 5.0 >>fiducials.asc\n"
65  "This example uses all alpha-carbons, assigns bins based on a distance cutoff of 5.0 angstroms\n"
66  "and writes the fiducials to fiducials.pdb and fiducials.dcd. A log of the selections\n"
67  "is stored in fiducials.asc\n"
68  "\n"
69  "SEE ALSO\n"
70  "\tsortfids\n";
71 
72  return(msg);
73 }
74 
75 
76 vector<uint> findFreeFrames(const vector<int>& map) {
77  vector<uint> indices;
78 
79  for (uint i=0; i<map.size(); ++i)
80  if (map[i] < 0)
81  indices.push_back(i);
82 
83  return(indices);
84 }
85 
86 
87 
88 int main(int argc, char *argv[]) {
89  string hdr = invocationHeader(argc, argv);
90 
91  if (argc < 7 || argc > 8) {
92  cerr << "Usage - fidpick model trajectory range|all selection output-name cutoff [seed]\n";
93  exit(-1);
94  }
95 
96  int opti = 1;
97  AtomicGroup model = createSystem(argv[opti++]);
98  model.clearBonds();
99 
100  pTraj traj = createTrajectory(argv[opti++], model);
101  string range(argv[opti++]);
102  string selection(argv[opti++]);
103  AtomicGroup subset = selectAtoms(model, selection);
104  string outname(argv[opti++]);
105  double cutoff = strtod(argv[opti++], 0);
106 
107  long seed;
108  if (opti < argc) {
109  seed = strtol(argv[opti++], 0, 10);
110  rng_singleton().seed(static_cast<unsigned int>(seed)); // This addresses a bug in earlier BOOSTs
111  } else
112  randomSeedRNG();
113 
114  vector<uint> frames;
115  if (range == "all")
116  for (uint i=0; i<traj->nframes(); ++i)
117  frames.push_back(i);
118  else
119  frames = parseRangeList<uint>(range);
120 
121 
122  boost::uniform_real<> rmap;
123  boost::variate_generator< base_generator_type&, boost::uniform_real<> > rng(rng_singleton(), rmap);
124 
125  vGroup fiducials;
126  vector<int> assignments(frames.size(), -1);
127 
128  cout << "Frames picked:\n";
129 
130  vector<uint> possible_frames = findFreeFrames(assignments);
131  while (! possible_frames.empty()) {
132  uint pick = possible_frames[static_cast<uint>(floor(possible_frames.size() * rng()))];
133 
134  traj->readFrame(frames[pick]);
135  traj->updateGroupCoords(model);
136 
137  AtomicGroup fiducial = subset.copy();
138  fiducial.centerAtOrigin();
139  uint myid = fiducials.size();
140  if (assignments[pick] >= 0) {
141  cerr << "INTERNAL ERROR - " << pick << " pick was already assigned to " << assignments[pick] << endl;
142  exit(-99);
143  }
144 
145  fiducials.push_back(fiducial);
146  assignments[pick] = myid;
147 
148  uint cluster_size = 0;
149  for (uint i = 0; i<assignments.size(); ++i) {
150  if (assignments[i] >= 0 || i == pick)
151  continue;
152  traj->readFrame(frames[i]);
153  traj->updateGroupCoords(model);
154  subset.centerAtOrigin();
155  subset.alignOnto(fiducial);
156  double d = subset.rmsd(fiducial);
157  if (d < cutoff) {
158  assignments[i] = myid;
159  ++cluster_size;
160  }
161  }
162 
163  cout << "\t" << frames[pick] << "\t" << cluster_size << endl;
164 
165  possible_frames = findFreeFrames(assignments);
166  }
167 
168  cerr << "Done!\nWrote " << fiducials.size() << " fiducials to " << outname << endl;
169 
170  DCDWriter dcd(outname + ".dcd", fiducials, hdr);
171 
172  PDB pdb = PDB::fromAtomicGroup(fiducials[0]);
173  pdb.remarks().add(hdr);
174  string fname = outname + ".pdb";
175  ofstream ofs(fname.c_str());
176  ofs << pdb;
177 
178 }
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
STL namespace.
greal rmsd(const AtomicGroup &)
Compute the RMSD between two groups.
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
A very lightweight class for writing simple DCDs.
Definition: dcdwriter.hpp:48
uint randomSeedRNG(void)
Randomly seeds the RNG.
GMatrix alignOnto(const AtomicGroup &)
Superimposes the current group onto the passed group.
Definition: AG_linalg.cpp:196
base_generator_type & rng_singleton(void)
Suite-wide random number generator singleton.
AtomicGroup copy(void) const
Creates a deep copy of this group.
Definition: AtomicGroup.cpp:56
void clearBonds(void)
Remove any bonding information present in contained atoms.
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.
GCoord centerAtOrigin(void)
Translates the group so that the centroid is at the origin.