Protein_representation

Authors: F. Cazals and T. Dreyfus and R. Tetley

Illustration of Proteus, a god of changing nature, by Andrea Alciato. Similarly to Proteus, proteins keep changing shape, which is of paramount importance for the regulation of their functions.

Introduction

Proteins: a mix of geometry, topology, biophysics, and biology. Biomolecules in general and polypeptide chains (PC in this package) in particular are complex objects. Their description indeed involves

topological information i.e. the covalent bonds – see also Molecular_covalent_structure,

geometric information i.e. the coordinates which may be Cartesian or internal – see also Molecular_coordinates ,

biophysical annotations inherently associated to the PCs: the hierarchical organization of PCs (atoms, amino-acids, whole chain), and also selected annotations (e.g. secondary structures). (See also ESBTL for data structures giving access to such pieces of information.)

It should be stressed that these three categories of information exist independently and may be used independently. For example:

The study of interfaces within protein complexes using space filling models (solvent accessible models) only require lists of atoms and their radii. See e.g. the applications Space_filling_model_interface_finder and Space_filling_model_interface.

The study of conformations, say in the context of the exploration of potential (and free) energy landscape requires geometric and topological information to compute energies See e.g. the packages Molecular_potential_energy and Landscape_explorer .

Functionalities offered. Because of these varying needs, this package provides functionalities to access all the information available at once. The main classes provided are:

Class SBL::CSB::T_Polypeptide_chain_representation : to access the aforementioned geometric, topological and biochemical pieces of information.

Class SBL::CSB::T_Protein_representation: to access the chains of a protein quaternary structure.

Class SBL::IO::T_Protein_representation_loader : to create instances of the previous two classes. Recall that a PDB file may define several geometric models, typically if the PDB file comes from NMR. If so, the loader forces the choice of one model. For this model, a SBL::CSB::T_Protein_representation containing instances of SBl::CSB::T_Polypeptide_chain_representation is returned.

Using Polypeptide chain and protein representations

The following piece of code illustrates how to access protein representations and their chains. For the latter, one gets access to

topological information encoded in the covalent structure,
geometric information encoded in the conformation,
biological/chemical information.

Loading proteins and polypeptide chains

Consider a PDB file containing one or several chains. The loader creates one SBL::CSB::T_Polypeptide_chain_representation for each chain, and these are stored within a SBL::CSB::T_Protein_representation .

Note that one can specify the chains to be loaded. In that case, SBL::CSB::T_Protein_representation contains a map mapping a chain id to the corresponding SBL::CSB::T_Polypeptide_chain_representation .

Note also that if there are missing residues in the input PDB file, the chain will be constructed with the missing residues. As a result, the corresponding SBL::CSB::T_Polypeptide_chain_representation will have several connected components. See the package Molecular_covalent_structure for further details.

In the sequel, we show how to access the various pieces of information. The following point should be stressed:

