[Chimera-users] problem with the scripting of rms deviations

[Eric Pettersen]

Hi Sebastian,
	Between Chimera and your script, there are a couple of problems with  
the RMSD calculation.  The executive summary of the solution is:   
explicitly list the atoms being matched in the RMSD calculation, i.e.:

rmsd #0 at pg@pb at o1b@o2b... #1.1 at p3@p2 at o6@o10...

	To figure out the names that Chimera assigned to your SDF molecules,  
you will have to open one and look at it.  The above approach  
explicitly ignores the hydrogens, so you can skip that in your  
script.  I would imagine that you might want to skip matching the  
oxygens of the terminal phosphate as well since that phosphate freely  
rotates and it's impossible to know which oxygen should match with  
which.

	So, the problems are:

1) Chimera doesn't do a connectivity analysis of the atoms in the rmsd  
command to figure out which to match with which.  In cases where the  
atoms aren't explicitly specified, it orders them by name (if the same  
name then by input order).  The names given to the atoms in the SDF  
file is completely different from those in the PDB file.

2) When you do "set operations" like intersect (&) and union (|),  
Chimera may lose track of the ordering of the atoms.  Chimera should  
do a better job of maintaining the ordering with intersection (and I  
will work on that) nonetheless commands like "rmsd #0 #1" or "rmsd  
#0:12 #1:1" will work (assuming the names of the atoms match) whereas  
more complicated commands involving intersection, union, or filtering  
by attribute values may not.

	I'll be opening a ticket in our Trac bug database for maintaining  
better ordering on intersection, with you on the cc list.

--Eric

                         Eric Pettersen
                         UCSF Computer Graphics Lab
                         http://www.cgl.ucsf.edu


On Aug 24, 2009, at 4:43 AM, Sebastian Kruggel wrote:

> Dear Chimera team and users,
>
> I tried to write a small chimera / python script to calculate rmsd  
> to the crystal structure for a bunch of docked ligands and transform  
> the values in a nice csv - because fred22 doesn't do this for me...
>
> The core function rms2ref is the one calculating the rmsd, than I  
> save the reply log and afterwards transform this text file.
>
> Things worked out fine BUT than I started the script several times,  
> and I got different orders in my list every time. Because I want to  
> assign rmsd values to scores etc it is essential to keep the order  
> provided in the sdf/mol2. I thought about problems with the format  
> but the problem exists with mol2 and sdf files. Than I thought,  
> chimera would change the order in opening the multimolecule files -  
> but saving the #1.1 molecule and opening it in a new session showed  
> - even more strange - that the molecule is exactly the same every  
> time. So there is only the possibility that either the calculation  
> of the rmsd or the order of the list (extraction of the  
> replyLog.txt) doesn't work - which I *really* don't understand at  
> all...
>
> Here is a small example where the rmsds-list has different orders  
> every time I start the script...
>
>
> ##############################
> from chimera import runCommand
> from chimera.tkgui import saveReplyLog
>
> def rms2ref(ref, pose, no):
>   runCommand('open %s; open %s' % (ref, pose))
>   i = 1
>   while i <= no:
>      # only hetatms are included --> ~@/element=h
>      runCommand('sel #1.%i & ~@/element=h' % i)
>      runCommand('rmsd #0 sel')
>      i += 1
>
> rms2ref('lig_cryst.pdb', 'poses.sdf', 10)
>
> # generate output from replyLog.txt
>
> saveReplyLog('replyLog.txt')
>
> rmsds = []
> for i in open('replyLog.txt').readlines():
>   if i.startswith('RMSD between'):
>      rmsds.append(i.split()[-2])
>
> print rmsds
> ##############################
>
>
> Maybe anybody can help me and tell me about the probably stupid  
> mistake I am making - would be really great! I attached the  
> lig_cryst.pdb and poses.sdf files just in case that anybody wants to  
> try out ;-)

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://plato.cgl.ucsf.edu/pipermail/chimera-users/attachments/20090824/b5649e73/attachment.html>