Skip to content

Home

Jump to bottom Edit New page
PyP0 edited this page Aug 21, 2015 · 1 revision

Structural superimposition project

Specifications:

Compute the optimal rigid superimposition between two sequence equivalent proteic structures.

Input:

  • Two proteic structures encoded as a PDB file.

Output:

  • Optimal superimposition of the two input structures (only C-alpha atoms), encoded in a PDB file;
  • Score metrics to evaluate the superimposition quality:
    • RMSD score;
    • TM-Score;
    • GDT-TS score;
    • MaxSub score.

A (brief) description of the work

Kabsch algorithm
The core of the project is the function which compute the optimal roto-translation that minimizes the RMSD (root mean squared deviation) between two paired sets of points. These two set of points are the 3D coordinates of the two input proteins’ C-alpha atoms. This function is implemented with the Kabsch algorithm, a popular way to resolve these kind of problems.
The algorithm’s input requires two sets of paired points (A, and B) expressed as Nx3 matrices, where each row is a 3D point and the columns represent the x,y,z coordinates.

The algorithm can be summarized in these steps:

  1. Find the centroids of both datasets A and B;

  2. Subtract the centroids values from each point of the two matrices, so that both centroids are at the origin;

  3. Compute the covariance matrix H of the two modified matrices;

  4. The optimal rotation matrix R can be easily calculated if Singular Value Decomposition (SVD) routines are available : [U,S,V]= svd(H); where U, S, V are matrices needed to compute R;

  5. Compute R as: R = V * transposed_U;

  6. Find the determinant’s sign of R to detect a possible “reflection” case, caused by the SVD method, and if it’s < 0 multiply the 3rd column of R by -1. This is the final optimal rotation matrix;

  7. Find the optimal translation vector t: t = - R * centroid_A + centroid_B.

A different and slightly more efficient way to implement the algorithm is by using quaternions.

Project classes
The project is composed by three classes:

  • RotoTranslation;
  • KabschAlgorithm;
  • StructuralSuperimposition.

The rotoTranslation class is an abstract class containing only a virtual method calculateRotoTranslation(), intended to be a base for multiple and different implementations of a roto-translation algorithm, like Kabsch, quaternion, etc.

The kabschAlgorithm, as the name suggests, is a concrete class which inherits from the abstract class, and implements the virtual method calculateRotoTranslation() with the Kabsch algorithm. The implemented version of Kabsch algorithm requires two Eigen nx3 matrices and, after computing the optimal rotation matrix and translation vector, stores them in a special Eigen object called Affine3d, which is a member of the KabschAlgorithm class.

The structuralSuperimposition class contains the iterative search algorithm used for each score, the functions to compute every score, the “interface methods” between Victor and Eigen (how to load the Spacer’s atoms coordinates into an Eigen matrix and vice-versa), and another Affine3d member object, which contains the optimal roto-translation matrix computed at the end of the search algorithm.

The iterative search algorithm
Every score is calculated using the MaxSub heuristic algorithm for finding the largest subset M = { (ai,bi), (ai+1, bi+1), …} such that for every pair (ai,bi) the Euclidean distance between ai and T(bi) is below the distance threshold d0, where T is the roto-translation that best superimposes the points of B over A (i.e. the roto-translation matrix obtained by executing Kabsch algorithm). The algorithm has been implemented following exactly the instructions given in [3]. Little implementation differences can be seen between different scores:

  • RMSD and MaxSub: the d0 threshold is not fixed, and can be changed by the user;
  • TM-Score: the d0 is automatically evaluated using the equation (5) in [1].
  • GDT-TS: following the considerations in paragraph 4.3 in [2], the score is implemented by calling four times the iterative algorithm with different d0 thresold (1, 2, 4, 8);

The Eigen library
From the Wikipedia page:

“Eigen is a high-level C++ library of template headers for linear algebra, matrix and vector operations, numerical solvers and related algorithms. Eigen is an open source library licensed under MPL2 starting from version 3.1.1. Earlier versions were licensed under LGPL3+. Eigen is often noted for its elegant API, versatile fixed and dynamic matrix capabilities and a range of dense and sparse solvers. To achieve high performance, Eigen utilizes explicit vectorization for the SSE 2/3/4, ARM NEON, and AltiVec instruction sets.”
Wikipedia: Eigen

