ex_peptide_hydrogen

ex_peptide_hydrogen reconstructs peptide hydrogen atoms using BioShell algorithm, where amide H is placed in reference to its N atom. Resulting coordinates are printed on the screen. The program also computes the amide-H positions using DSSP approach and calculates the average error (in Angstroms) between the two methods.

USAGE:
ex_peptide_hydrogen 5edw.pdb

Keywords:

Categories:

  • core::calc::structural::peptide_hydrogen()

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 <core/data/io/Pdb.hh>
#include <core/data/structural/structure_selectors.hh>
#include <core/calc/structural/BackboneHBondGeometry.hh>
#include <utils/exit.hh>

using namespace core::data::structural;
using namespace core::data::io;
using namespace core::data::basic;

std::string program_info = R"(

ex_peptide_hydrogen reconstructs peptide hydrogen atoms using BioShell algorithm, 
where amide H is placed in reference to its N atom. Resulting coordinates are printed
on the screen. The program also computes the amide-H positions using DSSP approach
and calculates the average error (in Angstroms) between the two methods.

USAGE:
    ex_peptide_hydrogen 5edw.pdb

)";

/** @brief Reconstructs peptide hydrogen atoms using two methods and compares the error between them.
 * CATEGORIES: core::calc::structural::peptide_hydrogen()
 * KEYWORDS:   PDB input; hydrogen reconstruction
 */
int main(const int argc, const char* argv[]) {
  if(argc < 2) utils::exit_OK_with_message(program_info); // --- complain about missing program parameter
  core::data::io::Pdb reader(argv[1], is_not_alternative); // --- Read in a PDB file
  core::data::structural::Structure_SP strctr = reader.create_structure(0); // --- create a Structure object from the first model

  core::real err = 0, n = 0;
  auto res_it = ++(strctr->first_residue()); // --- The residue being reconstructed
  auto prev_res_it = strctr->first_residue(); // --- preceding residue
  for (; res_it != strctr->last_residue(); ++res_it) {
    if (((*prev_res_it)->residue_type().parent_id > 20) || ((*res_it)->residue_type().parent_id > 20)) continue; // --- its not an amino acid
    try {
      // --- Rebuild the peptide hydrogen in a residue pointed by res_it iterator.
      // --- This method actually adds the newly created H atom to the residue
      core::calc::structural::peptide_hydrogen(*prev_res_it, *res_it);
      auto new_H = (*res_it)->find_atom_safe(" H  ");
      // --- Here we reconstruct amide H knowing the relevant atoms, but now according to the DSSP approach. Resulting H is not inserted
      auto prev_O = (*prev_res_it)->find_atom_safe(" O  ");
      auto prev_C = (*prev_res_it)->find_atom_safe(" C  ");
      auto this_N = (*res_it)->find_atom_safe(" N  ");
      PdbAtom other_H;
      core::calc::structural::peptide_hydrogen_dssp(*prev_C, *prev_O, *this_N, other_H);
      err += other_H.distance_to(*new_H);
      n++;
      for (const auto &atom : **res_it)
        std::cout << atom->to_pdb_line() << "\n"; //-- print all atoms in the current residue (in PDB format)
      ++prev_res_it; // --- advance one of the iterators by one residue; the other iterator is advanced by the loop
    } catch (utils::exceptions::AtomNotFound e) {
      std::cerr << e.what() << "\n";
      ++prev_res_it;
    }
  }
  std::cout << "# difference between two methods: " << err / n << "\n";
}

../_images/file_icon.png