PyMOL Wiki - User contributions [en]

Contact map visualizer

2021-11-17T09:34:20Z

Zhentg: /* Generate Contact Map */

{{Infobox script-repo
|type = plugin
|filename = plugins/contact_map_visualizer.py
|author = [[User:Venky|Venkatramanan Krishnamani]]
|license = -
}}

'''Enhanced version of this plugin is now available at [[CMPyMOL]]'''

The '''contact map visualizer''' plugin can link contact map images to the residues in PyMOL in a interactive way. Contact maps are pixel graphics where each protein residue corresponds to one line and one row of pixels. Thus for a 100 residue protein, such a image has 100x100 pixels. A common tool to generate such images is '''g_mdmat''' from the [http://www.gromacs.org gromacs] package.

== Usage ==

contact_map_visualizer [ image_file [, selection ]]

== Required Dependencies ==

* [http://pygame.org/ pygame]
* [http://wiki.python.org/moin/TkInter Tkinter] (optional and usually included with PyMOL)
* [http://www.pythonware.com/products/pil/ PIL] (optional, for automatically converting XPM images)

Example for installing all dependencies on a Ubuntu like system:

<source lang="bash">
sudo apt-get install python-tk python-imaging python-pygame
</source>

== Installation ==

* Navigate to plugins > install
* Locate the downloaded '''contact_map_visualizer.py''' file in the dialogbox and select 'OK'
* Quit and Restart 'pymol'

Alternative way: Just [[run]] the script, it will provide a command but no menu plugin entry.

== Generate Contact Map ==

[[File:Contact-Map-of-a-Trajectory.png|thumb|200px|right|Mean contact map of a protein trajectory generated from g_mdmat tool in the Gromacs analysis package.]]

Use the command [http://manual.gromacs.org/online/g_mdmat.html g_mdmat] from [http://www.gromacs.org Gromacs] analysis package. A typical contact map looks like the figure on the right.

To generate contact map of a single PDB. For example contact map for a PDB from RCSB, use the following command
g_mdmat -f <protein.pdb> -s <protein.pdb> -mean contact-map.xpm
To generate a mean contact map form a protein trajectory
g_mdmat -f <trajectory.pdb> -s <starting-frame.pdb> -mean contact-map.xpm

For Gromacs 2021
To generate contact map of a single PDB. For example contact map for a PDB from RCSB, use the following command
gmx g_mdmat -f <protein.pdb> -s <protein.pdb> -mean contact-map.xpm
To generate a mean contact map form a protein trajectory
gmx g_mdmat -f <trajectory.pdb> -s <starting-frame.pdb> -mean contact-map.xpm
To convert XPM to PNG format
convert contact-map.xpm contact-map.png

== Screenshot ==

[[File:Screenshot.png|thumb|800px|center]]

== Copyright ==
<source lang="text">
# Copyright Notice
# ================
#
# The PyMOL Plugin source code in this file is copyrighted, but you are
# free to use and copy it as long as you don't change or remove any of
# the copyright notices.
#
# -----------------------------------------------------------------------------------
# This PyMOL Plugin Contact Maps Visualizer is
# Copyright (C) 2012 by Venkatramanan Krishnamani <venks@andrew.cmu.edu>
#
# All Rights Reserved
#
# Permission to use, copy, modify, distribute, and distribute modified
# versions of this software and its documentation for any purpose and
# without fee is hereby granted, provided that the above copyright
# notice appear in all copies and that both the copyright notice and
# this permission notice appear in supporting documentation, and that
# the name(s) of the author(s) not be used in advertising or publicity
# pertaining to distribution of the software without specific, written
# prior permission.
#
# THE AUTHOR(S) DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
# INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN
# NO EVENT SHALL THE AUTHOR(S) BE LIABLE FOR ANY SPECIAL, INDIRECT OR
# CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF
# USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
# OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
# PERFORMANCE OF THIS SOFTWARE.
#
</source>

[[Category:Plugins]]
[[Category:Pymol-script-repo]]

ColorByDistance

2018-11-12T10:08:15Z

Zhentg: /* The Code */

{{Infobox script-repo
|type = script
|author = [[User:zhentg|Zhenting Gao]]
}}
== Introduction ==
Color the surface of a binding site by the distance to a specific atom.

== Usage ==
*Open PyMOL
*Load PDB files
*run this Python script inside PyMOL
*call the function
**colorByDistance anchor, bindingSite

== Required Arguments ==
Text

== Optional Arguments ==
Text

== Examples ==
Text
== The Code ==
<source>
def colorByDistance(anchor='anchor', site='site'):
#Based on: https://pymolwiki.org/index.php/Spectrum#Intermediate
#Author: Zhenting Gao (zhentgpicasa@gmail.com)
#Update: 11/12/2018
#Aim: Color atoms of a binding site based on their distance from a point
#May also refer to https://pymolwiki.org/index.php/Ramp_New#Ramp_.2B_Distance_Measure
# returns the length of the distance between atom A and atom B

diff_len = lambda x,y : math.sqrt((x[0]-y[0])*(x[0]-y[0]) + (x[1]-y[1])*(x[1]-y[1]) + (x[2]-y[2])*(x[2]-y[2]))

# fetch site from the PDB

#fetch site, async=0

# show it as surface

#as surface

# create the pseudoatom at the origin

cmd.pseudoatom("pOrig", anchor, label="origin")

# these are special PyMOL variables that will hold # the coordinates of
# the atoms and the pseudoatom

stored.origCoord = []
stored.distCoord = []

# copy the coordinates into those special variables

cmd.iterate_state(1, "pOrig", 'stored.origCoord.append((x,y,z))')
cmd.iterate_state(1, site, 'stored.distCoord.append((x,y,z))')

# extend origCoord to be the same length as the other

stored.origCoord *= len(stored.distCoord)

# calculate the distances

stored.newB = map(lambda x,y: diff_len(x,y), stored.distCoord, stored.origCoord)
#print(stored.newB)
# put them into the b-factor of the protein

cmd.alter( site, "b=stored.newB.pop(0)")

# color by rainbow_rev or any other
# palette listed in "help spectrum"

cmd.spectrum(expression="b", palette="rainbow", selection=site)
cmd.set("surface_color","-1",site) #color the surface of the binding site by corresponding atom colors
cmd.extend('colorByDistance', colorByDistance)
</source>

ColorByGroup

2018-11-12T09:54:36Z

Zhentg: /* The Code */

{{Infobox script-repo
|type = script
|author = [[User:zhentg|Zhenting Gao]]
}}
== Introduction ==
Color the objects by their groups.

== Usage ==
*Open PyMOL
*Load PDB files, create groups
*run this Python script inside PyMOL
*call the function
**colorByGroup

== Required Arguments ==
Text

== Optional Arguments ==
Text

== Examples ==
Text
== The Code ==
<source>
def colorByGroup(palette="rainbow"):
#Based on: https://pymolwiki.org/index.php/Get_Names#ARGUMENTS
#Author: Zhenting Gao (zhentgpicasa@gmail.com)
#Update: 11/12/2018
#Aim: Color the objects by their groups.
# returns the length of the distance between atom A and atom B
groups=cmd.get_names("group_objects")
numClusters=len(groups)
import math
for x in range(numClusters):
cmd.alter(groups[x],"b="+str(x))
cmd.spectrum( "b", palette, " ".join(groups))
#https://pymolwiki.org/index.php/Spectrum
cmd.extend('colorByGroup', colorByGroup)

</source>

ColorByGroup

2018-11-12T06:47:06Z

Zhentg: /* The Code */

{{Infobox script-repo
|type = script
|author = [[User:zhentg|Zhenting Gao]]
}}
== Introduction ==
Color the objects by their groups.

== Usage ==
*Open PyMOL
*Load PDB files, create groups
*run this Python script inside PyMOL
*call the function
**colorByGroup

== Required Arguments ==
Text

== Optional Arguments ==
Text

== Examples ==
Text
== The Code ==
<source>
def colorByGroup():
#Based on: https://pymolwiki.org/index.php/Get_Names#ARGUMENTS
#Author: Zhenting Gao (zhentgpicasa@gmail.com)
#Update: 11/12/2018
#Aim: Color the objects by their groups.
# returns the length of the distance between atom A and atom B
groups=cmd.get_names("group_objects")
numClusters=len(groups)
import math
for x in range(numClusters):
colorCode=x*int(math.floor(255/numClusters))
util.cba(colorCode,groups[x])
print(colorCode,groups[x])
cmd.extend('colorByGroup', colorByGroup)
</source>

ColorByGroup

2018-11-12T06:42:01Z

Zhentg: Created page with "{{Infobox script-repo |type = script |author = Zhenting Gao }} == Introduction == Color the objects by their groups. == Usage == *Open PyMOL *Load PD..."

{{Infobox script-repo
|type = script
|author = [[User:zhentg|Zhenting Gao]]
}}
== Introduction ==
Color the objects by their groups.

== Usage ==
*Open PyMOL
*Load PDB files, create groups
*run this Python script inside PyMOL
*call the function
**colorByGroup

== Required Arguments ==
Text

== Optional Arguments ==
Text

== Examples ==
Text
== The Code ==
<source>
def colorByGroup():
#Based on: https://pymolwiki.org/index.php/Get_Names#ARGUMENTS
#Author: Zhenting Gao (zhentgpicasa@gmail.com)
#Update: 11/12/2018
#Aim: Color atoms of a binding site based on their distance from a point
# returns the length of the distance between atom A and atom B
groups=cmd.get_names("group_objects")
numClusters=len(groups)
import math
for x in range(numClusters):
colorCode=x*int(math.floor(255/numClusters))
util.cba(colorCode,groups[x])
print(colorCode,groups[x])
cmd.extend('colorByGroup', colorByGroup)
</source>

ColorByDistance

2018-11-12T06:35:31Z

Zhentg:

{{Infobox script-repo
|type = script
|author = [[User:zhentg|Zhenting Gao]]
}}
== Introduction ==
Color the surface of a binding site by the distance to a specific atom.

== Usage ==
*Open PyMOL
*Load PDB files
*run this Python script inside PyMOL
*call the function
**colorByDistance anchor, bindingSite

== Required Arguments ==
Text

== Optional Arguments ==
Text

== Examples ==
Text
== The Code ==
<source>
def colorByDistance(anchor='anchor', site='site'):
#Based on: https://pymolwiki.org/index.php/Spectrum#Intermediate
#Author: Zhenting Gao (zhentgpicasa@gmail.com)
#Update: 11/12/2018
#Aim: Color atoms of a binding site based on their distance from a point
# returns the length of the distance between atom A and atom B

diff_len = lambda x,y : math.sqrt((x[0]-y[0])*(x[0]-y[0]) + (x[1]-y[1])*(x[1]-y[1]) + (x[2]-y[2])*(x[2]-y[2]))

# fetch site from the PDB

#fetch site, async=0

# show it as surface

#as surface

# create the pseudoatom at the origin

#pseudoatom pOrig, pos=(0,0,0), label=origin

# these are special PyMOL variables that will hold # the coordinates of
# the atoms and the pseudoatom

stored.origCoord = []
stored.distCoord = []

# copy the coordinates into those special variables

cmd.iterate_state(1, anchor, 'stored.origCoord.append((x,y,z))')
cmd.iterate_state(1, site, 'stored.distCoord.append((x,y,z))')

# extend origCoord to be the same length as the other

stored.origCoord *= len(stored.distCoord)

# calculate the distances

stored.newB = map(lambda x,y: diff_len(x,y), stored.distCoord, stored.origCoord)
print(stored.newB)
# put them into the b-factor of the protein

cmd.alter( "site", "b=stored.newB.pop(0)")

# color by rainbow_rev or any other
# palette listed in "help spectrum"

cmd.spectrum(expression="b", palette="rainbow", selection=site)
cmd.set("surface_color","-1",site) #color the surface of the binding site by corresponding atom colors
cmd.extend('colorByDistance', colorByDistance)

</source>

ColorByDistance

2018-11-12T06:33:08Z

Zhentg:

{{Infobox script-repo
|type = script
|author = [[User:zhentg|Zhenting Gao]]
}}

==Code==
<source>
def colorByDistance(anchor='anchor', site='site'):
#Based on: https://pymolwiki.org/index.php/Spectrum#Intermediate
#Author: Zhenting Gao (zhentgpicasa@gmail.com)
#Update: 11/12/2018
#Aim: Color atoms of a binding site based on their distance from a point
# returns the length of the distance between atom A and atom B

diff_len = lambda x,y : math.sqrt((x[0]-y[0])*(x[0]-y[0]) + (x[1]-y[1])*(x[1]-y[1]) + (x[2]-y[2])*(x[2]-y[2]))

# fetch site from the PDB

#fetch site, async=0

# show it as surface

#as surface

# create the pseudoatom at the origin

#pseudoatom pOrig, pos=(0,0,0), label=origin

# these are special PyMOL variables that will hold # the coordinates of
# the atoms and the pseudoatom

stored.origCoord = []
stored.distCoord = []

# copy the coordinates into those special variables

cmd.iterate_state(1, anchor, 'stored.origCoord.append((x,y,z))')
cmd.iterate_state(1, site, 'stored.distCoord.append((x,y,z))')

# extend origCoord to be the same length as the other

stored.origCoord *= len(stored.distCoord)

# calculate the distances

stored.newB = map(lambda x,y: diff_len(x,y), stored.distCoord, stored.origCoord)
print(stored.newB)
# put them into the b-factor of the protein

cmd.alter( "site", "b=stored.newB.pop(0)")

# color by rainbow_rev or any other
# palette listed in "help spectrum"

cmd.spectrum(expression="b", palette="rainbow", selection=site)
cmd.set("surface_color","-1",site) #color the surface of the binding site by corresponding atom colors
cmd.extend('colorByDistance', colorByDistance)

</source>

RmsdByResidue

2018-11-12T06:32:22Z

Zhentg: /* Introduction */

{{Infobox script-repo
|type = script
|author = [[User:zhentg|Zhenting Gao]]
}}

== Introduction ==
RMSD between two structures of the same protein. The concept is similar as RMSF between two structures.
*programming
**Platform
***PyMOL
****rms_cur
**Feature
***Output
****RMSD of all atoms of each residues pairs
****Least RMSD of all atoms of each residues pairs
*****symmetry of Phe, Tyr, His, Asp, Glu, Gln, Asn, Arg, Leu and Val needs to be considered
******switch the atom name and then calculate the RMSD again
******Selected least RMSD of a residue pair for report
*******RMSD of backbone atoms of each residues pairs
*******RMSD of C alpha atoms of each residues pairs
*****With defined residues pairs
******Residue pair can be limited to within binding site
**Workflow
***Read reference and target pdb files
****two structures should be superposed before using this function
*Note
**Python
**PyMOL
***Clean attributes
****otherwise rms_cur will fail
***How to get residue name?
****residue name, residue index and etc. can only be read from an atom
*Reference
**https://sourceforge.net/p/pymol/mailman/message/32710889/

== Usage ==
*Open PyMOL
*Load PDB files
*run this Python script inside PyMOL
*call the function
**rmsdByRes pdb1, resSelection, pdb2

== Required Arguments ==
Text

== Optional Arguments ==
Text

== Examples ==
Text
== The Code ==
<source lang="python">
#load library
from pymol import cmd, stored
import sys
import os

#function to judge if a file exists
def is_non_zero_file(fpath):
return os.path.isfile(fpath) and os.path.getsize(fpath) > 0

#function to switch atom names within a residue
def switchName(residueSelection, atomName1, atomName2):
"""
switch the name of two atoms
"""
cmd.alter(residueSelection+" and name "+atomName1, 'name="gzt"')
cmd.alter(residueSelection+" and name "+atomName2, 'name="'+atomName1+'"')
cmd.alter(residueSelection+" and name gzt", 'name="'+atomName2+'"')

#function to change atom names of some residues with symetric sidechain
def flipAtomName(targetResidueSelection):
"""
switch the atom names of specific residues
"""
# Create flipped residue
cmd.create("flippedRes",targetResidueSelection+" and not alt B")
targetResidueCa=cmd.get_model("flippedRes and name CA")
for g in targetResidueCa.atom:
# print g.resn
if g.resn=='ARG':
switchName("flippedRes", "NH1", "NH2")
elif g.resn=='HIS':
switchName("flippedRes", "ND1", "CD2")
switchName("flippedRes", "CE1", "NE2")
elif g.resn=='ASP':
switchName("flippedRes", "OD1", "OD2")
elif g.resn=='PHE':
switchName("flippedRes", "CD1", "CD2")
switchName("flippedRes", "CE1", "CE2")
elif g.resn=='GLN':
switchName("flippedRes", "OE1", "NE2")
elif g.resn=='GLU':
switchName("flippedRes", "OE1", "OE2")
elif g.resn=='LEU':
switchName("flippedRes", "CD1", "CD2")
elif g.resn=='ASN':
switchName("flippedRes", "OD1", "ND2")
elif g.resn=='TYR':
switchName("flippedRes", "CD1", "CD2")
switchName("flippedRes", "CE1", "CE2")
elif g.resn=='VAL':
switchName("flippedRes", "CG1", "CG2")
cmd.sort()
# cmd.label("flippedRes","name")
return "flippedRes"

#main function
def rmsdByRes(referenceProteinChain,sel, targetProteinChain):
"""
Update
Zhenting Gao on 7/28/2016

USAGE

rmsf referenceProteinChain, targetProteinChain, selection [,byres=0], [reference_state=1]

Calculate the RMSD for each residue pairs from two chains of the same protein from two crystal structures.

Workflow
Read reference and target pdb files
Align two structures
sel target, proA and chain A
#define target protein chain
sel refrence, proB and chain A
#define reference protein chain
align target, reference
#automatical alignment
Clean attributes
otherwise rms_cur will fail

"""
# Create temporary objects, exclude alternative conformation B
cmd.create("ref_gzt", referenceProteinChain+" and polymer and not alt B")
cmd.alter("ref_gzt", "chain='A'")
cmd.alter("ref_gzt", "segi=''")
cmd.create("target_gzt", targetProteinChain+" and polymer and not alt B")
cmd.alter("target_gzt", "chain='A'")
cmd.alter("target_gzt", "segi=''")
# cmd.align("target_gzt","ref_gzt",object="align")
# parameters
outputText=""
res2Check=['HIS','ASP','ARG','PHE','GLN','GLU','LEU','ASN','TYR','VAL']

# select alpha carbon of selected residues in reference structure
calpha=cmd.get_model(sel+" and name CA and not alt B")

for g in calpha.atom:
# print g.resi+g.resn
if cmd.count_atoms("ref_gzt and polymer and resi "+g.resi)==cmd.count_atoms("target_gzt and polymer and resi "+g.resi):
rmsdRes=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi,"target_gzt and polymer and resi "+g.resi)
rmsdResCa=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi+" and name ca","target_gzt and polymer and resi "+g.resi+" and name ca")
rmsdResBackbone=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi+" and name ca+n+c+o","target_gzt and polymer and resi "+g.resi+" and name ca+n+c+o")
# calculate minimum rmsd
rmsdResMin=rmsdRes
if g.resn in res2Check:
flippedRes=flipAtomName("target_gzt and polymer and resi "+g.resi)
rmsdFlippedRes=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi,flippedRes)
if rmsdFlippedRes<rmsdRes:
rmsdResMin=rmsdFlippedRes

# print cmd.count_atoms("ref_gzt and polymer and resi "+g.resi),cmd.count_atoms("target_gzt and polymer and resi "+g.resi)
outputText+="%s,%s,%s,%.3f,%.3f,%.3f,%.3f\n" % (targetProteinChain,g.resn,g.resi,rmsdRes,rmsdResCa,rmsdResBackbone,rmsdResMin)

print outputText
# Destroy temporary objects
cmd.delete("ref_gzt target_gzt align res_gzt "+flippedRes)

# Save data into csv
outputFile='rmsdByRes_'+sel+'.csv'
f=open(outputFile,'a')
if not is_non_zero_file(outputFile):
f.write("targe,residueName,residueId,allAtomRMSD,rmsdResCa,rmsdResBackbone,allAtomRMSDMin\n")
f.write(outputText)
f.close()
print "Results saved in "+outputFile

cmd.extend("rmsdByRes",rmsdByRes)

</source>

[[Category:Script Library]]
[[Category:Structural_Biology_Scripts]]
[[Category:Pymol-script-repo]]

ColorByDistance

2018-11-12T06:23:43Z

Zhentg:

===Code===
<source>
def colorByDistance(anchor='anchor', site='site'):
#Based on: https://pymolwiki.org/index.php/Spectrum#Intermediate
#Author: Zhenting Gao (zhentgpicasa@gmail.com)
#Update: 11/12/2018
#Aim: Color atoms of a binding site based on their distance from a point
# returns the length of the distance between atom A and atom B

diff_len = lambda x,y : math.sqrt((x[0]-y[0])*(x[0]-y[0]) + (x[1]-y[1])*(x[1]-y[1]) + (x[2]-y[2])*(x[2]-y[2]))

# fetch site from the PDB

#fetch site, async=0

# show it as surface

#as surface

# create the pseudoatom at the origin

#pseudoatom pOrig, pos=(0,0,0), label=origin

# these are special PyMOL variables that will hold # the coordinates of
# the atoms and the pseudoatom

stored.origCoord = []
stored.distCoord = []

# copy the coordinates into those special variables

cmd.iterate_state(1, anchor, 'stored.origCoord.append((x,y,z))')
cmd.iterate_state(1, site, 'stored.distCoord.append((x,y,z))')

# extend origCoord to be the same length as the other

stored.origCoord *= len(stored.distCoord)

# calculate the distances

stored.newB = map(lambda x,y: diff_len(x,y), stored.distCoord, stored.origCoord)
print(stored.newB)
# put them into the b-factor of the protein

cmd.alter( "site", "b=stored.newB.pop(0)")

# color by rainbow_rev or any other
# palette listed in "help spectrum"

cmd.spectrum(expression="b", palette="rainbow", selection=site)
cmd.set("surface_color","-1",site) #color the surface of the binding site by corresponding atom colors
cmd.extend('colorByDistance', colorByDistance)

</source>

ColorByDistance

2018-11-12T06:23:07Z

Zhentg: Created page with "===Code=== def colorByDistance(anchor='anchor', site='site'): #Based on: https://pymolwiki.org/index.php/Spectrum#Intermediate #Author: Zhenting Gao (zhentgpicasa@gmail.com)..."

===Code===
def colorByDistance(anchor='anchor', site='site'):
#Based on: https://pymolwiki.org/index.php/Spectrum#Intermediate
#Author: Zhenting Gao (zhentgpicasa@gmail.com)
#Update: 11/12/2018
#Aim: Color atoms of a binding site based on their distance from a point
# returns the length of the distance between atom A and atom B

diff_len = lambda x,y : math.sqrt((x[0]-y[0])*(x[0]-y[0]) + (x[1]-y[1])*(x[1]-y[1]) + (x[2]-y[2])*(x[2]-y[2]))

# fetch site from the PDB

#fetch site, async=0

# show it as surface

#as surface

# create the pseudoatom at the origin

#pseudoatom pOrig, pos=(0,0,0), label=origin

# these are special PyMOL variables that will hold # the coordinates of
# the atoms and the pseudoatom

stored.origCoord = []
stored.distCoord = []

# copy the coordinates into those special variables

cmd.iterate_state(1, anchor, 'stored.origCoord.append((x,y,z))')
cmd.iterate_state(1, site, 'stored.distCoord.append((x,y,z))')

# extend origCoord to be the same length as the other

stored.origCoord *= len(stored.distCoord)

# calculate the distances

stored.newB = map(lambda x,y: diff_len(x,y), stored.distCoord, stored.origCoord)
print(stored.newB)
# put them into the b-factor of the protein

cmd.alter( "site", "b=stored.newB.pop(0)")

# color by rainbow_rev or any other
# palette listed in "help spectrum"

cmd.spectrum(expression="b", palette="rainbow", selection=site)
cmd.set("surface_color","-1",site) #color the surface of the binding site by corresponding atom colors
cmd.extend('colorByDistance', colorByDistance)

</source>

RmsdByResidue

2016-12-09T08:56:47Z

Zhentg: /* Introduction */

{{Infobox script-repo
|type = script
|author = [[User:zhentg|Zhenting Gao]]
}}

== Introduction ==
RMSD between two structures of the same protein
*programming
**Platform
***PyMOL
****rms_cur
**Feature
***Output
****RMSD of all atoms of each residues pairs
****Least RMSD of all atoms of each residues pairs
*****symmetry of Phe, Tyr, His, Asp, Glu, Gln, Asn, Arg, Leu and Val needs to be considered
******switch the atom name and then calculate the RMSD again
******Selected least RMSD of a residue pair for report
*******RMSD of backbone atoms of each residues pairs
*******RMSD of C alpha atoms of each residues pairs
*****With defined residues pairs
******Residue pair can be limited to within binding site
**Workflow
***Read reference and target pdb files
****two structures should be superposed before using this function
*Note
**Python
**PyMOL
***Clean attributes
****otherwise rms_cur will fail
***How to get residue name?
****residue name, residue index and etc. can only be read from an atom
*Reference
**https://sourceforge.net/p/pymol/mailman/message/32710889/

== Usage ==
*Open PyMOL
*Load PDB files
*run this Python script inside PyMOL
*call the function
**rmsdByRes pdb1, resSelection, pdb2

== Required Arguments ==
Text

== Optional Arguments ==
Text

== Examples ==
Text
== The Code ==
<source lang="python">
#load library
from pymol import cmd, stored
import sys
import os

#function to judge if a file exists
def is_non_zero_file(fpath):
return os.path.isfile(fpath) and os.path.getsize(fpath) > 0

#function to switch atom names within a residue
def switchName(residueSelection, atomName1, atomName2):
"""
switch the name of two atoms
"""
cmd.alter(residueSelection+" and name "+atomName1, 'name="gzt"')
cmd.alter(residueSelection+" and name "+atomName2, 'name="'+atomName1+'"')
cmd.alter(residueSelection+" and name gzt", 'name="'+atomName2+'"')

#function to change atom names of some residues with symetric sidechain
def flipAtomName(targetResidueSelection):
"""
switch the atom names of specific residues
"""
# Create flipped residue
cmd.create("flippedRes",targetResidueSelection+" and not alt B")
targetResidueCa=cmd.get_model("flippedRes and name CA")
for g in targetResidueCa.atom:
# print g.resn
if g.resn=='ARG':
switchName("flippedRes", "NH1", "NH2")
elif g.resn=='HIS':
switchName("flippedRes", "ND1", "CD2")
switchName("flippedRes", "CE1", "NE2")
elif g.resn=='ASP':
switchName("flippedRes", "OD1", "OD2")
elif g.resn=='PHE':
switchName("flippedRes", "CD1", "CD2")
switchName("flippedRes", "CE1", "CE2")
elif g.resn=='GLN':
switchName("flippedRes", "OE1", "NE2")
elif g.resn=='GLU':
switchName("flippedRes", "OE1", "OE2")
elif g.resn=='LEU':
switchName("flippedRes", "CD1", "CD2")
elif g.resn=='ASN':
switchName("flippedRes", "OD1", "ND2")
elif g.resn=='TYR':
switchName("flippedRes", "CD1", "CD2")
switchName("flippedRes", "CE1", "CE2")
elif g.resn=='VAL':
switchName("flippedRes", "CG1", "CG2")
cmd.sort()
# cmd.label("flippedRes","name")
return "flippedRes"

#main function
def rmsdByRes(referenceProteinChain,sel, targetProteinChain):
"""
Update
Zhenting Gao on 7/28/2016

USAGE

rmsf referenceProteinChain, targetProteinChain, selection [,byres=0], [reference_state=1]

Calculate the RMSD for each residue pairs from two chains of the same protein from two crystal structures.

Workflow
Read reference and target pdb files
Align two structures
sel target, proA and chain A
#define target protein chain
sel refrence, proB and chain A
#define reference protein chain
align target, reference
#automatical alignment
Clean attributes
otherwise rms_cur will fail

"""
# Create temporary objects, exclude alternative conformation B
cmd.create("ref_gzt", referenceProteinChain+" and polymer and not alt B")
cmd.alter("ref_gzt", "chain='A'")
cmd.alter("ref_gzt", "segi=''")
cmd.create("target_gzt", targetProteinChain+" and polymer and not alt B")
cmd.alter("target_gzt", "chain='A'")
cmd.alter("target_gzt", "segi=''")
# cmd.align("target_gzt","ref_gzt",object="align")
# parameters
outputText=""
res2Check=['HIS','ASP','ARG','PHE','GLN','GLU','LEU','ASN','TYR','VAL']

# select alpha carbon of selected residues in reference structure
calpha=cmd.get_model(sel+" and name CA and not alt B")

for g in calpha.atom:
# print g.resi+g.resn
if cmd.count_atoms("ref_gzt and polymer and resi "+g.resi)==cmd.count_atoms("target_gzt and polymer and resi "+g.resi):
rmsdRes=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi,"target_gzt and polymer and resi "+g.resi)
rmsdResCa=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi+" and name ca","target_gzt and polymer and resi "+g.resi+" and name ca")
rmsdResBackbone=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi+" and name ca+n+c+o","target_gzt and polymer and resi "+g.resi+" and name ca+n+c+o")
# calculate minimum rmsd
rmsdResMin=rmsdRes
if g.resn in res2Check:
flippedRes=flipAtomName("target_gzt and polymer and resi "+g.resi)
rmsdFlippedRes=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi,flippedRes)
if rmsdFlippedRes<rmsdRes:
rmsdResMin=rmsdFlippedRes

