ap_ligand_contacts

ap_ligand_contacts finds contacts between a ligand molecule and a protein.

It reads a multi-model PDB file and detects contacts in everyone of them. The output provides the interacting residues (name and residueId) along with the number of observations for this contact

USAGE:
ap_ligand_contacts 5edw.pdb TTP 7.0

where 5edw.pdb id an input file, TTP the ligand code and 7.0 - contact distance in Angstroms

Categories:

  • core::data::io::Pdb

Input files:

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
#include <iostream>
#include <map>

#include <core/data/io/Pdb.hh>
#include <utils/string_utils.hh>
#include <utils/exit.hh>

std::string program_info = R"(

ap_ligand_contacts finds contacts between a ligand molecule and a protein.

It reads a multi-model PDB file and detects contacts in everyone of them. The output provides the interacting residues
(name and residueId) along with the number of observations for this contact

USAGE:
    ap_ligand_contacts 5edw.pdb TTP 7.0

where 5edw.pdb id an input file, TTP the ligand code and 7.0 - contact distance in Angstroms

)";

/** @brief Finds contacts between a ligand molecule and a protein.
 *
 * CATEGORIES: core::data::io::Pdb
 * KEYWORDS: PDB input; contact map; ligand
 * IMG: ap_ligand_contacts.png
 * IMG_ALT: Contacts found between 5EDW protein and its ligand TTP
 */
int main(const int argc, const char* argv[]) {

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

  core::data::io::Pdb reader(argv[1]); // --- file name (PDB format, may be gzip-ped)

  const std::string code(argv[2]);    // --- The ligand code is the second parameter of the program
  double cutoff = utils::from_string<double>(argv[3]); // The third parameter is the contact distance (in Angstroms)

  std::cout << " ---- ligand ---- | --------- partner -------- | distance\n";
  std::cout << "c  res  id atname |  c  res  id   type  atname | in Angstrom\n";

  for (size_t i = 0; i < reader.count_models(); ++i) { // --- Iterate over all models in the input file
    core::data::structural::Structure_SP strctr = reader.create_structure(i);

    // --- Here we use a standard <code>find_if</code> algorithm to find the ligand residue by its 3-letter code
    auto ligand = std::find_if(strctr->first_residue(), strctr->last_residue(), [&code](core::data::structural::Residue_SP res) {return (res->residue_type().code3==code);});
    if(ligand== strctr->last_residue()) { // --- If no ligand - print a message and take next structure
      std::cerr << "Model " << i << " of " << argv[1] << " has no " << argv[2] << " residue\n";
      continue;
    }
    if (reader.count_models() > 1) std::cout << "# Model " << i + 1 << "\n";
    for (auto it_resid = strctr->first_residue(); it_resid != strctr->last_residue(); ++it_resid) {
      if(*it_resid == *ligand) continue;
      double d = (*it_resid)->min_distance(*ligand);
      if (d < cutoff) { // --- if this is close enough,
        for(auto const & ligand_atom : **ligand) {
          for(auto const & other_atom : **it_resid) {
            if(ligand_atom->distance_to(*other_atom) <= cutoff)
              std::cout << utils::string_format("%c  %3s %4d %4s     %c  %3s %4d %6s %4s   %6.3f\n",(**ligand).owner()->id(),
                (**ligand).residue_type().code3.c_str(), (**ligand).id(), ligand_atom->atom_name().c_str(),
                (**it_resid).owner()->id(),
                (**it_resid).residue_type().code3.c_str(), (**it_resid).id(),
                core::chemical::monomer_type_name((**it_resid).residue_type()).c_str(), other_atom->atom_name().c_str(),
                ligand_atom->distance_to(*other_atom));
          }
        }
      }
    }
  }
}
../_images/ap_ligand_contacts.png