Why is this source code giving incorrect RMSD value?

Question

• Aligning PDB structures with Biopython

The following source code is obtained from the above link:

    import Bio.PDB    
    # Select what residues numbers you wish to align
    # and put them in a list
    start_id = 1
    end_id   = 70
    atoms_to_be_aligned = range(start_id, end_id + 1)
    
    # Start the parser
    pdb_parser = Bio.PDB.PDBParser(QUIET = True)
    
    # Get the structures
    ref_structure = pdb_parser.get_structure("reference", "c:/home/bbq_input_pdb/pdb1a6j.pdb")
    sample_structure = pdb_parser.get_structure("sample", "c:/home/bbq_output_pdb/1a6j.pdb")
    
    # Use the first model in the pdb-files for alignment
    # Change the number 0 if you want to align to another structure
    ref_model    = ref_structure[0]
    sample_model = sample_structure[0]
    
    # Make a list of the atoms (in the structures) you wish to align.
    # In this case we use CA atoms whose index is in the specified range
    ref_atoms = []
    sample_atoms = []
    
    # Iterate of all chains in the model in order to find all residues
    for ref_chain in ref_model:
      # Iterate of all residues in each model in order to find proper atoms
      for ref_res in ref_chain:
        # Check if residue number ( .get_id() ) is in the list
        if ref_res.get_id()[1] in atoms_to_be_aligned:
          # Append CA atom to list
          ref_atoms.append(ref_res['CA'])
    
    # Do the same for the sample structure
    for sample_chain in sample_model:
      for sample_res in sample_chain:
        if sample_res.get_id()[1] in atoms_to_be_aligned:
          sample_atoms.append(sample_res['CA'])
    
    # Now we initiate the superimposer:
    super_imposer = Bio.PDB.Superimposer()
    super_imposer.set_atoms(ref_atoms, sample_atoms)
    super_imposer.apply(sample_model.get_atoms())
    
    # Print RMSD:
    print("RMSD = ", super_imposer.rms)
    
    ## Save the aligned version of 1UBQ.pdb
    #io = Bio.PDB.PDBIO()
    #io.set_structure(sample_structure)
    #io.save("1UBQ_aligned.pdb")

This source code compares two PDB protein files and computes their RMSD values.

However, I observed that:

1. This tool gives different values than that of PyMOL with the same files
2. It shows values even if we compare a file with itself - in PyMOL, it is not possible to compare a file with itself

Can anyone tell me where the bug is, as I cannot locate it?

Answer 1

not an expert here just trying to figure out from code; see https://github.com/biopython/biopython/blob/master/Bio/PDB/Superimposer.py#L43 :

import numpy as np

from Bio.SVDSuperimposer import SVDSuperimposer -------------------------> HERE !!!!!!!!
from Bio.PDB.PDBExceptions import PDBException

class Superimposer:
    """Rotate/translate one set of atoms on top of another to minimize RMSD."""

    def __init__(self):
        """Initialize the class."""
        self.rotran = None
        self.rms = None

    def set_atoms(self, fixed, moving):
        """Prepare translation/rotation to minimize RMSD between atoms.
        Put (translate/rotate) the atoms in fixed on the atoms in
        moving, in such a way that the RMSD is minimized.
        :param fixed: list of (fixed) atoms
        :param moving: list of (moving) atoms
        :type fixed,moving: [L{Atom}, L{Atom},...]
        """
        if not len(fixed) == len(moving):
            raise PDBException("Fixed and moving atom lists differ in size")
        length = len(fixed)
        fixed_coord = np.zeros((length, 3))
        moving_coord = np.zeros((length, 3))
        for i in range(0, length):
            fixed_coord[i] = fixed[i].get_coord()
            moving_coord[i] = moving[i].get_coord()
------- sup = SVDSuperimposer()   -------------------------> HERE !!!!!!!!           
------- sup.set(fixed_coord, moving_coord)    --------------> HERE !!!!!!!!
        sup.run()
        self.rms = sup.get_rms()
        self.rotran = sup.get_rotran()

and from https://github.com/biopython/biopython/blob/master/Bio/SVDSuperimposer/init.py#L123

def set(self, reference_coords, coords):
        """Set the coordinates to be superimposed.
        coords will be put on top of reference_coords.
        - reference_coords: an NxDIM array
        - coords: an NxDIM array
        DIM is the dimension of the points, N is the number
        of points to be superimposed.
        """
        # clear everything from previous runs
        self._clear()
        # store coordinates
        self.reference_coords = reference_coords
-----   self.coords = coords    ---------------------------ADD HERE
        n = reference_coords.shape
        m = coords.shape
        if n != m or not (n[1] == m[1] == 3):
            raise Exception("Coordinate number/dimension mismatch.")
        self.n = n[0]

just add below the ADD HERE :

 if array_equal(self.reference_coords, self.coords):
            
            
            raise Exception("Sorry, you are aligning same structure")

needs from numpy import array_equal

Don't know, how the SVDSuperimposer (super_imposer.rotran) returns

the rotation and the translation from the coordinates fed to the

algorithm but could be you can get them even starting from same coords for

reference and sample ??? Not sure mine is just an hypothesis.

ADDENDUM:

from Biopython Documentation:

https://biopython.org/docs/1.76/api/Bio.SVDSuperimposer.html

as an example of how SVDSuperimposer works:

from Bio.SVDSuperimposer import SVDSuperimposer

from numpy import array, dot, set_printoptions


x = array([[51.65, -1.90, 50.07],

     [50.40, -1.23, 50.65],

     [50.68, -0.04, 51.54],

     [50.22, -0.02, 52.85]], 'f')

y = array([[51.30, -2.99, 46.54],

     [51.09, -1.88, 47.58],

     [52.36, -1.20, 48.03],

     [52.71, -1.18, 49.38]], 'f')


sup = SVDSuperimposer()

sup.set(x, y)

sup.run()

rms = sup.get_rms()

print('rms : ', rms)

rot, tran = sup.get_rotran()

print('rot, tran : ', rot, tran)

output:

rms :  0.0030430996905605216
rot :  [[ 0.68304944  0.53664494  0.4954349 ]
 [-0.5227743   0.83293146 -0.18147245]
 [-0.5100496  -0.13504598  0.8494775 ]]
tran :  [ 38.786064 -20.654562 -15.422249]

using y = x

output:

rms :  8.367921671944822e-08
rot :  [[ 1.0000000e+00 -1.9708821e-08 -5.0889071e-10]
 [-1.9708821e-08  1.0000000e+00  5.6950434e-08]
 [-5.0889071e-10  5.6950434e-08  1.0000000e+00]]
tran :  [ 0.0000000e+00 -1.9073486e-06  0.0000000e+00]