# print cmd.count_atoms("ref_gzt and polymer and resi "+g.resi),cmd.count_atoms("target_gzt and polymer and resi "+g.resi)
outputText+="%s,%s,%s,%.3f,%.3f,%.3f,%.3f\n" % (targetProteinChain,g.resn,g.resi,rmsdRes,rmsdResCa,rmsdResBackbone,rmsdResMin)

print outputText
# Destroy temporary objects
cmd.delete("ref_gzt target_gzt align res_gzt "+flippedRes)

# Save data into csv
outputFile='rmsdByRes_'+sel+'.csv'
f=open(outputFile,'a')
if not is_non_zero_file(outputFile):
f.write("targe,residueName,residueId,allAtomRMSD,rmsdResCa,rmsdResBackbone,allAtomRMSDMin\n")
f.write(outputText)
f.close()
print "Results saved in "+outputFile

cmd.extend("rmsdByRes",rmsdByRes)

</source>

[[Category:Script Library]]
[[Category:Structural_Biology_Scripts]]
[[Category:Pymol-script-repo]]

RmsdByResidue

2016-12-08T06:47:11Z

Zhentg:

{{Infobox script-repo
|type = script
|author = [[User:zhentg|Zhenting Gao]]
}}

== Introduction ==
RMSD between two structures of the same protein
*programming
**Platform
****PyMOL
*****rms_cur
**Feature
***Output
® RMSD of all atoms of each residues pairs
® Least RMSD of all atoms of each residues pairs
◊ symmetry of Phe, Tyr, His, Asp, Glu, Gln, Asn, Arg, Leu and Val needs to be considered
} switch the atom name and then calculate the RMSD again
◊ Selected least RMSD of a residue pair for report
® RMSD of backbone atoms of each residues pairs
® RMSD of C alpha atoms of each residues pairs
□ With defined residues pairs
® Residue pair can be limited to within binding site
§ Workflow
□ Read reference and target pdb files
® #two structures should be superposed before using this function
§ Note
□ Python
□ PyMOL
® Clean attributes
◊ otherwise rms_cur will fail
® How to get residue name?
◊ residue name, residue index and etc. can only be read from an atom

