#!/usr/bin/env python
# -*-coding:Utf-8 -*

#########################  Calcul des distances entre atomes (E. Jaspard - 2018) ########################
#Script inspire de https://warwick.ac.uk/fac/sci/moac/people/students/peter_cock/python/protein_carteDeContacts/ 

import Bio.PDB
from Bio.PDB.PDBParser import PDBParser
import numpy

codePDB = "1XI4"
nomDuFichier = "1XI4.pdb"

###############################   1ere partie    ############################### 
###   Definition d'une fonction "distanceCarbonesAlpha" qui renvoie la distance entre les carbones alpha de 2 residus 
def distanceCarbonesAlpha(residuUn, residuDeux) :
    vecteurDifference = residuUn["CA"].coord - residuDeux["CA"].coord   #["CA"] = selection des lignes du fichier PDB qui contiennent les coordonnees des carbones alpha 
    return numpy.sqrt(numpy.sum(vecteurDifference * vecteurDifference))

###   Definition d'une fonction qui utilise cette distance pour calculer une matrice de distance entiere
def matriceDistance(chaineUn, chaineDeux) :    # Matrice de distance des carbones alpha entre les deux chaines "chaineUn" et "chaineDeux" 
    resultat = numpy.zeros((len(chaineUn), len(chaineDeux)), numpy.float)
    for lignes, residuUn in enumerate(chaineUn) :
        for colonnes, residuDeux in enumerate(chaineDeux) :
            resultat[lignes, colonnes] = distanceCarbonesAlpha(residuUn, residuDeux)
    return resultat

###############################   2e partie    ###############################   
# Le module "Bio.PDB.PDBParser()" permet d'effectuer les opérations suivantes :
# Lecture du fichier PDB : les informations nécessaires sont stockees dans la variable "structure".  
# Cela fournit une liste des modeles (= les differentes chaines) contenus dans le fichier PDB.
decoupage = PDBParser()
structure = decoupage.get_structure(codePDB, nomDuFichier)   # "codePDB" et "nomDuFichier" = variables definies au debut
# Ces 2 lignes sont equivalentes a : structure = Bio.PDB.PDBParser().get_structure(codePDB, nomDuFichier)
modele = structure[0]

###############################   3e partie - matrice de distances    ###############################   
# Il y a 9 modeles de chaines lourdes et 9 modeles de chaines legeres.
# Le script calcule la matrice de distances entre la chaine lourde D et la legere chaine M car elles sont en contact étroit.
matriceDistance = matriceDistance(modele["D"], modele["M"])

# Puis transformation de la matrice de distances (sous forme d'un tableau de réels flottants)
# en carte de contacts (sous forme d'un tableau logique) en utilisant un seuil de 12 angstroms.
carteDeContacts = matriceDistance < 12.0

print "Distance minimale =", numpy.min(matriceDistance)   # Resultat : 4.62857341766
print "Distance maximale =", numpy.max(matriceDistance)   # Resultat : 201.617782593

# Remarque : les residus les plus proches sont Asn155 (chaine D) et Leu1504 (chaine M).

###############################   4e partie  -  affichage  ###############################   
import pylab
pylab.matshow(numpy.transpose(matriceDistance))
pylab.colorbar()
pylab.show()