LOOS  v2.3.2
blob_contact.cpp
1 /*
2  blob_contact.cpp
3 
4  Find residues within a given distance of a blob
5 
6 */
7 
8 /*
9  This file is part of LOOS.
10 
11  LOOS (Lightweight Object-Oriented Structure library)
12  Copyright (c) 2012, Tod D. Romo, Alan Grossfield
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 #include <loos.hpp>
32 
33 #include <DensityGrid.hpp>
34 #include <DensityTools.hpp>
35 
36 using namespace std;
37 using namespace loos;
38 using namespace loos::DensityTools;
39 
40 GCoord findMinCoord(const vector<GCoord>& coords) {
41  GCoord min(numeric_limits<double>::max(),
42  numeric_limits<double>::max(),
43  numeric_limits<double>::max());
44 
45  for (vector<GCoord>::const_iterator i = coords.begin(); i != coords.end(); ++i)
46  for (uint j=0; j<3; ++j)
47  if ((*i)[j] < min[j])
48  min[j] = (*i)[j];
49 
50  return(min);
51 }
52 
53 GCoord findMaxCoord(const vector<GCoord>& coords) {
54  GCoord max(numeric_limits<double>::min(),
55  numeric_limits<double>::min(),
56  numeric_limits<double>::min());
57 
58  for (vector<GCoord>::const_iterator i = coords.begin(); i != coords.end(); ++i)
59  for (uint j=0; j<3; ++j)
60  if ((*i)[j] > max[j])
61  max[j] = (*i)[j];
62 
63  return(max);
64 }
65 
66 
67 
68 vector<GCoord> findBlobCoords(DensityGrid<int>& grid, const int blobid) {
69  DensityGridpoint griddims = grid.gridDims();
70 
71  vector<GCoord> coords;
72  for (int k=0; k<griddims.z(); ++k)
73  for (int j=0; j<griddims.y(); ++j)
74  for (int i=0; i<griddims.x(); ++i) {
75  DensityGridpoint grid_coord(i, j, k);
76  if (grid(grid_coord) == blobid)
77  coords.push_back(grid.gridToWorld(grid_coord));
78  }
79 
80  return(coords);
81 }
82 
83 
84 vector<int> findResiduesNearBlob(const vector<GCoord>& blob, const vector<AtomicGroup>& residues, const double threshold) {
85 
86  double thresh2 = threshold * threshold;
87  vector<int> nearby(residues.size(), 0);
88 
89  for (uint k=0; k<residues.size(); ++k) {
90  bool flag = false;
91  for (uint j=0; j<residues[k].size() && !flag; ++j) {
92  GCoord c = residues[k][j]->coords();
93 
94  for (uint i=0; i<blob.size(); ++i) {
95  double d = c.distance2(blob[i]);
96  if (d <= thresh2) {
97  flag = true;
98  break;
99  }
100  }
101  }
102  nearby[k] = flag;
103  }
104 
105  return(nearby);
106 }
107 
108 
109 vector<double> calculatePercentageContacts(const RealMatrix& M) {
110  vector<long> sum(M.cols()-1, 0);
111 
112  for (uint i=1; i<M.cols(); ++i)
113  for (uint j=0; j<M.rows(); ++j)
114  sum[i-1] += M(j, i);
115 
116 
117  vector<double> occ(M.cols()-1, 0.0);
118  for (uint i=0; i<M.cols()-1; ++i)
119  occ[i] = static_cast<double>(sum[i]) / M.rows();
120 
121  return(occ);
122 }
123 
124 
125 
126 int main(int argc, char *argv[]) {
127  if (argc != 7) {
128  cerr << "Usage- blob_contact model traj selection skip blobid distance <grid >out.asc 2>report.txt\n";
129  cerr << "Note: requires an grid with blob ids (i.e. output from blobid)\n";
130  exit(-1);
131  }
132 
133  string hdr = invocationHeader(argc, argv);
134  int k = 1;
135  AtomicGroup model = createSystem(argv[k++]);
136  pTraj traj = createTrajectory(argv[k++], model);
137  AtomicGroup residue_subset = selectAtoms(model, argv[k++]);
138  vector<AtomicGroup> residues = residue_subset.splitByResidue();
139  uint skip = strtoul(argv[k++], 0, 10);
140  int blobid = strtol(argv[k++], 0, 10);
141  double distance = strtod(argv[k++], 0);
142 
143  DensityGrid<int> grid;
144  cin >> grid;
145  vector<GCoord> blob = findBlobCoords(grid, blobid);
146  if (blob.empty()) {
147  cerr << boost::format("ERROR- no voxels found for blobid %d\n") % blobid;
148  exit(-10);
149  }
150 
151  GCoord blobmin = findMinCoord(blob);
152  GCoord blobmax = findMaxCoord(blob);
153 
154 
155  ostringstream shdr;
156  shdr << hdr << endl;
157  shdr << boost::format("# Blob bounding box is %s x %s\n") % blobmin % blobmax;
158  shdr << boost::format("# Blob has %d voxels\n") % blob.size();
159  shdr << "# Residue list...\n";
160  for (uint i=0; i<residues.size(); ++i)
161  shdr << boost::format("# %d : %d %d %s %s\n")
162  % i
163  % residues[i][0]->id()
164  % residues[i][0]->resid()
165  % residues[i][0]->resname()
166  % residues[i][0]->segid();
167 
168  uint n = traj->nframes();
169  if (skip > 0) {
170  n -= skip;
171  traj->readFrame(skip-1);
172  }
173 
174  RealMatrix M(n, residues.size()+1);
175  uint t = 0;
176 
177  while (traj->readFrame()) {
178  traj->updateGroupCoords(model);
179  M(t, 0) = t + skip;
180  vector<int> nearby = findResiduesNearBlob(blob, residues, distance);
181  for (uint i=0; i<nearby.size(); ++i)
182  M(t, i+1) = nearby[i];
183  ++t;
184  }
185 
186 
187  writeAsciiMatrix(cout, M, shdr.str());
188 
189 
190  cerr << "# " << hdr << endl;
191  cerr << "# gnuplot-col\tresid\tatomid\tfractional contact\n";
192  vector<double> fraction = calculatePercentageContacts(M);
193  for (uint i=0; i<residues.size(); ++i) {
194  cerr << boost::format("%d\t%d\t%d\t%f\n")
195  % (i+2)
196  % residues[i][0]->resid()
197  % residues[i][0]->id()
198  % fraction[i];
199  }
200 }
std::vector< AtomicGroup > splitByResidue(void) const
Returns a vector of AtomicGroups, each comprising a single residue.
Basic 3-D coordinates class.
Definition: Coord.hpp:36
Simple matrix template class using policy classes to determine behavior.
Definition: MatrixImpl.hpp:53
STL namespace.
double distance2(const Coord< T > &o) const
Distance squared between two coordinates.
Definition: Coord.hpp:429
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
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
loos::GCoord gridToWorld(const DensityGridpoint &v) const
Converts grid coords to real-space (world) coords.
Class for handling groups of Atoms (pAtoms, actually)
Definition: AtomicGroup.hpp:87
Namespace for most things not already encapsulated within a class.