Python：计算两个蛋白或小分子之间的RMSD

Python脚本：计算两个蛋白或小分子之间的RMSD

用法：

python rmsd.py protein1.pdb protein2.pdb

rmsd.py

  

   

    

     

    
    

     

      # Root-mean-square deviation (RMSD) for proteins and/or ligands
     
    
   

    

     

    
    

     

      # in PDB files.
     
    
   

    

     

    
    

     

      #
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      class Pdb(object):
     
    
   

    

     

    
    

     
 """ Object that allows operations with protein files in PDB format. """
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 def __init__(self, file_cont = [], pdb_code = ""):
     
    
   

    

     

    
    

     

       self.cont = []
     
    
   

    

     

    
    

     

       self.atom = []
     
    
   

    

     

    
    

     

       self.hetatm = []
     
    
   

    

     

    
    

     

       self.fileloc = ""
     
    
   

    

     

    
    

     
 if isinstance(file_cont, list):
     
    
   

    

     

    
    

     

       self.cont = file_cont[:]
     
    
   

    

     

    
    

     
 elif isinstance(file_cont, str):
     
    
   

    

     

    
    

     
 try:
     
    
   

    

     

    
    

     
 with open(file_cont, 'r') as pdb_file:
     
    
   

    

     

    
    

     

       self.cont = [row.strip() for row in pdb_file.read().split('\n') if row.strip()]
     
    
   

    

     

    
    

     
 except FileNotFoundError as err:
     
    
   

    

     

    
    

     

       print(err)
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 if self.cont:
     
    
   

    

     

    
    

     

       self.atom = [row for row in self.cont if row.startswith('ATOM')]
     
    
   

    

     

    
    

     

       self.hetatm = [row for row in self.cont if row.startswith('HETATM')]
     
    
   

    

     

    
    

     

       self.conect = [row for row in self.cont if row.startswith('CONECT')]
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 def rmsd(self, sec_molecule, ligand=False, atoms="no_h"):
     
    
   

    

     

    
    

     
 """
     
    
   

    

     

    
    

     

       Calculates the Root Mean Square Deviation (RMSD) between two
     
    
   

    

     

    
    

     

       protein or ligand molecules in PDB format.
     
    
   

    

     

    
    

     

       Requires that both molecules have the same number of atoms in the
     
    
   

    

     

    
    

     

       same numerical order.
     
    
   

    

     

    
    

     

       Keyword arguments:
     
    
   

    

     

    
    

     

       sec_molecule (PdbObj): the second molecule as PdbObj object.
     
    
   

    

     

    
    

     

       ligand (bool): If true, calculates the RMSD between two
     
    
   

    

     

    
    

     

       ligand molecules (based on HETATM entries), else RMSD
     
    
   

    

     

    
    

     

       between two protein molecules (ATOM entries) is calculated.
     
    
   

    

     

    
    

     

       hydrogen (bool): If True, hydrogen atoms will be included in the
     
    
   

    

     

    
    

     

       RMSD calculation.
     
    
   

    

     

    
    

     

       atoms (string) [all/c/no_h/ca]: "all" includes all atoms in the RMSD calculation,
     
    
   

    

     

    
    

     

       "c" only considers carbon atoms, "no_h" considers all but hydrogen atoms,
     
    
   

    

     

    
    

     

       and "ca" compares only C-alpha protein atoms.
     
    
   

    

     

    
    

     

       Returns:
     
    
   

    

     

    
    

     

       Calculated RMSD value as float or None if RMSD not be
     
    
   

    

     

    
    

     

       calculated.
     
    
   

    

     

    
    

     

       """
     
    
   

    

     

    
    

     

       rmsd = None
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 if not ligand:
     
    
   

    

     

    
    

     

       coords1, coords2 = self.atom, sec_molecule.atom
     
    
   

    

     

    
    

     
 else:
     
    
   

    

     

    
    

     

       coords1, coords2 = self.hetatm, sec_molecule.hetatm
     
    
   

    

     

    
    

     
 if atoms == "c":
     
    
   

    

     

    
    

     

       coords1 = [row for row in coords1 if row[77:].startswith('C')]
     
    
   

    

     

    
    

     

       coords2 = [row for row in coords2 if row[77:].startswith('C')]
     
    
   

    

     

    
    

     
 elif atoms == "no_h":
     
    
   

    

     

    
    

     

       coords1 = [row for row in coords1 if not row[77:].startswith('H')]
     
    
   

    

     

    
    

     

       coords2 = [row for row in coords2 if not row[77:].startswith('H')]
     
    
   

    

     

    
    

     
 elif atoms == "ca":
     
    
   

    

     

    
    

     

       coords1 = self.calpha()
     
    
   

    

     

    
    

     

       coords2 = sec_molecule.calpha()
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 if all((coords1, coords2, len(coords1) == len(coords2))):
     
    
   

    

     

    
    

     

       total = 0
     
    
   

    

     

    
    

     
 for (i, j) in zip(coords1, coords2):
     
    
   

    

     

    
    

     

       total += ( float(i[30:38]) - float(j[30:38]) )**2 +\
     
    
   

    

     

    
    

     

       ( float(i[38:46]) - float(j[38:46]) )**2 +\
     
    
   

    

     

    
    

     

       ( float(i[46:54]) - float(j[46:54]) )**2 
     
    
   

    

     

    
    

     

       rmsd = round(( total / len(coords1) )**0.5, 4)
     
    
   

    

     

    
    

     
 return rmsd
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      if __name__ == '__main__':
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 import argparse
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

       parser = argparse.ArgumentParser(
     
    
   

    

     

    
    

     

       description='The RMSD measures the average distance between atoms \n'\
     
    
   

    

     

    
    

     
 'of 2 protein or ligand structures.\n'\
     
    
   

    

     

    
    

     
 'By default, all atoms but hydrogen atoms of the protein are included in the RMSD calculation.\n'\
     
    
   

    

     

    
    

     
 'NOTE: Both structures must contain the same number of atoms in similar order.',
     
    
   

    

     

    
    

     

       epilog='Example:\n'\
     
    
   

    

     

    
    

     
 'rmsd.py ~/Desktop/pdb1.pdb ~/Desktop/pdb2.pdb\n'\
     
    
   

    

     

    
    

     
 '0.7377',
     
    
   

    

     

    
    

     

       formatter_class=argparse.RawTextHelpFormatter
     
    
   

    

     

    
    

     

       )
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

       parser.add_argument('PDBfile1')
     
    
   

    

     

    
    

     

       parser.add_argument('PDBfile2')
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

       parser.add_argument('-l', '--ligand', action='store_true', help='Calculates RMSD between ligand (HETATM) atoms.')
     
    
   

    

     

    
    

     

       parser.add_argument('-c', '--carbon', action='store_true', help='Calculates the RMSD between carbon atoms only.')
     
    
   

    

     

    
    

     

       parser.add_argument('-ca', '--calpha', action='store_true', help='Calculates the RMSD between alpha-carbon atoms only.')
     
    
   

    

     

    
    

     

       parser.add_argument('-v', '--version', action='version', version='rmsd v. 1.0')
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

       args = parser.parse_args()
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

       pdb1 = Pdb(args.PDBfile1)
     
    
   

    

     

    
    

     

       pdb2 = Pdb(args.PDBfile2)
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 if args.carbon and args.calpha:
     
    
   

    

     

    
    

     

       print('\nERROR: Please provide EITHER -c OR -ca, not both.\n')
     
    
   

    

     

    
    

     

       parser.print_help()
     
    
   

    

     

    
    

     

       quit()
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 if args.ligand and args.carbon:
     
    
   

    

     

    
    

     

       print(pdb1.rmsd(sec_molecule=pdb2, ligand=True, atoms="c"))
     
    
   

    

     

    
    

     
 elif args.ligand and args.calpha:
     
    
   

    

     

    
    

     

       print(pdb1.calpha())
     
    
   

    

     

    
    

     

       print(pdb1.rmsd(sec_molecule=pdb2, ligand=True, atoms="ca"))
     
    
   

    

     

    
    

     
 elif args.ligand:
     
    
   

    

     

    
    

     

       print(pdb1.rmsd(sec_molecule=pdb2, ligand=True, atoms="no_h"))
     
    
   

    

     

    
    

     
 elif args.calpha:
     
    
   

    

     

    
    

     

       print(pdb1.rmsd(sec_molecule=pdb2, ligand=False, atoms="ca"))
     
    
   

    

     

    
    

     
 elif args.carbon:
     
    
   

    

     

    
    

     

       print(pdb1.rmsd(sec_molecule=pdb2, ligand=False, atoms="c"))
     
    
   

    

     

    
    

     
 else:
     
    
   

    

     

    
    

     

       print(pdb1.rmsd(sec_molecule=pdb2, ligand=False, atoms="no_h"))
     
    
  
 

文章来源: drugai.blog.csdn.net，作者：DrugAI，版权归原作者所有，如需转载，请联系作者。

原文链接：drugai.blog.csdn.net/article/details/81008659