ex_chi_correlation

Calculates Chi dihedral angles of amino acid side chains measured in two different protein structures.

The program reads two homologous protein structures. The proteins are assumed to be already aligned and Chi angles are calculated for every pair of identical positions.

USAGE:
ex_chi_correlation ./test_inputs/2fd2.pdb ./test_inputs/5fd1.pdb

Keywords:

Categories:

  • core/calc/structural/evaluate_chi

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
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
#include <iostream>
#include <iomanip>
#include <core/data/io/Pdb.hh>

#include <core/calc/structural/protein_angles.hh>
#include <core/chemical/ChiAnglesDefinition.hh>
#include <utils/exit.hh>
#include <core/calc/structural/angles.hh>


std::string program_info = R"(

Calculates Chi dihedral angles of amino acid side chains measured in two different protein structures.

The program reads two homologous protein structures. The proteins are assumed to be already aligned and  Chi angles
are calculated for every pair of identical positions.

USAGE:
    ex_chi_correlation ./test_inputs/2fd2.pdb ./test_inputs/5fd1.pdb

)";

/** @brief Calculates correlation between Chi dihedral angles measured in two different protein structures.
 *
 * CATEGORIES: core/calc/structural/evaluate_chi;
 * KEYWORDS:   PDB input; chi angles; rotamers
 * IMG: ex_chi_correlation_plot.png
 * IMG_ALT: Example correlation of chi angles between two homologus structures
 */
int main(const int argc, const char* argv[]) {

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

  using namespace core::data::structural;
  core::data::io::Pdb readerA(argv[1]); // file name (PDB format, may be gzip-ped)
  Structure_SP proteinA = readerA.create_structure(0);
  core::data::io::Pdb readerB(argv[2]);
  Structure_SP proteinB = readerB.create_structure(0);

  std::vector<std::vector<core::real>> chiA, chiB;
  std::vector<std::string> labels; // for nice output on a screen
  std::vector<std::string> mpt_annotationsA; // MPT string - one per residue only (for other lines of a given residue empty strings are inserted)
  std::vector<std::string> mpt_annotationsB;

  if (proteinA->count_residues() != proteinB->count_residues()) {
    std::cerr << "The two input PDB files should contain only the aligned parts of input proteins and be equal in length!\n";
    return 0;
  }
  auto a_res_it = proteinA->first_const_residue();
  auto b_res_it = proteinB->first_const_residue();
  while(a_res_it!=proteinA->last_const_residue()) {

    if ((*a_res_it)->residue_type() == (*b_res_it)->residue_type()) { // check chi angles

      std::vector<core::real> ca, cb;
      for (core::index2 k = 1; k <= core::chemical::ChiAnglesDefinition::count_chi_angles((*a_res_it)->residue_type()); ++k) {
        try {
          double aA = core::calc::structural::evaluate_chi((**a_res_it), k);
          double aB = core::calc::structural::evaluate_chi((**b_res_it), k);
          ca.push_back(aA);
          cb.push_back(aB);
          labels.push_back(utils::string_format("%4d%c %4d%c %c %1d", (*a_res_it)->id(), (*a_res_it)->icode(),
            (*b_res_it)->id(), (*b_res_it)->icode(), (*a_res_it)->residue_type().code1, k));
        } catch (utils::exceptions::AtomNotFound e) {
          std::cerr << e.what() << "\n";
          std::cerr << "Can't define chi angle for residue " << (**a_res_it) << "\n";
        }
      }
      if (ca.size() == cb.size() && ca.size() > 0) {
        chiA.push_back(ca);
        chiB.push_back(cb);
        mpt_annotationsA.push_back(core::calc::structural::define_rotamer(ca));
        mpt_annotationsB.push_back(core::calc::structural::define_rotamer(cb));
      }
    }
    ++a_res_it;
    ++b_res_it;
  } // end of while loop over residues

  std::cout << "#ires jres aa k  ichi_k   jchi_k delta(chi) irot jrot\n";
  double err = 0.0;
  core::index2 n_reoriented = 0;

  size_t ilabel=0;
  for(size_t ires=0;ires<chiA.size();++ires) {
    for (size_t i = 0; i < chiA[ires].size(); ++i) {
      double e = fabs(chiA[ires][i] - chiB[ires][i]);
      if (e > M_PI) e = 2.0 * M_PI - e;
      if (e > 0.523) ++n_reoriented; // --- i.e. 30 degrees of error
      std::cout << labels[ilabel] << utils::string_format("%8.2f %8.2f %8.2f",
        core::calc::structural::to_degrees(chiA[ires][i]), core::calc::structural::to_degrees(chiB[ires][i]),
        core::calc::structural::to_degrees(e));
      if (i == 0)
        std::cout << std::setw(5) << mpt_annotationsA[ires] << std::setw(5) << mpt_annotationsB[ires] << "\n";
      else std::cout << "\n";
      err += e;
      ++ilabel;
    }
  }
  std::cout <<"# avg_err, n_diff, n_res: " << err / double(chiA.size()) << " "<<n_reoriented<<" "<<chiA.size()<<"\n";
}
../_images/ex_chi_correlation_plot.png