Residue contact map

Question

How to use Biopython to calculate the protein residue contact map from a PDBfile,and the final contact map is storage in Numpy array as a .npy file.

if you can, pls give me some explicit guidence or examples . thanx for your help.

Answer 1

Ok , my attempt , given as Contact Map (from Contact maps) :

A protein contact map represents the distances between all possible aminoacid resid pairs in the form of a matrix. Usually contact maps are defined in a binary manner: two residues are in contact if their distances is below a certain threshold, and not in contact otherwise..

and calc_dist_matrix from Protein Contact Maps using Biopython

input 1rei.pdb , 1REI :

Code :

import Bio

print('Biopython version : ', Bio.__version__)


from Bio import PDB

import numpy as np 



def calc_dist_matrix(chain_one, chain_two , cutoff) :
    """Returns a matrix of C-alpha distances between two chains"""
    answer = np.zeros((len(chain_one), len(chain_two)), np.intc)
    for row, residue_one in enumerate(chain_one) :
        for col, residue_two in enumerate(chain_two) :
            
            if (residue_one['CA'] - residue_two['CA']) <= cutoff :
                
                distance = 1
                
            else :
                
                distance = 0

            answer[row, col] = distance
            
    return answer



parser = PDB.PDBParser(QUIET = True)

pdb1 ='1rei.pdb' 

print(pdb1.split('.pdb')[0])


structure = parser.get_structure(pdb1.split('.pdb')[0], pdb1)


resi1 = []

for chain in structure[0]:
    
    resi1.extend([res for res in chain.get_residues() if res.id[0] == ' '])


print(resi1, len(resi1))


resi2 = resi1.copy()

cutoff = 20.0

contact_map = calc_dist_matrix(resi1 , resi2 , cutoff)

print('contact map : ', contact_map.shape, contact_map.size,contact_map.dtype,
      
      np.min(contact_map), np.max(contact_map), type(contact_map) )

from PIL import Image, ImageOps

import time

img = Image.fromarray(contact_map* 255, mode='L').convert('1')


img.show()

time.sleep(4)

img_inv = ImageOps.invert(img)

img_inv.show()

using a cutoff of 20.0 A and Cα-Cα as distance between residues

here is the image output with PIL , my map is :

`contact map :  (214, 214) 45796 int32 0 1 <class 'numpy.ndarray'>` ;

grayscale (0 or 255)

inverted