00001 """
00002 Interface to HB+ command line program
00003
00004 Author: Marco Biasini
00005 """
00006 import tempfile
00007 from ost import settings
00008 import subprocess
00009 import re
00010 import os
00011 from ost import io
00012 from ost import mol
00013
00014 def _LocateHBPlus(hbplus_bin):
00015 return settings.Locate('hbplus', explicit_file_name=hbplus_bin,
00016 env_name="HBPLUS")
00017
00018 class HBond:
00019 def __init__(self, donor, acceptor):
00020 self.donor=donor
00021 self.acceptor=acceptor
00022
00023 def _ParseOutput(ent, output):
00024 hbonds=[]
00025
00026 lines=list(output)
00027 for index, line in enumerate(lines):
00028 if line.startswith('<---DONOR---> <-ACCEPTOR-->'):
00029 lines=lines[index+4:]
00030 break
00031 for line in lines:
00032 if len(line.strip())==0:
00033 continue
00034 don_chain=line[0]
00035 don_rnum=int(line[1:5])
00036 don_ins_c=line[5]
00037 don_atom=line[9:13]
00038 if don_chain=='-':
00039 don_chain=' '
00040 if don_ins_c=='-':
00041 don_ins_c='\0'
00042 acc_chain=line[14]
00043 acc_rnum=int(line[15:19])
00044 acc_ins_c=line[19]
00045 acc_atom=line[24:28]
00046 if acc_chain=='-':
00047 acc_chain=' '
00048 if acc_ins_c=='-':
00049 acc_ins_c='\0'
00050 donor=ent.FindAtom(don_chain, mol.ResNum(don_rnum, don_ins_c),
00051 don_atom.strip())
00052 acc=ent.FindAtom(acc_chain, mol.ResNum(acc_rnum, acc_ins_c),
00053 acc_atom.strip())
00054 assert donor.IsValid()
00055 assert acc.IsValid()
00056 hbonds.append(HBond(donor, acc))
00057 return hbonds
00058
00059 def HBondList(ent, hbplus_bin=None):
00060 """
00061 returns a list of HBonds found in the given entity (handle or view)
00062 """
00063 full_bin=_LocateHBPlus(hbplus_bin)
00064 temp_d=tempfile.mkdtemp(prefix='hbplus_')
00065 hb_proc=subprocess.Popen(full_bin, shell=True, stdout=subprocess.PIPE,
00066 stdin=subprocess.PIPE)
00067 file_name=os.path.join(temp_d, 'ent.pdb')
00068 io.SaveEntity(ent, file_name)
00069 hb_proc.stdin.write('%s\n' % temp_d)
00070 hb_proc.stdin.write('%s\n\n' % file_name)
00071 for line in hb_proc.stdout:
00072 match=re.match(r'Configured for (\d+) atoms and\s+(\d+) residues\.', line)
00073 if match:
00074 assert ent.atom_count<int(match.group(1))
00075 assert ent.residue_count<int(match.group(2))
00076 hb_out=open(os.path.join(temp_d, 'ent.hb2'), 'r')
00077 hbonds=_ParseOutput(ent, hb_out)
00078 os.system('rm -rf "%s"' % temp_d)
00079 return hbonds
00080
00081 def HBondScore(ent1, ent2, hbplus_bin=None):
00082 """
00083 returns percentage of hydrogen bonds in ent1 that are also present in ent2 in
00084 the range of 0 to 1.
00085
00086 this function is slow as hell, and could be improved drastically by first
00087 sorting the hydrogen bond lists by donor residue number.
00088 """
00089
00090 hbonds_a=HBondList(ent1, hbplus_bin=hbplus_bin)
00091 hbonds_b=HBondList(ent2, hbplus_bin=hbplus_bin)
00092 matching=0
00093 for a in hbonds_a:
00094 found=False
00095 for b in hbonds_b:
00096 acc_names_match=b.acceptor.name==a.acceptor.name
00097 don_names_match=b.donor.name==a.donor.name
00098 don_num=b.donor.residue.number==a.donor.residue.number
00099 acc_num=b.acceptor.residue.number==a.acceptor.residue.number
00100 if acc_names_match and don_names_match and acc_num and don_num:
00101 found=True
00102 break
00103 if found:
00104 matching+=1
00105 if len(hbonds_a)>0:
00106 return float(matching)/len(hbonds_a)
00107 else:
00108 return 0.0
00109
00110 if __name__=='__main__':
00111 print 'HBond Score:', HBondScore(io.LoadPDB(sys.argv[1]),
00112 io.LoadPDB(sys.argv[2]))
