ap_ligand_trajectory

ap_ligand_trajectory finds contacts between a ligand molecule and a protein.

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

USAGE:
ap_ligand_trajectory 2kwi.pdb GNP 3.5

where 2kwi.pdb id an input file, GNP the ligand code and 3.5 - 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
#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_trajectory finds contacts between a ligand molecule and a protein.

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

USAGE:
    ap_ligand_trajectory 2kwi.pdb GNP 3.5

where 2kwi.pdb id an input file, GNP the ligand code and 3.5 - contact distance in Angstroms

)";

/** @brief Finds contacts between a ligand molecule and a multimodel-protein.
 *
 * CATEGORIES: core::data::io::Pdb
 * KEYWORDS: PDB input; contact map; ligand
 */
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::map<std::string, int> results; // --- Map used to store results
  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 << " has no " << argv[2] << " residue\n";
      continue;
    }
    for (auto it_resid = strctr->first_residue(); it_resid != strctr->last_residue(); ++it_resid) {
      double d = (*it_resid)->min_distance(*ligand);
      if (d < cutoff) { // --- if this is close enough,
        std::string key = utils::string_format("%3s  %c  %4d",(*it_resid)->residue_type().code3.c_str(),
          (*it_resid)->owner()->id(), (*it_resid)->id());
        if (results.find(key) == results.end()) results[key] = 1;
        else results[key]++;
      }
    }
  }

  // --- print results
  std::cout <<"#resn chain resid counts\n";
  for(const auto & p:results)
    std::cout << p.first<<" "<<utils::string_format("%5d",p.second)<<"\n";
}
../_images/file_icon.png