LOOS  v2.3.2
near_blobs.cpp
1 /*
2  near_blobs.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 
41 typedef vector<GCoord> vCoords;
42 typedef vector<vCoords> vvCoords;
43 typedef vector<AtomicGroup> vGroup;
44 
45 
46 
47 vvCoords separateBlobs(const DensityGrid<int>& grid) {
48 
49  int max_blobid = 0;
50 
51  for (long i = 0; i < grid.size(); ++i)
52  if (grid(i) > max_blobid)
53  max_blobid = grid(i);
54 
55  vvCoords blobs(max_blobid);
56 
57  DensityGridpoint dims = grid.gridDims();
58  for (int k=0; k<dims.z(); ++k)
59  for (int j=0; j<dims.y(); ++j)
60  for (int i=0; i<dims.x(); ++i) {
61  DensityGridpoint p(i, j, k);
62  if (grid(p)) {
63  int id = grid(p);
64  GCoord c = grid.gridToWorld(p);
65  blobs[id-1].push_back(c);
66  }
67  }
68 
69  return(blobs);
70 }
71 
72 
73 
74 vector<uint> findBlobsNearResidue(const vvCoords& blobs, const AtomicGroup& residue, const double dist) {
75  vector<uint> blobids;
76 
77  double d2 = dist * dist;
78  for (uint k=0; k<blobs.size(); ++k) {
79  bool flag = true;
80 
81  for (uint j=0; j<residue.size() && flag; ++j) {
82  GCoord c = residue[j]->coords();
83 
84  for (uint i=0; i<blobs[k].size() && flag; ++i)
85  if (c.distance2(blobs[k][i]) <= d2)
86  flag = false;
87 
88  }
89 
90  if (!flag)
91  blobids.push_back(k);
92 
93  }
94 
95  return(blobids);
96 }
97 
98 
99 
100 
101 int main(int argc, char *argv[]) {
102 
103  if (argc != 4) {
104  cerr << "Usage- near_blobs model selection distance <blobs.grid\n";
105  cerr << "NOTE: grid must have IDs (i.e. output from blobid)\n";
106  exit(-1);
107  }
108 
109  string hdr = invocationHeader(argc, argv);
110 
111  int k = 1;
112  AtomicGroup model = createSystem(argv[k++]);
113  AtomicGroup subset = selectAtoms(model, argv[k++]);
114  double distance = strtod(argv[k++], 0);
115 
116  DensityGrid<int> the_grid;
117  cin >> the_grid;
118 
119  vvCoords blobs = separateBlobs(the_grid);
120  vGroup residues = subset.splitByResidue();
121 
122  cout << "# " << hdr << endl;
123  cout << "# Atomid Resid Resname Segid Bloblist...\n";
124  for (uint i=0; i<residues.size(); ++i) {
125  vector<uint> ids = findBlobsNearResidue(blobs, residues[i], distance);
126  if (ids.size() == 0)
127  continue;
128  cout << boost::format("%d\t%d\t%s\t%s\t")
129  % residues[i][0]->id()
130  % residues[i][0]->resid()
131  % residues[i][0]->resname()
132  % residues[i][0]->segid();
133  for (uint j=0; j<ids.size(); ++j)
134  cout << ids[j]+1 << (j == ids.size()-1 ? "" : ",");
135  cout << endl;
136  }
137 }
std::vector< AtomicGroup > splitByResidue(void) const
Returns a vector of AtomicGroups, each comprising a single residue.
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
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.