import ost as _ost
import ost.geom as _geom
"""
Module written by Niklaus Johner (niklaus.johner@a3.epfl.ch) 2012
This module is a flexible rewrite of HBPlus, allowing to calculate hbond
conservation between different structures or over a trajectory.
It uses customizable dictionaries to define donors and acceptors and can
account for equivalent HBonds such as involving for example either oxygen atom
of an Aspartic Acid.
"""
__all__=('HBondDonor','HBondAcceptor','BuildCHARMMHBondDonorAcceptorDict','AreHBonded','GetHbondDonorAcceptorList',\
'GetHbondListFromDonorAcceptorLists','GetHbondListFromView','GetHbondListFromTraj','GetHbondListBetweenViews'\
'CalculateHBondScore','AnalyzeHBondScore')
class HBondableAtoms:
def __init__(self,donors=[],acceptors=[]):
self.donors=donors
self.acceptors=acceptors
def BuildCHARMMHBondDonorAcceptorDict():
hb_da_dict={}
bb_donors=[['N','HN']]
bb_acceptors=[['O',['C']]]
#hydrophobic
hb_da_dict['ALA']=HBondableAtoms(bb_donors,bb_acceptors)
hb_da_dict['GLY']=HBondableAtoms(bb_donors,bb_acceptors)
hb_da_dict['ILE']=HBondableAtoms(bb_donors,bb_acceptors)
hb_da_dict['LEU']=HBondableAtoms(bb_donors,bb_acceptors)
hb_da_dict['PHE']=HBondableAtoms(bb_donors,bb_acceptors)
hb_da_dict['MET']=HBondableAtoms(bb_donors,bb_acceptors)
hb_da_dict['VAL']=HBondableAtoms(bb_donors,bb_acceptors)
#special case
hb_da_dict['PRO']=HBondableAtoms([],bb_acceptors)
#sidechain donors
ne=[['NE','HE']]
cz=[['NH1','HH11'],['NH1','HH12'],['NH2','HH21'],['NH2','HH22']]
hb_da_dict['ARG']=HBondableAtoms(bb_donors+ne+cz,bb_acceptors)
nz=[['NZ','HZ1'],['NZ','HZ2'],['NZ','HZ3']]
hb_da_dict['LYS']=HBondableAtoms(bb_donors+nz,bb_acceptors)
ne1=[['NE1','HE1']]
hb_da_dict['TRP']=HBondableAtoms(bb_donors+ne1,bb_acceptors)
sg=[['SG','HG1']]
hb_da_dict['CYS']=HBondableAtoms(bb_donors+sg,bb_acceptors)
#sidechain acceptors
od12=[['OD1',['CG']],['OD2',['CG']]]
hb_da_dict['ASP']=HBondableAtoms(bb_donors,bb_acceptors+od12)
oe12=[['OE1',['CD']],['OE2',['CD']]]
hb_da_dict['GLU']=HBondableAtoms(bb_donors,bb_acceptors+oe12)
#sidechain donor and acceptor
od1=[['OD1',['CG']]]
nd2=[['ND2','HD21'],['ND2','HD22']]
hb_da_dict['ASN']=HBondableAtoms(bb_donors+nd2,bb_acceptors+od1)
ne2=[['NE2','HE21'],['NE2','HE22']]
oe1=[['OE1',['CD']]]
hb_da_dict['GLN']=HBondableAtoms(bb_donors+ne2,bb_acceptors+oe1)
og_d=[['OG','HG1']]
og_a=[['OG',['CB']]]
hb_da_dict['SER']=HBondableAtoms(bb_donors+og_d,bb_acceptors+og_a)
og1_d=[['OG1','HG1']]
og1_a=[['OG1',['CB']]]
hb_da_dict['THR']=HBondableAtoms(bb_donors+og1_d,bb_acceptors+og1_a)
oh_d=[['OH','HH']]
oh_a=[['OH',['CZ']]]
hb_da_dict['TYR']=HBondableAtoms(bb_donors+oh_d,bb_acceptors+oh_a)
#histidine
nd1_d=[['ND1','HD1']]
ne2_a=[['NE2',['CD2','CE1']]]
hb_da_dict['HSD']=HBondableAtoms(bb_donors+nd1_d,bb_acceptors+ne2_a)
ne2_d=[['NE2','HE2']]
nd1_a=[['ND1',['CG','CE1']]]
hb_da_dict['HSE']=HBondableAtoms(bb_donors+ne2_d,bb_acceptors+nd1_a)
hb_da_dict['HSP']=HBondableAtoms(bb_donors+nd1_d+ne2_d,bb_acceptors)
#non-standard protonation state:
oe12=[['OE1',['CD']],['OE2',['CD']]]
oe2=[['OE2','HE2']]
hb_da_dict['GLUP']=HBondableAtoms(bb_donors+oe2,bb_acceptors+oe12)
od12=[['OD1',['CG']],['OD2',['CG']]]
od2=[['OD2','HD2']]
hb_da_dict['ASPP']=HBondableAtoms(bb_donors+od2,bb_acceptors+od12)
return hb_da_dict
def BuildCHARMMHBondDonorEquivalenceDict():
donor_swap_dict={}
donor_swap_dict['ARG']=[[['NH1','HH11'],['NH1','HH12'],['NH2','HH21'],['NH2','HH22']]]
donor_swap_dict['ASN']=[[['ND2','HD21'],['ND2','HD22']]]
donor_swap_dict['GLN']=[[['NE2','HE21'],['NE2','HE22']]]
donor_swap_dict['LYS']=[[['NZ','HZ1'],['NZ','HZ2'],['NZ','HZ3']]]
return donor_swap_dict
def BuildCHARMMHBondAcceptorEquivalenceDict():
acceptor_swap_dict={}
acceptor_swap_dict['ASP']=[[['OD1',['CG']],['OD2',['CG']]]]
acceptor_swap_dict['GLU']=[[['OE1',['CD']],['OE2',['CD']]]]
return acceptor_swap_dict
def ListEquivalentDonors(donor,donor_swap_dict):
if not donor.heavy_atom.residue.name in donor_swap_dict:return [donor]
donor_list=[donor]
donor_atom_names=[donor.heavy_atom.name,donor.hydrogen.name]
res=donor.heavy_atom.residue
for equivalence_list in donor_swap_dict[donor.heavy_atom.residue.name]:
if not donor_atom_names in equivalence_list:continue
for atom_names in equivalence_list:
if atom_names==donor_atom_names:continue
else:donor_list.append(HBondDonor.FromResidue(res,atom_names[0],atom_names[1]))
return donor_list
def ListEquivalentAcceptors(acceptor,acceptor_swap_dict):
if not acceptor.atom.residue.name in acceptor_swap_dict:return [acceptor]
acceptor_list=[acceptor]
acceptor_atom_names=[acceptor.atom.name,[a.name for a in acceptor.antecedent_list]]
res=acceptor.atom.residue
for equivalence_list in acceptor_swap_dict[acceptor.atom.residue.name]:
if not acceptor_atom_names in equivalence_list:continue
for atom_names in equivalence_list:
if atom_names==acceptor_atom_names:continue
else:acceptor_list.append(HBondAcceptor.FromResidue(res,atom_names[0],atom_names[1]))
return acceptor_list
class HBondDonor:
def __init__(self,donor,hydrogen):
self.heavy_atom=donor
self.hydrogen=hydrogen
def __eq__(self,hbd):
if not isinstance(hbd,HBondDonor):return False
else:return self.heavy_atom==hbd.heavy_atom and self.hydrogen==hbd.hydrogen
def __hash__(self):
return hash((self.heavy_atom,self.hydrogen))
@classmethod
def FromResidue(cls,res,donor_name,hydrogen_name,verbose=True):
_donor=res.FindAtom(donor_name).handle
_hydrogen=res.FindAtom(hydrogen_name).handle
if not _donor.IsValid():
if verbose:print 'Could not find '+donor_name+' in residue '+str(res)
return
if not _hydrogen.IsValid():
if verbose:print 'Could not find '+hydrogen_name+' in residue '+str(res)
return
return cls(_donor,_hydrogen)
def IsHBondedTo(self,acceptor):
return AreHBonded(self,acceptor)
class HBondAcceptor:
def __init__(self,acceptor,antecedent_list):
self.atom=acceptor
self.antecedent_list=antecedent_list
def __eq__(self,hba):
if not isinstance(hba,HBondAcceptor):return False
else:return self.atom==hba.atom
def __hash__(self):
return hash(self.atom)
@classmethod
def FromResidue(cls,res,acceptor_name,antecedent_name_list,verbose=True):
_acceptor=res.FindAtom(acceptor_name).handle
_antecedent_list=_ost.mol.AtomHandleList([res.FindAtom(name).handle for name in antecedent_name_list])
if not _acceptor.IsValid():
if verbose:print 'Could not find '+acceptor_name+' in residue '+str(res)
return
for i,a in enumerate(_antecedent_list):
if not a.IsValid():
if verbose:print 'Could not find '+antecedent_name_list[i]+' in residue '+str(res)
return
return cls(_acceptor,_antecedent_list)
def IsHBondedTo(self,donor):
return AreHBonded(donor,self)
class HBond:
def __init__(self,donor,acceptor):
self.donor=donor
self.acceptor=acceptor
def __eq__(self,hb):
if not isinstance(hb,HBond):return False
else:return self.acceptor==hb.acceptor and self.donor==hb.donor
def __hash__(self):
return hash((self.donor,self.acceptor))
def IsFormed(self):return AreHBonded(self.donor,self.acceptor)
def DonorAcceptorDistance(self):return _geom.Distance(self.donor.heavy_atom.pos,self.acceptor.atom.pos)
def HydrogenAcceptorDistance(self):return _geom.Distance(self.donor.hydrogen.pos,self.acceptor.atom.pos)
def DonorHydrogenAcceptorAngle(self):
dp=self.donor.heavy_atom.pos;ap=self.acceptor.atom.pos;hp=self.donor.hydrogen.pos
return _geom.Angle(dp-hp,ap-hp)
def DonorAcceptorAcceptorAntecedentAngle(self):
dp=self.donor.heavy_atom.pos;ap=self.acceptor.atom.pos
a=min([_geom.Angle(aa.pos-ap,dp-ap) for aa in self.acceptor.antecedent_list])
return a
def HydrogenAcceptorAcceptorAntecedentAngle(self):
hp=self.donor.hydrogen.pos;ap=self.acceptor.atom.pos
a=min([_geom.Angle(aa.pos-ap,hp-ap) for aa in acceptor.antecedent_list])
return a
def AreHBonded(donor,acceptor,da_dist=3.9,ha_dist=2.5,dha_angle=1.57,daaa_angle=1.57,haaa_angle=1.57):
"""
determines if a donor/acceptor pair is hydrogen bonded or not
"""
dp=donor.heavy_atom.pos
ap=acceptor.atom.pos
if _geom.Distance(dp,ap)>da_dist:return False
hp=donor.hydrogen.pos
if _geom.Distance(hp,ap)>ha_dist:return False
if _geom.Angle(dp-hp,ap-hp)<dha_angle:return False
a=min([_geom.Angle(aa.pos-ap,dp-ap) for aa in acceptor.antecedent_list])
if a<daaa_angle:return False
a=min([_geom.Angle(aa.pos-ap,hp-ap) for aa in acceptor.antecedent_list])
if a<haaa_angle:return False
return True
def GetHbondDonorAcceptorList(eh,hbond_donor_acceptor_dict={},verbose=True):
"""
returns a list of hydrogen-bond donors and acceptors from an Entity or EntityView.
It relies on atom names to determine the list of H-bond donors and acceptors.
These names are given in a dictionary, which defaults to CHARMM.
"""
if not hbond_donor_acceptor_dict:
print 'Using default CHARMM atom namings'
hbond_donor_acceptor_dict=BuildCHARMMHBondDonorAcceptorDict()
donor_list=[]
acceptor_list=[]
for r in eh.residues:
if not r.name in hbond_donor_acceptor_dict:
print 'donors and acceptors for',r,'are not defined in the dictionary and will not be included'
continue
res_da_dict=hbond_donor_acceptor_dict[r.name]
for acceptor in res_da_dict.acceptors:
a=HBondAcceptor.FromResidue(r,acceptor[0],acceptor[1],verbose)
if not a==None:acceptor_list.append(a)
for donor in res_da_dict.donors:
d=HBondDonor.FromResidue(r,donor[0],donor[1],verbose)
if not d==None:donor_list.append(d)
return [donor_list,acceptor_list]
def GetHbondListFromDonorAcceptorLists(donor_list,acceptor_list):
"""
return a list of hydrogen bonds between donors and acceptors from
a list of donors and a list of acceptors.
"""
hbond_list=[]
for donor in donor_list:
for acceptor in acceptor_list:
if AreHBonded(donor,acceptor):hbond_list.append(HBond(donor,acceptor))
return hbond_list
def GetHbondListFromView(eh,hbond_donor_acceptor_dict={},verbose=True):
"""
return a list of hydrogen bonds from an Entity or EntityView.
if no dictionary for the hbond donors and acceptors is specified
it will use the standard CHARMM names to determine them
"""
[donor_list,acceptor_list]=GetHbondDonorAcceptorList(eh,hbond_donor_acceptor_dict,verbose)
hbond_list=[]
for donor in donor_list:
for acceptor in acceptor_list:
if AreHBonded(donor,acceptor):hbond_list.append(HBond(donor,acceptor))
return hbond_list
def GetHbondListFromTraj(t,eh,cutoff=0.7,stride=1,swap=False,donor_swap_dict={},acceptor_swap_dict={},hbond_donor_acceptor_dict={},verbose=True):
"""
return a list of hydrogen bonds from an Entity or EntityView that are
present in a fraction of the frames larger than *cutoff*.
if no dictionary for the hbond donors and acceptors is specified
it will use the standard CHARMM names to determine them
"""
if swap:
if not donor_swap_dict:
print 'use of standard CHARMM HBond donor swap dictionary'
donor_swap_dict=BuildCHARMMHBondDonorEquivalenceDict()
if not acceptor_swap_dict:
print 'use of standard CHARMM HBond acceptor swap dictionary'
acceptor_swap_dict=BuildCHARMMHBondAcceptorEquivalenceDict()
[donor_list,acceptor_list]=GetHbondDonorAcceptorList(eh,hbond_donor_acceptor_dict,verbose)
hb_list=[]
hb_score=[]
nframes=0
for i in range(0,t.GetFrameCount(),stride):
t.CopyFrame(i)
nframes+=1
hbond_list=GetHbondListFromDonorAcceptorLists(donor_list,acceptor_list)
if swap:
hbond_list=GetEquivalentHBonds(hbond_list,eh,swap,donor_swap_dict,acceptor_swap_dict,verbose)
for hb in hbond_list:
if hb in hbond_list[hbond_list.index(hb)+1:]:hbond_list.pop(hbond_list.index(hb))
for hb in hbond_list:
if hb in hb_list:hb_score[hb_list.index(hb)]+=1
else:
hb_score.append(1)
hb_list.append(hb)
hbond_list=[]
hbond_score=[]
for hb,s in zip(hb_list,hb_score):
if s/float(nframes)>=cutoff:
if swap:hbond_list.append(list(hb)[0])
else:hbond_list.append(hb)
hbond_score.append(s/float(nframes))
return (hbond_list,hbond_score)
def GetHbondListBetweenViews(eh1,eh2,hbond_donor_acceptor_dict={},verbose=True):
"""
return the list of hydrogen bonds formed between two Entity or EntityView.
if no dictionary for the hbond donors and acceptors is specified
it will use the standard CHARMM names to determine them
"""
[donor_list1,acceptor_list1]=GetHbondDonorAcceptorList(eh1,hbond_donor_acceptor_dict,verbose)
[donor_list2,acceptor_list2]=GetHbondDonorAcceptorList(eh2,hbond_donor_acceptor_dict,verbose)
hbond_list=[]
for donor in donor_list1:
for acceptor in acceptor_list2:
if AreHBonded(donor,acceptor):
hb=HBond(donor,acceptor)
if hb in hbond_list:continue
hbond_list.append(hb)
for donor in donor_list2:
for acceptor in acceptor_list1:
if AreHBonded(donor,acceptor):
hb=HBond(donor,acceptor)
if hb in hbond_list:continue
hbond_list.append(hb)
return hbond_list
def GetEquivalentHBonds(ref_hbond_list,eh,swap=False,donor_swap_dict={},acceptor_swap_dict={},verbose=True):
hbond_list=[]
if not swap:
for hbond in ref_hbond_list:
res1=eh.FindResidue(hbond.donor.heavy_atom.chain.name,hbond.donor.heavy_atom.residue.number.num)
res2=eh.FindResidue(hbond.acceptor.atom.chain.name,hbond.acceptor.atom.residue.number.num)
donor=HBondDonor.FromResidue(res1,hbond.donor.heavy_atom.name,hbond.donor.hydrogen.name,verbose)
acceptor=HBondAcceptor.FromResidue(res2,hbond.acceptor.atom.name,[a.name for a in hbond.acceptor.antecedent_list],verbose)
hbond_list.append(set([HBond(donor,acceptor)]))
else:
if not donor_swap_dict:
print 'use of standard CHARMM HBond donor swap dictionary'
donor_swap_dict=BuildCHARMMHBondDonorEquivalenceDict()
if not acceptor_swap_dict:
print 'use of standard CHARMM HBond acceptor swap dictionary'
acceptor_swap_dict=BuildCHARMMHBondAcceptorEquivalenceDict()
for hbond in ref_hbond_list:
res1=eh.FindResidue(hbond.donor.heavy_atom.chain.name,hbond.donor.heavy_atom.residue.number.num)
res2=eh.FindResidue(hbond.acceptor.atom.chain.name,hbond.acceptor.atom.residue.number.num)
donor=HBondDonor.FromResidue(res1,hbond.donor.heavy_atom.name,hbond.donor.hydrogen.name,verbose)
acceptor=HBondAcceptor.FromResidue(res2,hbond.acceptor.atom.name,[a.name for a in hbond.acceptor.antecedent_list],verbose)
donor_list=ListEquivalentDonors(donor,donor_swap_dict)
acceptor_list=ListEquivalentAcceptors(acceptor,acceptor_swap_dict)
hbond_list.append(set([HBond(d,a) for d in donor_list for a in acceptor_list]))
return hbond_list
def CalculateHBondScore(ref_eh,eh2,ref_eh2=None,hbond_donor_acceptor_dict={},swap=False,donor_swap_dict={},acceptor_swap_dict={},verbose=True):
"""
Returns the fraction of H-bonds from ref_eh that are also present in eh2.
If ref_eh2 is specified, it uses as reference the Hbonds between ref_eh and ref_eh2.
This allows to look at H-bonds between specific parts of proteins or so.
Alternatively ref_eh can be a list of H-bonds.
This function relies on atom names to determine the list of H-bond donors and acceptors.
These names are given in a dictionary, which defaults to CHARMM.
If swap is set to True, a dictionary for equivalent donors and one for equivalent acceptors
(defaults to CHARMM) is used to check for equivalent HBonds in eh2.
If swap is set to True, if two equivalent hydrogen bonds are present in the reference entity
(for example both oxygens of ASP H-bonding the same atom), it suffices that on of these bonds is
present in eh2 for both of them to be counted as present in eh2.
"""
if ref_eh2:hbond_list1=GetHbondListBetweenViews(ref_eh,ref_eh2,hbond_donor_acceptor_dict,verbose)
elif type(ref_eh)==list:hbond_list1=ref_eh
else:hbond_list1=GetHbondListFromView(ref_eh,hbond_donor_acceptor_dict,verbose)
nbonds=float(len(hbond_list1))
if nbonds==0:
print 'No HBonds in reference view'
return None
hbond_list2=GetEquivalentHBonds(hbond_list1,eh2,swap,donor_swap_dict,acceptor_swap_dict,verbose)
c=0
for hl in hbond_list2:
for hbond in hl:
if HBond(hbond.donor,hbond.acceptor).IsFormed():
c+=1
break
return c/float(len(hbond_list1))
def AnalyzeHBondScore(ref_eh,t,eh2,ref_eh2=None,hbond_donor_acceptor_dict={},swap=False,donor_swap_dict={},acceptor_swap_dict={},first=0,last=-1,stride=1,verbose=True):
"""
Returns the same score as CalculateHBondScore, but for a trajectory.
"""
if swap:
if not donor_swap_dict:
print 'use of standard CHARMM HBond donor swap dictionary'
donor_swap_dict=BuildCHARMMHBondDonorEquivalenceDict()
if not acceptor_swap_dict:
print 'use of standard CHARMM HBond acceptor swap dictionary'
acceptor_swap_dict=BuildCHARMMHBondAcceptorEquivalenceDict()
if ref_eh2:hbond_list1=GetHbondListBetweenViews(ref_eh,ref_eh2,hbond_donor_acceptor_dict,verbose)
elif type(ref_eh)==list:hbond_list1=ref_eh
else:hbond_list1=GetHbondListFromView(ref_eh,hbond_donor_acceptor_dict,verbose)
nbonds=float(len(hbond_list1))
if nbonds==0:
print 'No HBonds in reference view'
return None
print 'number of hbonds in ref_eh:',nbonds
hbond_list2=GetEquivalentHBonds(hbond_list1,eh2,swap,donor_swap_dict,acceptor_swap_dict,verbose)
if last==-1:last=t.GetFrameCount()
score=FloatList()
for f in range(first,last,stride):
t.CopyFrame(f)
c=0
for hl in hbond_list2:
for hbond in hl:
if hbond.IsFormed():
c+=1
break
score.append(c/nbonds)
return score
def GetHBondListIntersection(ref_hbond_list,ref_eh,hbond_list,swap=False,donor_swap_dict={},acceptor_swap_dict={}):
if swap:
if not donor_swap_dict:
print 'use of standard CHARMM HBond donor swap dictionary'
donor_swap_dict=BuildCHARMMHBondDonorEquivalenceDict()
if not acceptor_swap_dict:
print 'use of standard CHARMM HBond acceptor swap dictionary'
acceptor_swap_dict=BuildCHARMMHBondAcceptorEquivalenceDict()
hbond_list=GetEquivalentHBonds(hbond_list,ref_eh,swap,donor_swap_dict,acceptor_swap_dict)
ref_hbond_list=GetEquivalentHBonds(ref_hbond_list,ref_eh,swap,donor_swap_dict,acceptor_swap_dict)
out_hbond_list=[]
for hb in ref_hbond_list:
if hb in hbond_list:
out_hbond_list.append(hb)
return out_hbond_list