3

I have written a python script to plot the 'Ramachandran Plot' of Ubiquitin protein. I am using biopython. I am working with pdb files. My script is as below :

import Bio.PDB
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

phi_psi = ([0,0])
phi_psi = np.array(phi_psi)
pdb1 ='/home/devanandt/Documents/VMD/1UBQ.pdb'

for model in Bio.PDB.PDBParser().get_structure('1UBQ',pdb1) :
    for chain in model :
        polypeptides = Bio.PDB.PPBuilder().build_peptides(chain)
        for poly_index, poly in enumerate(polypeptides) :
            print "Model %s Chain %s" % (str(model.id), str(chain.id)),
            print "(part %i of %i)" % (poly_index+1, len(polypeptides)),
            print "length %i" % (len(poly)),
            print "from %s%i" % (poly[0].resname, poly[0].id[1]),
            print "to %s%i" % (poly[-1].resname, poly[-1].id[1])
            phi_psi = poly.get_phi_psi_list()
            for res_index, residue in enumerate(poly) :
                #res_name = "%s%i" % (residue.resname, residue.id[1])
                #print res_name, phi_psi[res_index]
                phi_psi = np.vstack([phi_psi \
                ,np.asarray(phi_psi[res_index])]).astype(np.float) 
                #np.float - conversion to float array from object

phi, psi = np.transpose(phi_psi)

phi = np.degrees(phi)
psi = np.degrees(psi)

phi = phi[~np.isnan(phi)] # avoiding nan
psi = psi[~np.isnan(psi)]

f,ax = plt.subplots(1)
plt.title('Ramachandran Plot for Ubiquitin')
plt.xlabel('$\phi^o$', size=20,fontsize=15)
plt.ylabel('$\psi^o$ ', size=20,fontsize=15)

h=ax.hexbin(phi, psi,  extent=[-180,180,-180,180],cmap=plt.cm.Blues)
#h=ax.hexbin(phi, psi,gridsize=35,  extent=[-180,180,-180,180],cmap=plt.cm.Blues)

f.colorbar(h)
plt.grid()
plt.show()

I would like to modify this code so as to neglect the GLYCINE amino acid and then plot Ramachandran plot. My output is as below: enter image description here

dexterdev
  • 537
  • 4
  • 22

1 Answers1

5

You can remove them after indexing the GLYs:

for poly_index, poly in enumerate(polypeptides):
    gly_index = [i for i, res in enumerate(poly) if res.get_resname() == "GLY"]

After the main loop and phy/psi calculation, delete the points from the array:

new_phi_psi = np.delete(phi_psi, gly_index, 0)
phi, psi = np.transpose(new_phi_psi)

Remove the step where you get rid of the NaNs. Now plot the points to get something like this:

h=ax.hexbin(phi, psi, extent=[-180,180,-180,180],cmap=plt.cm.Blues)
h=ax.hexbin(n_phi, n_psi, extent=[-180,180,-180,180],cmap=plt.cm.Reds)

The darker points comes from the full phi/psi with GLY residues

xbello
  • 7,223
  • 3
  • 28
  • 41
  • I got the plot, I would like to verify it using VMD or so? But using the inbuilt ramaplot in VMD can I remove the GLY and see the plot? – dexterdev Oct 06 '14 at 17:21
  • 1
    I'm not sure of what you want. Read http://www2.warwick.ac.uk/fac/sci/moac/people/students/peter_cock/python/ramachandran/drawing/ – xbello Oct 06 '14 at 17:56
  • 1
    I did this some time ago, using a external programs. From Molprobity I got the "chiropraxis.jar" library. Launch the jar with `java -cp chiropraxis.jar chiropraxis.rotarama.Ramalyze -pdf YourPDB output.pdf`. Chiropraxis download link: http://kinemage.biochem.duke.edu/software/dangle.php – xbello Oct 06 '14 at 18:04
  • What I want is compare Ramachandran plots of ubiquitin with and without glycine. The without glycine image is here : http://imgur.com/aKqPFwc . It seems that more points are there in histogram now. But how? Any way I forgot to say thanks for your effort. – dexterdev Oct 06 '14 at 18:20
  • 1
    I guess you can generate two sets of points to graph with matplotlib, and assign a different color to each one. Here I just overlapped the two graphs (the original in your question and the new one you just uploaded to imgur): http://imgur.com/mlD29IN. – xbello Oct 06 '14 at 18:30
  • I deleted my answer, because as I was writing, my solution and xbello's will alter the phi and psi data. – Roberto Oct 06 '14 at 21:59
  • @Roberto : So you are saying that both of your methods alter phi and psi data. Now I am confused. – dexterdev Oct 07 '14 at 02:51
  • @xbello : I am getting those plots, but what about Roberto's thoughts. – dexterdev Oct 07 '14 at 02:52
  • 1
    I'm updating the answer. – xbello Oct 07 '14 at 06:46