ap_stacking_interactions

Finds stacking interactions in a given PDB file.

The program prints relative orientation (three Euler angles and the distance) between any two aromatic rings found in amino acid side chains that are closer than 5.0. The rings are assumed to be flat rigid moieties.

USAGE:
ap_stacking_interactions 5edw.pdb

Keywords:

Categories:

  • core::data::io::Pdb; core::calc::structural::transformations::Rototranslation

Output files:

Program source:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
#include <iostream>
#include <numeric>
#include <core/data/structural/structure_selectors.hh>
#include <core/data/io/Pdb.hh>
#include <core/calc/structural/transformations/transformation_utils.hh>
#include <utils/exit.hh>

std::string program_info = R"(

Finds stacking interactions in a given PDB file.

The program prints relative orientation (three Euler angles and the distance) between any two aromatic rings
found in amino acid side chains that are closer than 5.0. The rings are assumed to be flat rigid moieties.
USAGE:
    ap_stacking_interactions 5edw.pdb

)";

using namespace core::data::io; // --- Pdb reader and PdbLineFilter lives there
using core::chemical::Monomer;
using namespace core::calc::structural::transformations; // --- transformations

const float DISTANCE_CUTOFF = 5.0;

/// ---------- The map defines three atoms to build a local coordinate system for each ring
std::map <std::string, std::vector<std::string>> reference_atoms {
  { "HIS", {" CG "," CE1", " NE2"} },
  { "TRP", {" CG "," CZ2", " CE3"} },
  { "TYR", {" CG "," CE1", " CE2"} },
  { "PHE", {" CG "," CE1", " CE2"} }
};

/** @brief Finds stacking interactions in a given PDB file.
 *
 * CATEGORIES: core::data::io::Pdb; core::calc::structural::transformations::Rototranslation
 * KEYWORDS:   PDB input; PDB line filter; stacking interactions; rototranslation
 */
int main(const int argc, const char *argv[]) {

  if(argc < 2) utils::exit_OK_with_message(program_info); // --- complain about missing program parameter

  // --- Read a PDB file given as an argument to this program
  Pdb reader(argv[1], // --- input PDB file
    all_true(is_not_water, is_not_alternative, is_not_hydrogen,
      invert_filter(is_bb)), // --- Inverted backbone selector  reads only side chains
    true); // --- yes, read header

  core::data::structural::IsAromaticAA if_aromatic;
  core::data::structural::Structure_SP s = reader.create_structure(0);

  std::cout << "#alpha   beta     gamma    r    c   id aa  c   id  aa  PDB\n";

  for (auto r_it_a = s->first_residue(); r_it_a != s->last_residue(); ++r_it_a) {
    if (!if_aromatic(**r_it_a)) continue;
    Vec3 cm_a;
    for (auto vi:(**r_it_a)) cm_a += *vi;
    cm_a /= (*r_it_a)->size();
    for (auto r_it_b = s->first_residue(); r_it_b != r_it_a; ++r_it_b) {
      if (!if_aromatic(**r_it_b)) continue;
      Vec3 cm_b;
      for (auto vi:(**r_it_b)) cm_b += *vi;
      cm_b /= (*r_it_b)->size();
      if (cm_a.distance_to(cm_b) < DISTANCE_CUTOFF) {
        try {
          Rototranslation_SP rt_a = local_coordinates_three_atoms(**r_it_a, reference_atoms[(*r_it_a)->residue_type().code3]);
          Rototranslation_SP rt_b = local_coordinates_three_atoms(**r_it_a, reference_atoms[(*r_it_b)->residue_type().code3]);
          Vec3 e = euler_angles(*rt_a, *rt_b);
          e *= 180.0 / 3.14159;
          std::cout << utils::string_format("%7.2f %7.2f %7.2f %5.3f  %c %4d %3s %c %4d %3s  %s\n",e.x,e.y,e.z,cm_a.distance_to(cm_b),
            (*r_it_a)->owner()->id(), (*r_it_a)->id(),  (*r_it_a)->residue_type().code3.c_str(),
            (*r_it_a)->owner()->id(), (*r_it_a)->id(),  (*r_it_a)->residue_type().code3.c_str(), s->code().c_str());
        } catch (utils::exceptions::AtomNotFound e) {
          std::cerr << "Missing atom: " << e.missing_atom_name << " in " << e.source_residue_name << "\n";
        }
      }
    }
  }
}
../_images/file_icon.png