SBL::IO::T_Protein_representation_loader::get_protein (0) : returns the SBL::CSB::T_Protein_representation associated with the first PDB file loaded (if the loader was fed with several PDB files).

  //Loads the covalent structure and the conformation from a PDB file.
  Protein_representation_loader loader;
  loader.add_input_pdb_file(argv[1]);
  loader.set_allow_incomplete_chains(true);
  loader.set_occupancy_factor(3);
  //loader.set_loaded_chains("C"); 
  bool b = loader.load(true, std::cout);
  if(!b)
    return -1;
  Protein_representation& prot = loader.get_protein(0);
 
  for(Protein_representation::Polypeptide_chains_iterator it = prot.chains_begin(); it != prot.chains_end(); it++)
    {
      Polypeptide_chain_representation& P = it->second;

Covalent Structure File Loader statistics:
Number of loaded covalent structures: 1
Details for each covalent structure :
-- structure 1:
-- -- Number of loaded atoms: 3408
-- -- Number of particles: 5406
-- -- Number of modeled particles: 3408
-- -- Number of loaded conformations: 1
-- -- Number of bonds: 5467
-- -- Number of modeled bonds: 3469
-- -- Number of built disulfide bonds: 4 / 4

A common difficulty with PDB files is the presence of non-standard molecules / residues. Since this package is concerned with protein chains only, upon facing a residue which is non standard for a protein, the loader issues a message and halts the construction of the molecule. In other words, for a file containing a mixture of molecules (proteins, nucleic acids, drugs, etc), the loader only returns the protein chains found.

Enumerating atoms, residues, and the associated information

In the sequel, we focus on SBL::CSB::T_Polypeptide_chain_representation , and show how to access information associated with atoms. As an illustration, the first part of the snippet shows how to count elements of each type via a map; the second one collects the temperature factors of all atoms.

In terms of data structures, dereferencing the iterator on atoms via (*it) gives access to the ESBTL::Molecular_atom, data structure.

      //Count the number of elements of each type :
      //Compute min / max / average of temperature factors :
      double min_B = 0, max_B = 0, ave_B = 0;
      std::cout << std::endl << "Number of elements of each type : " << std::endl;
      std::map<std::string, std::size_t> elements;
      for(Polypeptide_chain_representation::Atoms_iterator it = P.atoms_begin(); it != P.atoms_end(); it++)
        {
          //counting elements
          std::string element = (*it).element();
          std::map<std::string, std::size_t>::iterator find = elements.find(element);
          if(find == elements.end())
            find = elements.insert(std::make_pair(element, 0)).first;
          find->second++;       
 
          //temp factor
          double B = (*it).temperature_factor();
          if(it == P.atoms_begin())
            {min_B = B;max_B = B;}
          else
            {
              if(max_B < B)max_B = B;
              if(B < min_B)min_B = B;
            }
          ave_B += B;
 
        }
      ave_B /= P.get_number_of_embedded_particles();
      for(std::map<std::string, std::size_t>::const_iterator it = elements.begin(); it != elements.end(); it++)
        std::cout << it->first << " : " << it->second << std::endl;
      std::cout << std::endl << "min / max / average of temperature factors : " << min_B << " " << max_B << " " << ave_B << std::endl;

Number of elements of each type :
C : 613
H : 264
N : 193
O : 185
S : 10

min / max / average of temperature factors : 0 84.09 28.3898

In a number of protein chains, residues are accompanied by insertion codes. Resid and insertion codes are accessed directly from the residue as follows :

residue.residue_sequence_number(); residue.insertion_code();

If one or more residues are missing in a chain, that chain contains several connected components. Several methods apply both to the whole chain, and to connected components. These methods are parameterized by an integer: by default, this number is negative and the method applied to the whole chain irrespective of connected components. If there are say n connected components, the index in 0..n-1 indicates which c.c. should be processed.

Iterating on the backbone

In the following, an iterator on the backbone is used to store all backbone atoms into three containers, respectively for Calpha, C, N. As previously, backbone atoms are returned as ESBTL::Molecular_atom .

      //Collect the CA, C and N in the backbone
      std::vector<Polypeptide_chain_representation::Atom> cas, cs, ns;
      for(Polypeptide_chain_representation::Backbone_iterator it = P.backbone_begin(); it != P.backbone_end(); it++)
        {
          std::string atom_name = (*it).atom_name();
          if(atom_name.compare("CA") == 0)    cas.push_back(*it);
          else if(atom_name.compare("C") == 0) cs.push_back(*it); 
          else if(atom_name.compare("N") == 0) ns.push_back(*it);
        }
      std::cout << std::endl  << "Found in the backbone : " << cas.size() << " CA, " << cs.size() << " C, " << ns.size() << " N." << std::endl;

Found in the backbone : 129 CA, 129 C, 129 N.

We only iterate on those atoms which are provided in the PDB. That is, for a structure presenting gaps, backbone atoms in the gaps are skipped. Since the primary sequence information is easily accessed, users can identify gaps based on resids. It is also possible to iterate over the backbone for each connected component of the chain : the begin and end iterators are parameterized by the number of the connected component (if none is specified, the whole chain is considered). See the code snippet below for an illustration.

      //Collect the CA, C and N in the backbone
      for(std::size_t i = 0; i < P.get_number_of_components(); i++)
        {
          cas.clear();cs.clear();ns.clear();
          for(Polypeptide_chain_representation::Backbone_iterator it = P.backbone_begin(i); it != P.backbone_end(i); it++)
            {
              std::string atom_name = (*it).atom_name();
              if(atom_name.compare("CA") == 0)    cas.push_back(*it);
              else if(atom_name.compare("C") == 0) cs.push_back(*it); 
              else if(atom_name.compare("N") == 0) ns.push_back(*it);
            }
          std::cout << std::endl  << "Found in the backbone between residues " << P.get_first_residue(i).residue_sequence_number() << P.get_first_residue(i).insertion_code() << " and " <<  P.get_last_residue(i).residue_sequence_number() << P.get_last_residue(i).insertion_code()<<  ": " << cas.size() << " CA, " << cs.size() << " C, " << ns.size() << " N." << std::endl;
        }

Counting residues by type

Molecular residues are accessed via the class ESBTL::Molecular_residue from ESBTL. Collecting residues of a given type is straightforward:

      //Count the residues of each type
      std::cout << std::endl << "Number of residues of each type : " << std::endl;
      std::map<std::string, std::size_t> res_names;
 
      for(Polypeptide_chain_representation::Residues_iterator it = P.residues_begin(); it != P.residues_end(); it++)
        {
          std::string res_name = it->residue_name();
          std::map<std::string, std::size_t>::iterator find = res_names.find(res_name);
          if(find == res_names.end())
            find = res_names.insert(std::make_pair(res_name, 0)).first;
          find->second++;       
        }
 
      for(std::map<std::string, std::size_t>::const_iterator it = res_names.begin(); it != res_names.end(); it++)
        std::cout << it->first << " : " << it->second << std::endl;

Number of residues of each type :
ALA : 12
ARG : 11
ASN : 14
ASP : 7
CYS : 8
GLN : 3
GLU : 2
GLY : 12
HIS : 1
ILE : 6
LEU : 8
LYS : 6
MET : 2
PHE : 3
PRO : 2
SER : 10
THR : 7
TRP : 6
TYR : 3
VAL : 6

Accessing the Cartesian atomic coordinates

In the following, we show how to compute the center of mass of Calpha carbons, in a pedestrian way.

The example also calls the function SBL::CSB::T_Polypeptide_chain_representation::compute_heavy_atoms_center_of_mass() , whose name is self-explanatory.

      //Compute the center of mass of the CA
      double x = 0, y = 0, z = 0, num_calphas = 0;
      for(Polypeptide_chain_representation::Calpha_iterator it = P.calphas_begin(); it != P.calphas_end(); it++)
        {
          x += P.get_x(*it); y += P.get_y(*it); z += P.get_z(*it);
          num_calphas++;
        }
      x /= num_calphas; y /= num_calphas; z /= num_calphas;
      Polypeptide_chain_representation::Point_3 com = P.compute_heavy_atoms_center_of_mass();
      std::cout  << std::endl << "Center of mass of CA : (" << x << ", " << y << ", " << z << ")" << std::endl;
      std::cout << std::endl << "Center of mass of Heavy atoms: (" << com.x() << ", " << com.y() << ", " << com.z() << ")" << std::endl;

Center of mass of CA : (53.2452, -17.1442, 7.57498)

Center of mass of Heavy atoms: (53.1526, -17.1834, 7.31916)

(Advanced) Note that the functions get_x() get_y() get_z() retrieve cartesian coordinates from the conformation. As mentioned in Introduction, the Cartesian coordinates are available via two channels: the data structures for the information contained in the PDB file (the ESBTL data structures) and the data structure for the conformation. Both are coherent if the coordinates are not changed; nevertheless, functions get_x() get_y() get_z() use the conformation.

Changing the Cartesian atomic coordinates

The following example shows how to change Cartesian coordinates:

      //Translate the protein two angstroms in the x direction
      for(Polypeptide_chain_representation::Atoms_iterator it = P.atoms_begin(); it != P.atoms_end(); it++)
        P.get_x(*it) += 2;
 
      Polypeptide_chain_representation::Point_3 new_com = P.compute_heavy_atoms_center_of_mass();
      double distance = CGAL::sqrt(CGAL::squared_distance(com, new_com));
      std::cout << "Distance between new and old center of mass: " << distance << std::endl;

Distance between new and old center of mass: 2

(Advanced) Following the remark above, upon changing the coordinates, Cartesian coordinates must be accessed via the conformation and not the ESBTL data structures. This is naturally taken care of by the get_x() get_y() get_z() functions. (NB: to put it sharply, the user should not mine the data structures to access the hidden ESBTL data structures!)

Accessing internal coordinates

In the following, we show how to access internal coordinates, namely bond lengths, valence angles, and dihedral angles. Note in particular that the latter can be used to produce the so-called ramachandran plot. See also Fig. dihedral-angles-backbone.

This first snippet illustrates iterators returning all internal coordinates, which are accessed via dedicated iterators:

      //Compute internal coordinates :
      double mean_length = 0, num_bonds = 0;
      for(Polypeptide_chain_representation::Bond_iterator it = P.bonds_begin(); it != P.bonds_end(); it++)
        {
          mean_length += (*it).get_bond_length();
          num_bonds++;
        }
      std::cout << "Mean bond length: " << (mean_length/num_bonds) << std::endl;
 
      double mean_valence_angle = 0, num_angles = 0;
      for(Polypeptide_chain_representation::Bond_angle_iterator it = P.bond_angles_begin(); it != P.bond_angles_end(); it++)
        {
          mean_valence_angle +=  (*it).get_valence_angle();
          num_angles++;
        }
      std::cout << "Mean valence angle: " << (mean_valence_angle/num_angles) << std::endl;
 
      double mean_dihedral_angle = 0;
      num_angles = 0;
      for(Polypeptide_chain_representation::Dihedral_angle_iterator it = P.dihedral_angles_begin(); it != P.dihedral_angles_end(); it++)
        {
          mean_dihedral_angle +=  (*it).get_dihedral_angle();
          num_angles++;
        }
      std::cout << "Mean dihedral angle: " << (mean_dihedral_angle/num_angles) << std::endl;
  
      double mean_phi_angle = 0;
      num_angles = 0;
      for(Polypeptide_chain_representation::Phi_iterator it = P.phi_angles_begin(); it != P.phi_angles_end(); it++)
        {
          mean_phi_angle +=  (*it).get_dihedral_angle();
          num_angles++;
        }
      std::cout << "Mean phi angle: " << (mean_phi_angle/num_angles) << std::endl;
 
      double mean_psi_angle = 0;
      num_angles = 0;
      for(Polypeptide_chain_representation::Psi_iterator it = P.psi_angles_begin(); it != P.psi_angles_end(); it++)
        {
          mean_psi_angle +=  (*it).get_dihedral_angle();
          num_angles++;
        }
      std::cout << "Mean psi angle: " << (mean_psi_angle/num_angles) << std::endl;
 
      double mean_omega_angle = 0;
      num_angles = 0;
      for(Polypeptide_chain_representation::Omega_iterator it = P.omega_angles_begin(); it != P.omega_angles_end(); it++)
        {
          mean_omega_angle +=  (*it).get_dihedral_angle();
          num_angles++;
        }
 
      std::cout << "Mean omega angle: " << (mean_omega_angle/num_angles) << std::endl;

Mean bond length: 2.0449
Mean valence angle: 1.63797
Mean dihedral angle: 0.222753
Mean phi angle: -1.05698
Mean psi angle: 0.52554
Mean omega angle: 0.425752

This second snippet focuses on the dihedral angles associated with a given residue:

      std::cout << "Residue 33" << " has Phi angle: " << P.get_phi_angle(33)
                << ", Psi angle: " << P.get_psi_angle(33)
                << " and Omega angle: " << P.get_omega_angle(33)
                << std::endl;    

Residue 2 has Phi angle: -1.87983, Psi angle: 2.15193 and Omega angle: 3.11996
Residue 33 has Phi angle: -1.04222, Psi angle: -0.895227 and Omega angle: -3.08817

Dihedral angles along the backbone of a polypeptide chain. By convention, the three angles display, namely $\Phi, \Psi, \omega$ , are associated with the i-th amino-acid.

Implementation and functionalities

Data structures and classes

As mentioned in Introduction, PC come with topological, geometric, and biophysical information. This package provides three main classes:

SBL::IO::T_Protein_representation_loader : the loader , which processes PDB files. It follows the loader concept as defined in Module_base .

SBL::CSB::T_Polypeptide_chain_representation: data structure to access the topological, geometric and biophysical information associated with a PC.

SBL::CSB::T_Protein_representation: data structure mapping chain ids to their SBL::CSB::T_Polypeptide_chain_representation .

Internally, the topological, geometric, and biophysical pieces of information are stored into the following DS:

topological information: SBL::CSB::T_Molecular_covalent_structure, from package Molecular_covalent_structure
geometric information: ConformationTraits, namely any data structure accommodating a conformation i.e. a d-dimensional point. The simplest such data structure is a vector of doubles, but more involved options are made available in package ConformationTraits .
biophysical annotations: in ESBTL, recall that a molecular system refers to a collection of molecular models – as found in a typical PDB file. Biophysical annotations, including the hierarchy associated with a PC (chain, a.a., atoms), are stored within

Options offered by these classes

The loader SBL::IO::T_Protein_representation_loader has a number of options listed below, that can be set either from the command line using the Module_base framework, or directly using appropriate methods :

–pdb-file string: Input PDB file
(SBL::IO::T_Protein_representation_loader::add_input_pdb_file)
–water: loads water molecules
(SBL::IO::T_Protein_representation_loader::set_loaded_water)
–hetatoms: loads hetero-atoms
(SBL::IO::T_Protein_representation_loader::set_loaded_hetatoms)
–coarse-level int: coarse level for the covalent structure (0 : none, 1 : heavy atoms, 2 : residues)
(SBL::IO::T_Protein_representation_loader::set_coarse_level)
–no-ss-bonds: does not perform disulfide bonds calculations
(SBL::IO::T_Protein_representation_loader::set_without_disulfide_bonds)
–occupancy-policy: Selection policy for atoms with no alternate location and occupancy not equal to 1. Allowed values are: 1 (all), 2 (forbbiden, default), 3 (none), 4 (max), 5 (min)
(SBL::IO::T_Protein_representation_loader::set_occupancy_factor)
–alternate string: Alternate coordinates: alternative to be used (char)
(SBL::IO::T_Protein_representation_loader::set_alternate_coordinate)
–model-number int: ID of the model to be loaded frome the PDB file (default is 1)
(SBL::IO::T_Protein_representation_loader::set_PDB_model_number)
–chains string: Subset of chains to be loaded.
(SBL::IO::T_Protein_representation_loader::set_loaded_chains)
–allow-incomplete-chains: Allows to load residues even if there are missing ones.
(SBL::IO::T_Protein_representation_loader::set_allow_incomplete_chains)

Example

The following example is the tutorial example presented in section Using Polypeptide chain and protein representations snippet by snippet :

#include <iostream>
#include <sstream>
#include <iterator>
#include <algorithm>
 
#include <SBL/Models/Atom_with_hierarchical_info_traits.hpp>
#include <SBL/CSB/Particle_info_for_proteins.hpp>
#include <SBL/CSB/Molecular_covalent_structure.hpp>
#include <SBL/CSB/Molecular_covalent_structure_builder_for_proteins.hpp>
#include <SBL/IO/Protein_representation_loader.hpp>
#include <SBL/CSB/Molecular_primitive_internal_coordinates.hpp>
#include <SBL/CSB/Polypeptide_chain_representation.hpp>
#include <SBL/CSB/Protein_representation.hpp>
 
typedef SBL::Models::T_Atom_with_hierarchical_info_traits<>                               Particle_traits;
typedef typename Particle_traits::Molecular_system                                        Molecular_system;
typedef SBL::CSB::T_Particle_info_for_proteins<Particle_traits>                           Particle_info;
typedef SBL::CSB::T_Particle_info_traits<Particle_info>                                   Particle_info_traits;
 
typedef SBL::CSB::T_Molecular_covalent_structure<Particle_info>                           Covalent_structure;
typedef SBL::CSB::T_Molecular_covalent_structure_builder_for_proteins<Covalent_structure> Covalent_structure_builder;
 
typedef std::vector<double>                                                               Conformation;
 
typedef SBL::CSB::T_Molecular_primitive_internal_coordinates
<Conformation, Covalent_structure>                                                        Internal_coordinates;
 
typedef SBL::CSB::T_Polypeptide_chain_representation <Particle_traits, Covalent_structure, Conformation>   Polypeptide_chain_representation;
 
typedef SBL::CSB::T_Protein_representation<Polypeptide_chain_representation> Protein_representation;
 
typedef SBL::IO::T_Protein_representation_loader
<Protein_representation, Covalent_structure_builder>           Protein_representation_loader;
 
int main(int argc, char *argv[])
{
  if(argc < 2)
    return -1;
  
  //Loads the covalent structure and the conformation from a PDB file.
  Protein_representation_loader loader;
  loader.add_input_pdb_file(argv[1]);
  loader.set_allow_incomplete_chains(true);
  loader.set_occupancy_factor(3);
  //loader.set_loaded_chains("C"); 
  bool b = loader.load(true, std::cout);
  if(!b)
    return -1;
  Protein_representation& prot = loader.get_protein(0);
 
  for(Protein_representation::Polypeptide_chains_iterator it = prot.chains_begin(); it != prot.chains_end(); it++)
    {
      Polypeptide_chain_representation& P = it->second;
 
      //Count the number of elements of each type :
      //Compute min / max / average of temperature factors :
      double min_B = 0, max_B = 0, ave_B = 0;
      std::cout << std::endl << "Number of elements of each type : " << std::endl;
      std::map<std::string, std::size_t> elements;
      for(Polypeptide_chain_representation::Atoms_iterator it = P.atoms_begin(); it != P.atoms_end(); it++)
        {
          //counting elements
          std::string element = (*it).element();
          std::map<std::string, std::size_t>::iterator find = elements.find(element);
          if(find == elements.end())
            find = elements.insert(std::make_pair(element, 0)).first;
          find->second++;       
 
          //temp factor
          double B = (*it).temperature_factor();
          if(it == P.atoms_begin())
            {min_B = B;max_B = B;}
          else
            {
              if(max_B < B)max_B = B;
              if(B < min_B)min_B = B;
            }
          ave_B += B;
 
        }
      ave_B /= P.get_number_of_embedded_particles();
      for(std::map<std::string, std::size_t>::const_iterator it = elements.begin(); it != elements.end(); it++)
        std::cout << it->first << " : " << it->second << std::endl;
      std::cout << std::endl << "min / max / average of temperature factors : " << min_B << " " << max_B << " " << ave_B << std::endl;
 
      //Collect the CA, C and N in the backbone
      std::vector<Polypeptide_chain_representation::Atom> cas, cs, ns;
      for(Polypeptide_chain_representation::Backbone_iterator it = P.backbone_begin(); it != P.backbone_end(); it++)
        {
          std::string atom_name = (*it).atom_name();
          if(atom_name.compare("CA") == 0)    cas.push_back(*it);
          else if(atom_name.compare("C") == 0) cs.push_back(*it); 
          else if(atom_name.compare("N") == 0) ns.push_back(*it);
        }
      std::cout << std::endl  << "Found in the backbone : " << cas.size() << " CA, " << cs.size() << " C, " << ns.size() << " N." << std::endl;
 
      //Collect the CA, C and N in the backbone
      for(std::size_t i = 0; i < P.get_number_of_components(); i++)
        {
          cas.clear();cs.clear();ns.clear();
          for(Polypeptide_chain_representation::Backbone_iterator it = P.backbone_begin(i); it != P.backbone_end(i); it++)
            {
              std::string atom_name = (*it).atom_name();
              if(atom_name.compare("CA") == 0)    cas.push_back(*it);
              else if(atom_name.compare("C") == 0) cs.push_back(*it); 
              else if(atom_name.compare("N") == 0) ns.push_back(*it);
            }
          std::cout << std::endl  << "Found in the backbone between residues " << P.get_first_residue(i).residue_sequence_number() << P.get_first_residue(i).insertion_code() << " and " <<  P.get_last_residue(i).residue_sequence_number() << P.get_last_residue(i).insertion_code()<<  ": " << cas.size() << " CA, " << cs.size() << " C, " << ns.size() << " N." << std::endl;
        }
 
      //Count the residues of each type
      std::cout << std::endl << "Number of residues of each type : " << std::endl;
      std::map<std::string, std::size_t> res_names;
 
      for(Polypeptide_chain_representation::Residues_iterator it = P.residues_begin(); it != P.residues_end(); it++)
        {
          std::string res_name = it->residue_name();
          std::map<std::string, std::size_t>::iterator find = res_names.find(res_name);
          if(find == res_names.end())
            find = res_names.insert(std::make_pair(res_name, 0)).first;
          find->second++;       
        }
 
      for(std::map<std::string, std::size_t>::const_iterator it = res_names.begin(); it != res_names.end(); it++)
        std::cout << it->first << " : " << it->second << std::endl;
 
      //Compute the center of mass of the CA
      double x = 0, y = 0, z = 0, num_calphas = 0;
      for(Polypeptide_chain_representation::Calpha_iterator it = P.calphas_begin(); it != P.calphas_end(); it++)
        {
          x += P.get_x(*it); y += P.get_y(*it); z += P.get_z(*it);
          num_calphas++;
        }
      x /= num_calphas; y /= num_calphas; z /= num_calphas;
      Polypeptide_chain_representation::Point_3 com = P.compute_heavy_atoms_center_of_mass();
      std::cout  << std::endl << "Center of mass of CA : (" << x << ", " << y << ", " << z << ")" << std::endl;
      std::cout << std::endl << "Center of mass of Heavy atoms: (" << com.x() << ", " << com.y() << ", " << com.z() << ")" << std::endl;
 
      //Translate the protein two angstroms in the x direction
      for(Polypeptide_chain_representation::Atoms_iterator it = P.atoms_begin(); it != P.atoms_end(); it++)
        P.get_x(*it) += 2;
 
      Polypeptide_chain_representation::Point_3 new_com = P.compute_heavy_atoms_center_of_mass();
      double distance = CGAL::sqrt(CGAL::squared_distance(com, new_com));
      std::cout << "Distance between new and old center of mass: " << distance << std::endl;
 
      //Compute internal coordinates :
      double mean_length = 0, num_bonds = 0;
      for(Polypeptide_chain_representation::Bond_iterator it = P.bonds_begin(); it != P.bonds_end(); it++)
        {
          mean_length += (*it).get_bond_length();
          num_bonds++;
        }
      std::cout << "Mean bond length: " << (mean_length/num_bonds) << std::endl;
 
      double mean_valence_angle = 0, num_angles = 0;
      for(Polypeptide_chain_representation::Bond_angle_iterator it = P.bond_angles_begin(); it != P.bond_angles_end(); it++)
        {
          mean_valence_angle +=  (*it).get_valence_angle();
          num_angles++;
        }
      std::cout << "Mean valence angle: " << (mean_valence_angle/num_angles) << std::endl;
 
      double mean_dihedral_angle = 0;
      num_angles = 0;
      for(Polypeptide_chain_representation::Dihedral_angle_iterator it = P.dihedral_angles_begin(); it != P.dihedral_angles_end(); it++)
        {
          mean_dihedral_angle +=  (*it).get_dihedral_angle();
          num_angles++;
        }
      std::cout << "Mean dihedral angle: " << (mean_dihedral_angle/num_angles) << std::endl;
  
      double mean_phi_angle = 0;
      num_angles = 0;
      for(Polypeptide_chain_representation::Phi_iterator it = P.phi_angles_begin(); it != P.phi_angles_end(); it++)
        {
          mean_phi_angle +=  (*it).get_dihedral_angle();
          num_angles++;
        }
      std::cout << "Mean phi angle: " << (mean_phi_angle/num_angles) << std::endl;
 
      double mean_psi_angle = 0;
      num_angles = 0;
      for(Polypeptide_chain_representation::Psi_iterator it = P.psi_angles_begin(); it != P.psi_angles_end(); it++)
        {
          mean_psi_angle +=  (*it).get_dihedral_angle();
          num_angles++;
        }
      std::cout << "Mean psi angle: " << (mean_psi_angle/num_angles) << std::endl;
 
      double mean_omega_angle = 0;
      num_angles = 0;
      for(Polypeptide_chain_representation::Omega_iterator it = P.omega_angles_begin(); it != P.omega_angles_end(); it++)
        {
          mean_omega_angle +=  (*it).get_dihedral_angle();
          num_angles++;
        }
 
      std::cout << "Mean omega angle: " << (mean_omega_angle/num_angles) << std::endl;
 
      std::cout << "Residue 33" << " has Phi angle: " << P.get_phi_angle(33)
                << ", Psi angle: " << P.get_psi_angle(33)
                << " and Omega angle: " << P.get_omega_angle(33)
                << std::endl;    
    }
  return 0;
 
}

Applications

This package also offers several useful programs to inspect properties of proteins / their conformations:

sbl-protein-info.exe: loads protein chain(s) form a PDB file, and parses it to deliver a variety of statistics (number of a.a. of each type, center of mass of the structure, calculation of internal coordinates and average statistics, etc. Note that the corresponding code of interest for developers that wish to use this package to develop novel applications.

sbl-protein-ramachandran.exe: computes the Ramachandran plot of a protein; alternatively, for conformations of the same protein, computes the Ramachandran plot of specified amino-acids.

sbl-protein-pdb-cleaner.py : a python script cleaning PDB files to make sure the contain the required information for ESBTL loaders (in particular chemical element names), and builders of the covalent structure (correct atom naming required to build the molecular topology).

Table of Contents