[PyMOL] how to get RMSD from align command

Discussion:

Thomas Evangelidis

2009-10-29 02:48:36 UTC

Simple question, it must have been answered before but couldn't find
it so far:

how can I get the RMSD value from the align command in a python
script? cmd.align() returnes a tuple of 8 numbers and none of them is
the actual RMSD value I get when I align these 2 structures manually.

thanks,Tom

Warren DeLano

2009-10-29 03:29:32 UTC

Permalink

Tom,

The complication with cmd.align() is that it is doing a whole lot more than a simple alignment. The first number is in fact the RMS, but it covers only the subset of the input atoms remaining after refinement is completed. The count of aligned atoms is the second field.

If you're looking for exact RMS fit values over a well-defined set of atoms, try using cmd.pair_fit(sele1, sele2) instead, for example:

load $TUT/1hpv.pdb

create loopA, A/46-55/

create loopB, B/46-55/

show sticks, loop*

print cmd.pair_fit("loopA////CA", "loopB////CA")

Cheers,
Warren

-----Original Message-----
From: Thomas Evangelidis [mailto:***@mbg.duth.gr]
Sent: Wed 10/28/2009 8:09 PM
To: pymol-***@lists.sourceforge.net
Subject: [PyMOL] how to get RMSD from align command

Simple question, it must have been answered before but couldn't find
it so far:

how can I get the RMSD value from the align command in a python
script? cmd.align() returnes a tuple of 8 numbers and none of them is
the actual RMSD value I get when I align these 2 structures manually.

thanks,Tom

------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
PyMOL-users mailing list (PyMOL-***@lists.sourceforge.net)
Info Page: https://lists.sourceforge.net/lists/listinfo/pymol-users
Archives: http://www.mail-archive.com/pymol-***@lists.sourceforge.net

Thomas Evangelidis

2009-10-29 12:14:57 UTC

Permalink

Hi Warren,

Now I am more confused. I used align command to measure the overall
RMSD between homologous structures. Apparently align is not
appropriate if not sufficient sequence similarity is present. I later
came across CEalign plugin, which does structure-based
superimposition. This command also returns an RMSD value, which
happens to be higher for the set of homologous structures I analysed.
Based on this I concluded that align command is more appropriate to
measure measure structural similarity, but I'm not sure now.

My objective is to select between a small set of homologous structures
(co-crystallized with ligands), the one with the highest structural
similarity and use it to place a dummy ligand into one where it's
missing. I.e. "1ebh" lacks a ligand, so I compared it with "1ebg",
"1els", "1one", "2one", "5enl", "6enl", "7enl" which have
co-crystallized ligands. According to both align and cealign, 5ENL is
the most similar with overall RMSD 0.246 and 0.420907 respectively.
The question is which command produces the best superimposition
provided that these proteins belong to the enolase family and thus
have high sequence and structural similarity? Based on the RMSD
values, align does, but as you said this value pertains to only a
subset of atoms after refinement.

Apparently pair_fit command cannot be used for my purpose as it
requires identical selections.

Tom

Post by Warren DeLano
Tom,
The complication with cmd.align() is that it is doing a whole lot
more than a simple alignment. The first number is in fact the RMS,
but it covers only the subset of the input atoms remaining after
refinement is completed. The count of aligned atoms is the second
field.
If you're looking for exact RMS fit values over a well-defined set
load $TUT/1hpv.pdb
create loopA, A/46-55/
create loopB, B/46-55/
show sticks, loop*
print cmd.pair_fit("loopA////CA", "loopB////CA")
Cheers,
Warren
-----Original Message-----
Sent: Wed 10/28/2009 8:09 PM
Subject: [PyMOL] how to get RMSD from align command
Simple question, it must have been answered before but couldn't find
how can I get the RMSD value from the align command in a python
script? cmd.align() returnes a tuple of 8 numbers and none of them is
the actual RMSD value I get when I align these 2 structures manually.
thanks,Tom
------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Info Page: https://lists.sourceforge.net/lists/listinfo/pymol-users

----- End message from ***@delsci.com -----

Andreas Forster

2009-10-29 09:42:05 UTC

Permalink

Hey Thomas,

if you want to use align anyway, make sure to use the quite=0 option.
The "quiet" option (if present) is set to zero by default for parsed
PyMOL commands, but is not set for Python API calls.

align
is nearly equal to
cmd.align(quiet=0)

Thus, if you want to get rmsd output, include quiet=0 and run your
script with output redirection.

Andreas

Post by Thomas Evangelidis
Simple question, it must have been answered before but couldn't find
how can I get the RMSD value from the align command in a python
script? cmd.align() returnes a tuple of 8 numbers and none of them is
the actual RMSD value I get when I align these 2 structures manually.
thanks,Tom
------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Info Page: https://lists.sourceforge.net/lists/listinfo/pymol-users

Thomas Evangelidis

2009-10-29 11:40:56 UTC

Permalink

Thank you Adnreas.

I'm straggling with cealign now. When I run it in batch mode I get
"Selector-Error: Invalid selection name" for certain structures, which
in fact are very close homologs, whereas the number of RMSD values I
get varies (sometimes 3, 5, etc). Here's the code:

query_template_chains = {
"1ebh" : ["1ebg", "1els", "1one", "2one",
"5enl", "6enl", "7enl"]
}

for query in query_template_chains.keys():
for template in query_template_chains[query]:
cmd.fetch(query)
cmd.fetch(template)
print "Superimposing ", query, " onto ", template
cmd.do("cealign "+query+" and c. A, "+template+" and c. A")
cmd.delete(query)
cmd.delete(template)

Do you know what's wrong?

Thomas

Post by Andreas Forster
Hey Thomas,
if you want to use align anyway, make sure to use the quite=0 option.
The "quiet" option (if present) is set to zero by default for parsed
PyMOL commands, but is not set for Python API calls.
align
is nearly equal to
cmd.align(quiet=0)
Thus, if you want to get rmsd output, include quiet=0 and run your
script with output redirection.
Andreas
On Thu, Oct 29, 2009 at 2:48 AM, Thomas Evangelidis

----- End message from ***@gmail.com -----

Warren DeLano

2009-10-29 13:55:59 UTC

Permalink

Tom,

Apologies for the confusion, but there really isn't any way to measure overall structural similarity in PyMOL other than pair_fit, but you have to do the work of explicitly stating which pairs are to be aligned.

The pair_fit selections don't have to be identical, but they do have to correspond in a pairwise fashion. In other words, you can compare any number of C-alpha positions so long as the same number of atomic positions are specified in each structure.

Cheers,
Warren

-----Original Message-----
From: Thomas Evangelidis [mailto:***@mbg.duth.gr]
Sent: Thu 10/29/2009 5:24 AM
To: Warren DeLano
Cc: pymol-***@lists.sourceforge.net
Subject: Re: [PyMOL] how to get RMSD from align command

Hi Warren,

Now I am more confused. I used align command to measure the overall
RMSD between homologous structures. Apparently align is not
appropriate if not sufficient sequence similarity is present. I later
came across CEalign plugin, which does structure-based
superimposition. This command also returns an RMSD value, which
happens to be higher for the set of homologous structures I analysed.
Based on this I concluded that align command is more appropriate to
measure measure structural similarity, but I'm not sure now.

My objective is to select between a small set of homologous structures
(co-crystallized with ligands), the one with the highest structural
similarity and use it to place a dummy ligand into one where it's
missing. I.e. "1ebh" lacks a ligand, so I compared it with "1ebg",
"1els", "1one", "2one", "5enl", "6enl", "7enl" which have
co-crystallized ligands. According to both align and cealign, 5ENL is
the most similar with overall RMSD 0.246 and 0.420907 respectively.
The question is which command produces the best superimposition
provided that these proteins belong to the enolase family and thus
have high sequence and structural similarity? Based on the RMSD
values, align does, but as you said this value pertains to only a
subset of atoms after refinement.

Apparently pair_fit command cannot be used for my purpose as it
requires identical selections.

Tom

----- End message from ***@delsci.com -----

------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
PyMOL-users mailing list (PyMOL-***@lists.sourceforge.net)
Info Page: https://lists.sourceforge.net/lists/listinfo/pymol-users
Archives: http://www.mail-archive.com/pymol-***@lists.sourceforge.net

Thomas Evangelidis

2009-10-29 23:54:19 UTC

Permalink

Warren,

before I write something incorrect into Methods section, do you think
the overall RMSD returned by align command is a valid criterion to
select which structure of the family is most similar?

Tom

Post by Warren DeLano
Tom,
Apologies for the confusion, but there really isn't any way to
measure overall structural similarity in PyMOL other than pair_fit,
but you have to do the work of explicitly stating which pairs are to
be aligned.
The pair_fit selections don't have to be identical, but they do have
to correspond in a pairwise fashion. In other words, you can
compare any number of C-alpha positions so long as the same number
of atomic positions are specified in each structure.
Cheers,
Warren
-----Original Message-----
Sent: Thu 10/29/2009 5:24 AM
To: Warren DeLano
Subject: Re: [PyMOL] how to get RMSD from align command
Hi Warren,
Now I am more confused. I used align command to measure the overall
RMSD between homologous structures. Apparently align is not
appropriate if not sufficient sequence similarity is present. I later
came across CEalign plugin, which does structure-based
superimposition. This command also returns an RMSD value, which
happens to be higher for the set of homologous structures I analysed.
Based on this I concluded that align command is more appropriate to
measure measure structural similarity, but I'm not sure now.
My objective is to select between a small set of homologous structures
(co-crystallized with ligands), the one with the highest structural
similarity and use it to place a dummy ligand into one where it's
missing. I.e. "1ebh" lacks a ligand, so I compared it with "1ebg",
"1els", "1one", "2one", "5enl", "6enl", "7enl" which have
co-crystallized ligands. According to both align and cealign, 5ENL is
the most similar with overall RMSD 0.246 and 0.420907 respectively.
The question is which command produces the best superimposition
provided that these proteins belong to the enolase family and thus
have high sequence and structural similarity? Based on the RMSD
values, align does, but as you said this value pertains to only a
subset of atoms after refinement.
Apparently pair_fit command cannot be used for my purpose as it
requires identical selections.
Tom