Eigen is free software, doesn't have any dependencies other than the C++ standard library and is standard C++98. Eigen is being successfully used with the following compilers:

  • GCC, version 4.1 and newer. Very good performance with GCC 4.2 and newer.
  • MSVC (Visual Studio), 2008 and newer (the old 2.x versions of Eigen support MSVC 2005, but without vectorization).
  • Intel C++ compiler. Very good performance.
  • LLVM/CLang++ (2.8 and newer).
  • MinGW, recent versions. Very good performance when using GCC 4.
  • QNX's QCC compiler.

Adding the Eigen library (eigen.tuxfamily.org) to the Victor project was felt useful to take advantage of the SVD routines and the matrices manipulation methods (avoiding to re-implement them from scratch). It also implements the quaternions, so it will be useful if, sometimes in the future, the Kabsch method will be re-implemented using them.

The used modules are:

  • Eigen/Core: Matrix and Array classes, basic linear algebra, array manipulation;
  • Eigen/Geometry: Transform (Affine3d), Translation, Scaling, Rotation2D and 3D rotations (Quaternion, AngleAxis);
  • Eigen/SVD: SVD decomposition with least-squares solver (JacobiSVD).

Version
The Eigen version used in the project is the 3.2.5 (the latest stable release, at the moment), but to accommodate the Victor’s GNU GPL it was also tested on the last LGPL3+ version, the 3.1.0 (June 24, 2012).

Installation
In order to use Eigen, no installations are required. The header files in the Eigen subdirectory are the only files required to compile programs using Eigen, and they are the same for all platforms. It is not necessary to use CMake or install anything. The latest Eigen's source code can be downloaded from Eigen download page, while the old LGPL3+ version can be downloaded from here. In the #include directives will be placed the path where the Eigen directory will be saved.

Notice: Probably for a library specific reason, during the "make" of the Test folder an error about Eigen's max template depth reached in Core module can occur:

error: template instantiation depth exceeds maximum of 36 (use -ftemplate-depth= to increase the maximum)

A quick workaround can be made raising the -ftemplate-depth flag in Makefile.global from 36 to 40. This error seems related only to the file Testkabsch.h added in that folder, but following the "stack" trace seems related to the Core library module.

How to use the application
The application requires a minimum number of mandatory parameters, in this order:

  1. A PDB file path for the first protein structure (the “model”, usually referred as A);

  2. A PDB file path for the second protein structure (the “experimental structure”, usually referred as B);

  3. An output PDB file name that will contain the optimally roto-translated B protein structure;

  4. An integer L (>= 4) representing the length of the starting seeds in the Max Sub iterative search algorithm; obviously it can’t be bigger than the protein size;

  5. A letter representing which score to compute:

    • –t : TM-Score;
    • –m : MaxSub;
    • –g : GDT-TS score;
    • –r : RMSD score;
    • –a : all of the above;

  6. Optional parameter: if “-m”, “-a” or “-r” is chose, a custom d0 distance threshold (in Angstrom) can be selected. Default value is 3.5 A.

Example:
Given the two protein PDB files:

  • A.PDB;
  • B.PDB;

To view every metric score with L = 4 and default d0 distance:
./PdbSuperimposition A.PDB B.PDB output.PDB 4 –a

To compute MaxSub score with L=50 and d0 = 5.0 A:
./PdbSuperimposition A.PDB B.PDB output1.PDB 4 –m 5.0

Tests results
A first test has been made with the proteins:

  • 15C8_H_input.pdb
  • 25C8_H_input.pdb

The program scores will be tested against those considered "exacts", that is the ones taken from Zhanglab.

The “exact” scores output was:

  • 217 residues found;
  • RMSD of the common residues= 1.604
  • TM-score = 0.9350 (d0= 5.48)
  • MaxSub-score= 0.8572 (d0= 3.50)
  • GDT-TS-score= 0.8733 %(d<1)=0.6452 %(d<2)=0.8710 %(d<4)=0.9770 %(d<8)=1.0000

