BbPlane: Difference between revisions

From PyMOLWiki
Jump to navigation Jump to search
No edit summary
(support multiple chains and proline)
Line 28: Line 28:
from pymol.cgo import *    # get constants
from pymol.cgo import *    # get constants
from pymol import cmd, stored
from pymol import cmd, stored
from chempy import cpv


def bbPlane(objSel, r=1.0, g=1.0, b=1.0, transp=0.0):
def bbPlane(objSel='(all)', color='white', transp=0.0):
     """
     """
     Draws a plane across the backbone for a selection
     Draws a plane across the backbone for a selection
Line 37: Line 38:
         (object or selection) protein object or selection
         (object or selection) protein object or selection


     r, g, b, transp
     color = string: color name or number {default: white}
        (float) red, blue, green and transparency, components (0.0--1.0)
 
    transp = float: transparency component (0.0--1.0) {default: 0.0}


     RETURNS
     RETURNS
Line 51: Line 53:
     """
     """
     # format input
     # format input
     r, g, b, transp = float(r), float(g), float(b), float(transp)
     transp = float(transp)
     stored.AAs = []
     stored.AAs = []


Line 58: Line 60:
      
      
     # get the list of residue ids
     # get the list of residue ids
     cmd.iterate( "(" + objSel + ") and n. CA" , "stored.AAs.append(resi)")
     for obj in cmd.get_object_list(objSel):
        cmd.iterate( obj + " and (" + objSel + ") and n. CA" , "stored.AAs.append('/"+obj+"/%s/%s/%s' % (segi,chain,resi))")


     # need at least two amino acids
     # need at least two amino acids
Line 68: Line 71:
     obj = [
     obj = [
         BEGIN, TRIANGLES,
         BEGIN, TRIANGLES,
         COLOR, r, g, b,
         COLOR,
         ]
         ]
    obj.extend(cmd.get_color_tuple(color))


     for res in range(0, len(stored.AAs)-1):
     for res in range(0, len(stored.AAs)-1):
         curIdx, nextIdx = str(stored.AAs[res]), str(stored.AAs[res+1])
         curIdx, nextIdx = str(stored.AAs[res]), str(stored.AAs[res+1])


         # ignore prolines for now
        corners = [
         if cmd.count_atoms( "i. %s and resn PRO" % nextIdx ):
            curIdx + "/CA",
             continue
            curIdx + "/O",
            nextIdx + "/N extend 1 and e. H",
            nextIdx + "/CA"
            ]
 
         # prolines
         if cmd.count_atoms( "%s and resn PRO" % nextIdx ):
             corners[2] = nextIdx + "/CD"
 
         # if the data are incomplete for any residues, ignore
         # if the data are incomplete for any residues, ignore
         if cmd.count_atoms( "i. " + curIdx + " and n. CA") == 0 or \
         if map(cmd.count_atoms, corners) != [1,1,1,1]:
                cmd.count_atoms( "i. " + curIdx + " and n. O") == 0 or \
            print 'peptide bond %s -> %s incomplete' % (curIdx, nextIdx)
                cmd.count_atoms("((i. " + nextIdx + " and n. N) extend 1) and (e. H)")==0 or \
                cmd.count_atoms("i. " + nextIdx + " and n. CA")== 0:
             continue
             continue
             
 
         # populate the position array
         # populate the position array
         pos = [ cmd.get_atom_coords( "i. " + curIdx + " and n. CA"),
         pos = map(cmd.get_atom_coords, corners)
                cmd.get_atom_coords( "i. " + curIdx + " and n. O"),
 
                cmd.get_atom_coords( "((i. " + nextIdx + " and n. N) extend 1) and (e. H)"),
        if cpv.distance(pos[0], pos[3]) > 4.0:
                cmd.get_atom_coords( "i. " + nextIdx + " and n. CA")  
            print '%s and %s not adjacent' % (curIdx, nextIdx)
                ]
            continue


         # fill in the vertex data for the triangles;  
         # fill in the vertex data for the triangles;  

