Menu
▾
▴
[PyMOL] aforementioned NCS snapshot script
|
From: Seth H. <sh...@ms...> - 2003-04-16 23:40:40
|
I was wondering how the attachment would fare on the email list.
Apparently, "not well" so I just include it in line below, but watch
out for proper spacing/line returns in the process.
### Purpose: Get snapshots of each of several NCS-related molecules
in identical position
### Usage: Line up your view on molecule 1 (i.e. chain A) and type "run NCS.py"
### Notes: Turn on any maps, cartoons, colorings, whatever it is that
you want to examine
### Make sure you have NCS operators defined in a CNS-format file and that
### this file is appropriately called in this script (current default
is ncs_matrices.list)
### Be careful that these operators are the ones to transform chain A
### onto chain X, rather than the inverse which CNS lists in the bottom half of
### its NCS matrices files
### Then type "run NCS.py", and Bob's your mother's brother, or at
least a second cousin
### Output: The script writes .png snapshots for each NCS molecule,
called molX.png where
### X is some number. Also, it will output each "set_view" command
that it came up
### with so that if you want to go back and interactively look at
mol4, for example, you
### can just type @mol4 in pymol (while the molecule is already
loaded, etc.) and it
### will take you there by appropriately setting the view.
### If you have ImageMagick installed, you can uncomment the penultimate line
### about os.system call to montage. The montage command takes the preceding
### mol*.png files produced of the views and makes a single png file with all
### of the panels represented. The man page for montage has the options fully
### explained, but briefly the way I use it has -geometry with
options for x and y
### pixels for each panel plus how wide of a buffer or border you
want between panels
### (e.g. +5+5) and the -tiel 4x2 describes how you want the panels
laid out, in this
### case 4 columns and 2 rows. These options are followed by the
input file names
### (mol*.png) and the output (montage.png). You can, of course,
also run it manually
### for only those cases that you want since it takes a little longer
than the main
### part of the script. Also, if you want to ray trace each frame
you need to uncomment
### the cmd.ray() line in the outview function.
### Caveats: Um, basically, in python, I have no idea what I am doing.
### This is the first, and so far only, Python script I have ever written.
### I'm sure there are much, much better ways to do a lot of this, but I have
### just plowed ahead and spelled things out in a basic klunky style.
### Feel free to send me improvements or suggestions, like better,
more flexible ways
### to get the NCS operators out of the file so it won't be so
stringently tied to
### the CNS file format. It would probably be best to have pymol
extract the NCS
### relationshiops itself, for example. Also, I have only used this
on one case
### (8 NCS copies, though, which is what motivated this script) but I
think it will work
### generally.
### Overview Method:
### Get starting view for mol1
### Get NCS matrices from file for chain A->B, A->C, A->D, ... A->X
### Multiply rotation matrix of starting view (first 9) by NCS
rotation matrix for each mol
### Get origin position from starting view (elements 12,13,14 of starting view
### recalling that array begins at elemnt 0)
### Multiply by origin x,y,z by ncs rotation matrix as with view
### Additionally, add the x,y,z translations to the result to get the
new origin position
### Note: Translations for inverse matrices sometimes, but not
always, match those
### for the original matrices. I don't fully understand
this, but that's
### probably just me. In any case, be
careful about getting the right matrices.
### The new view is written out as a "set_view" command, saved in a
file, also applied
### and a png snapshot taken before going on to the next molecule.
### PyMOL's get/set_view matrix:
### Note that what I learned of pymol's set_view is that camera
position (elements 9-11 of
### the set view matrix) is relative to the origin. Since we're
moving the origin to the
### appropriate new location, we can just leave the relative camera
position as it was.
### Same seems to go for the clipping planes (elements 15,16). The
final element is the
### orthoscopic flag, and this too is simply copied over from the first view.
### This script has been hacked out (and I do mean hacked) by Seth Harris.
### sh...@ms... for any suggestions/feedback/gratitude
### Thanks to Warren DeLano for PyMOL!
from pymol import cmd
from Numeric import *
from LinearAlgebra import *
import re
### define function to put out set_view command for a given molecule
def outview (mol, n):
of=open('mol'+repr(mol),'w')
of.write('### for molecule '+repr(mol)+'\n')
of.write('set_view (\\\n')
of.write(repr(n[0])+', '+repr(n[1])+',
'+repr(n[2])+',\\\n')
of.write(repr(n[3])+', '+repr(n[4])+',
'+repr(n[5])+',\\\n')
of.write(repr(n[6])+', '+repr(n[7])+',
'+repr(n[8])+',\\\n')
of.write(repr(n[9])+', '+repr(n[10])+',
'+repr(n[11])+',\\\n')
of.write(repr(n[12])+', '+repr(n[13])+',
'+repr(n[14])+',\\\n')
of.write(repr(n[15])+', '+repr(n[16])+',
'+repr(n[17])+')\n')
cmd.set_view(n)
#uncomment following if you want each
snapshot ray traced
#cmd.ray()
cmd.png('mol'+repr(mol)+'.png')
#filename=input("File name for montage (e.g. residue name&number)
[montage.png]: ")
# get starting view and put that out for molecule 1 unchanged
m = cmd.get_view(0)
mol=1
outview(mol, m)
start=array(((m[0],m[1],m[2]),(m[3],m[4],m[5]),(m[6],m[7],m[8])), Float)
origin=array(((m[12],m[13],m[14])), Float)
count=0
# search for matrix lines in ncs_matrices.list file
# note that these are all matrices for moving a->x, where x is the NCS copy
# this rotation matrix will be applied to the original view and to the
# original origin. The resulting origin position will then be modified by
# adding the translation shift to give the final new origin
# Strangely, note that while the inverse operator works fine on the rotation
# matrix, the inverse translation that CNS comes up with is not simply the
# opposite sign of the original translation values. Hmmm.
# Anyway, it almost killed me, but finally this seems to work correctly.
p=re.compile('.*matrix.*')
for rot in open('ncs_matrices.list'):
p=re.compile('.*matrix.*')
x = p.match(rot)
if x or count>0:
print 'Match found: ', rot
# search CNS-style file for three numbers separated
by one or more spaces
# between parentheses in each of these lines. Once
found use count to
# get 4 lines total for the whole matrix
nums=re.compile(r'.*\(\s+([0-9.-]*)\s+([0-9.-]*)\s+([0-9.-]*).*\).*')
y = nums.match(rot)
index=count*3
ncs=list((y.group(1), y.group(2), y.group(3)))
if count==0:
nl=ncs
else:
new=nl
nl[-1:]=nl[-1:]+ncs
nl=new
count=count+1
if count==4:
mol=mol+1
# get rotation matrix from ncslist (first
nine elements)
rotmat=
array(((float(nl[0]),float(nl[1]),float(nl[2])),\
(float(nl[3]),float(nl[4]),float(nl[5])),\
(float(nl[6]),float(nl[7]),float(nl[8]))), Float)
#rotinv = inverse(rotmat)
print '\nstarting position matrix is:'
print start
print '\norigin is:'
print origin
print '\nncs rotation matrix is:'
print rotmat
p = matrixmultiply(rotmat,start)
print '\nproduct of matrixmultiply rotmat, start is'
print p
no = matrixmultiply(rotmat,origin)
print '\nproduct of matrixmultiply for
origin, prior to translation is:'
print no
print '\ntranslation shift will be:'
print nl[9],nl[10],nl[11]
nopos =
array(((no[0]+float(nl[9]),no[1]+float(nl[10]),no[2]+float(nl[11]))),
Float)
n =
array((p[0][0],p[0][1],p[0][2],p[1][0],p[1][1],p[1][2],\
p[2][0],p[2][1],p[2][2],m[9],m[10],m[11],\
nopos[0],nopos[1],nopos[2],m[15],m[16],m[17]))
print 'Outputting view for molecule '+repr(mol)
outview (mol,n)
count=0
#os.system ("montage -geometry 640x480+5+5 -tile 4x2 mol*.png montage.png")
print 'All done'
--
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Seth F. Harris, PhD
Agard Laboratory
University of California, San Francisco
415-502-2930
|