RMSD

How can I map two structures onto each other and compute their RMSD?

Mapping two systems onto each other can be done using BALL's StructureMapper. After successful mapping, the RMSD can be easily computed calling the function calculateRMSD().

In the following we assume, that the two structures, e.g. proteins, are already "moved" onto each other, i.e. a transformation has been applied. How to do this is described here.

C++

#include <BALL/FORMAT/PDBFile.h>
#include <BALL/KERNEL/system.h>
#include <BALL/STRUCTURE/structureMapper.h>

using namespace BALL;
using namespace std;


  // read in two systems, of which the RMSD should be determined
  PDBFile pdbfile1("test_complex_1.pdb", std::ios::in);
  PDBFile pdbfile2("test_complex_2.pdb", std::ios::in);

  System S1;
  System S2;

  pdbfile1 >> S1;
  pdbfile2 >> S2;

  float rmsd = 0.0;

  // use the StructureMapper to map both systems onto each other
  StructureMapper mapper(S1, S2);
  mapper.calculateDefaultBijection();
  
  // determine the RMSD
  rmsd = mapper.calculateRMSD();
  cout << " The RMSD of the given structures is: " << rmsd << endl;
  
  pdbfile1.close();
  pdbfile2.close();

If you want to limit the number of atoms, for which you want to compute the rmsd, you can create a AtomBijection and pass it to the StructureMapper.

  AtomBijection atom_bijection;

  // create an AtomBijection for all backbone C alpha atoms
  atom_bijection.assignCAlphaAtoms(S1, S2);
  // or
  atom_bijection.assignBackboneAtoms(S1, S2);
 
  // pass the bijection to the StructureMapper
  rmsd = mapper.calculateRMSD(atom_bijection);

For cases where the atom set limitations are more elaborate, e.g. a binding pocket, the method assignAtomsByProperty(S1, S2) creates an AtomBijection for all matching atoms where at least one atoms is marked with a property "ATOMBIJECTION_RMSD_SELECTION".

An elegant way to create such a limitation is BALL's selector mechanism:

  System sys;
  PDBFile pdb("mypdb.pdb"));
  pdb.read(sys);

  // create a selection, e.g. using BALL's selector class
  Selector sel;
  sys.deselect();
  Expression e("residue(ALA)");
  sel.setExpression(e);
  sys.apply(sel);
  list<Atom*> selected_atoms = sel.getSelectedAtoms();
  list<BALL::Atom*>::iterator it = selected_atoms.begin();
  for ( ; it != selected_atoms.end(); ++it)
  {
    // set the property
    (*it)->setProperty("ATOMBIJECTION_RMSD_SELECTION", true);
  }

  // create a bijection for all matching atoms where
  // one atom has the property "ATOMBIJECTION_RMSD_SELECTION"
  atom_bijection.assignAtomsByProperty(sys, S2);
  rmsd = mapper.calculateRMSD(atom_bijection);

In the context of docking, computing each pairwise AtomBijection becomes far to laborius. Fortunately, the properties can be passed forward to our class ConformationSet, that is e.g. used by our class PoseClustering.

  // read a "reference" structure
  System sys;
  PDBFile pdb("mypdb.pdb"));
  pdb.read(sys);

  // set the properties as shown above  

  // read the conformation set and the reference structure
  ConformationSet cs2;
  cs2.setup(sys);
  cs2.readDCDFile("my_pdb_set.dcd"));
  cs2.resetScoring();

  // create a PoseClustering
  PoseClustering pc;
  pc.options.set(PoseClustering::Option::RMSD_LEVEL_OF_DETAIL, PoseClustering::PROPERTY_BASED_ATOM_BIJECTION);
  pc.setConformationSet(&cs2);

  // cluster the structures
  pc.compute();

  // get the number of clusters
  cout << pc.getNumberOfClusters() << endl;
  
  // get a representative structure for cluster 5
  boost::shared_ptr<System> sys5 = pc.getClusterRepresentative(5);
  
  // get a ConformationSet containing one structure per final cluster 
  boost::shared_ptr<ConformationSet> cs4 = pc.getReducedConformationSet();

  // get all conformations of cluster 6 
  boost::shared_ptr<ConformationSet> cs3 = pc.getClusterConformationSet(6);

Note that the PoseClustering assumes, that all necessary transformations have been performed previously!

You may want to have a look onto our tool APPLICATIONS/TOOLS/DockPoseClustering.C.

Python

 
from BALL import *

# get two PDBFiles
pdbfile1 = PDBFile("/path/to/a/StructureMapper_test.pdb")
pdbfile2 = PDBFile("/path/to/a/StructureMapper_test.pdb")

S1 = System()
S2 = System()

pdbfile1 >> S1
pdbfile2 >> S2

# change the atom positions a little bit
for at in atoms(S1):
  at.setPosition(at.getPosition() + Vector3(1, 1, 1))
rmsd = float(0.0)

# use the StructureMapper to map both systems onto each other
mapper = StructureMapper(S1, S2)
mapper.calculateDefaultBijection()
  
# determine the RMSD
rmsd = mapper.calculateRMSD()
print " The RMSD of the given structures is: ", rmsd
  
#done
pdbfile1.close()
pdbfile2.close()