[MMTK] Comparing normal modes for different proteins

[Konrad Hinsen]

On 18.03.2007, at 21:32, Ramon Crehuet wrote:

>    I'd like to compare the normal modes of a set of proteins which are
> quite similar (in RMS and sequence identity). I know this can be
> tricky because it involves comparing things of different universes (I
> can't see a way of putting all the proteins in the same universe).

That is not difficult (just put them at a large enough distance that  
there are no interactions), but it wouldn't really help in your  
situation.

>    I need to compare residue-by-residue so that, for a given mode, I'd
> like to know to which residue and which atom does mode[i] refer to. I
> know the usual way to do things is the opposite, letting MMTK do the
> transformation as in positions = universe.configuration().

Indeed the recommended way is to index ParticleVector objects (of  
which Mode objects are a subclass) with atom objects, not numbers.

The link between the two is quite simple:  mode[atom] is equivalent  
to mode[atom.index]. So if you want to know which atom mode[i] refers  
to, you can just make up an appropriately ordered list of atoms:

	atoms = universe.numberOfAtoms*[None]
	for atom in universe.atomList():
		atoms[atom.index] = atom

The reason why this is not recommended is that atom numbering is not  
a documented feature of the MMTP API. I make no promises as to the  
conditions under which atom.index is available and valid, nor do I  
promise that there will be an attribute atom.index in future MMTK  
releases. In all releases up to now, you can count on atom.index  
being valid after an energy calculation.

However, I have yet to see a situation in which referring to  
atom.index simplifies any problem. For your case, I would create a  
map between the atoms of the two universes. For example, this piece  
of code maps the C-alphas on each other, assuming sequences of  
identical length:

	atom_map_1_2 = {}
	atom_map_2_1 = {}
	for i_chain in range(len(protein_1)):
		chain = protein_1[i_chain]
		for i_residue in range(len(chain)):
			residue = chain[i_residue]
			ca_1 = residue.peptide.C_alpha
			ca_2 = protein_2[i_chain][i_residue].peptide.C_alpha
			atom_map_1_2[ca_1] = ca_2
			atom_map_2_1[ca_2] = ca_1

Then mode_1[atom] corresponds to mode_2[atom_map_1_2[atom]], which is  
all you need to know.

Konrad.
--
---------------------------------------------------------------------
Konrad Hinsen
Centre de Biophysique Moléculaire, CNRS Orléans
Synchrotron Soleil - Division Expériences
Saint Aubin - BP 48
91192 Gif sur Yvette Cedex, France
Tel. +33-1 69 35 97 15
E-Mail: hinsen at cnrs-orleans.fr
---------------------------------------------------------------------