The scores computed by the program, with L = 4 and d0 = 3.5 are:

  • 211 residues found;
  • RMSD score: 1.58475
  • TM-score: 0.932606 (with automatically calculated d0 threshold: 5.40289)
  • MaxSub score: 0.857753
  • GDT-TS score: 0.767476; (d<1A): 0.526608;(d<2A): 0.6955;(d<4A): 0.881262;(d<8A): 0.966536

Increasing the L value we can see a really small change in the values:

L = 50:

  • RMSD score: 1.58475
  • TM-score: 0.932606
  • MaxSub score: 0.857753
  • GDT-TS score: 0.771526

L = 150:

  • RMSD score: 1.58475
  • TM-score: 0.932606
  • MaxSub score: 0.857753
  • GDT-TS score: 0.718056

jmolView

More tests have been made using CASP11 predictions.

A couple of examples took from CASP predictioncenter:

T0759TS006_1.pdb T0759TS006_2.pdb

  • “exact” score:

    • Structure1: A538104 Length= 109
    • Structure2: B538104 Length= 109 (by which all scores are normalized)
    • Number of residues in common= 109
    • RMSD of the common residues= 15.425
    • TM-score = 0.1983 (d0= 3.84)
    • MaxSub-score= 0.1226 (d0= 3.50)
    • GDT-TS-score= 0.1766 %(d<1)=0.0734 %(d<2)=0.1284 %(d<4)=0.1743 %(d<8)=0.3303

  • PdbSuperimposition score (L = 4, d0 = 3.5):

    • numero AA nello spacer: 109
    • numero AA nello spacer: 109
    • RMSD score: 34.8957
    • TM-score (d0 = 3.83808): 0.108219
    • MaxSub score: 0.10984
    • gdt-ts score (d<1A): 0.0512168
    • gdt-ts score (d<2A): 0.0742871
    • gdt-ts score (d<4A): 0.119353
    • gdt-ts score (d<8A): 0.24389
    • GDT-TS score: 0.122187

T0759TS006_3.pdb T0759TS006_4.pdb

  • “exact” score:

    • Structure1: A677926 Length= 109
    • Structure2: B677926 Length= 109 (by which all scores are normalized)
    • Number of residues in common= 109
    • RMSD of the common residues= 15.755
    • TM-score = 0.2147 (d0= 3.84)
    • MaxSub-score= 0.1232 (d0= 3.50)
    • GDT-TS-score= 0.2064 %(d<1)=0.1101 %(d<2)=0.1284 %(d<4)=0.1927 %(d<8)=0.3945

  • PdbSuperimposition score (L = 4, d0 = 3.5):

    • numero AA nello spacer: 109
    • numero AA nello spacer: 109
    • RMSD score: 28.2708
    • TM-score (d0 = 3.83808): 0.124511
    • MaxSub score: 0.122333
    • gdt-ts score (d<1A): 0.0604032
    • gdt-ts score (d<2A): 0.102734
    • gdt-ts score (d<4A): 0.125398
    • gdt-ts score (d<8A): 0.217372
    • GDT-TS score: 0.126477

Problems still open, autocriticism and possible future work:

  • More testing: from the small set of tests made it's hard to evaluate the precision of the scores. A really good precision can be seen in the MaxSub evaluation, and also a quite good one in TM-Score. The precision of GDT-TS is constrained by his specific implementation, as said in [2]. The RMSD score precision, if compared with the "exact" one, is really bad.
  • "Missing Aminoacids": looking at the "exact" score, we can see a difference between the AA found by the two programs. This can be one of the reasons why the scores are not very accurates.
  • Calibrations: a possible way to improve the scores would be a calibration of the parameters L and d0.
  • Kabsch algorithm: re-implementing the algorithm with quaternions (and maybe without Eigen) would create a more numerical-stable algorithm.

References:

[1] Scoring function for automated assessment of protein structure template quality.
Zhang Y, Skolnick J.
Proteins. 2004 Dec 1;57(4):702-10.

[2] http://predictioncenter.org/casp/casp7/public/doc/LCS_GDT.README
Manuale ClusCo (file pdf)

[3] MaxSub: an automated measure for the assessment of protein structure prediction quality.
Siew N, Elofsson A, Rychlewski L, Fischer D.
Bioinformatics. 2000 Sep;16(9):776-85.

Add a custom footer
Add a custom sidebar

Clone this wiki locally