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

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

Note: It shows values even if we compare a file with itself - in PyMOL, it is not possible to compare a file with itself

Answer 1

Edit: CA does not refer to Calcium atoms. They refer to alpha carbon atoms as defined in the PDB format.

To get an accurate RMSD value between two molecules they would need to be first re-centered and then rotated unto each other to get the true minimal RMSD.

The code that you posted compares only aligns alpha carbon atoms. This can be seen in the following lines of the code:

# Select what residues numbers you wish to align
# and put them in a list
start_id = 1
end_id   = 70

########
# 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 = []

#########

ref_atoms.append(ref_res['CA'])
sample_atoms.append(sample_res['CA'])

To get an accurate RMSD value all of the atoms would need to re-centered and rotated accordingly. Not just alpha carbon atoms. Comparing only the alpha atoms could be a bug in case you were after the true minimal RMSD.

I used the script posted above with the same PDB file of https://www.rcsb.org/structure/1a6j and it printed out the RMSD to be

RMSD =  1.2115634513933915e-14

When I used the rmsd library, it gives me a better (true) RMSD value:

RMSD =  6.8963066943027425e-15