4
$\begingroup$

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

$\endgroup$
3
  • 3
    $\begingroup$ Normally, there are three ways to compute the RMSD: considering the structures "as it" (i.e. keep them in their current coordinates), centering the structures and rotating/aligning the structures (example). Even when you align two structures, the alignment should be exactly the same. $\endgroup$
    – Camps
    Oct 21, 2022 at 17:13
  • 1
    $\begingroup$ Another possible source of difference would be if one software uses raw cartesian coordinates and the other superimposes in mass-weighted coordinates. Not sure which people prefer in this context though. $\endgroup$
    – jheindel
    Nov 3, 2022 at 4:54
  • 1
    $\begingroup$ Just to confirm, how different are the RMSD values you are seeing? If it's on the order of what is shown in Vandan Revanur's answer below, this shouldn't make any functional difference. Note the script you are using also seems to be at least partially tailored to specific PDB files. For example, the range of residue numbers to consider is hard coded. $\endgroup$
    – Tyberius
    Feb 20, 2023 at 13:31

1 Answer 1

3
$\begingroup$

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
$\endgroup$
2
  • 1
    $\begingroup$ 1. Both the RMSDs shown are very small, nearly on the order of machine precision. The different arise from different methodologies, but that small a difference in RMSD is almost certainly not meaningful. 2. CA does not refer to calcium, but rather to an alpha carbon, as defined by the PDB format $\endgroup$
    – Tyberius
    Feb 19, 2023 at 22:50
  • $\begingroup$ @Tyberius thanks for pointing it out the mistake with CA atoms. In any case it would be better to use rmsd library instead of the above code right ? $\endgroup$ Feb 20, 2023 at 6:30

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .