Averaging

In this practical, we will try to make use of the non-crystallographic symmetry (NCS) in the P2 myelin structure. We will use the experimental electron density that was based on 2 SIRAS identical atom sites.

For averaging you need a map (obviously), a set of operators that describe how one of the NCS-units is related to the others, and a volume that defines the NCS-unit (the ‘mask’). The mask and map are usually on the same grid.

As before, make a subdirectory (if it does not exist) with the name of the Linux box you are using. If the box is #3, there should be a subdirectory of this name so do the following

% cd 3
% mkdir averaging
% cd averaging
% cp /home/alwyn/o/averaging/* .
% ono
Mono enabled only
Gamepad disabled
O > This version of O is compiled with the Absoft Co. compiler
O > To allow runtime distribution, I am required to state:
O > "Portions if this program include material copyrighted ? by
O > Absoft Corporation 1988-1998"
O > Use of this program, implies acceptance of the present
O > set of conditions of use.
O > Contact alwyn@xray.bmc.uu.se if in doubt.
O > O version 9.0.7 , Build 040316
O > Define an O file (terminate with blank):
O > Menu names are not defined.
O > Enter file name [/Users/alwyn/o/data/menu.odb]:
O > Startup file was never loaded
O > Enter file name [/Users/alwyn/o/data/startup.odb]:
O > Stereochemistry file was never loaded
O > Enter file name [/Users/alwyn/o/data/stereo_chem.odb]:
O > Connectivity used is : all
O > Maximum linkage distance = 2.00
O > There were 47 residues.
O > 329 atoms.
O > ODAT : /Users/alwyn/o/data/
O > Making visibility data structures.
O > Making visibility data structures.
Fm>
Fm> Map type is DSN6
Fm>
Fm> Parameters as read from the map file:
Fm> Grid ................. 76 82 48
Fm> Origin ............... 0 0 0
Fm> Extent ............... 77 83 49
Fm> Fast, medium, slow.... X Y Z
Fm> Cell axes ............ 91.80 99.50 56.50
Fm> Cell angles .......... 90.00 90.00 90.00
Fm> No reslicing of map necessary
Fm> Prod ................. 20.00
Fm> Plus ................. 0
Fm> Min, max, sigma ...... 0.00000 8.35000 0.95545
Fm> Scale ................ 1.000
Fm> Symmetry information may be used
Fm> O-style symm-ops for spacegroup P212121
Fm> 4 symmetry operators
New> Number of mouse buttons 1


The electron density map ano1.map has been read into O and is available under the Density slider as map MIR. From the above output, we can see that the separation between grid points is ~91.8/76 = 1.2 Å. This is a bit coarse for averaging.

We will start off with the same skeleton as in the p2map practical, so type @p2_start and remind yourself how many molecules there are in the asymmetric unit. I call them A, B and C. Then click on @p2_1 to generate a main-chain skeleton of the ‘A’ molecule. This will be the molecule from which we will try to derive the transformation operators to take us to the NCS related molecules.

The commands that are related to averaging all start with NCS_*
> ncs
As2> NCS is not a unique keyword.
As2> NCS_mask_mak is a possibility.
As2> NCS_mask_wri is a possibility.
As2> NCS_mask_edi is a possibility.
As2> NCS_operator is a possibility.
As2> NCS_average is a possibility.
As2> NCS_mask_sph is a possibility.
As2> NCS_mask_lay is a possibility.
As2> NCS_map_writ is a possibility.
As2> NCS is not a visible command.



Initial operators
First, we need to get the NCS-operators, or at least approximations to them. If we have molecular models, this is easy to do with the Lsq_* commands in O. If we have multiple heavy atom binding sites per NCS-molecule and if we could decide on which are related to which, we could also use the Lsq_* commands. We don’t have this kind of information, so we’ll do it another way use NCS_operator command. For this we need a skeleton object of the NCS origin molecule and then a rough idea where the NCS related molecule is positioned.

> ncs_op
New> What map? [MIR]:
New> Object name [A]:
New> Step search width [3]:
New> Number of atoms in object 369
New> Belongs to molecule START
New> Pivot coordinates in object [ 49.58 62.36 32.40]:
New> This takes some time ...
New> Starting 1 of 3
New> Starting 2 of 3
New> Starting 3 of 3
New> Best score: 1.161


After the fourth prompt, you need to ID an atom roughly in the centre of what I call the ‘B’ molecule. A translational grid search is carried out around this identified atom, and the atoms in the skeleton object are pivoted in a full 3-D search. The default search width is +/3 steps, each 1 Å wide around the ID’d atom. The score is the average density at each point in the skeleton for the best fitting translation/orientation.

The command has generated a new entry in the O DB with the operator to take you from the A to the B molecule:

> dir .lsq_*rt*
Heap> .LSQ_RT_IDENT R W 12
Heap> .LSQ_RT_BEST R W 12


The first operator is the identity operation, taking A to itself. Now we can see how good the A molecule fits on the B molecule using the Lsq_obj command. This command takes a rotation-translation operator and applies it to an atomic object. The coordinates in the object are not affected by this command.

> lsq
Heap> LSQ is not a unique keyword.
Heap> Lsq_explicit is a possibility.
Heap> Lsq_improve is a possibility.
Heap> Lsq_object is a possibility.
Heap> Lsq_molecule is a possibility.
Heap> Lsq_Paired_a is a possibility.
Heap> LSQ is not a visible command.
> lsq_obj
Lsq > Apply a transformation to an existing object.
Lsq > There are these transformations in the database
Lsq > IDENT BEST
Lsq > Which alignment [<CR>=restore a transformed object] ? BEST
Lsq > There is an object called P2ALL
Lsq > There is an object called A
Lsq > Which object ? a

If this is good enough, you can jump forward to the ‘finding the operator to the C molecule’ section, after you write out the result to a file atob.odb

write .LSQ_RT_BEST atob.odb ;

Lets first try increasing the search range:

> ncs_op
New> What map? [MIR]:
New> Object name [A]:
New> Step search width [3]: 5
New> Number of atoms in object 369
New> Belongs to molecule START
New> Pivot coordinates in object [ 49.58 62.36 32.40]:
New> This takes some time ...
New> Starting 1 of 5
New> Starting 2 of 5
New> Starting 3 of 5
New> Starting 4 of 5
New> Starting 5 of 5
New> Best score: 2.066
> lsq_obj best a


If it still looks very bad, try identifying a different atom as the pivot point for the search. In the above example, I got a much better score and the overlaid objects looked good, so I write them out with:

write .LSQ_RT_BEST atob.odb ;

If you still have problems at this stage, centre on (77., 50., 21.) and identify atom Bones atom 4764 as the pivot point.

Now we have to repeat the same procedure for the third NCS molecule, C. This is one of the pair that is sitting on a crystallographic 2-fold screw axis along the Z-axis. Best of luck with ID’ing an atom, here is my attempt:

> ncs_op
New> What map? [MIR]:
New> Object name [A]:
New> Step search width [3]:
New> Number of atoms in object 369
New> Belongs to molecule START
New> Pivot coordinates in object [ 49.58 62.36 32.40]:
New> This takes some time ...
New> Starting 1 of 3
New> Starting 2 of 3
New> Starting 3 of 3
New> Best score: 1.604
lsq_obj best a


If at first you don’t succeed, try, try again. Bones 2796 is a good place to start.
Write out your result:

write .LSQ_RT_BEST atoc.odb ;


Spherical mask

These first operators may not be good enough to use for averaging. To improve them, we need a spherical mask that is centred on the NCS-asymmetric unit molecule that has a sufficiently large radius to cover the whole of this molecule. To make it, we’ll first make a small skeleton object that is positioned at the centre of the A-molecule, so centre there SVP.

> b_set
Bone> Skeleton molecule name [START ]:
Bone> Object name [A ]: 1
Bone> Sphere radius [ 22.0 ]: 2
Bone> What bone levels? [ 1 3 ]:
> b_dr


That makes an object call ‘1’ that has just a few atoms in it. Now we’ll place large spheres at each of these atoms to make the simple mask, and write it out to the file system:

> NCS_mask_sph
New> What map (for cell, grid units only)? [MIR]:
New> Object name [1]:
New> Atomic radii [3.0]: 22.
New> Mask name [MASK]: sphere
New> Number of atoms in object 3
New> Belongs to molecule START
> NCS_mask_wri
Mask> What map? [MIR]: sphere
Mask> Mask file name? [sphere.mask]:

I have used a radius of 22 Å around each of the skeleton atoms; it should be big enough. If you use too large a radius, you’ll pull in density from abutting molecules. If you use too little a radius, you’ll reduce the signal to noise in the next step where we improve the operators.


Improving the operators.

This step is done outside of O, so fire up a new terminal, set the directory to the one we are working in and activate the program we will use that is called ‘imp’. This program will prompt for the map (ano1.E), the mask (sphere.mask if you used my name), a file containing the crystal symmetry (p212121.sym in this case), and then a file containing the NCS operator (atob.odb or atoc.odb if you use my names). The program gives you a number of options for improving the operator:

Select one of the following options:
Q = save current operator and quit
L = list settings
N = save current operator and read in a new one
T = do a 3D translation search
R = do a 3D rotation search
A = do an automatic R/T-search
P = do a pseudo-6D R/T-search
6 = do a complete 6D R/T-search


The easiest one (and best) is ‘A’. Just take the defaults. The program then runs until it finds a local maximum for what it calls a correlation coefficient. This is the CC of the map points in the mask containing the A molecule, with the same points after the application of the NCS-operator. The latter operation does not usually transform a map point in A onto another point in B, say. Instead, the operation is to a non-integral point in B. Therefore, an interpolation has to carried out.

You need to do this for both the B and C molecules, and maybe you have a few different starting estimates for each operator. Here is what happened with my operators:
A-> B imp
I made one operator with the default step ise of 3 in NCS_operator in O. The correlation coefficient rose from 0.11 to 0.13, i.e. it did not converge to a good solution. The 5-step operator, however, produced a CC of 0.52 and was improved to 0.53

A->C imp
The 3-step operator that I generated started off with a CC of 0.39 and was improved to 0.45.

My operators either worked easily or did not work at all. Therfore, I’ve put another operator from Gerard K. that is not very good, but will converge to something much better. This is in the file atob_wrong.odb try it.

You need to write out your improved operators with the ‘Q’ command in imp. Use names that make some sort of sense, e.g. atob_imp.odb and atoc_imp.odb

We will read these new operators into O but first we have to modify their internal O names. This is what my improved A->B operator file looked like:

.LSQ_RT_BEST R 12 (3f15.8)
0.39147559 0.08517487 0.91624177
0.11143136 0.98398775 -0.13908407
-0.91341609 0.15654117 0.37571535
64.83197784 -34.40064240 -5.91784286


The text .LSQ_RT_BEST is the operator’s name in O. Both the A->B and A-> operators will have these names, so we need to change them to something like

.LSQ_RT_ATOB_imp R 12 (3f15.8)
0.39147559 0.08517487 0.91624177
0.11143136 0.98398775 -0.13908407
-0.91341609 0.15654117 0.37571535
64.83197784 -34.40064240 -5.91784286

with a text editor. Them we can read them into O and use them.


> read atob_imp.odb
> read atoc_imp.odb
> dir .lsq*rt*
Heap> .LSQ_RT_IDENT R W 12
Heap> .LSQ_RT_BEST R W 12
Heap> .LSQ_RT_ATOB_IMP R W 12
Heap> .LSQ_RT_ATOC_IMP R W 12


These operators can now be tested with the A object to see how well they fit. Try it.

lsq_obj ATOB_IMP a

and

lsq_obj ATOC_IMP a

Read in Gerard’s operator to see how far off it was from the B-molecule


Making a Better Mask.

The first step is to use the spherical mask to carry out an averaging step with the improved operators.

> nc_av
New> What map? [MIR]: mir
number SPG ops 4
New> What mask? [?]: sphere
New> New averaged map: av1
New> Define an NCS operator .lsq_rt_ident
New> Define an NCS operator .lsq_rt_ATOC_IMP
New> Define an NCS operator .lsq_rt_ATOB_IMP
New> Define an NCS operator .lsq_rt_
New> Average.
134830 28304 0 0
New> Min, max, sigma ...... -0.30110 5.68663 0.38489

Note that when you are prompted for the operators, you just define the last part of the name! We have now generated a new map (I called it AV1) in O’s FastMap system. You can open a slider for it, and take a look at it. Note that the density only exists within the simple spherical mask. There is also other density from neighbouring molecules. To get rid of this, we need a better mask and we’ll make it by skeletonizing the new map, make an object of main-chain atoms, edit it a bit, and make a new mask. Sounds complicated bit it is NOT. The skeleton step first:

> skel
Qm> What map? [AV1]:
Qm> All map [Y]/N?
Qm> Base level [1.25]: 3.5
Qm> Skeleton name [SKL]: av


I used a level of 3.5 the RMS value of the averaged map to make the skeleton. I called the new skeleton AV. Now use the Bones pull-down to make a main-chain object of it: define the skeleton molecule, the radius, then make the object.

You should see that there are a few bonds from skeleton atoms outside the A-molecule. Break these bonds, and when you are happy make a new mask, carry on:

> ncs_ma_s
New> What map (for cell, grid units only)? [AV1]: mir
New> Object name [MCAV]:
New> Atomic radii [3.0]:
New> Mask name [MASK]: mask1
New> Number of atoms in object 375
New> Belongs to molecule AV


The new mask is called MASK1 in the FastMap system. Admire your handy work. With a 3 Å radius around each atom, the mask is not good enough since it has holes in it. These we will rempove by a series of mask expansions and contractions. Here is what I did:

> NCS_mask_lay
New> What mask ? [MIR]: mask1
New> Expand layer (+), or peel it (-)? [+]: +
> NCS_mask_lay mask1 +
> NCS_mask_lay mask1 -
> NCS_mask_lay mask1 -
> NCS_mask_lay mask1 +


Time to save things, O DB and the mask:

> save
> NCS_mask_wri
Mask> What map? [MASK1]:
Mask> Mask file name? [MASK1.mask]:


As I mentioned earlier, the MIR map was calculated on a ~1.2 Å grid, which is a bit too rough for using with cyclic averaging at the edge resolution of the P2 myelin study (some date to 2.8  Å). Therefore, we’ll transfer this mask to a finer grid, using ‘mama’. Start the program and use these inputs:

re m1 MASK1.mask
new grid 100 110 64
new cell 91.800 99.500 56.500 90.000 90.000 90.000
new unit m2 m1
wr m2 finer.mask


A new mask gets made on a ~0.9 Å grid. Read it into O and look at it (fm_file, no symmetry to be used). How does sit compare to the earlier mask, MASK1?

Before cyclic averaging, we could also check for overlap, but my mask had almost none so we can leave that step out.


Cyclic Averaging

This is done outside O. The steps are:

• Average the map using the current mask
• Expand the averaged NCA-unit to fill the asymmetric unit or cell
• FFT the averaged map
• Use these phases and scaled 2|Fo|-|Fc| amplitudes to generate a new map
• Go back to step 1

The result of 5 such cycles is the av_5.map that we looked at in the p2map practical.