ap_contact_map

ap_contact_map calculates a contact map for a given protein structure

If a multi-model PDB file was given, the program prints contact count observed in all models

USAGE:
ap_contact_map  CA 2kwi.pdb.pdb 4.5

where 2kwi.pdb is the input file and 4.5 the contact distance in Angstroms. CA defines the contact map type; allowed options are: CA CB and SC for Calpha, C-beta and all atom side chain, respectively

Categories:

  • core::calc::structural::ContactMap

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

#include <core/index.hh>
#include <core/data/io/Pdb.hh>
#include <core/calc/structural/ContactMap.hh>
#include <core/data/structural/structure_selectors.hh>
#include <utils/exit.hh>

std::string program_info = R"(

ap_contact_map calculates a contact map for a given protein structure

If a multi-model PDB file was given, the program prints contact count observed in all models

USAGE:
    ap_contact_map  CA 2kwi.pdb.pdb 4.5

where 2kwi.pdb is the input file and 4.5 the contact distance in Angstroms. CA defines the contact map type;
allowed options are: CA CB and SC for Calpha, C-beta and all atom side chain, respectively

)";

/** @brief Calculates a contact map for a given protein structure
 *
 * CATEGORIES: core::calc::structural::ContactMap
 * KEYWORDS: PDB input; contact map
 * IMG: ap_contact_map.png
 * IMG_ALT: Contact map calculated for 2KWI protein structure solved by NMR. The 2KWI deposit holds 51 models, the color scale shows how popular is a given contact among the models
 */
int main(const int argc, const char* argv[]) {

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

  core::data::structural::AtomSelector_SP selector = std::make_shared<core::data::structural::IsSC>();
  core::data::io::PdbLineFilter filter = core::data::io::is_not_water;
  if (std::strcmp(argv[1],"CA")==0 ) {
    selector = std::make_shared<core::data::structural::IsNamedAtom>(" CA ");
    core::data::io::PdbLineFilter filter = core::data::io::is_ca;
  }
  if (std::strcmp(argv[1],"CB")==0) {
    core::data::io::PdbLineFilter filter = core::data::io::is_cb;
    selector = std::make_shared<core::data::structural::IsNamedAtom>(" CB ");
  }

  double cutoff = utils::from_string<double>(argv[3]); // The third parameter is the contact distance (in Angstroms)
  core::data::io::Pdb reader(argv[2],filter); // --- file name (PDB format, may be gzip-ped)

  core::data::structural::Structure_SP structure = reader.create_structure(0);
  core::calc::structural::ContactMap cmap(*structure,cutoff,selector);
  for (int i_model = 1; i_model < reader.count_models(); ++i_model) {
    reader.fill_structure(i_model,*structure);
    cmap.add(*structure);
  }

  std::vector<std::pair<core::index2,core::index2>> contacts;
  cmap.nonempty_indexes(contacts);

  for(const std::pair<core::index2,core::index2> ij : contacts) {
    core::index2 i_res = ij.first;
    core::index2 j_res = ij.second;
    std::cout << utils::string_format("%4d %c %4d%c %4d %c %4d%c %d\n", i_res, cmap.residue_index(i_res).chain_id,
      cmap.residue_index(i_res).residue_id, cmap.residue_index(i_res).i_code,
      j_res, cmap.residue_index(j_res).chain_id,
      cmap.residue_index(j_res).residue_id, cmap.residue_index(j_res).i_code,
      cmap.at(i_res, j_res, 0));

  }
}
../_images/ap_contact_map.png