------------------------------------------------------------------------------
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart your
developing skills, take BlackBerry mobile applications to market and stay
ahead of the curve. Join us from November 9 - 12, 2009. Register now!
http://p.sf.net/sfu/devconference
_______________________________________________
Info Page: https://lists.sourceforge.net/lists/listinfo/pymol-users

----- End message from ***@delsci.com -----

Warren DeLano

2009-10-30 00:14:48 UTC

Permalink

Tom,

Because the number of atoms in the alignment can vary with each pairwise
structure comparison, the RMS number returned by align cannot be
considered valid criterion except perhaps if presented along with the
number of atoms in each comparison. E.g.:

RMS over X atoms

0.5 1451
0.3 1273
0.4 1343

But remember, each final RMS value will involve a different set of
atoms...and that's the root of the problem...there's no common reference
point.

Therefore, I must recommand against using 'align' for this purpose --
it's just not the right tool for the job.

Cheers,
Warren

Post by Warren DeLano
-----Original Message-----
Sent: Thursday, October 29, 2009 4:56 PM
To: Warren DeLano
Subject: RE: [PyMOL] how to get RMSD from align command
Warren,
before I write something incorrect into Methods section, do you think
the overall RMSD returned by align command is a valid criterion to
select which structure of the family is most similar?
Tom

later

Post by Warren DeLano

Post by Warren DeLano
came across CEalign plugin, which does structure-based
superimposition. This command also returns an RMSD value, which
happens to be higher for the set of homologous structures I

analysed.

Post by Warren DeLano

Post by Warren DeLano
Based on this I concluded that align command is more appropriate to
measure measure structural similarity, but I'm not sure now.
My objective is to select between a small set of homologous

structures

Post by Warren DeLano

Post by Warren DeLano
(co-crystallized with ligands), the one with the highest structural
similarity and use it to place a dummy ligand into one where it's
missing. I.e. "1ebh" lacks a ligand, so I compared it with "1ebg",
"1els", "1one", "2one", "5enl", "6enl", "7enl" which have
co-crystallized ligands. According to both align and cealign, 5ENL

Post by Warren DeLano

Post by Warren DeLano
the most similar with overall RMSD 0.246 and 0.420907 respectively.
The question is which command produces the best superimposition
provided that these proteins belong to the enolase family and thus
have high sequence and structural similarity? Based on the RMSD
values, align does, but as you said this value pertains to only a
subset of atoms after refinement.
Apparently pair_fit command cannot be used for my purpose as it
requires identical selections.
Tom

find

Post by Warren DeLano

Post by Warren DeLano
how can I get the RMSD value from the align command in a python
script? cmd.align() returnes a tuple of 8 numbers and none of them

Post by Warren DeLano

Post by Warren DeLano
the actual RMSD value I get when I align these 2 structures

manually.

Post by Warren DeLano

Post by Warren DeLano
thanks,Tom

-----------------------------------------------------------------------

Post by Warren DeLano
-------

Post by Warren DeLano

Post by Warren DeLano
Come build with us! The BlackBerry(R) Developer Conference in SF,

Post by Warren DeLano

Post by Warren DeLano
is the only developer event you need to attend this year. Jumpstart

your

Post by Warren DeLano

Post by Warren DeLano
developing skills, take BlackBerry mobile applications to market

and

Post by Warren DeLano
stay

Post by Warren DeLano

Post by Warren DeLano
ahead of the curve. Join us from November 9 - 12, 2009. Register

now!

Post by Warren DeLano

Post by Warren DeLano
http://p.sf.net/sfu/devconference
_______________________________________________
Info Page: https://lists.sourceforge.net/lists/listinfo/pymol-users

------------------------------------------------------------------------

Post by Warren DeLano
------

Post by Warren DeLano
Come build with us! The BlackBerry(R) Developer Conference in SF, CA
is the only developer event you need to attend this year. Jumpstart

your

Post by Warren DeLano

Post by Warren DeLano
developing skills, take BlackBerry mobile applications to market and

stay

Post by Warren DeLano
ahead of the curve. Join us from November 9 - 12, 2009. Register

