crmsd_on_ligands.py

Calculates crmsd value on all atoms of a ligand conformations.

This scripts reads one or more PDB files, extracts all ligands that maches a given three-letter code and calculates crmsd between them on all atoms

USAGE:
python3 crmsd_on_ligands.py HEM 5ofq.pdb 4rm4.pdb

Categories:

  • core/calc/structural/transformations/CrmsdonVec3

Program source:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
import sys, math

sys.path.append('../../../../../bin/')
from pybioshell.core.data.io import Pdb
from pybioshell.std import vector_core_data_basic_Vec3
from pybioshell.core.calc.structural.transformations import *
from pybioshell.utils import LogManager
LogManager.INFO()

if len(sys.argv) < 3 :
    print("""

    Calculates crmsd value on all atoms of a ligand conformations.
    
    This scripts reads one or more PDB files, extracts all ligands that maches a given 
    three-letter code and calculates crmsd between them on all atoms

USAGE:
    python3 crmsd_on_ligands.py HEM 5ofq.pdb 4rm4.pdb

    CATEGORIES: core/calc/structural/transformations/CrmsdonVec3
    KEYWORDS:   PDB input; ligand; crmsd

  """)
    sys.exit()

rms = CrmsdOnVec3()

ligands = []
resids = []
pdb_codes = []
structures = []
for pdb_code in sys.argv[2:] :
  pdb = Pdb(pdb_code,"",False)
  structure = pdb.create_structure(0)
  structures.append(structure)
  for ic in range(structure.count_chains()) :
    chain = structure[ic]
    for ir in range(chain.terminal_residue_index() + 1,chain.size()) :
        resid = chain[ir]
        code3 = resid.residue_type().code3
        if code3 != sys.argv[1] : continue # Skip other ligands
        resids.append(resid)
        pdb_codes.append(pdb_code)
        atoms = vector_core_data_basic_Vec3()
        for ia in range(resid.count_atoms()):
            atoms.append(resid[ia])
        ligands.append(atoms)

for i_ligand in range(0, len(ligands)) :
      ir = resids[i_ligand]
      for j_ligand in range(i_ligand):
          jr = resids[j_ligand]
          crmsd_val = rms.crmsd(ligands[i_ligand],ligands[j_ligand], len(ligands[j_ligand]) )
          print("%s %4d %3s %c - %s %4d %3s %c : %7.3f" % 
            (pdb_codes[i_ligand], ir.id(), ir.residue_type().code3, ir.owner().id(), 
            pdb_codes[j_ligand], jr.id(), jr.residue_type().code3, jr.owner().id(), crmsd_val))
../_images/file_icon.png