contactMap

From p3d
Jump to: navigation, search

ContactMapMini.py

A contact map shows the distance and residues between two peptide chains.

The contact map script creates, as the name already might suggest, the data to plot a contact map. In order to keep it short, the plotting was not included and, in this case was done using ProFit.

usage

contactSurfaceMini.py <pdb file> <peptide Chain 1> <peptide Chain 2> <max distance>

peptide chain 1 and peptide chain 2 are the two chains for which the contact map data should be generated and max distance is the maximal distance that should be taken into account.

The output is e.g.

138 	 113 	 2.86066390896 
113 	 138 	 2.58556550874
136 	 117 	 2.75362397578
117 	 136 	 2.84670985525
249 	 122 	 2.8287354065
117 	 138 	 2.87679491796
134 	 247 	 2.89939062563
247 	 134 	 2.73693185885
122 	 249 	 2.85885186745
114 	 138 	 2.7297146371

where the first column is residue number for chain 1, the second column is residue number for chain 2 and the last column is the minimum distance between those residues in column 1 and 2.

example

Violaxanthine deepoxidase dimeric and active form was used to create this example.

maxDistance = 7Å
maxDistance = 14Å
maxDistance = 21Å

code

#!/usr/bin/env python
'''contact surface by Ch. Fufezan 2009

usage: contactSurface.py <pdb file> <chain1> <chain2> <maxDistance to take into account>'''

import sys, os
from p3d import protein as protein
from collections import defaultdict as ddict

def estimateMinContacts(Protein=None,ResidueAtom=None,targetChain=None,DistanceMatrix=None,maxDistance=None):
	residueAtoms = Protein.query('protein and resid {0} and chain {1}'.format(ResidueAtom.resid,ResidueAtom.chain))
	surounding = Protein.query('protein and within {0} of'.format(maxDistance),residueAtoms,'and chain {0}'.format(targetChain))
	for residAtom in residueAtoms:
		for atom in surounding:
			d = atom.distanceTo(residAtom)
			ordered_tuple = (residAtom.resid,atom.resid) if targetChain == 'B' else (atom.resid,residAtom.resid)
			if DistanceMatrix[ordered_tuple] > d or DistanceMatrix[ordered_tuple] == 0:
				DistanceMatrix[ordered_tuple] = d
	return DistanceMatrix

if (__name__ == '__main__'):
	if (len(sys.argv) != 5):
		print (__doc__)
		sys.exit(1)
	pdb = protein.Protein(sys.argv[1],DunbrackNaming=False)
	chain_A , chain_B, maxDistance = sys.argv[2:5]
	DistanceMatrix = ddict(float)
	''' query Set alphas to have one atom per residue from both chains '''
	chain_A_alphas = pdb.query('chain {0} and alpha'.format(chain_A))
	chain_B_alphas = pdb.query('chain {0} and alpha'.format(chain_B))
	for alpha in chain_A_alphas:
		DistanceMatrix = estimateMinContacts(Protein=pdb,ResidueAtom=alpha,targetChain=chain_B,DistanceMatrix=DistanceMatrix,maxDistance=maxDistance)
	for alpha in chain_B_alphas:
		DistanceMatrix = estimateMinContacts(Protein=pdb,ResidueAtom=alpha,targetChain=chain_A,DistanceMatrix=DistanceMatrix,maxDistance=maxDistance)
	for coords,value in DistanceMatrix.items():
		print(coords[0],'\t',coords[1],'\t',value)