== Usage ==
*Open PyMOL
*Load PDB files
*run this Python script inside PyMOL
*call the function
**rmsdByRes pdb1, resSelection, pdb2

== Required Arguments ==
Text

== Optional Arguments ==
Text

== Examples ==
Text
== The Code ==
<source lang="python">
#load library
from pymol import cmd, stored
import sys
import os

#function to judge if a file exists
def is_non_zero_file(fpath):
return os.path.isfile(fpath) and os.path.getsize(fpath) > 0

#function to switch atom names within a residue
def switchName(residueSelection, atomName1, atomName2):
"""
switch the name of two atoms
"""
cmd.alter(residueSelection+" and name "+atomName1, 'name="gzt"')
cmd.alter(residueSelection+" and name "+atomName2, 'name="'+atomName1+'"')
cmd.alter(residueSelection+" and name gzt", 'name="'+atomName2+'"')

#function to change atom names of some residues with symetric sidechain
def flipAtomName(targetResidueSelection):
"""
switch the atom names of specific residues
"""
# Create flipped residue
cmd.create("flippedRes",targetResidueSelection+" and not alt B")
targetResidueCa=cmd.get_model("flippedRes and name CA")
for g in targetResidueCa.atom:
# print g.resn
if g.resn=='ARG':
switchName("flippedRes", "NH1", "NH2")
elif g.resn=='HIS':
switchName("flippedRes", "ND1", "CD2")
switchName("flippedRes", "CE1", "NE2")
elif g.resn=='ASP':
switchName("flippedRes", "OD1", "OD2")
elif g.resn=='PHE':
switchName("flippedRes", "CD1", "CD2")
switchName("flippedRes", "CE1", "CE2")
elif g.resn=='GLN':
switchName("flippedRes", "OE1", "NE2")
elif g.resn=='GLU':
switchName("flippedRes", "OE1", "OE2")
elif g.resn=='LEU':
switchName("flippedRes", "CD1", "CD2")
elif g.resn=='ASN':
switchName("flippedRes", "OD1", "ND2")
elif g.resn=='TYR':
switchName("flippedRes", "CD1", "CD2")
switchName("flippedRes", "CE1", "CE2")
elif g.resn=='VAL':
switchName("flippedRes", "CG1", "CG2")
cmd.sort()
# cmd.label("flippedRes","name")
return "flippedRes"

#main function
def rmsdByRes(referenceProteinChain,sel, targetProteinChain):
"""
Update
Zhenting Gao on 7/28/2016

USAGE

rmsf referenceProteinChain, targetProteinChain, selection [,byres=0], [reference_state=1]

Calculate the RMSD for each residue pairs from two chains of the same protein from two crystal structures.

Workflow
Read reference and target pdb files
Align two structures
sel target, proA and chain A
#define target protein chain
sel refrence, proB and chain A
#define reference protein chain
align target, reference
#automatical alignment
Clean attributes
otherwise rms_cur will fail

"""
# Create temporary objects, exclude alternative conformation B
cmd.create("ref_gzt", referenceProteinChain+" and polymer and not alt B")
cmd.alter("ref_gzt", "chain='A'")
cmd.alter("ref_gzt", "segi=''")
cmd.create("target_gzt", targetProteinChain+" and polymer and not alt B")
cmd.alter("target_gzt", "chain='A'")
cmd.alter("target_gzt", "segi=''")
# cmd.align("target_gzt","ref_gzt",object="align")
# parameters
outputText=""
res2Check=['HIS','ASP','ARG','PHE','GLN','GLU','LEU','ASN','TYR','VAL']

# select alpha carbon of selected residues in reference structure
calpha=cmd.get_model(sel+" and name CA and not alt B")

for g in calpha.atom:
# print g.resi+g.resn
if cmd.count_atoms("ref_gzt and polymer and resi "+g.resi)==cmd.count_atoms("target_gzt and polymer and resi "+g.resi):
rmsdRes=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi,"target_gzt and polymer and resi "+g.resi)
rmsdResCa=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi+" and name ca","target_gzt and polymer and resi "+g.resi+" and name ca")
rmsdResBackbone=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi+" and name ca+n+c+o","target_gzt and polymer and resi "+g.resi+" and name ca+n+c+o")
# calculate minimum rmsd
rmsdResMin=rmsdRes
if g.resn in res2Check:
flippedRes=flipAtomName("target_gzt and polymer and resi "+g.resi)
rmsdFlippedRes=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi,flippedRes)
if rmsdFlippedRes<rmsdRes:
rmsdResMin=rmsdFlippedRes

# print cmd.count_atoms("ref_gzt and polymer and resi "+g.resi),cmd.count_atoms("target_gzt and polymer and resi "+g.resi)
outputText+="%s,%s,%s,%.3f,%.3f,%.3f,%.3f\n" % (targetProteinChain,g.resn,g.resi,rmsdRes,rmsdResCa,rmsdResBackbone,rmsdResMin)

print outputText
# Destroy temporary objects
cmd.delete("ref_gzt target_gzt align res_gzt "+flippedRes)

# Save data into csv
outputFile='rmsdByRes_'+sel+'.csv'
f=open(outputFile,'a')
if not is_non_zero_file(outputFile):
f.write("targe,residueName,residueId,allAtomRMSD,rmsdResCa,rmsdResBackbone,allAtomRMSDMin\n")
f.write(outputText)
f.close()
print "Results saved in "+outputFile

cmd.extend("rmsdByRes",rmsdByRes)

</source>

[[Category:Script Library]]
[[Category:Structural_Biology_Scripts]]
[[Category:Pymol-script-repo]]

RmsdByResidue

2016-12-08T06:46:20Z

Zhentg:

{{Infobox script-repo
|type = script
|author = [[User:zhentg|Zhenting Gao]]
}}

== Introduction ==
RMSD between two structures of the same protein
*programming
**Platform
****PyMOL
*****rms_cur
**Feature
***Output
® RMSD of all atoms of each residues pairs
® Least RMSD of all atoms of each residues pairs
◊ symmetry of Phe, Tyr, His, Asp, Glu, Gln, Asn, Arg, Leu and Val needs to be considered
} switch the atom name and then calculate the RMSD again
◊ Selected least RMSD of a residue pair for report
® RMSD of backbone atoms of each residues pairs
® RMSD of C alpha atoms of each residues pairs
□ With defined residues pairs
® Residue pair can be limited to within binding site
§ Workflow
□ Read reference and target pdb files
® #two structures should be superposed before using this function
§ Note
□ Python
□ PyMOL
® Clean attributes
◊ otherwise rms_cur will fail
® How to get residue name?
◊ residue name, residue index and etc. can only be read from an atom

== Usage ==
*Open PyMOL
*Load PDB files
*run this Python script inside PyMOL
*call the function
**rmsdByRes pdb1, resSelection, pdb2

== Required Arguments ==
* '''sel1''' = first selection
* '''sel2''' = second selection
* '''max_dist''' = max distance in Angstroms

== Optional Arguments ==
* '''output''' = accepts Screen/Print/None (default N)
* '''sidechain''' = limits (Y) results to sidechain atoms (default N)
* '''show''' = shows (Y) individual distances in pymol menu (default=N)

== Examples ==
'''example #1'''
<source lang="python">
PyMOL>pairwise_dist 1efa and chain D, 1efa and chain B, 3, output=S, show=Y

1efa/D/DC/13/OP1 to 1efa/B/TYR/47/OH: 2.765
1efa/D/DC/13/OP2 to 1efa/B/LEU/6/N: 2.983
1efa/D/DC/13/OP2 to 1efa/B/LEU/6/CB: 2.928
1efa/D/DT/14/O4' to 1efa/B/ALA/57/CB: 2.827
1efa/D/DT/14/OP1 to 1efa/B/ASN/25/OD1: 2.858
1efa/D/DT/14/OP1 to 1efa/B/GLN/54/NE2: 2.996
1efa/D/DT/14/OP2 to 1efa/B/SER/21/OG: 2.517
1efa/D/DC/15/N4 to 1efa/B/GLN/18/NE2: 2.723
1efa/D/DA/16/N6 to 1efa/B/GLN/18/NE2: 2.931

Number of distances calculated: 9
</source>
== The Code ==
<source lang="python">
#load library
from pymol import cmd, stored
import sys
import os

#function to judge if a file exists
def is_non_zero_file(fpath):
return os.path.isfile(fpath) and os.path.getsize(fpath) > 0

#function to switch atom names within a residue
def switchName(residueSelection, atomName1, atomName2):
"""
switch the name of two atoms
"""
cmd.alter(residueSelection+" and name "+atomName1, 'name="gzt"')
cmd.alter(residueSelection+" and name "+atomName2, 'name="'+atomName1+'"')
cmd.alter(residueSelection+" and name gzt", 'name="'+atomName2+'"')

#function to change atom names of some residues with symetric sidechain
def flipAtomName(targetResidueSelection):
"""
switch the atom names of specific residues
"""
# Create flipped residue
cmd.create("flippedRes",targetResidueSelection+" and not alt B")
targetResidueCa=cmd.get_model("flippedRes and name CA")
for g in targetResidueCa.atom:
# print g.resn
if g.resn=='ARG':
switchName("flippedRes", "NH1", "NH2")
elif g.resn=='HIS':
switchName("flippedRes", "ND1", "CD2")
switchName("flippedRes", "CE1", "NE2")
elif g.resn=='ASP':
switchName("flippedRes", "OD1", "OD2")
elif g.resn=='PHE':
switchName("flippedRes", "CD1", "CD2")
switchName("flippedRes", "CE1", "CE2")
elif g.resn=='GLN':
switchName("flippedRes", "OE1", "NE2")
elif g.resn=='GLU':
switchName("flippedRes", "OE1", "OE2")
elif g.resn=='LEU':
switchName("flippedRes", "CD1", "CD2")
elif g.resn=='ASN':
switchName("flippedRes", "OD1", "ND2")
elif g.resn=='TYR':
switchName("flippedRes", "CD1", "CD2")
switchName("flippedRes", "CE1", "CE2")
elif g.resn=='VAL':
switchName("flippedRes", "CG1", "CG2")
cmd.sort()
# cmd.label("flippedRes","name")
return "flippedRes"

#main function
def rmsdByRes(referenceProteinChain,sel, targetProteinChain):
"""
Update
Zhenting Gao on 7/28/2016

USAGE

rmsf referenceProteinChain, targetProteinChain, selection [,byres=0], [reference_state=1]

Calculate the RMSD for each residue pairs from two chains of the same protein from two crystal structures.

Workflow
Read reference and target pdb files
Align two structures
sel target, proA and chain A
#define target protein chain
sel refrence, proB and chain A
#define reference protein chain
align target, reference
#automatical alignment
Clean attributes
otherwise rms_cur will fail

"""
# Create temporary objects, exclude alternative conformation B
cmd.create("ref_gzt", referenceProteinChain+" and polymer and not alt B")
cmd.alter("ref_gzt", "chain='A'")
cmd.alter("ref_gzt", "segi=''")
cmd.create("target_gzt", targetProteinChain+" and polymer and not alt B")
cmd.alter("target_gzt", "chain='A'")
cmd.alter("target_gzt", "segi=''")
# cmd.align("target_gzt","ref_gzt",object="align")
# parameters
outputText=""
res2Check=['HIS','ASP','ARG','PHE','GLN','GLU','LEU','ASN','TYR','VAL']

# select alpha carbon of selected residues in reference structure
calpha=cmd.get_model(sel+" and name CA and not alt B")

for g in calpha.atom:
# print g.resi+g.resn
if cmd.count_atoms("ref_gzt and polymer and resi "+g.resi)==cmd.count_atoms("target_gzt and polymer and resi "+g.resi):
rmsdRes=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi,"target_gzt and polymer and resi "+g.resi)
rmsdResCa=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi+" and name ca","target_gzt and polymer and resi "+g.resi+" and name ca")
rmsdResBackbone=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi+" and name ca+n+c+o","target_gzt and polymer and resi "+g.resi+" and name ca+n+c+o")
# calculate minimum rmsd
rmsdResMin=rmsdRes
if g.resn in res2Check:
flippedRes=flipAtomName("target_gzt and polymer and resi "+g.resi)
rmsdFlippedRes=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi,flippedRes)
if rmsdFlippedRes<rmsdRes:
rmsdResMin=rmsdFlippedRes

# print cmd.count_atoms("ref_gzt and polymer and resi "+g.resi),cmd.count_atoms("target_gzt and polymer and resi "+g.resi)
outputText+="%s,%s,%s,%.3f,%.3f,%.3f,%.3f\n" % (targetProteinChain,g.resn,g.resi,rmsdRes,rmsdResCa,rmsdResBackbone,rmsdResMin)

print outputText
# Destroy temporary objects
cmd.delete("ref_gzt target_gzt align res_gzt "+flippedRes)

