OpenStructure
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
hbplus.py
Go to the documentation of this file.
1 """
2 Interface to HB+ command line program
3 
4 Author: Marco Biasini
5 """
6 import tempfile
7 from ost import settings
8 import subprocess
9 import re
10 import os
11 import shutil
12 from ost import io
13 from ost import mol
14 
15 def _LocateHBPlus(hbplus_bin):
16  return settings.Locate('hbplus', explicit_file_name=hbplus_bin,
17  env_name="HBPLUS")
18 
19 class HBond:
20  def __init__(self, donor, acceptor):
21  self.donor=donor
22  self.acceptor=acceptor
23 
24 def _ParseOutput(ent, output):
25  hbonds=[]
26  # skip header
27  lines=list(output)
28  for index, line in enumerate(lines):
29  if line.startswith('<---DONOR---> <-ACCEPTOR-->'):
30  lines=lines[index+4:]
31  break
32  for line in lines:
33  if len(line.strip())==0:
34  continue
35  don_chain=line[0]
36  don_rnum=int(line[1:5])
37  don_ins_c=line[5]
38  don_atom=line[9:13]
39  if don_chain=='-':
40  don_chain=' '
41  if don_ins_c=='-':
42  don_ins_c='\0'
43  acc_chain=line[14]
44  acc_rnum=int(line[15:19])
45  acc_ins_c=line[19]
46  acc_atom=line[24:28]
47  if acc_chain=='-':
48  acc_chain=' '
49  if acc_ins_c=='-':
50  acc_ins_c='\0'
51  donor=ent.FindAtom(don_chain, mol.ResNum(don_rnum, don_ins_c),
52  don_atom.strip())
53  acc=ent.FindAtom(acc_chain, mol.ResNum(acc_rnum, acc_ins_c),
54  acc_atom.strip())
55  assert donor.IsValid()
56  assert acc.IsValid()
57  hbonds.append(HBond(donor, acc))
58  return hbonds
59 
60 def HBondList(ent, hbplus_bin=None):
61  """
62  returns a list of HBonds found in the given entity (handle or view)
63  """
64  full_bin=_LocateHBPlus(hbplus_bin)
65  temp_d=tempfile.mkdtemp(prefix='hbplus_')
66  file_name=os.path.join(temp_d, 'ent.pdb')
67  io.SaveEntity(ent, file_name)
68  hb_proc=subprocess.Popen(full_bin, shell=True, stdout=subprocess.PIPE,
69  stdin=subprocess.PIPE)
70  hb_proc.stdin.write(('%s\n' % temp_d).encode())
71  hb_proc.stdin.write(('%s\n\n' % file_name).encode())
72  stdout,_ = hb_proc.communicate()
73 
74  for line in stdout.decode().splitlines():
75  match=re.match(r'Configured for (\d+) atoms and\s+(\d+) residues\.', line)
76  if match:
77  assert ent.atom_count<int(match.group(1))
78  assert ent.residue_count<int(match.group(2))
79  hb_out=open(os.path.join(temp_d, 'ent.hb2'), 'r')
80  hbonds=_ParseOutput(ent, hb_out)
81  hb_out.close()
82  shutil.rmtree(temp_d)
83  return hbonds
84 
85 def HBondScore(ent1, ent2, hbplus_bin=None):
86  """
87  returns percentage of hydrogen bonds in ent1 that are also present in ent2 in
88  the range of 0 to 1.
89 
90  this function is slow as hell, and could be improved drastically by first
91  sorting the hydrogen bond lists by donor residue number.
92  """
93 
94  hbonds_a=HBondList(ent1, hbplus_bin=hbplus_bin)
95  hbonds_b=HBondList(ent2, hbplus_bin=hbplus_bin)
96  matching=0
97  for a in hbonds_a:
98  found=False
99  for b in hbonds_b:
100  acc_names_match=b.acceptor.name==a.acceptor.name
101  don_names_match=b.donor.name==a.donor.name
102  don_num=b.donor.residue.number==a.donor.residue.number
103  acc_num=b.acceptor.residue.number==a.acceptor.residue.number
104  if acc_names_match and don_names_match and acc_num and don_num:
105  found=True
106  break
107  if found:
108  matching+=1
109  if len(hbonds_a)>0:
110  return float(matching)/len(hbonds_a)
111  else:
112  return 0.0
113 
114 if __name__=='__main__':
115  print('HBond Score:', HBondScore(io.LoadPDB(sys.argv[1]),
116  io.LoadPDB(sys.argv[2])))