DistancesRH

Introduction

Given an object, a reference amino acid and a radius, this scripts calculates pairwise distances between atom(s) in the reference and all atoms that fall within the given range.

The script is intended to provide a range of useful tools while requiring minimal coding skills (i.e. it is not optimised for efficiency).

By changing input parameters it is possible to:

• choose the output (none, on screen, on file)
• measure distances from one single atom
• display distance objects in pymol
• restrict distances between hydrogen-bond forming atoms only

The script also runs a number of controls to verify that inputs and outputs exists/are suitable.

Usage

distancesRH obj, ref, dist, [chain, [output S/P/N, [show Y/N, [Hbonds Y/N, [aname]]]]]

Note #1: the selected molecule must be visible. Note #2: soon a version that manages multiple states.

Required Arguments

• obj = object name
• ref = reference amino acid id
• dist = radius in Angstroms

Optional Arguments

• chain = selected chain
• output = accepts S (screen), P (print) or N (none), default=S
• show = shows distances in pymol window, default=N
• Hbonds = restricts results to atoms that can form Hbonds, default=N
• aname = selects a specific atom in reference object and to calculate distances from one individual atom

Examples

example #1

distancesRH 2WHH, 25, 3.2

Will show distances between all atoms of residue 25 in 2WHH and any atom within 3.2 Å and return the result on screen.

example #2

distancesRH 1LVM, 151, 3, B, P, Y, Y, N

Will show distances between CB of residue 151 in chain B of 1LVM and all atoms within 3 Å. Distances are limited to atoms that can form hydrogen bonds and will be printed in "dist_rH_151_N.txt". Dist objects will be shown in pymol.

The Code

from pymol import cmd, stored, math
	
def distancesRH (obj, ref, dist, chain=0, output="S", show="N", Hbonds="N", aname="*"):
	"""	
	usage: distancesRH obj, ref, dist, [chain, [output S/P/N, [show Y/N, [Hbonds Y/N, [aname]]]]]
	
	obj: object name, ref:reference aa, dist: radius in A, chain: selected chain
	output: accepts S (screen), P (print) or N (none), default=S
	show: shows distances in pymol window, default=N
	Hbonds: restricts results to atoms that can form Hbonds, default=N
	aname: name of a specific atom in ref (to calculate distances from one atom only)
		
	example: distancesRH 2WHH, 25, 3.2, [B, [S, [N, [Y, [CB]]]]]
	"""
	print ""
	#delete previous selections and displayed distances
	cmd.delete ("dist*")
	cmd.delete ("Atoms")
	cmd.delete ("Residues")
	cmd.delete ("Reference")
	cmd.delete ("Hbonds")
	cmd.delete ("RefAtom")
	#creates and names selections	
	if chain==0: cen=obj
	else: cen=obj+" and chain "+chain
	refA=" and resi "+ref+" and name "+aname
	cmd.select("Reference", "byres "+cen+refA)
	m1=cmd.get_model ("Reference")
	if aname!="*":
		cmd.select("RefAtom", cen+refA)
		m1=cmd.get_model ("RefAtom")
	cmd.select("Atoms", cen+refA+" around "+str(dist)+" and not resi "+ref+" and "+cen)
	cmd.show("nb_spheres", "(Atoms or Reference) and (resn HOH or HET)")
	cmd.select("Residues", "byres Atoms")
	m2=cmd.get_model("Atoms and visible")	
	if Hbonds=="Y":
		cmd.select("__Hbonds1", "Reference and not symbol c")
		if aname!="*":
			cmd.select("__Hbonds1", "RefAtom and not symbol c")
		cmd.select("__Hbonds2", "Atoms and not symbol c")
		cmd.select("Hbonds", "__Hbonds1 or __Hbonds2")
		cmd.show("lines", "Hbonds")
		m1=cmd.get_model ("__Hbonds1")
		m2=cmd.get_model ("__Hbonds2 and visible")	
	
	#controllers	
	if (not (output=="S" or output=="N" or output=="P")) or (not (show=="Y" or show=="N")) or (not (Hbonds=="Y" or Hbonds=="N")):
		print "Malformed input, please try again"
		return
	abort=1
	abort2=0
	mc=cmd.get_model(cen)
	for c0 in range(len(mc.atom)):
		if mc.atom[c0].resi==ref: abort=0
	mc2=cmd.get_model(cen+refA)
	if aname!="*":
		abort2=1
		for c00 in range (len(mc2.atom)):	
			if mc2.atom[c00].name==aname: abort2=0
	if Hbonds=="Y":
		if aname[0]=="C":
			print "Warning: no atoms in the selection can form H-bonds"
			return
	if abort==1: print "Warning: no such residue in the reference object"
	if abort2==1: 
		if abort==0: print "Warning: no such atom in the reference object"
	if abort==1 or abort2==1: return

	#measures distances
	distance2=[]
	distances=[]
	screenOP=""	
	fileOP=""
	for c1 in range(len(m1.atom)):
		for c2 in range(len(m2.atom)):
			if math.sqrt(sum(map(lambda f: (f[0]-f[1])**2, zip(m1.atom[c1].coord,m2.atom[c2].coord))))<float(dist):
				distance=cmd.distance (obj+"//"+m1.atom[c1].chain+"/"+m1.atom[c1].resi+"/"+m1.atom[c1].name, obj+"//"+m2.atom[c2].chain+"/"+m2.atom[c2].resi+"/"+m2.atom[c2].name)
				distances.append(distance)
				screenOP+="%s/%s/%s/%s to %s/%s/%s/%s: %.3f\n" % (m1.atom[c1].chain,m1.atom[c1].resn,m1.atom[c1].resi,m1.atom[c1].name,m2.atom[c2].chain,m2.atom[c2].resn,m2.atom[c2].resi,m2.atom[c2].name, distance)
				fileOP+="%s/%s/%s/%s\t\t\t%s/%s/%s/%s\t\t\t\t%.5f\n"%  (m1.atom[c1].chain,m1.atom[c1].resn,m1.atom[c1].resi,m1.atom[c1].name,m2.atom[c2].chain,m2.atom[c2].resn,m2.atom[c2].resi,m2.atom[c2].name, distance)	
	#last controller
	if len(distances)==0:
		if abort2==0:
			print "Warning: no atoms within the selected range"
		return

	#data outputs
	if output=="N": print "Data recorded, but not printed"
	if output=="S":
		print "Distance of atoms in "+obj+"/"+ref+"/"+aname+" to all atoms within "+dist+" Angstroms"
		if Hbonds=="Y": print "Distances restricted to potential H-bond pairs"
		print screenOP	
	if output=="P":
		print 'data printed in "distances_rH_'+ref+'-'+aname+'.txt" file'
		f=open('distances_rH_'+ref+'-'+aname+'.txt','w')
		f.write("Pairwise distances betweet atoms in %s and all atoms within %s Angstroms\n" % (obj+"/"+ref+"/"+aname, dist))
		if Hbonds=="Y": f.write("Distances restricted to potential H-bond pairs\n")
		f.write("Atom in reference\tAtom in destination\t\tdistance\n")
		f.write(fileOP)
		f.close()
	print "All done! Good luck with your data"

	#graphical output
	if show=="N": cmd.delete ("dist*")
	if show=="Y":
		cmd.show("lines", "Reference")
		cmd.show("lines", "Residues")
				
cmd.extend("distancesRH", distancesRH);