# Save data into csv
outputFile='rmsdByRes_'+sel+'.csv'
f=open(outputFile,'a')
if not is_non_zero_file(outputFile):
f.write("targe,residueName,residueId,allAtomRMSD,rmsdResCa,rmsdResBackbone,allAtomRMSDMin\n")
f.write(outputText)
f.close()
print "Results saved in "+outputFile

cmd.extend("rmsdByRes",rmsdByRes)

</source>

[[Category:Script Library]]
[[Category:Structural_Biology_Scripts]]
[[Category:Pymol-script-repo]]

User:Zhentg

2016-12-08T06:36:52Z

Zhentg: Created page with "Zhenting Gao [http://www.linkedin.com/in/zhentg >>LinkedIn]"

Zhenting Gao
[http://www.linkedin.com/in/zhentg >>LinkedIn]

RmsdByResidue

2016-08-10T08:08:35Z

Zhentg:

<source lang="python">
#load library
from pymol import cmd, stored
import sys
import os

#function to judge if a file exists
def is_non_zero_file(fpath):
return os.path.isfile(fpath) and os.path.getsize(fpath) > 0

#function to switch atom names within a residue
def switchName(residueSelection, atomName1, atomName2):
"""
switch the name of two atoms
"""
cmd.alter(residueSelection+" and name "+atomName1, 'name="gzt"')
cmd.alter(residueSelection+" and name "+atomName2, 'name="'+atomName1+'"')
cmd.alter(residueSelection+" and name gzt", 'name="'+atomName2+'"')

#function to change atom names of some residues with symetric sidechain
def flipAtomName(targetResidueSelection):
"""
switch the atom names of specific residues
"""
# Create flipped residue
cmd.create("flippedRes",targetResidueSelection+" and not alt B")
targetResidueCa=cmd.get_model("flippedRes and name CA")
for g in targetResidueCa.atom:
# print g.resn
if g.resn=='ARG':
switchName("flippedRes", "NH1", "NH2")
elif g.resn=='HIS':
switchName("flippedRes", "ND1", "CD2")
switchName("flippedRes", "CE1", "NE2")
elif g.resn=='ASP':
switchName("flippedRes", "OD1", "OD2")
elif g.resn=='PHE':
switchName("flippedRes", "CD1", "CD2")
switchName("flippedRes", "CE1", "CE2")
elif g.resn=='GLN':
switchName("flippedRes", "OE1", "NE2")
elif g.resn=='GLU':
switchName("flippedRes", "OE1", "OE2")
elif g.resn=='LEU':
switchName("flippedRes", "CD1", "CD2")
elif g.resn=='ASN':
switchName("flippedRes", "OD1", "ND2")
elif g.resn=='TYR':
switchName("flippedRes", "CD1", "CD2")
switchName("flippedRes", "CE1", "CE2")
elif g.resn=='VAL':
switchName("flippedRes", "CG1", "CG2")
cmd.sort()
# cmd.label("flippedRes","name")
return "flippedRes"

#main function
def rmsdByRes(referenceProteinChain,sel, targetProteinChain):
"""
Update
Zhenting Gao on 7/28/2016

USAGE

rmsf referenceProteinChain, targetProteinChain, selection [,byres=0], [reference_state=1]

Calculate the RMSD for each residue pairs from two chains of the same protein from two crystal structures.

Workflow
Read reference and target pdb files
Align two structures
sel target, proA and chain A
#define target protein chain
sel refrence, proB and chain A
#define reference protein chain
align target, reference
#automatical alignment
Clean attributes
otherwise rms_cur will fail

"""
# Create temporary objects, exclude alternative conformation B
cmd.create("ref_gzt", referenceProteinChain+" and polymer and not alt B")
cmd.alter("ref_gzt", "chain='A'")
cmd.alter("ref_gzt", "segi=''")
cmd.create("target_gzt", targetProteinChain+" and polymer and not alt B")
cmd.alter("target_gzt", "chain='A'")
cmd.alter("target_gzt", "segi=''")
# cmd.align("target_gzt","ref_gzt",object="align")
# parameters
outputText=""
res2Check=['HIS','ASP','ARG','PHE','GLN','GLU','LEU','ASN','TYR','VAL']

# select alpha carbon of selected residues in reference structure
calpha=cmd.get_model(sel+" and name CA and not alt B")

for g in calpha.atom:
# print g.resi+g.resn
if cmd.count_atoms("ref_gzt and polymer and resi "+g.resi)==cmd.count_atoms("target_gzt and polymer and resi "+g.resi):
rmsdRes=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi,"target_gzt and polymer and resi "+g.resi)
rmsdResCa=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi+" and name ca","target_gzt and polymer and resi "+g.resi+" and name ca")
rmsdResBackbone=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi+" and name ca+n+c+o","target_gzt and polymer and resi "+g.resi+" and name ca+n+c+o")
# calculate minimum rmsd
rmsdResMin=rmsdRes
if g.resn in res2Check:
flippedRes=flipAtomName("target_gzt and polymer and resi "+g.resi)
rmsdFlippedRes=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi,flippedRes)
if rmsdFlippedRes<rmsdRes:
rmsdResMin=rmsdFlippedRes

# print cmd.count_atoms("ref_gzt and polymer and resi "+g.resi),cmd.count_atoms("target_gzt and polymer and resi "+g.resi)
outputText+="%s,%s,%s,%.3f,%.3f,%.3f,%.3f\n" % (targetProteinChain,g.resn,g.resi,rmsdRes,rmsdResCa,rmsdResBackbone,rmsdResMin)

print outputText
# Destroy temporary objects
cmd.delete("ref_gzt target_gzt align res_gzt "+flippedRes)

# Save data into csv
outputFile='rmsdByRes_'+sel+'.csv'
f=open(outputFile,'a')
if not is_non_zero_file(outputFile):
f.write("targe,residueName,residueId,allAtomRMSD,rmsdResCa,rmsdResBackbone,allAtomRMSDMin\n")
f.write(outputText)
f.close()
print "Results saved in "+outputFile

cmd.extend("rmsdByRes",rmsdByRes)

</source>

[[Category:Script Library]]
[[Category:Structural_Biology_Scripts]]
[[Category:Pymol-script-repo]]

RmsdByResidue

2016-08-10T08:06:41Z

Zhentg: Created page with "<source lang="python"> #load library from pymol import cmd, stored import sys import os #function to judge if a file exists def is_non_zero_file(fpath): return os.path.is..."

BiologicalUnit

2011-09-26T05:06:30Z

Zhentg: /* Usage */