now!

Post by Warren DeLano

Post by Warren DeLano
http://p.sf.net/sfu/devconference
_______________________________________________
Info Page: https://lists.sourceforge.net/lists/listinfo/pymol-users

Jason Vertrees

2009-10-30 21:50:39 UTC

Permalink

Message: 6
Date: Thu, 29 Oct 2009 13:40:56 +0200
Subject: Re: [PyMOL] how to get RMSD from CEalign command
Content-Type: text/plain; charset=ISO-8859-7; DelSp="Yes";
format="flowed"
Thank you Adnreas.
I'm straggling with cealign now. When I run it in batch mode I get
"Selector-Error: Invalid selection name" for certain structures, which
in fact are very close homologs, whereas the number of RMSD values I
query_template_chains = {
"1ebh" : ["1ebg", "1els", "1one", "2one",
"5enl", "6enl", "7enl"]
}
cmd.fetch(query)
cmd.fetch(template)
print "Superimposing ", query, " onto ", template
cmd.do("cealign "+query+" and c. A, "+template+" and c. A")
cmd.delete(query)
cmd.delete(template)
Do you know what's wrong?
Thomas

Thomas,

The PyMOL UI runs asynchronously. I think you are deleting your PDB
files before they can be aligned. I made a couple changes in your
script. Instead of calling cmd.do just call cealign directly. See
below,

query_template_chains = {
"1ebh" : ["1ebg", "1els", "1one", "2one",
"5enl", "6enl", "7enl"]
}

for query in query_template_chains.keys():
for template in query_template_chains[query]:
cmd.fetch(query, async=0)
cmd.fetch(template, async=0)
print "Superimposing ", query, " onto ", template
cealign( query+" and c. A", template+" and c. A")
cmd.delete(query)
cmd.delete(template)

Hope this helps,

-- Jason

Jason Vertrees, PhD

PyMOLWiki -- http://www.pymolwiki.org

Thomas Evangelidis

2009-10-31 14:11:29 UTC

Permalink

Post by Andreas Forster
Thomas,
The PyMOL UI runs asynchronously. I think you are deleting your PDB
files before they can be aligned. I made a couple changes in your
script. Instead of calling cmd.do just call cealign directly. See
below,
query_template_chains = {
"1ebh" : ["1ebg", "1els", "1one", "2one",
"5enl", "6enl", "7enl"]
}
cmd.fetch(query, async=0)
cmd.fetch(template, async=0)
print "Superimposing ", query, " onto ", template
cealign( query+" and c. A", template+" and c. A")
cmd.delete(query)
cmd.delete(template)
Hope this helps,
-- Jason
Jason Vertrees, PhD

Jason,

I can't run cealign() from a script:

NameError: global name 'cealign' is not defined

I installed pymol from Fedora repositories and can't find .pymolrc,
therefore I load CEalig in my script like this:

cmd.do("run /home/thomas/Documents/cealign-0.9/cealign.py") # load
CEalign plugin
cmd.do("run /home/thomas/Documents/cealign-0.9/qkabsch.py")

The only way to run CEalign is using cmd.do().

Thomas

And of

Jason Vertrees

2009-11-01 04:23:16 UTC

Permalink

Thomas,

The command

cmd.do(".../cealign.py")

runs cealign.py, which runs cmd.extend. So you should have the cealign
command in your namespace after the cmd.do. Have you tried running

cmd.do("your/path/to/cealign")
cealign(protA, protB)?

I just tried this on my system and it worked for me. Let me know if it
works for you.

-- J

Jason Vertrees, PhD

PyMOLWiki -- http://www.pymolwiki.org

Post by Thomas Evangelidis

Jason,
NameError: global name 'cealign' is not defined
I installed pymol from Fedora repositories and can't find .pymolrc,
cmd.do("run /home/thomas/Documents/cealign-0.9/cealign.py") # load CEalign
plugin
cmd.do("run /home/thomas/Documents/cealign-0.9/qkabsch.py")
The only way to run CEalign is using cmd.do().
Thomas
And of

Thomas Evangelidis

2009-11-01 22:41:02 UTC

Permalink

No it doesn't, both in Fedora 11 and Ubuntu 9.04. I get:

NameError: global name 'cealign' is not defined

Thomas

Post by Andreas Forster
Thomas,
The command
cmd.do(".../cealign.py")
runs cealign.py, which runs cmd.extend. So you should have the cealign
command in your namespace after the cmd.do. Have you tried running
cmd.do("your/path/to/cealign")
cealign(protA, protB)?
I just tried this on my system and it worked for me. Let me know if it
works for you.
-- J
Jason Vertrees, PhD
PyMOLWiki -- http://www.pymolwiki.org

Post by Thomas Evangelidis

----- End message from ***@gmail.com -----