PDBUtils.py

import math

def is_number(s):
	try:
		float(s)
		return True
	except ValueError:
		return False

class Atom:
	def __init__(self,typ,coords,line):
		self.typ = typ
		self.coords = coords
		self.pdb_line = line

def atom_distance(a1,a2):
		dist = euclidean_distance(a1.coords[0]-a2.coords[0],a1.coords[1]-a2.coords[1],a1.coords[2]-a2.coords[2])
		return dist
#Euclidean distance when the difference in coordinates has already been submitted
#x= (x_0 - x_1) etc...
def euclidean_distance(x,y,z):
		
		result = math.sqrt(x*x+y*y+z*z)
		return result

#Class to hold a single chain from the PDB
class Residue:
	def __init__(self,atoms,res_type):
		self.atoms = atoms
		self.res_type = res_type
	
				
	def distance(self,res1,res2):
		mindist = 99999999
		for a1 in res1.atoms:
			for a2 in res2.atoms:
				dist = self.euclid(a1.coords[0]-a2.coords[0],a1.coords[1]-a2.coords[1],a1.coords[2]-a2.coords[2])
				if dist<mindist:
					mindist = dist
					
		return mindist
	
	def get_coords(self,aname):
		mindist = 99999999
		for a1 in self.atoms:
			if a1.typ!=aname:
				continue
			else:
				return (a1.coords[0],a1.coords[1],a1.coords[2])
	def euclid(self,x,y,z):
		result = math.sqrt(x*x+y*y+z*z)
		
		return result

class PDBchain:
	
	#Get the contact mappings between the two structures
	def contact_map(self,str1,str2):
		mapping = dict()
		types = dict()
		for chain_1 in str1:
			for chain_2 in str2:
				for res_id_1 in str1[chain_1].residues:
					if len(res_id_1)==0:
						continue
					r_1 = str1[chain_1].residues[res_id_1]
					index_id_1 = chain_1+str(res_id_1[0])+res_id_1[1]
					if index_id_1 not in mapping:
						mapping[index_id_1] = dict()
					types[index_id_1] = r_1.res_type
					for res_id_2 in str2[chain_2].residues:
						if len(res_id_2)==0:
							continue
						index_id_2 = chain_2+str(res_id_2[0])+res_id_2[1]
						r_2 = str2[chain_2].residues[res_id_2]
						types[index_id_2] = r_2.res_type
						d = r_1.distance(r_1,r_2)
						if d<11.5:
							mapping[index_id_1][index_id_2] = d
	
		return mapping,types
	#Change the b-factor on the pdb line
	def color_line(self,line,bf):
		while(len(bf))<6:
			bf=" "+bf
		line = line[0:60]+bf+line[66:len(line)]
			
		return line
	
	def __init__(self,pdb_file,which_chain):
		
		try:
			residues = dict()
			atoms = []
			curr_id = ""
			res_type = ""
			for line in pdb_file:
				

				if line[0:4]!='ATOM':
					continue
				
				line = line.strip()
				#print "Chain:",line[21]
				chain = line[21]
				if chain!=which_chain:
					continue
				#print "Sid:",line[23:28].replace(" ","")
				sid = line[22:28].replace(" ","")
				if  is_number(sid):
					sid = (int(sid),'')
				else:
					sid = (int(sid[0:len(sid)-1]),sid[len(sid)-1])
				if len(curr_id) ==0:
					curr_id = sid
				#print "atom: ",line[12:17].replace(" ","")
				atom_type = 	line[12:17].replace(" ","")			
				
				#print "x:",float(line[30:38].replace(" ",""))
				x = float(line[30:38].replace(" ",""))
				#print "y:",float(line[38:46].replace(" ",""))
				y = float(line[38:46].replace(" ",""))
				#print "z:",float(line[46:54].replace(" ",""))
				z = float(line[46:54].replace(" ",""))
				atom = Atom(atom_type,(x,y,z),line)				
				#Push the previous residue				
				if sid!=curr_id:
					#push the residue
					#print "Pushing residue",sid
					res = Residue(atoms,res_type)
					residues[curr_id] = res
					curr_id = sid
					atoms = [atom]
				else:
					atoms.append(atom)
				#print "Res: ",line[17:21].replace(" ","")
				res_type = line[17:21].replace(" ","")

			#Last residue			
			res = Residue(atoms,res_type)
			residues[curr_id] = res
			#print "Initialized with ",len(residues)," residues"
			self.residues = residues		
		except IOError:
			print "IOError: ",path_to_pdb,"does not exist"
			quit()





if __name__=='__main__':
	print 'running main'
	path_to_pdb = '/databases/cpdb/database/1L7F.pdb'
	pdbchain = PDBchain(path_to_pdb,'A')
	ligands = LigandAssembly(path_to_pdb)
	str1 = {'A':pdbchain}
	str2 = {'L':ligands}
	maps,typs= pdbchain.contact_map(str1,str2)
	print maps