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