This script can be used to re-create biological units for proteins. (This was created as a workaround of PyMOL's semi-functioning [[Symexp]] command.) It's also a fun script to play with for learning about symmetry.

<gallery widths=250px heights=250px>
Image:before.png|single unit, before running this program
Image:after.png|after the expansion of the 60 units
Image:1rmv fiber.png|Example creating the biological unit for PDB 1RMV created from fiber diffraction.
</gallery>

= Usage =
<source lang="python">
load /path/to/some/pdbFile.pdb
symMat = readSymmetry("/path/to/some/pdbFile.pdb","pdbFile")
biologicalUnit("mates", "pdbFile", symMat)
</source>

= The Code =
<source lang="python">
#
# Jason Vertrees <Jason-dot-Vertrees-at-schrodinger_dot_com>, 2010.
#
import pymol
from pymol import cmd

def readSymmetry(inFile, verbose=None):
"""
This function will read "inFile" and glean the
symmetry operations, if any, from it.

PARAMS
inFile
(string) path to PDB file

verbose
(boolean) if verbose is not None, print more

RETURNS
matrix
Array of lists. One 16-element list per symmetry operation. Feed this matrix
into manualSymExp in order to make the other symmetry mates in the biological unit
"""
# a remark-350 lines has:
# REMARK 350 BIOMTn TAGn X Y Z Tx
REM, TAG, BIOMT, OPNO, X, Y, Z, TX = range(8)

thePDB = open(inFile, 'rb').readlines()

matrices = []
curTrans = -1

# The transformation is,
# output = U*input + Tx

for l in thePDB:
tokens = l.split()
if len(tokens)!=8:
continue
if tokens[REM]=="REMARK" and tokens[TAG]=="350" and tokens[BIOMT].startswith("BIOMT"):
if tokens[OPNO]!=curTrans:
# new transformation matrix
matrices.append([])

matrices[-1].append( map( lambda s: float(s), tokens[X:]))
curTrans = tokens[OPNO]

if verbose!=None:
print "Found %s symmetry operators in %s." % (len(matrices), inFile)
return matrices

def biologicalUnit(prefix, objSel, matrices ):
"""
Manually expands the object in "objSel" by the symmetry operations provided in "matrices" and
prefixes the new objects with "prefix".

PARAMS
prefix
(string) prefix name for new objects

objSel
(string) name of object to expand

matrices
(list of 16-element lists) array of matrices from readSymmetry

RETUNRS
None

SIDE EFFECTS
Creates N new obects each rotated and translated according to the symmetry operators, where N
equals len(matrices).
"""
for m in matrices:
n = cmd.get_unused_name(prefix)
cmd.create(n, objSel)
s1 = "%s + (x*%s + y*%s + z*%s)" % (m[0][3], m[0][0], m[0][1], m[0][2])
s2 = "%s + (x*%s + y*%s + z*%s)" % (m[1][3], m[1][0], m[1][1], m[1][2])
s3 = "%s + (x*%s + y*%s + z*%s)" % (m[2][3], m[2][0], m[2][1], m[2][2])
cmd.alter_state(1, n, "(x,y,z) = (%s, %s, %s)" % (s1, s2, s3) )
</source>

= Notes =
This is slow compared to [[Symexp]]; use the above for learning, playing and when [[Symexp]] doesn't work as advertised.

= See Also =
*[[BiologicalUnit/Quat]] (alternative implementation)
*[[Symexp]]
*[[SuperSym]]
*[http://pdbbeta.rcsb.org/pdb/static.do?p=education_discussion/Looking-at-Structures/bioassembly_tutorial.html PDB Tutorial Biol. Units]
*[http://en.wikipedia.org/wiki/Fiber_diffraction Wikipedia article]

[[Category:Script_Library]]
[[Category:Math_Scripts]]
[[Category:Structural_Biology_Scripts]]

Talk:BiologicalUnit

2010-05-20T06:07:47Z

Zhentg:

[[User:Zhentg|Dr. Zhenting Gao]] 06:07, 20 May 2010 (UTC)

The script does not work for me. Can anyone help me out?

My procedure is
#Open PyMOL
#run E:/readSymmetry.py
#symMat = readSymmetry("E:/pdb1k4c.pdb","pdbFile")
#biologicalUnit("mates", "pdbFile", symMat)
#Then the following error occurs
<pre>
Traceback (most recent call last):
File "d:\cadd\pymol/modules\pymol\parser.py", line 455, in parse
exec(layer.com2+"\n",self.pymol_names,self.pymol_names)
File "<string>", line 1, in ?
File "E:/readSymmetry.py", line 77, in biologicalUnit
n = cmd.get_unused_name(prefix)
AttributeError: 'module' object has no attribute 'get_unused_name'
</pre>

Talk:BiologicalUnit

2010-05-20T06:07:19Z

Zhentg: Created page with '~~~~ The script does not work for me. My procedure is #Open PyMOL #run E:/readSymmetry.py #symMat = readSymmetry("E:/pdb1k4c.pdb","pdbFile") #biologicalUnit("mates", "pdbFile",…'

[[User:Zhentg|Dr. Zhenting Gao]] 06:07, 20 May 2010 (UTC)

The script does not work for me.

My procedure is
#Open PyMOL
#run E:/readSymmetry.py
#symMat = readSymmetry("E:/pdb1k4c.pdb","pdbFile")
#biologicalUnit("mates", "pdbFile", symMat)
#Then the following error occurs
<pre>
Traceback (most recent call last):
File "d:\cadd\pymol/modules\pymol\parser.py", line 455, in parse
exec(layer.com2+"\n",self.pymol_names,self.pymol_names)
File "<string>", line 1, in ?
File "E:/readSymmetry.py", line 77, in biologicalUnit
n = cmd.get_unused_name(prefix)
AttributeError: 'module' object has no attribute 'get_unused_name'
</pre>

Cealign plugin

2009-06-08T07:43:28Z

Zhentg: /* Requirements */

== Introduction ==

'''Go directly to [[Cealign#Version_0.8-RBS|DOWNLOAD]]'''

This page is the home page of the open-source CEAlign PyMOL plugin. The CE algorithm is a fast and accurate protein structure alignment algorithm, pioneered by Drs. Shindyalov and Bourne (See
References). There are a few changes from the original CE publication (See Notes).

The source code is implemented in C with the rotations finally done by Numpy in Python. Because the computationally complex portion of the code is written in C, it's quick. That is, on my machines --- relatively fast 64-bit machines --- I can align two 400+ amino acid structures in about 0.300 s with the C++ implementation.

This plugs into PyMol very easily. See [[Cealign#The_Code|the code]] and [[Cealign#Examples|examples]] for installation and usage.

== Comparison to PyMol ==
'''Why should you use this?'''

PyMOL's structure alignment algorithm is fast and robust. However, its first step is to perform a sequence alignment of the two selections. Thus, proteins in the '''twilight zone''' or those having a low sequence identity, may not align well. Because CE is a structure-based alignment, this is not a problem. Consider the following example. The image at LEFT was the result of CE-aligning two proteins (1C0M chain B to 1BCO). The result is '''152''' aligned (alpha carbons) residues (not atoms) at '''4.96 Angstroms'''. The image on the RIGHT shows the results from PyMol's align command: an alignment of '''221 atoms''' (not residues) at an RMSD of '''15.7 Angstroms'''.

<gallery>
Image:cealign_ex1.png|Cealign's results (152 aligned; 4.96 Ang.)
Image:pymol_align.png|PyMol's results (763 atoms; 18.4 Ang. )
</gallery>

=== Fit vs. optAlign ===
====Take Home messages====
* [[fit]] and [[optAlign]] perform nearly equally as well
* if you need an algorithm with an appropriate reference, use [[optAlign]] (references at bottom of page).
* [[fit]] is faster -- if you're aligning many structures, use it over [[optAlign]]

====Discussion====
[[optAlign]] is a function within the [[Cealign]] package that performs the optimal superposition of two objects of equal length. [[optAlign]] follows the Kabsch algorithm which is a closed form, and provably optimal solution to the problem. [[fit]] on the other hand uses the Jacobi rotations to iteratively arrive at the solution of optimal superposition. The difference in error between [[optAilgn]] and [[fit]] seems to be a non-issue (see below) as they both arrive at equivalent solutions for the rotation matrix. The two algorithms are undertake different approaches to orthogonally diagonalizing the correlation matrix.

PyMOL's [[fit]] is fast and works well. If you have to use something with a known reference then check out the "optAlign" function from the qkabsch.py file that comes with this [[Calign]] package. If not, you can just use [[fit]] and avoid installing new software. :-)

optAlign is slower than fit. I just tested both on a sample NMR ensemble; and, while not an extensive validation of "fit" it shows that (1) fit is faster; and (2) fit gets the same exact RMSD as "optAlign" (when optAlign is told to use all atoms, not just CA). To make optAlign use all atoms and not just the alpha-carbon backbones, comment out (that is, put a "#" at the start of lines 183 and 184 in qkabsch.py, where it says "CUT HERE").

<source lang="python">
fetch 1nmr
split_states 1nmr
delete 1nmr

# compare fit and optAlign RMSDs
for x in cmd.get_names(): print cmd.fit("1nmr_0001", x)
for x in cmd.get_names(): optAlign(x, "1nmr_0001")
</source>
<source lang="bash">
# results from fit
0.0
4.50344991684
5.33588504791
5.78613853455
7.25597000122
6.67145586014
3.25131297112
3.36766290665
6.74802017212
5.1579709053
5.96959495544
6.68093347549
4.13217163086
5.51539039612
6.24266338348
6.03838825226
5.01363992691
5.33336305618
6.87617444992
7.797062397

#results from optAlign
RMSD=0.000000
RMSD=4.503450
RMSD=5.335886
RMSD=5.786138
RMSD=7.255970
RMSD=6.671456
RMSD=3.251313
RMSD=3.367663
RMSD=6.748021
RMSD=5.157971
RMSD=5.969595
RMSD=6.680934
RMSD=4.132172
RMSD=5.515390
RMSD=6.242664
RMSD=6.038388
RMSD=5.013640
RMSD=5.333363
RMSD=6.876174
RMSD=7.797062
</source>

== Examples ==
=== Usage ===
==== Syntax ====

CEAlign has the semantic, and syntactic formalism of
<source lang="python">
cealign MASTER, TARGET
</source>
where a post-condition of the algorithm is that the coordinates of the '''MASTER''' protein are unchanged. This allows for easier multi-protein alignments. For example,
<source lang="python">
cealign 1AUE, 1BZ4
cealign 1AUE, 1B68
cealign 1AUE, 1A7V
cealign 1AUE, 1CPR
</source>
will superimpose all the TARGETS onto the MASTER.

=====Examples=====
<source lang="python">
cealign 1cll and i. 42-55, 1ggz and c. A
cealign 1kao, 1ctq
cealign 1fao, 1eaz
</source>

=====Multiple Structure Alignments=====
Use the '''alignto''' command, now provided with cealign. Just type,
<source lang="python">
alignto PROT
</source>
to align all your proteins in PyMOL to the one called, '''PROT'''.

=== Results ===
See '''Changes''' for updates. But, overall, the results here are great.

<gallery>
Image:v7_1fao_1eaz.png|EASY: 1FAO vs. 1EAZ; 96 residues, 1.28 Ang
Image:v7_1cbs_1hmt.png|EASY: 1CBS vs. 1HMT; 128 residues, 2.01 Ang
Image:v7_1a15_1b50.png|MODERATE: 1A15 vs 1B50; 56 residues, 2.54 Ang.
Image:v7_1oan_1s6n.png|EASY: 1OAN vs. 1S6N (state 1); 96 residues aligned to 3.83 Ang. RMSD.
Image:v7_1rlw_1byn.png|HARD: 1RLW to 1BYN; 104 residues; 2.21 Ang.
Image:v7_1ten_3hhr.png|HARD: 1TEN vs. 3HHR; 80 residues, 2.91 Ang.
Image:v7_2sim_1nsb.png|HARD: 2SIM vs. 1NSB; 272 residues, 4.93 Ang.
Image:v7_1cew_1mol.png|HARD: 1CEW vs. 1MOL; 80 residues, 4.03 Ang.
</gallery>

== Installation ==

===Mac OS X (10.5)===
[[Image:Cealign mac os x.png|300px|thumb|center|CEAlign running on Mac OS X (10.5)]]
* Install PyMOL under fink.
* Install Numpy for fink:
<source lang="bash">
/sw/bin/fink install scipy-core-py25
</source>
* Download and install cealign (download instructions below)
<source lang="bash">
sudo /sw/bin/python setup.py install
</source>
* In PyMOL, run the two scripts needed for cealing: "cealign.py" and "qkabsch.py". These are located in the cealign directory you previously downloaded.
* Voila!
* Note that the above python version must match the same version that is used by PyMOL. If you are using the pre-compiled version of MacPyMOL, the above instructions won't work.

===Windows systems===
This is a quick and dirty method to get it working on Win32 right now, more details coming soon.
====Requirements====
* Latest PyMol, installed on your system
* Numpy for python 2.4 -- quick download of just what's needed: http://users.umassmed.edu/Shivender.Shandilya/pymol/numpy.zip[Note: the numpy.zip here are corrupt, users may download the latest 'Numpy for python 2.4' directly from SourceForge.net <[[User:Zhentg|Dr. Zhenting Gao]] 07:43, 8 June 2009 (UTC)>]
* Pre-compiled ccealign.pyd python module: http://users.umassmed.edu/Shivender.Shandilya/pymol/ccealign.zip
* Modified pymolrc: http://users.umassmed.edu/Shivender.Shandilya/pymol/pymolrc
* cealign.py and qkabsch.py from the Cealign-0.8-RBS package: download below

====Directions====
# Unzip the numpy.zip file, which will give you a folder named '''numpy'''
# Move this entire folder to: C:\Program Files\DeLano Scientific\PyMOL\modules\ (or the corresponding location on your system)
# Unzip ccealign.zip, which will give you a file called '''ccealign.pyd'''
# Move this pyd file to: C:\Program Files\DeLano Scientific\PyMOL\py24\DLLs\ (or the corresponding location on your system)
# Copy the downloaded '''pymolrc''' file to: C:\Program Files\DeLano Scientific\PyMOL\ (or the corresponding location on your system)
# Extract and copy the files cealign.py and qkabsch.py from the Cealign-0.8-RBS package to: C:\Program Files\DeLano Scientific\PyMOL\py24\Lib\ (or the corresponding location on your system)
# Run PyMol and load some molecules
# Run this command in Pymol: '''cealign molecule1, molecule2'''
# Enjoy!

===Gentoo Linux===
Add the science overlay via
layman -a sci
and emerge the cealign plugin
emerge pymol-plugins-cealign

===*nix systems===
====Requirements====
* C compiler
* Python 2.4+ with distutils
* Numpy
** for User-compiled PyMOL: <source lang="python">python setup.py install</source>
** for the precompiled version of PyMOL <source lang="python">python setup.py install --prefix "" --root /DIR_TO/pymol/ext/</source>

====Directions====
# uncompress the distribution file '''cealign-VERSION.tgz'''
# cd cealign-VERSION
# sudo python setup.py install # if you installed by PyMOL by hand
## python setup.py install --prefix "" --root /DIR/TO/pymol/ext/ # if you are using the precompiled binary download
# insert "run DIR_TO_CEALIGN/cealign.py" and "run DIR_TO_CEALIGN/qkabsch.py" into your '''.pymolrc''' file, or just run the two Python scripts by hand.
# load some molecules
# run, '''cealign molecule1, molecule2'''
# enjoy

=====Pre-compiled Hackish Install=====
For those people that prefer to use the pre-compiled version of PyMOL, here are the basics for your install. '''This is a poor method of installing Cealign. I suggest users compile and install their own PyMOL.''' The final goal is to get
# '''ccealign.so''' module into '''PYMOL/ext/lib/python2.4/site-packages'''
# numpy installed (get the numpy directory into (or linked into) '''PYMOL/ext/lib/python2.4/site-packages'''
# and be able to run cealign.py and qkabsch.py from PyMOL.
If you can do the above three steps, '''cealign''' should run from the pre-compiled PyMOL.

In more detail, on a completely fictitious machine --- that is, I created the following commands from a fake machine and I don't expect a copy/paste of this to work '''anywhere''', but the commands should be helpful enough to those who need it:
<source lang="python">
# NOTES:
# This is fake code: don't copy/paste it.
#
# PYMOL='dir to precompiled PyMOL install'
# CEALIGN='dir where you will unpack cealign'
# replace lib with lib64 for x86-64
# install numpy
apt-get install numpy

# link numpy to PyMOL
ln -s /usr/local/lib/python2.4/site-packages/numpy PYMOL/ext/lib/python2.4/site-packages

# download and install Cealign
wget http://www.pymolwiki.org/images/e/ed/Cealign-0.6.tar.bz2
tar -jxvf Cealign-0.6.tar.bz2
cd cealign-0.6
sudo python setup.py build
cp build/lib-XYZ-linux/ccealign.so PYMOL/ext/lib/python2.4/site-packages

# run pymol and try it out
pymol
run CEALIGN/cealign.py
run CEALIGN/qkabsch.py
fetch 1cew 1mol, async=0
cealign 1c, 1m
</source>

== The Code ==
Please unpack and read the documentation. All comments/questions should be directed to Jason Vertrees (javertre _at_ utmb ...dot... edu).

'''LATEST IS v0.8-RBS'''. (Dedicated to Bryan Sutton for allowing me to use his computer for testing.)

=== Version 0.8-RBS ===
* '''Download: [[Media:Cealign-0.8-RBS.tar.bz2|CE Align v0.8-RBS]] (bz2)'''
* '''Download: [[Media:Cealign-0.8-RBS.zip|CE Align v0.8-RBS]] (zip)'''

=== Beta Version 0.9 ===
Use at your own peril. Please report any problems or inconsistent alignments to this discussion page, or to me directly; my email address all over this page.

'''Improvements/Changes''':
* All C++
** So, faster
** comes with the dependencies built in
* No numpy

''' Download: [[Media:Cealign-0.9.zip|CE Align v0.9]] (zip)'''

== Coming Soon ==
* Windows binary
* Linux Binaries (32bit, x86-64)
* Better instructions for precompiled distributions
* Optimization

== Updates ==

=== 2008-03-25 ===
Pure C++ code released. See the beta version above.

=== 2007-04-14 ===
v0.8-RBS source updated. Found the bug that had been plaguing 32-bit machines. This should be the last release for a little while.

Also, I provide the option of aligning based solely upon RMSD or upon the better CE-Score. See the '''References''' for information on the '''CE Score'''.

== Troubleshooting ==

Post your problems/solutions here.

=== Unicode Issues in Python/Numpy ===
'''Problem''': Running/Installing cealign gives
<source lang="python">
Traceback (most recent call last):
File "/home/byron/software/pymol_1.00b17/pymol/modules/pymol/parser.py",
line 308, in parse
File "/home/byron/software/pymol_1.00b17/pymol/modules/pymol/parsing.py",
line 410, in run_file
File "qkabsch.py", line 86, in ?
import numpy
File "/usr/lib/python2.4/site-packages/numpy/__init__.py", line 36, in ?
import core
File "/usr/lib/python2.4/site-packages/numpy/core/__init__.py", line 5, in ?
import multiarray
ImportError: /home/byron/software/pymol/ext/lib/python2.4/site-packages/numpy/core/multiarray.so:
undefined symbol: _PyUnicodeUCS4_IsWhitespace
</source>
where the important line is
<source lang="python">
undefined symbol: _PyUnicodeUCS4_IsWhitespace
</source>

This problem indicates that your Numpy Unicode is using a different byte-size for unicode characters than is the Python distribution your PyMOL is running from. For example, this can happen if you use the pre-built PyMOL and some other pre-built Numpy package.

'''Solution''': Hand-install Numpy.

=== LinAlg Module Not Found ===
'''Problem''': Running CE Align gives the following error message:
<source lang="python">
run qkabsch.py
Traceback (most recent call last):
File "/usr/lib/python2.4/site-packages/pymol/parser.py", line 285, in parse
parsing.run_file(exp_path(args[nest][0]),pymol_names,pymol_names)
File "/usr/lib/python2.4/site-packages/pymol/parsing.py", line 407, in run_file
execfile(file,global_ns,local_ns)
File "qkabsch.py", line 86, in ?
import numpy
File "/usr/lib/python2.4/site-packages/numpy/__init__.py", line 40, in ?
import linalg
ImportError: No module named linalg
</source>

'''Solution''': You do not have the linear algebra module installed (or Python can't find it) on your machine. One workaround is to install [http://www.scipy.org/ Scientific Python]. (on debian/ubuntu this can be done by: sudo apt-get install python-scipy) Another is to reinstall the Numpy package from source, ensuring that you have the necessary requirements for the linear algebra module (linpack, lapack, fft, etc.).

=== CCEAlign & NumPy Modules Not Found ===
'''Problem''': Running CE Align gives the following error message:
<source lang="python">
PyMOL>run cealign.py
Traceback (most recent call last):
File "/home/local/warren/MacPyMOL060530/build/Deployment/MacPyMOL.app/pymol/modules/pymol/parser.py", line 297, in parse
File "/home/local/warren/MacPyMOL060530/build/Deployment/MacPyMOL.app/pymol/modules/pymol/parsing.py", line 408, in run_file
File "/usr/local/pymol/scripts/cealign-0.1/cealign.py", line 59, in ?
from ccealign import ccealign
ImportError: No module named ccealign
run qkabsch.py
Traceback (most recent call last):
File "/home/local/warren/MacPyMOL060530/build/Deployment/MacPyMOL.app/pymol/modules/pymol/parser.py", line 297, in parse
File "/home/local/warren/MacPyMOL060530/build/Deployment/MacPyMOL.app/pymol/modules/pymol/parsing.py", line 408, in run_file
File "qkabsch.py", line 86, in ?
import numpy
ImportError: No module named numpy
</source>

'''Solution''': This problem occurs under [http://www.apple.com/macosx Apple Mac OS X] if (a) the Apple's python executable on your machine (/usr/bin/python, currently version 2.3.5) is superseded by [http://fink.sourceforge.net/ Fink]'s python executable (/sw/bin/python, currently version 2.5) and (b) you are using [http://delsci.com/rel/099/#MacOSX precompiled versions of PyMOL] (MacPyMOL, PyMOLX11Hybrid or PyMOL for Mac OS X/X11). These executables ignore Fink's python and instead use Apple's - so, in order to run CE Align, one must install NumPy (as well as CE Align itself) using Apple's python. To do so, first download the [http://sourceforge.net/project/showfiles.php?group_id=1369&package_id=175103 Numpy source code archive] (currently version 1.0.1), unpack it, change directory to numpy-1.0.1 and specify the full path to Apple's python executable during installation: <tt>sudo /usr/bin/python setup.py install | tee install.log</tt>. Then, donwload the [http://www.pymolwiki.org/index.php/Cealign#The_Code CE Align source code archive] (currently version 0.2), unpack it, change directory to cealign-0.2 and finally install CE Align as follows: <tt>sudo /usr/bin/python setup.py install | tee install.log</tt>.
[[User:Lucajovine|Luca Jovine]] 05:11, 25 January 2007 (CST).

=== The Function SimpAlign() is not found ===
'''Problem''': Running CE Align gives the following error message:
<source lang="python">
PyMOL>cealign 1CLL,1GGZ
Traceback (most recent call last):
File "C:\Program Files (x86)\DeLano Scientific\PyMOL/modules\pymol\parser.py", line 203, in parse
result=apply(kw[nest][0],args[nest],kw_args[nest])
File "py24/Lib/cealign.py", line 177, in cealign
curScore = simpAlign( matA, matB, mol1, mol2, stored.mol1, stored.mol2, align=0, L=len(matA) )
NameError: global name 'simpAlign' is not defined
</source>
I am running PyMOL v. 0.99rc6 on Win XP Professional x64 edition version 2003 sp2 and have followed the windows install procedure as described above.

'''Answer''': This simply means that PyMOL couldn't find the simplAlign function. To let PyMOL know about this, you must run the following commands before running [[cealign]]:
<source lang="python">
run /your/path/to/cealign/qkabsch.py
run /your/path/to/cealign/cealign.py
</source>
but most people that use cealign would just put these two lines in their '''.pymolrc''' file.

=== Short Alignments Don't Work ===
If you are trying to align fewer than 16 residues then use [[align]], [[super]], or [[optAlign]]. CE uses a window size of 8; and to build a path of more than one window, you need 2*8=16 residues. I will insert some code to re-route small alignments to one of the aforementioned alignment algorithms.

=== It Worked A Second Ago! ===
[[Image:Rewind.png|thumb|right|Showing the rewind button to rewind to state 1.]]

If you were using cealign (or alignto) and now the commands don't work -- that is, they return an RMSD, but don't actually superimpose the objects, then you have a simple problem dealing with states. Most likely the cause of this oddness was (1) when you issued "cealign prot1, prot2" one of them was actually an ensemble of states or (2) you are trying to align to proteins with only one state, but are not looking at state one (because the last protein you were considering had more than one state and you quit editing that protein on a state that's not state 1). To fix this, use the rewind button to get the proteins back into state 1 & reissue the cealign/alignto command.

== References ==
Text taken from PubMed and formatted for the wiki. The first reference is the most important for this code.

# Shindyalov IN, Bourne PE. '''Protein structure alignment by incremental combinatorial extension (CE) of the optimal path.''' ''Protein Eng.'' 1998 Sep;11(9):739-47. PMID: 9796821 [PubMed - indexed for MEDLINE]
# Jia Y, Dewey TG, Shindyalov IN, Bourne PE. '''A new scoring function and associated statistical significance for structure alignment by CE.''' ''J Comput Biol.'' 2004;11(5):787-99. PMID: 15700402 [PubMed - indexed for MEDLINE]
# Pekurovsky D, Shindyalov IN, Bourne PE. '''A case study of high-throughput biological data processing on parallel platforms.''' ''Bioinformatics.'' 2004 Aug 12;20(12):1940-7. Epub 2004 Mar 25. PMID: 15044237 [PubMed - indexed for MEDLINE]
# Shindyalov IN, Bourne PE. '''An alternative view of protein fold space.''' ''Proteins.'' 2000 Feb 15;38(3):247-60. PMID: 10713986 [PubMed - indexed for MEDLINE]

== License ==
The CEAlign and all its subprograms that I wrote, are released under the open source Free BSD License (BSDL).

[[Category:Script_Library]]
[[Category:Structure_Alignment|Cealign]]

Cealign

2009-06-08T07:43:28Z

Zhentg: /* Requirements */

Cealign plugin

2009-06-08T07:41:45Z

Zhentg: /* Requirements */

== Introduction ==

'''Go directly to [[Cealign#Version_0.8-RBS|DOWNLOAD]]'''

This page is the home page of the open-source CEAlign PyMOL plugin. The CE algorithm is a fast and accurate protein structure alignment algorithm, pioneered by Drs. Shindyalov and Bourne (See
References). There are a few changes from the original CE publication (See Notes).

The source code is implemented in C with the rotations finally done by Numpy in Python. Because the computationally complex portion of the code is written in C, it's quick. That is, on my machines --- relatively fast 64-bit machines --- I can align two 400+ amino acid structures in about 0.300 s with the C++ implementation.

This plugs into PyMol very easily. See [[Cealign#The_Code|the code]] and [[Cealign#Examples|examples]] for installation and usage.

== Comparison to PyMol ==
'''Why should you use this?'''

PyMOL's structure alignment algorithm is fast and robust. However, its first step is to perform a sequence alignment of the two selections. Thus, proteins in the '''twilight zone''' or those having a low sequence identity, may not align well. Because CE is a structure-based alignment, this is not a problem. Consider the following example. The image at LEFT was the result of CE-aligning two proteins (1C0M chain B to 1BCO). The result is '''152''' aligned (alpha carbons) residues (not atoms) at '''4.96 Angstroms'''. The image on the RIGHT shows the results from PyMol's align command: an alignment of '''221 atoms''' (not residues) at an RMSD of '''15.7 Angstroms'''.

<gallery>
Image:cealign_ex1.png|Cealign's results (152 aligned; 4.96 Ang.)
Image:pymol_align.png|PyMol's results (763 atoms; 18.4 Ang. )
</gallery>

=== Fit vs. optAlign ===
====Take Home messages====
* [[fit]] and [[optAlign]] perform nearly equally as well
* if you need an algorithm with an appropriate reference, use [[optAlign]] (references at bottom of page).
* [[fit]] is faster -- if you're aligning many structures, use it over [[optAlign]]

====Discussion====
[[optAlign]] is a function within the [[Cealign]] package that performs the optimal superposition of two objects of equal length. [[optAlign]] follows the Kabsch algorithm which is a closed form, and provably optimal solution to the problem. [[fit]] on the other hand uses the Jacobi rotations to iteratively arrive at the solution of optimal superposition. The difference in error between [[optAilgn]] and [[fit]] seems to be a non-issue (see below) as they both arrive at equivalent solutions for the rotation matrix. The two algorithms are undertake different approaches to orthogonally diagonalizing the correlation matrix.

PyMOL's [[fit]] is fast and works well. If you have to use something with a known reference then check out the "optAlign" function from the qkabsch.py file that comes with this [[Calign]] package. If not, you can just use [[fit]] and avoid installing new software. :-)

optAlign is slower than fit. I just tested both on a sample NMR ensemble; and, while not an extensive validation of "fit" it shows that (1) fit is faster; and (2) fit gets the same exact RMSD as "optAlign" (when optAlign is told to use all atoms, not just CA). To make optAlign use all atoms and not just the alpha-carbon backbones, comment out (that is, put a "#" at the start of lines 183 and 184 in qkabsch.py, where it says "CUT HERE").

<source lang="python">
fetch 1nmr
split_states 1nmr
delete 1nmr

# compare fit and optAlign RMSDs
for x in cmd.get_names(): print cmd.fit("1nmr_0001", x)
for x in cmd.get_names(): optAlign(x, "1nmr_0001")
</source>
<source lang="bash">
# results from fit
0.0
4.50344991684
5.33588504791
5.78613853455
7.25597000122
6.67145586014
3.25131297112
3.36766290665
6.74802017212
5.1579709053
5.96959495544
6.68093347549
4.13217163086
5.51539039612
6.24266338348
6.03838825226
5.01363992691
5.33336305618
6.87617444992
7.797062397

#results from optAlign
RMSD=0.000000
RMSD=4.503450
RMSD=5.335886
RMSD=5.786138
RMSD=7.255970
RMSD=6.671456
RMSD=3.251313
RMSD=3.367663
RMSD=6.748021
RMSD=5.157971
RMSD=5.969595
RMSD=6.680934
RMSD=4.132172
RMSD=5.515390
RMSD=6.242664
RMSD=6.038388
RMSD=5.013640
RMSD=5.333363
RMSD=6.876174
RMSD=7.797062
</source>

== Examples ==
=== Usage ===
==== Syntax ====

CEAlign has the semantic, and syntactic formalism of
<source lang="python">
cealign MASTER, TARGET
</source>
where a post-condition of the algorithm is that the coordinates of the '''MASTER''' protein are unchanged. This allows for easier multi-protein alignments. For example,
<source lang="python">
cealign 1AUE, 1BZ4
cealign 1AUE, 1B68
cealign 1AUE, 1A7V
cealign 1AUE, 1CPR
</source>
will superimpose all the TARGETS onto the MASTER.

=====Examples=====
<source lang="python">
cealign 1cll and i. 42-55, 1ggz and c. A
cealign 1kao, 1ctq
cealign 1fao, 1eaz
</source>

=====Multiple Structure Alignments=====
Use the '''alignto''' command, now provided with cealign. Just type,
<source lang="python">
alignto PROT
</source>
to align all your proteins in PyMOL to the one called, '''PROT'''.

=== Results ===
See '''Changes''' for updates. But, overall, the results here are great.

<gallery>
Image:v7_1fao_1eaz.png|EASY: 1FAO vs. 1EAZ; 96 residues, 1.28 Ang
Image:v7_1cbs_1hmt.png|EASY: 1CBS vs. 1HMT; 128 residues, 2.01 Ang
Image:v7_1a15_1b50.png|MODERATE: 1A15 vs 1B50; 56 residues, 2.54 Ang.
Image:v7_1oan_1s6n.png|EASY: 1OAN vs. 1S6N (state 1); 96 residues aligned to 3.83 Ang. RMSD.
Image:v7_1rlw_1byn.png|HARD: 1RLW to 1BYN; 104 residues; 2.21 Ang.
Image:v7_1ten_3hhr.png|HARD: 1TEN vs. 3HHR; 80 residues, 2.91 Ang.
Image:v7_2sim_1nsb.png|HARD: 2SIM vs. 1NSB; 272 residues, 4.93 Ang.
Image:v7_1cew_1mol.png|HARD: 1CEW vs. 1MOL; 80 residues, 4.03 Ang.
</gallery>

== Installation ==

===Mac OS X (10.5)===
[[Image:Cealign mac os x.png|300px|thumb|center|CEAlign running on Mac OS X (10.5)]]
* Install PyMOL under fink.
* Install Numpy for fink:
<source lang="bash">
/sw/bin/fink install scipy-core-py25
</source>
* Download and install cealign (download instructions below)
<source lang="bash">
sudo /sw/bin/python setup.py install
</source>
* In PyMOL, run the two scripts needed for cealing: "cealign.py" and "qkabsch.py". These are located in the cealign directory you previously downloaded.
* Voila!
* Note that the above python version must match the same version that is used by PyMOL. If you are using the pre-compiled version of MacPyMOL, the above instructions won't work.

===Windows systems===
This is a quick and dirty method to get it working on Win32 right now, more details coming soon.
====Requirements====
* Latest PyMol, installed on your system
* Numpy for python 2.4 -- quick download of just what's needed: http://users.umassmed.edu/Shivender.Shandilya/pymol/numpy.zip[Note: the numpy.zip here are corrupt, users may download the latest 'Numpy for python 2.4' directly from SourceForge.net]
* Pre-compiled ccealign.pyd python module: http://users.umassmed.edu/Shivender.Shandilya/pymol/ccealign.zip
* Modified pymolrc: http://users.umassmed.edu/Shivender.Shandilya/pymol/pymolrc
* cealign.py and qkabsch.py from the Cealign-0.8-RBS package: download below

====Directions====
# Unzip the numpy.zip file, which will give you a folder named '''numpy'''
# Move this entire folder to: C:\Program Files\DeLano Scientific\PyMOL\modules\ (or the corresponding location on your system)
# Unzip ccealign.zip, which will give you a file called '''ccealign.pyd'''
# Move this pyd file to: C:\Program Files\DeLano Scientific\PyMOL\py24\DLLs\ (or the corresponding location on your system)
# Copy the downloaded '''pymolrc''' file to: C:\Program Files\DeLano Scientific\PyMOL\ (or the corresponding location on your system)
# Extract and copy the files cealign.py and qkabsch.py from the Cealign-0.8-RBS package to: C:\Program Files\DeLano Scientific\PyMOL\py24\Lib\ (or the corresponding location on your system)
# Run PyMol and load some molecules
# Run this command in Pymol: '''cealign molecule1, molecule2'''
# Enjoy!

===Gentoo Linux===
Add the science overlay via
layman -a sci
and emerge the cealign plugin
emerge pymol-plugins-cealign

===*nix systems===
====Requirements====
* C compiler
* Python 2.4+ with distutils
* Numpy
** for User-compiled PyMOL: <source lang="python">python setup.py install</source>
** for the precompiled version of PyMOL <source lang="python">python setup.py install --prefix "" --root /DIR_TO/pymol/ext/</source>

====Directions====
# uncompress the distribution file '''cealign-VERSION.tgz'''
# cd cealign-VERSION
# sudo python setup.py install # if you installed by PyMOL by hand
## python setup.py install --prefix "" --root /DIR/TO/pymol/ext/ # if you are using the precompiled binary download
# insert "run DIR_TO_CEALIGN/cealign.py" and "run DIR_TO_CEALIGN/qkabsch.py" into your '''.pymolrc''' file, or just run the two Python scripts by hand.
# load some molecules
# run, '''cealign molecule1, molecule2'''
# enjoy

=====Pre-compiled Hackish Install=====
For those people that prefer to use the pre-compiled version of PyMOL, here are the basics for your install. '''This is a poor method of installing Cealign. I suggest users compile and install their own PyMOL.''' The final goal is to get
# '''ccealign.so''' module into '''PYMOL/ext/lib/python2.4/site-packages'''
# numpy installed (get the numpy directory into (or linked into) '''PYMOL/ext/lib/python2.4/site-packages'''
# and be able to run cealign.py and qkabsch.py from PyMOL.
If you can do the above three steps, '''cealign''' should run from the pre-compiled PyMOL.

In more detail, on a completely fictitious machine --- that is, I created the following commands from a fake machine and I don't expect a copy/paste of this to work '''anywhere''', but the commands should be helpful enough to those who need it:
<source lang="python">
# NOTES:
# This is fake code: don't copy/paste it.
#
# PYMOL='dir to precompiled PyMOL install'
# CEALIGN='dir where you will unpack cealign'
# replace lib with lib64 for x86-64
# install numpy
apt-get install numpy

# link numpy to PyMOL
ln -s /usr/local/lib/python2.4/site-packages/numpy PYMOL/ext/lib/python2.4/site-packages

# download and install Cealign
wget http://www.pymolwiki.org/images/e/ed/Cealign-0.6.tar.bz2
tar -jxvf Cealign-0.6.tar.bz2
cd cealign-0.6
sudo python setup.py build
cp build/lib-XYZ-linux/ccealign.so PYMOL/ext/lib/python2.4/site-packages

# run pymol and try it out
pymol
run CEALIGN/cealign.py
run CEALIGN/qkabsch.py
fetch 1cew 1mol, async=0
cealign 1c, 1m
</source>

== The Code ==
Please unpack and read the documentation. All comments/questions should be directed to Jason Vertrees (javertre _at_ utmb ...dot... edu).

'''LATEST IS v0.8-RBS'''. (Dedicated to Bryan Sutton for allowing me to use his computer for testing.)

=== Version 0.8-RBS ===
* '''Download: [[Media:Cealign-0.8-RBS.tar.bz2|CE Align v0.8-RBS]] (bz2)'''
* '''Download: [[Media:Cealign-0.8-RBS.zip|CE Align v0.8-RBS]] (zip)'''

=== Beta Version 0.9 ===
Use at your own peril. Please report any problems or inconsistent alignments to this discussion page, or to me directly; my email address all over this page.

'''Improvements/Changes''':
* All C++
** So, faster
** comes with the dependencies built in
* No numpy

''' Download: [[Media:Cealign-0.9.zip|CE Align v0.9]] (zip)'''

== Coming Soon ==
* Windows binary
* Linux Binaries (32bit, x86-64)
* Better instructions for precompiled distributions
* Optimization

== Updates ==

=== 2008-03-25 ===
Pure C++ code released. See the beta version above.

=== 2007-04-14 ===
v0.8-RBS source updated. Found the bug that had been plaguing 32-bit machines. This should be the last release for a little while.

Also, I provide the option of aligning based solely upon RMSD or upon the better CE-Score. See the '''References''' for information on the '''CE Score'''.

== Troubleshooting ==

Post your problems/solutions here.

=== Unicode Issues in Python/Numpy ===
'''Problem''': Running/Installing cealign gives
<source lang="python">
Traceback (most recent call last):
File "/home/byron/software/pymol_1.00b17/pymol/modules/pymol/parser.py",
line 308, in parse
File "/home/byron/software/pymol_1.00b17/pymol/modules/pymol/parsing.py",
line 410, in run_file
File "qkabsch.py", line 86, in ?
import numpy
File "/usr/lib/python2.4/site-packages/numpy/__init__.py", line 36, in ?
import core
File "/usr/lib/python2.4/site-packages/numpy/core/__init__.py", line 5, in ?
import multiarray
ImportError: /home/byron/software/pymol/ext/lib/python2.4/site-packages/numpy/core/multiarray.so:
undefined symbol: _PyUnicodeUCS4_IsWhitespace
</source>
where the important line is
<source lang="python">
undefined symbol: _PyUnicodeUCS4_IsWhitespace
</source>

This problem indicates that your Numpy Unicode is using a different byte-size for unicode characters than is the Python distribution your PyMOL is running from. For example, this can happen if you use the pre-built PyMOL and some other pre-built Numpy package.

'''Solution''': Hand-install Numpy.

=== LinAlg Module Not Found ===
'''Problem''': Running CE Align gives the following error message:
<source lang="python">
run qkabsch.py
Traceback (most recent call last):
File "/usr/lib/python2.4/site-packages/pymol/parser.py", line 285, in parse
parsing.run_file(exp_path(args[nest][0]),pymol_names,pymol_names)
File "/usr/lib/python2.4/site-packages/pymol/parsing.py", line 407, in run_file
execfile(file,global_ns,local_ns)
File "qkabsch.py", line 86, in ?
import numpy
File "/usr/lib/python2.4/site-packages/numpy/__init__.py", line 40, in ?
import linalg
ImportError: No module named linalg
</source>

'''Solution''': You do not have the linear algebra module installed (or Python can't find it) on your machine. One workaround is to install [http://www.scipy.org/ Scientific Python]. (on debian/ubuntu this can be done by: sudo apt-get install python-scipy) Another is to reinstall the Numpy package from source, ensuring that you have the necessary requirements for the linear algebra module (linpack, lapack, fft, etc.).

=== CCEAlign & NumPy Modules Not Found ===
'''Problem''': Running CE Align gives the following error message:
<source lang="python">
PyMOL>run cealign.py
Traceback (most recent call last):
File "/home/local/warren/MacPyMOL060530/build/Deployment/MacPyMOL.app/pymol/modules/pymol/parser.py", line 297, in parse
File "/home/local/warren/MacPyMOL060530/build/Deployment/MacPyMOL.app/pymol/modules/pymol/parsing.py", line 408, in run_file
File "/usr/local/pymol/scripts/cealign-0.1/cealign.py", line 59, in ?
from ccealign import ccealign
ImportError: No module named ccealign
run qkabsch.py
Traceback (most recent call last):
File "/home/local/warren/MacPyMOL060530/build/Deployment/MacPyMOL.app/pymol/modules/pymol/parser.py", line 297, in parse
File "/home/local/warren/MacPyMOL060530/build/Deployment/MacPyMOL.app/pymol/modules/pymol/parsing.py", line 408, in run_file
File "qkabsch.py", line 86, in ?
import numpy
ImportError: No module named numpy
</source>

'''Solution''': This problem occurs under [http://www.apple.com/macosx Apple Mac OS X] if (a) the Apple's python executable on your machine (/usr/bin/python, currently version 2.3.5) is superseded by [http://fink.sourceforge.net/ Fink]'s python executable (/sw/bin/python, currently version 2.5) and (b) you are using [http://delsci.com/rel/099/#MacOSX precompiled versions of PyMOL] (MacPyMOL, PyMOLX11Hybrid or PyMOL for Mac OS X/X11). These executables ignore Fink's python and instead use Apple's - so, in order to run CE Align, one must install NumPy (as well as CE Align itself) using Apple's python. To do so, first download the [http://sourceforge.net/project/showfiles.php?group_id=1369&package_id=175103 Numpy source code archive] (currently version 1.0.1), unpack it, change directory to numpy-1.0.1 and specify the full path to Apple's python executable during installation: <tt>sudo /usr/bin/python setup.py install | tee install.log</tt>. Then, donwload the [http://www.pymolwiki.org/index.php/Cealign#The_Code CE Align source code archive] (currently version 0.2), unpack it, change directory to cealign-0.2 and finally install CE Align as follows: <tt>sudo /usr/bin/python setup.py install | tee install.log</tt>.
[[User:Lucajovine|Luca Jovine]] 05:11, 25 January 2007 (CST).

=== The Function SimpAlign() is not found ===
'''Problem''': Running CE Align gives the following error message:
<source lang="python">
PyMOL>cealign 1CLL,1GGZ
Traceback (most recent call last):
File "C:\Program Files (x86)\DeLano Scientific\PyMOL/modules\pymol\parser.py", line 203, in parse
result=apply(kw[nest][0],args[nest],kw_args[nest])
File "py24/Lib/cealign.py", line 177, in cealign
curScore = simpAlign( matA, matB, mol1, mol2, stored.mol1, stored.mol2, align=0, L=len(matA) )
NameError: global name 'simpAlign' is not defined
</source>
I am running PyMOL v. 0.99rc6 on Win XP Professional x64 edition version 2003 sp2 and have followed the windows install procedure as described above.

'''Answer''': This simply means that PyMOL couldn't find the simplAlign function. To let PyMOL know about this, you must run the following commands before running [[cealign]]:
<source lang="python">
run /your/path/to/cealign/qkabsch.py
run /your/path/to/cealign/cealign.py
</source>
but most people that use cealign would just put these two lines in their '''.pymolrc''' file.

=== Short Alignments Don't Work ===
If you are trying to align fewer than 16 residues then use [[align]], [[super]], or [[optAlign]]. CE uses a window size of 8; and to build a path of more than one window, you need 2*8=16 residues. I will insert some code to re-route small alignments to one of the aforementioned alignment algorithms.

=== It Worked A Second Ago! ===
[[Image:Rewind.png|thumb|right|Showing the rewind button to rewind to state 1.]]

If you were using cealign (or alignto) and now the commands don't work -- that is, they return an RMSD, but don't actually superimpose the objects, then you have a simple problem dealing with states. Most likely the cause of this oddness was (1) when you issued "cealign prot1, prot2" one of them was actually an ensemble of states or (2) you are trying to align to proteins with only one state, but are not looking at state one (because the last protein you were considering had more than one state and you quit editing that protein on a state that's not state 1). To fix this, use the rewind button to get the proteins back into state 1 & reissue the cealign/alignto command.

== References ==
Text taken from PubMed and formatted for the wiki. The first reference is the most important for this code.

# Shindyalov IN, Bourne PE. '''Protein structure alignment by incremental combinatorial extension (CE) of the optimal path.''' ''Protein Eng.'' 1998 Sep;11(9):739-47. PMID: 9796821 [PubMed - indexed for MEDLINE]
# Jia Y, Dewey TG, Shindyalov IN, Bourne PE. '''A new scoring function and associated statistical significance for structure alignment by CE.''' ''J Comput Biol.'' 2004;11(5):787-99. PMID: 15700402 [PubMed - indexed for MEDLINE]
# Pekurovsky D, Shindyalov IN, Bourne PE. '''A case study of high-throughput biological data processing on parallel platforms.''' ''Bioinformatics.'' 2004 Aug 12;20(12):1940-7. Epub 2004 Mar 25. PMID: 15044237 [PubMed - indexed for MEDLINE]
# Shindyalov IN, Bourne PE. '''An alternative view of protein fold space.''' ''Proteins.'' 2000 Feb 15;38(3):247-60. PMID: 10713986 [PubMed - indexed for MEDLINE]

== License ==
The CEAlign and all its subprograms that I wrote, are released under the open source Free BSD License (BSDL).

[[Category:Script_Library]]
[[Category:Structure_Alignment|Cealign]]

Cealign

2009-06-08T07:41:45Z

Zhentg: /* Requirements */

== Introduction ==

'''Go directly to [[Cealign#Version_0.8-RBS|DOWNLOAD]]'''

This page is the home page of the open-source CEAlign PyMOL plugin. The CE algorithm is a fast and accurate protein structure alignment algorithm, pioneered by Drs. Shindyalov and Bourne (See
References). There are a few changes from the original CE publication (See Notes).

The source code is implemented in C with the rotations finally done by Numpy in Python. Because the computationally complex portion of the code is written in C, it's quick. That is, on my machines --- relatively fast 64-bit machines --- I can align two 400+ amino acid structures in about 0.300 s with the C++ implementation.

This plugs into PyMol very easily. See [[Cealign#The_Code|the code]] and [[Cealign#Examples|examples]] for installation and usage.

== Comparison to PyMol ==
'''Why should you use this?'''

PyMOL's structure alignment algorithm is fast and robust. However, its first step is to perform a sequence alignment of the two selections. Thus, proteins in the '''twilight zone''' or those having a low sequence identity, may not align well. Because CE is a structure-based alignment, this is not a problem. Consider the following example. The image at LEFT was the result of CE-aligning two proteins (1C0M chain B to 1BCO). The result is '''152''' aligned (alpha carbons) residues (not atoms) at '''4.96 Angstroms'''. The image on the RIGHT shows the results from PyMol's align command: an alignment of '''221 atoms''' (not residues) at an RMSD of '''15.7 Angstroms'''.

<gallery>
Image:cealign_ex1.png|Cealign's results (152 aligned; 4.96 Ang.)
Image:pymol_align.png|PyMol's results (763 atoms; 18.4 Ang. )
</gallery>

=== Fit vs. optAlign ===
====Take Home messages====
* [[fit]] and [[optAlign]] perform nearly equally as well
* if you need an algorithm with an appropriate reference, use [[optAlign]] (references at bottom of page).
* [[fit]] is faster -- if you're aligning many structures, use it over [[optAlign]]

====Discussion====
[[optAlign]] is a function within the [[Cealign]] package that performs the optimal superposition of two objects of equal length. [[optAlign]] follows the Kabsch algorithm which is a closed form, and provably optimal solution to the problem. [[fit]] on the other hand uses the Jacobi rotations to iteratively arrive at the solution of optimal superposition. The difference in error between [[optAilgn]] and [[fit]] seems to be a non-issue (see below) as they both arrive at equivalent solutions for the rotation matrix. The two algorithms are undertake different approaches to orthogonally diagonalizing the correlation matrix.

PyMOL's [[fit]] is fast and works well. If you have to use something with a known reference then check out the "optAlign" function from the qkabsch.py file that comes with this [[Calign]] package. If not, you can just use [[fit]] and avoid installing new software. :-)

optAlign is slower than fit. I just tested both on a sample NMR ensemble; and, while not an extensive validation of "fit" it shows that (1) fit is faster; and (2) fit gets the same exact RMSD as "optAlign" (when optAlign is told to use all atoms, not just CA). To make optAlign use all atoms and not just the alpha-carbon backbones, comment out (that is, put a "#" at the start of lines 183 and 184 in qkabsch.py, where it says "CUT HERE").

<source lang="python">
fetch 1nmr
split_states 1nmr
delete 1nmr

# compare fit and optAlign RMSDs
for x in cmd.get_names(): print cmd.fit("1nmr_0001", x)
for x in cmd.get_names(): optAlign(x, "1nmr_0001")
</source>
<source lang="bash">
# results from fit
0.0
4.50344991684
5.33588504791
5.78613853455
7.25597000122
6.67145586014
3.25131297112
3.36766290665
6.74802017212
5.1579709053
5.96959495544
6.68093347549
4.13217163086
5.51539039612
6.24266338348
6.03838825226
5.01363992691
5.33336305618
6.87617444992
7.797062397

#results from optAlign
RMSD=0.000000
RMSD=4.503450
RMSD=5.335886
RMSD=5.786138
RMSD=7.255970
RMSD=6.671456
RMSD=3.251313
RMSD=3.367663
RMSD=6.748021
RMSD=5.157971
RMSD=5.969595
RMSD=6.680934
RMSD=4.132172
RMSD=5.515390
RMSD=6.242664
RMSD=6.038388
RMSD=5.013640
RMSD=5.333363
RMSD=6.876174
RMSD=7.797062
</source>

== Examples ==
=== Usage ===
==== Syntax ====

CEAlign has the semantic, and syntactic formalism of
<source lang="python">
cealign MASTER, TARGET
</source>
where a post-condition of the algorithm is that the coordinates of the '''MASTER''' protein are unchanged. This allows for easier multi-protein alignments. For example,
<source lang="python">
cealign 1AUE, 1BZ4
cealign 1AUE, 1B68
cealign 1AUE, 1A7V
cealign 1AUE, 1CPR
</source>
will superimpose all the TARGETS onto the MASTER.

=====Examples=====
<source lang="python">
cealign 1cll and i. 42-55, 1ggz and c. A
cealign 1kao, 1ctq
cealign 1fao, 1eaz
</source>

=====Multiple Structure Alignments=====
Use the '''alignto''' command, now provided with cealign. Just type,
<source lang="python">
alignto PROT
</source>
to align all your proteins in PyMOL to the one called, '''PROT'''.

=== Results ===
See '''Changes''' for updates. But, overall, the results here are great.

<gallery>
Image:v7_1fao_1eaz.png|EASY: 1FAO vs. 1EAZ; 96 residues, 1.28 Ang
Image:v7_1cbs_1hmt.png|EASY: 1CBS vs. 1HMT; 128 residues, 2.01 Ang
Image:v7_1a15_1b50.png|MODERATE: 1A15 vs 1B50; 56 residues, 2.54 Ang.
Image:v7_1oan_1s6n.png|EASY: 1OAN vs. 1S6N (state 1); 96 residues aligned to 3.83 Ang. RMSD.
Image:v7_1rlw_1byn.png|HARD: 1RLW to 1BYN; 104 residues; 2.21 Ang.
Image:v7_1ten_3hhr.png|HARD: 1TEN vs. 3HHR; 80 residues, 2.91 Ang.
Image:v7_2sim_1nsb.png|HARD: 2SIM vs. 1NSB; 272 residues, 4.93 Ang.
Image:v7_1cew_1mol.png|HARD: 1CEW vs. 1MOL; 80 residues, 4.03 Ang.
</gallery>

== Installation ==

===Mac OS X (10.5)===
[[Image:Cealign mac os x.png|300px|thumb|center|CEAlign running on Mac OS X (10.5)]]
* Install PyMOL under fink.
* Install Numpy for fink:
<source lang="bash">
/sw/bin/fink install scipy-core-py25
</source>
* Download and install cealign (download instructions below)
<source lang="bash">
sudo /sw/bin/python setup.py install
</source>
* In PyMOL, run the two scripts needed for cealing: "cealign.py" and "qkabsch.py". These are located in the cealign directory you previously downloaded.
* Voila!
* Note that the above python version must match the same version that is used by PyMOL. If you are using the pre-compiled version of MacPyMOL, the above instructions won't work.

===Windows systems===
This is a quick and dirty method to get it working on Win32 right now, more details coming soon.
====Requirements====
* Latest PyMol, installed on your system
* Numpy for python 2.4 -- quick download of just what's needed: http://users.umassmed.edu/Shivender.Shandilya/pymol/numpy.zip[Note: the numpy.zip here are corrupt, users may download the latest 'Numpy for python 2.4' directly from SourceForge.net]
* Pre-compiled ccealign.pyd python module: http://users.umassmed.edu/Shivender.Shandilya/pymol/ccealign.zip
* Modified pymolrc: http://users.umassmed.edu/Shivender.Shandilya/pymol/pymolrc
* cealign.py and qkabsch.py from the Cealign-0.8-RBS package: download below

====Directions====
# Unzip the numpy.zip file, which will give you a folder named '''numpy'''
# Move this entire folder to: C:\Program Files\DeLano Scientific\PyMOL\modules\ (or the corresponding location on your system)
# Unzip ccealign.zip, which will give you a file called '''ccealign.pyd'''
# Move this pyd file to: C:\Program Files\DeLano Scientific\PyMOL\py24\DLLs\ (or the corresponding location on your system)
# Copy the downloaded '''pymolrc''' file to: C:\Program Files\DeLano Scientific\PyMOL\ (or the corresponding location on your system)
# Extract and copy the files cealign.py and qkabsch.py from the Cealign-0.8-RBS package to: C:\Program Files\DeLano Scientific\PyMOL\py24\Lib\ (or the corresponding location on your system)
# Run PyMol and load some molecules
# Run this command in Pymol: '''cealign molecule1, molecule2'''
# Enjoy!

===Gentoo Linux===
Add the science overlay via
layman -a sci
and emerge the cealign plugin
emerge pymol-plugins-cealign

===*nix systems===
====Requirements====
* C compiler
* Python 2.4+ with distutils
* Numpy
** for User-compiled PyMOL: <source lang="python">python setup.py install</source>
** for the precompiled version of PyMOL <source lang="python">python setup.py install --prefix "" --root /DIR_TO/pymol/ext/</source>

====Directions====
# uncompress the distribution file '''cealign-VERSION.tgz'''
# cd cealign-VERSION
# sudo python setup.py install # if you installed by PyMOL by hand
## python setup.py install --prefix "" --root /DIR/TO/pymol/ext/ # if you are using the precompiled binary download
# insert "run DIR_TO_CEALIGN/cealign.py" and "run DIR_TO_CEALIGN/qkabsch.py" into your '''.pymolrc''' file, or just run the two Python scripts by hand.
# load some molecules
# run, '''cealign molecule1, molecule2'''
# enjoy

=====Pre-compiled Hackish Install=====
For those people that prefer to use the pre-compiled version of PyMOL, here are the basics for your install. '''This is a poor method of installing Cealign. I suggest users compile and install their own PyMOL.''' The final goal is to get
# '''ccealign.so''' module into '''PYMOL/ext/lib/python2.4/site-packages'''
# numpy installed (get the numpy directory into (or linked into) '''PYMOL/ext/lib/python2.4/site-packages'''
# and be able to run cealign.py and qkabsch.py from PyMOL.
If you can do the above three steps, '''cealign''' should run from the pre-compiled PyMOL.

In more detail, on a completely fictitious machine --- that is, I created the following commands from a fake machine and I don't expect a copy/paste of this to work '''anywhere''', but the commands should be helpful enough to those who need it:
<source lang="python">
# NOTES:
# This is fake code: don't copy/paste it.
#
# PYMOL='dir to precompiled PyMOL install'
# CEALIGN='dir where you will unpack cealign'
# replace lib with lib64 for x86-64
# install numpy
apt-get install numpy

# link numpy to PyMOL
ln -s /usr/local/lib/python2.4/site-packages/numpy PYMOL/ext/lib/python2.4/site-packages

# download and install Cealign
wget http://www.pymolwiki.org/images/e/ed/Cealign-0.6.tar.bz2
tar -jxvf Cealign-0.6.tar.bz2
cd cealign-0.6
sudo python setup.py build
cp build/lib-XYZ-linux/ccealign.so PYMOL/ext/lib/python2.4/site-packages

# run pymol and try it out
pymol
run CEALIGN/cealign.py
run CEALIGN/qkabsch.py
fetch 1cew 1mol, async=0
cealign 1c, 1m
</source>

== The Code ==
Please unpack and read the documentation. All comments/questions should be directed to Jason Vertrees (javertre _at_ utmb ...dot... edu).

'''LATEST IS v0.8-RBS'''. (Dedicated to Bryan Sutton for allowing me to use his computer for testing.)

=== Version 0.8-RBS ===
* '''Download: [[Media:Cealign-0.8-RBS.tar.bz2|CE Align v0.8-RBS]] (bz2)'''
* '''Download: [[Media:Cealign-0.8-RBS.zip|CE Align v0.8-RBS]] (zip)'''

=== Beta Version 0.9 ===
Use at your own peril. Please report any problems or inconsistent alignments to this discussion page, or to me directly; my email address all over this page.

'''Improvements/Changes''':
* All C++
** So, faster
** comes with the dependencies built in
* No numpy

''' Download: [[Media:Cealign-0.9.zip|CE Align v0.9]] (zip)'''

== Coming Soon ==
* Windows binary
* Linux Binaries (32bit, x86-64)
* Better instructions for precompiled distributions
* Optimization

== Updates ==

=== 2008-03-25 ===
Pure C++ code released. See the beta version above.

=== 2007-04-14 ===
v0.8-RBS source updated. Found the bug that had been plaguing 32-bit machines. This should be the last release for a little while.

Also, I provide the option of aligning based solely upon RMSD or upon the better CE-Score. See the '''References''' for information on the '''CE Score'''.

== Troubleshooting ==

Post your problems/solutions here.

=== Unicode Issues in Python/Numpy ===
'''Problem''': Running/Installing cealign gives
<source lang="python">
Traceback (most recent call last):
File "/home/byron/software/pymol_1.00b17/pymol/modules/pymol/parser.py",
line 308, in parse
File "/home/byron/software/pymol_1.00b17/pymol/modules/pymol/parsing.py",
line 410, in run_file
File "qkabsch.py", line 86, in ?
import numpy
File "/usr/lib/python2.4/site-packages/numpy/__init__.py", line 36, in ?
import core
File "/usr/lib/python2.4/site-packages/numpy/core/__init__.py", line 5, in ?
import multiarray
ImportError: /home/byron/software/pymol/ext/lib/python2.4/site-packages/numpy/core/multiarray.so:
undefined symbol: _PyUnicodeUCS4_IsWhitespace
</source>
where the important line is
<source lang="python">
undefined symbol: _PyUnicodeUCS4_IsWhitespace
</source>

This problem indicates that your Numpy Unicode is using a different byte-size for unicode characters than is the Python distribution your PyMOL is running from. For example, this can happen if you use the pre-built PyMOL and some other pre-built Numpy package.

'''Solution''': Hand-install Numpy.

=== LinAlg Module Not Found ===
'''Problem''': Running CE Align gives the following error message:
<source lang="python">
run qkabsch.py
Traceback (most recent call last):
File "/usr/lib/python2.4/site-packages/pymol/parser.py", line 285, in parse
parsing.run_file(exp_path(args[nest][0]),pymol_names,pymol_names)
File "/usr/lib/python2.4/site-packages/pymol/parsing.py", line 407, in run_file
execfile(file,global_ns,local_ns)
File "qkabsch.py", line 86, in ?
import numpy
File "/usr/lib/python2.4/site-packages/numpy/__init__.py", line 40, in ?
import linalg
ImportError: No module named linalg
</source>

'''Solution''': You do not have the linear algebra module installed (or Python can't find it) on your machine. One workaround is to install [http://www.scipy.org/ Scientific Python]. (on debian/ubuntu this can be done by: sudo apt-get install python-scipy) Another is to reinstall the Numpy package from source, ensuring that you have the necessary requirements for the linear algebra module (linpack, lapack, fft, etc.).

=== CCEAlign & NumPy Modules Not Found ===
'''Problem''': Running CE Align gives the following error message:
<source lang="python">
PyMOL>run cealign.py
Traceback (most recent call last):
File "/home/local/warren/MacPyMOL060530/build/Deployment/MacPyMOL.app/pymol/modules/pymol/parser.py", line 297, in parse
File "/home/local/warren/MacPyMOL060530/build/Deployment/MacPyMOL.app/pymol/modules/pymol/parsing.py", line 408, in run_file
File "/usr/local/pymol/scripts/cealign-0.1/cealign.py", line 59, in ?
from ccealign import ccealign
ImportError: No module named ccealign
run qkabsch.py
Traceback (most recent call last):
File "/home/local/warren/MacPyMOL060530/build/Deployment/MacPyMOL.app/pymol/modules/pymol/parser.py", line 297, in parse
File "/home/local/warren/MacPyMOL060530/build/Deployment/MacPyMOL.app/pymol/modules/pymol/parsing.py", line 408, in run_file
File "qkabsch.py", line 86, in ?
import numpy
ImportError: No module named numpy
</source>

'''Solution''': This problem occurs under [http://www.apple.com/macosx Apple Mac OS X] if (a) the Apple's python executable on your machine (/usr/bin/python, currently version 2.3.5) is superseded by [http://fink.sourceforge.net/ Fink]'s python executable (/sw/bin/python, currently version 2.5) and (b) you are using [http://delsci.com/rel/099/#MacOSX precompiled versions of PyMOL] (MacPyMOL, PyMOLX11Hybrid or PyMOL for Mac OS X/X11). These executables ignore Fink's python and instead use Apple's - so, in order to run CE Align, one must install NumPy (as well as CE Align itself) using Apple's python. To do so, first download the [http://sourceforge.net/project/showfiles.php?group_id=1369&package_id=175103 Numpy source code archive] (currently version 1.0.1), unpack it, change directory to numpy-1.0.1 and specify the full path to Apple's python executable during installation: <tt>sudo /usr/bin/python setup.py install | tee install.log</tt>. Then, donwload the [http://www.pymolwiki.org/index.php/Cealign#The_Code CE Align source code archive] (currently version 0.2), unpack it, change directory to cealign-0.2 and finally install CE Align as follows: <tt>sudo /usr/bin/python setup.py install | tee install.log</tt>.
[[User:Lucajovine|Luca Jovine]] 05:11, 25 January 2007 (CST).

=== The Function SimpAlign() is not found ===
'''Problem''': Running CE Align gives the following error message:
<source lang="python">
PyMOL>cealign 1CLL,1GGZ
Traceback (most recent call last):
File "C:\Program Files (x86)\DeLano Scientific\PyMOL/modules\pymol\parser.py", line 203, in parse
result=apply(kw[nest][0],args[nest],kw_args[nest])
File "py24/Lib/cealign.py", line 177, in cealign
curScore = simpAlign( matA, matB, mol1, mol2, stored.mol1, stored.mol2, align=0, L=len(matA) )
NameError: global name 'simpAlign' is not defined
</source>
I am running PyMOL v. 0.99rc6 on Win XP Professional x64 edition version 2003 sp2 and have followed the windows install procedure as described above.

'''Answer''': This simply means that PyMOL couldn't find the simplAlign function. To let PyMOL know about this, you must run the following commands before running [[cealign]]:
<source lang="python">
run /your/path/to/cealign/qkabsch.py
run /your/path/to/cealign/cealign.py
</source>
but most people that use cealign would just put these two lines in their '''.pymolrc''' file.

=== Short Alignments Don't Work ===
If you are trying to align fewer than 16 residues then use [[align]], [[super]], or [[optAlign]]. CE uses a window size of 8; and to build a path of more than one window, you need 2*8=16 residues. I will insert some code to re-route small alignments to one of the aforementioned alignment algorithms.

=== It Worked A Second Ago! ===
[[Image:Rewind.png|thumb|right|Showing the rewind button to rewind to state 1.]]

If you were using cealign (or alignto) and now the commands don't work -- that is, they return an RMSD, but don't actually superimpose the objects, then you have a simple problem dealing with states. Most likely the cause of this oddness was (1) when you issued "cealign prot1, prot2" one of them was actually an ensemble of states or (2) you are trying to align to proteins with only one state, but are not looking at state one (because the last protein you were considering had more than one state and you quit editing that protein on a state that's not state 1). To fix this, use the rewind button to get the proteins back into state 1 & reissue the cealign/alignto command.

== References ==
Text taken from PubMed and formatted for the wiki. The first reference is the most important for this code.

# Shindyalov IN, Bourne PE. '''Protein structure alignment by incremental combinatorial extension (CE) of the optimal path.''' ''Protein Eng.'' 1998 Sep;11(9):739-47. PMID: 9796821 [PubMed - indexed for MEDLINE]
# Jia Y, Dewey TG, Shindyalov IN, Bourne PE. '''A new scoring function and associated statistical significance for structure alignment by CE.''' ''J Comput Biol.'' 2004;11(5):787-99. PMID: 15700402 [PubMed - indexed for MEDLINE]
# Pekurovsky D, Shindyalov IN, Bourne PE. '''A case study of high-throughput biological data processing on parallel platforms.''' ''Bioinformatics.'' 2004 Aug 12;20(12):1940-7. Epub 2004 Mar 25. PMID: 15044237 [PubMed - indexed for MEDLINE]
# Shindyalov IN, Bourne PE. '''An alternative view of protein fold space.''' ''Proteins.'' 2000 Feb 15;38(3):247-60. PMID: 10713986 [PubMed - indexed for MEDLINE]

== License ==
The CEAlign and all its subprograms that I wrote, are released under the open source Free BSD License (BSDL).

[[Category:Script_Library]]
[[Category:Structure_Alignment|Cealign]]