Revision as of 04:29, 30 June 2010

This script will draw a CGO plane between the backbone atoms of two neighboring residues. This is to show the planarity of the atoms. The image style this is meant to represent can be found many places, like "Introduction to Protein Structure" by Branden and Tooze (2nd ed. pp. 8).

Examples

# download the source and save as bbPlane.py
run bbPlane.py
fetch 1cll
# make planes for residues 4-9
bbPlane i. 4-10

The Source

#
# -- bbPLane.py - draws a CGO plane across the backbone atoms of
#                 neighboring amino acids
# 
# Author: Jason Vertrees, 06/2010
# Copyright (C) Schrodinger
# Open Source License: MIT
#
from pymol.cgo import *    # get constants
from pymol import cmd, stored
from chempy import cpv

def bbPlane(objSel='(all)', color='white', transp=0.0):
    """
    Draws a plane across the backbone for a selection

    PARAMS
    objSel,  
        (object or selection) protein object or selection

    color = string: color name or number {default: white}

    transp = float: transparency component (0.0--1.0) {default: 0.0}

    RETURNS
    A new object representing the planes

    NOTES
    * You need to pass in an object or selection with at least two
    amino acids.  The plane spans CA_i, O_i, N-H_(i+1), and CA_(i+1)

    * Avoid two-sided lighting until I rearrange the vertices

    """
    # format input
    transp = float(transp)
    stored.AAs = []

    # need hydrogens
    cmd.h_add(objSel)
    
    # get the list of residue ids
    for obj in cmd.get_object_list(objSel):
        cmd.iterate( obj + " and (" + objSel + ") and n. CA" , "stored.AAs.append('/"+obj+"/%s/%s/%s' % (segi,chain,resi))")

    # need at least two amino acids
    if len(stored.AAs) <= 1:
        print "ERROR: Please provide at least two amino acids, the alpha-carbon on the 2nd is needed."
        return

    # prepare the cgo
    obj = [
        BEGIN, TRIANGLES,
        COLOR,
        ]
    obj.extend(cmd.get_color_tuple(color))

    for res in range(0, len(stored.AAs)-1):
        curIdx, nextIdx = str(stored.AAs[res]), str(stored.AAs[res+1])

        corners = [
            curIdx + "/CA",
            curIdx + "/O",
            nextIdx + "/N extend 1 and e. H",
            nextIdx + "/CA"
            ]

        # prolines
        if cmd.count_atoms( "%s and resn PRO" % nextIdx ):
            corners[2] = nextIdx + "/CD"

        # if the data are incomplete for any residues, ignore
        if map(cmd.count_atoms, corners) != [1,1,1,1]:
            print 'peptide bond %s -> %s incomplete' % (curIdx, nextIdx)
            continue

        # populate the position array
        pos = map(cmd.get_atom_coords, corners)

        if cpv.distance(pos[0], pos[3]) > 4.0:
            print '%s and %s not adjacent' % (curIdx, nextIdx)
            continue

        # fill in the vertex data for the triangles; 
        # two-sided lighting is an issue
        obj.extend( [
            VERTEX, pos[0][0], pos[0][1], pos[0][2],
            VERTEX, pos[1][0], pos[1][1], pos[1][2],
            VERTEX, pos[2][0], pos[2][1], pos[2][2],

            VERTEX, pos[1][0], pos[1][1], pos[1][2],
            VERTEX, pos[2][0], pos[2][1], pos[2][2],            
            VERTEX, pos[3][0], pos[3][1], pos[3][2],
            ])

    # finish the CGO
    obj.append(END)

    # update the UI
    newName =  cmd.get_unused_name("backbonePlane")
    cmd.load_cgo(obj, newName)
    cmd.set("cgo_transparency", transp, newName)


cmd.extend("bbPlane", bbPlane)