ap_contact_map_overlap

ap_contact_map_overlap calculates overlap between contact maps calculated for two (or more) structures

USAGE:
ap_contact_map_overlap CA native.pdb models.pdb cutoff 8.0
In this case the program evaluates contact map overlap (measured by Jaccard coefficient) between the native structure
and every model found in models.pdb. 8.0 is the contact distance in Angstroms and 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

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
#include <vector>
#include <algorithm>

#include <core/data/io/Pdb.hh>
#include <core/calc/structural/ContactMap.hh>

#include <utils/exit.hh>

std::string program_info = R"(

ap_contact_map_overlap calculates overlap between contact maps calculated for two (or more) structures

USAGE:
    ap_contact_map_overlap CA native.pdb models.pdb cutoff 8.0
In this case the program evaluates contact map overlap (measured by Jaccard coefficient) between the native structure
and every model found in models.pdb. 8.0 is the contact distance in Angstroms and 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 overlap between contact maps calculated for two (or more) structures
 *
 * CATEGORIES: core::calc::structural::ContactMap
 * KEYWORDS: PDB input; contact map
 */
int main(const int argc, const char *argv[]) {

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

  using namespace core::data::io; // PDB is from this namespace
  using namespace core::data::structural;
  using namespace core::calc::structural;

  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[(argc==5) ? 4 : 3]); // The third/fourth parameter is the contact distance (in Angstroms)

  // --- This is the case when user gave a reference structure (e.g. the native one)
  if (argc == 5) {
    Pdb pdb_native = Pdb(argv[2], filter);
    core::data::structural::Structure_SP reference_structure = pdb_native.create_structure(0);
    ContactMap_SP reference_map = std::make_shared<ContactMap>(*reference_structure, cutoff, selector);

    Pdb models_pdb = Pdb(argv[3], filter);
    ContactMap cmap(*reference_structure, cutoff, selector);
    std::vector<std::pair<core::index2, core::index2>> contacts;
    for (core::index4 i = 0; i < models_pdb.count_models(); ++i) {
      models_pdb.fill_structure(i,*reference_structure);
      ContactMap cmap(*reference_structure, cutoff, selector);
      std::cout << i << " " << reference_map->jaccard_overlap_coefficient(cmap) << "\n";
    }
  } else {
    Pdb models_pdb = Pdb(argv[2], filter);
    core::data::structural::Structure_SP structure = models_pdb.create_structure(0);
    std::vector<std::vector<std::pair<core::index2, core::index2>>> models(models_pdb.count_models());
    for (int i = 0; i < models_pdb.count_models(); ++i) {
      models_pdb.fill_structure(i,*structure);
      ContactMap cmap(*models_pdb.create_structure(i), cutoff, selector);
      cmap.nonempty_indexes(models[i]);
    }

    for (core::index4 i = 1; i < models.size(); ++i) {
      for (core::index4 j = 0; j < i; ++j)
        std::cout << i << " " << j << " " << core::calc::structural::jaccard_overlap_coefficient(models[i], models[j]) << "\n";
    }
  }
}
../_images/file_icon.png