[Chimera-users] Massive superposing...

[Eric Pettersen]

On Jan 21, 2010, at 3:15 PM, Yasser Almeida Hernández wrote:

> Hi...
> I'm performing a massive superposing of aminoacids in different  
> conformations. The first model is a arginine "template" reference  
> structure that is the fixed model in all matching iterations. I've  
> loaded 236 models and i want to match the arginine guanidine group  
> on the template guanidine. In each iteration i execute two matching,  
> with different atoms selections. The first selection is the "right"  
> selection: CD,NE,CZ,NH1,NH2 and the second is the "inverse"  
> selection: CD,NE,CZ,NH2,NH1. I run the two superpositions, and i  
> want to show that with the minimal RMSD between the guanidine  
> groups. This is the code:

Well, I can see a variety of problems.  Starting from the top...

> import chimera
> from chimera import runCommand
>
> right = "@CD,NE,CZ,NH1,NH2"
> inverse = "@CD,NE,CZ,NH2,NH1"
>
> for i in range(2,328):

'i' will range from 2 to 327, yet you say you've loaded 236 models.   
I'm guessing you actually want 'i' to vary from 1 to 235 (or 236 if  
the "236 models" you mention don't include the template model).   
Therefore you should be using "range(1, 236)".

> 	print "Superposing..."
> 	runCommand("match #"+str(i)+right+" #0"+right+"; close 1")

What is with all the "close 1" on the end of all your commands?  Do  
you mean "wait 1" (to draw a new frame)?  "Close 1" will close model  
1.  It makes no sense to do it over and over.  Also, since you are not  
saving any images in this script you don't need "wait" commands either.

> 	right_rmsd = runCommand("rmsd #"+str(i)+right+" #0"+right+"; close  
> 1")

'runCommand' doesn't return a value.  Actually it does, but it's  
whether or not a 'wait' is needed after the command string to update  
the display, not a value related to the commands' functions.  It  
doesn't return a command-specific value pretty much for the reason  
exemplified by the string you used:  it is two commands strung  
together, so which command's value should be returned?

At any rate, to get the RMSD you will need to call the underlying rmsd  
function directly.  It is in the Midas module.  So this:

	from Midas import rmsd
	right_rmsd = rmsd("#"+str(i)+right, "#0"+right)

>        runCommand("match #"+str(i)+right+" #0"+right+"; close 1")
> 	inverse = runCommand("rmsd #"+str(i)+inverse+" #0"+right+"; close 1")

You need to assign the value (once you change to the rmsd command) to  
inverse_rmsd, not inverse.  After the above statement, inverse no  
longer has the value "@CD,NE,CZ,NH2,NH1", so a lot of things will  
begin to go badly wrong at this point.

>
> 	if right < inverse:

if right_rmsd < inverse_rmsd:

> 		print "Minimum RMSD =",right,"\n\n"

right_rmsd, not 'right'.

> 		runCommand("match #"+str(i)+right+" #0"+right+"; close 1")
>
> 	elif right > inverse:

Again, need "_rmsd" on these.  So the script does nothing if the RMSDs  
are exactly equal (unlikely, I know)?  It seems like you could just  
use a simple "else:" here instead of the elif.

> 		print "Minimum RMSD =",inverse_rmsd,"\n\n"
> 		runCommand("match #"+str(i)+right+" #0"+right+"; close 1")
>
> In principle this script should show the matching with the minimum  
> RMSD, but it doesn't work.
>
> Please, help...
> Thanks

--Eric

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





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