Learn how to work with a part of a map. In particular, this involves keeping track of coordinate shifts and symmetry.
It is often convenient to work with just a part of a map. For example, if your structure contains identical 24 copies of a molecule, you might want to work with just one copy and then duplicate it to create the other 23 later. In this section we are going to create a new map that represents part of the map_data map we have been working with. For this example, we are going to just look at the raw boxed map that results so that we can see how it works. In a real situation you would not only cut out the box but shift its origin using the map_model_manager (see further below).
Cutting out part of a map is referred to here as boxing. The key step in cutting out a part of a map is specifying the part of the map to cut out. The grid point with the smallest indices in each direction in the cut-out map is called the “lower bounds” and the grid point with the highest indices is the “upper bounds”. If the lower_bounds is (la,lb,lc) and the upper_bounds is (ua,ub,uc) then the resulting box will have (1+ua-la, 1+ub-lb, 1+uc-lc) grid points along a, b and c. This is the value of all for the boxed map. The box will start at lower_bounds, which will be the value of origin for the boxed map.
First, lets set up a map:
from iotbx.map_model_manager import map_model_manager # load map_model_manager
mmm=map_model_manager() # get initialized instance of the map_model_manager
mmm.generate_map() # get a model and calculate a map for it
map_data = mmm.map_data() # our 3D map object
Let’s set some lower and upper bounds:
lower_bounds = (10,10,10) # lower bounds for boxed map
upper_bounds = (21,31,21) # upper bounds for boxed map
Now let’s cut out a part of our map. In this example we are going to use maptbx.copy() to do this. Normally you will instead use the map_model_manager but using maptbx.copy() shows exactly what is happening:
from cctbx import maptbx # import maptbx
box_map_data = maptbx.copy(map_data, lower_bounds, upper_bounds) # box the map
Let’s look at the origin and dimensions of this boxed map. We expect the origin to be at our lower_bounds (10,10,10) and the dimensions to be (12, 22, 12), and the last available point is (21,31,21):
print( box_map_data.origin()) # prints (10, 10, 10)
print( box_map_data.all()) # prints (12, 22, 12)
print( box_map_data.last(False)) # prints (21, 31, 21)
This boxed map is defined for any grid points from lower_bounds to upper_bounds, inclusive. We can get the value at grid point (11,12,13) just like we did before and we will get the same answer:
print( box_map_data[11,12,13]) # prints 0.0416163499881
The boxed map is not defined at grid point (0,0,0) or any other points outside of the bounds. The boxed map is also a different size than the original (which had a size of 38400):
print( box_map_data.size() ) # prints 3168
One final point: the boxed map is a copy, not a pointer to part of the original map. So if you change a value in your boxed map, the original map is unchanged.
Notice that indices of grid points in the boxed map we just created are still in the same reference frame as our original map (i.e., the grid point (11, 12, 13) has the same value in the original map as it has in the boxed map. However the origin of the boxed map (10,10,10) is different from the original origin (0,0,0).
In cctbx, almost all work with maps is carried out after shifting the origin of the map to (0,0,0). If there is a model associated with a map, the coordinates of atoms in the model are also shifted. Normally this is done for you with the map_model_manager. In this section we are going to do an origin shift directly so that we can see how it works.
Let’s shift the origin of our boxed map (box_map_data) to (0,0,0) with the shift_origin() function to see what happens:
shifted_box_map_data = box_map_data.shift_origin() # shift origin to (0,0,0)
Note that the shift_origin() function does not change the original object but rather returns a new object with shifted origin but the same data. Note also that if you change data in the shifted map you are changing the data in the unshifted map as well.
Let’s look at the origin and dimensions of our boxed map after shifting the origin. The size should not change but the origin should be at (0,0,0) and the last available point should be (11,21,11) instead of (21,31,21):
print(shifted_box_map_data.origin()) # prints (0, 0, 0)
print(shifted_box_map_data.all()) # prints (12, 22, 12)
print(shifted_box_map_data.last(False)) # prints (11, 21, 11)
This shifted version of our boxed map data is defined from (0,0,0) to (11,21,11). There has been an origin shift of (10,10,10). If we want to see the value at the grid point that used to be called (11,12,13), we now have to look at the grid point (1,2,3):
print(shifted_box_map_data[1,2,3]) # prints 0.0416163499881
The dimensions of our boxed map (12,22,12) in grid units and in Angstrom units are different from our original map. That means that the crystal_symmetry of our boxed map has changed. We need a way to keep track of several things:
Here is the convention for keeping track of origins, gridding and dimensions:
The number of grid points and crystal_symmetry corresponding to the full original map (the full map from cryo-EM or the crystallographic unit cell map for crystallography) are referred to as the unit_cell_grid and unit_cell_crystal_symmetry. Here unit_cell is intended to mean the full unit cell of the map.
The origin shift between original and boxed maps is recorded in grid units as origin_shift_grid_units, where the origin shift is the shift needed to move the origin back to its original position. That means that if the origin of the boxed map was at (10,10,10) and it was moved to (0,0,0), the origin shift was (10,10,10).
At the same time the origin shift is recorded as a shift in Cartesian coordinates called shift_cart. The value of shift_cart is the shift that has been applied to any coordinates. It has the opposite sign of origin_shift_grid_units.
The crystal_symmetry of the boxed map is referred to as crystal_symmetry. This is the crystal_symmetry that is used for any calculations that are done with a boxed map and model.
Now that we have seen how boxing and shifting origins works, we can appreciate how the map_model_manager (mmm) has methods to take care of everything for us. Let’s start from the beginning and box our original map and the model that goes with it using the map_model_manager. Let’s use the same lower and upper bounds we did before and create a new map_model_manager that is boxed and origin-shifted, but that knows all about how to put everything back where it was:
boxed_mmm = mmm.extract_all_maps_with_bounds( # create box
lower_bounds = lower_bounds, # lower bounds
upper_bounds = upper_bounds) # upper bounds
If we look at the map_data object from this boxed map_model_manager we will find that it looks just like the boxed, origin-shifted one we created a moment ago:
new_shifted_box_map_data = boxed_mmm.map_manager().map_data() #
print(new_shifted_box_map_data.origin()) # prints (0, 0, 0)
print(new_shifted_box_map_data.all()) # prints (12, 22, 12)
print(new_shifted_box_map_data.last(False)) # prints (11, 21, 11)
print(new_shifted_box_map_data[1,2,3]) # prints 0.0416163499881
However the map_model_manager knows how to write out this boxed map so that it superimposes on the original:
boxed_mmm.write_map('boxed_map.ccp4') # superimposes on orig
Similarly, if you write out your map with the DataManager it will superimpose on the original.
The map_model_manager shifts the origin of your map and it also shifts the origin of any model that is part of the map_model_manager. That means that the model in our new map_model_manager (boxed_mmm) matches the map in that map_model_manager. The way this happens is that the coordinates of atoms in the model were all shifted by the value of shift_cart() (calculated from the origin_shift_grid_units) that was set when the origin of the map was shifted. When you write out your model, the reverse shift is automatically applied by the DataManager.
When your working map has been extracted from a larger map with a map_model_manager, the manager will shift the origin of all maps, models, and associated reconstruction symmetry (ncs_spec) objects to (0,0,0) before working with them.
Then you do work with the origin at (0,0,0). You automatically get coordinates from your models that are relative to this origin if the model came from the map_model_manager where the origin was automatically shifted. For example you might say:
working_sites_cart = boxed_mmm.model().get_sites_cart() # sites
Note that when you are working with your boxed model, the values of coordinates that you get from your model with the function get_sites_cart() will not be the coordinates in the original reference frame of your model. They were shifted to match your boxed, shifted map.
Finally when writing out maps, models, or ncs objects, the shift origin is shifted back. This shift back is automatic so you do not need to do anything.
If you have created amap_model_manager with one map and model, and you have a new map that you want to match, just shift the origin of the new map with shift_origin(). If the maps are compatible, they will now match. You can check that with the is_similar() function.
If for some reason you just create a map_manager and want to work with it and do not want to combine it with a model, you should shift the origin of your map_manager to (0,0,0) using the shift_origin() method before working with it.
The map_model_manager has two sets of methods for creating boxed versions of itself. One set extracts a new copy that is boxed, and the other changes itself by boxing. These two methods have different properties. The extracted version is a completely new object. In contrast, boxing a map_model_manager in place changes the map_model_manager itself. Its map is boxed,creating a new (differently-sized) map, but its model is changed in place, with new, shifted coordinates but keeping the same model object.
To extract a boxed version of a map_model_manager, you use a method such as the extract_all_maps_with_bounds() method that we used above.
boxed_mmm = mmm.extract_all_maps_with_bounds( # create box
lower_bounds = lower_bounds, # lower bounds
upper_bounds = upper_bounds) # upper bounds
The extracted map and model are deep copies, not pointers to the originals. This means that if you change the extracted map or model, you do not change the original. Let’s do this just to see. Here are the current coordinates of a site in the model and a position in the map in the map_model_manager:
print (mmm.model().get_sites_cart()[0]) # 14.476000000000003, 10.57, 8.342)
print (mmm.map_manager().map_data()[11,12,13]) # prints 0.0416163499881
The coordinates for the extracted map_model_manager are shifted:
print(boxed_mmm.model().get_sites_cart()[0])# (7.005666666666668, 3.339250000000002, 0.967625000000001)
The indices of the map in the extracted map_model_manager are offset by (10,10,10) so now the index (1,2,3) corresponds to the index (11,12,13) in the original map:
print (boxed_mmm.map_manager().map_data()[1,2,3]) # prints 0.0416163499881
Now we change the value of a coordinate or a value in the map for our extracted map_model_manager. (Note that the procedure for changing coordinates in a model has three steps: get the sites_cart, change something in sites_cart, set the sites cart in the model with the changed sites_cart):
boxed_sites_cart = boxed_mmm.model().get_sites_cart() # get boxed sites
boxed_sites_cart[0] = (10,10,10) # set value of one coordinate in boxed sites
boxed_mmm.model().set_sites_cart(boxed_sites_cart) # set coordinates in model
boxed_mmm.map_manager().map_data()[1,2,3] = 77. # change map value
Nothing happens to the original because the extracted map_model_manager was a deep copy:
print (mmm.model().get_sites_cart()[0]) # 14.476000000000003, 10.57, 8.342)
print (mmm.map_manager().map_data()[11,12,13]) # prints 0.0416163499881
On the other hand, to box a map_model_manager in place, you use a method such as the box_all_maps_with_bounds_and_shift_origin method. Let’s make references to the model and map in our map_model_manager so that we can track what happens to them when we box the map_model_manager in place:
mmm_model_ref = mmm.model() # reference to model in mmm
mmm_map_manager_ref = mmm.map_manager() # reference to model in mmm
Now let’s instead box the map_model_manager in place. This is going to produce a new map but will alter the model and it will not create a new model:
mmm.box_all_maps_with_bounds_and_shift_origin( # change mmm in place
lower_bounds = lower_bounds, # lower bounds
upper_bounds = upper_bounds) # upper bounds
The map_model_manager is changed in place( nothing is returned). The new map will be a new copy (its size has changed and changing it does not change the original), but the new model is the same model as the original, with shifted coordinates. If you change the model in this map_model_manager, any references to the original model also change.
Let’s do this just to see. The current coordinates of the first site in the model are now shifted:
print (mmm.model().get_sites_cart()[0]) # (7.005666666666668, 3.339250000000002, 0.967625000000001)
and the indices of the map in the map_model_manager are offset by (10,10,10) so now the index (1,2,3) corresponds to the index (11,12,13) in the original map:
print (mmm.map_manager().map_data()[1,2,3]) # prints 0.0416163499881
Now we change the value of a coordinate or a value in the map for our (now-shifted) map_model_manager:
sites_cart = mmm.model().get_sites_cart() # get boxed sites
sites_cart[0] = (20,20,20) # set value of one coordinate in sites_cart
mmm.model().set_sites_cart(sites_cart) # set coordinates in model
mmm.map_manager().map_data()[1,2,3] = 222. # change map value
The value in our original model now changes because the model in our map_model_manager was shifted in place:
print (mmm_model_ref.get_sites_cart()[0]) # (20.0, 20.000000000000004, 20.0)
map_model_manager made a fresh copy of the map:
print (mmm_map_manager_ref.map_data()[11,12,13]) # prints 0.0416163499881
If you have a map_manager and you have a new model that you want to match, use the function shift_model_to_match_map(new_model) which will shift the new model (in place, changing new_model) to match this map.
shift_any_model_to_match and when to use set_model_symmetries_and_shift_cart_to_match_map When you work with a model and a map_model_manager, you may want to adjust the origin for the model (the reference frame that it is in) to match the origin of the map_model_manager.
If you load the model into the map_model_manager, this will happen automatically. Suppose we take the model from our boxed map_model_manager:
m = boxed_mmm.model() # get the model
print(m) # prints info about the model including origin shift
Now let's shift this model back to its original location (as we would have when we read in a model with the DataManager):
shift_cart = m.shift_cart() # current origin shift
shift_to_apply = tuple([-x for x in shift_cart]) # opposite shift (to apply()
m.shift_model_and_set_crystal_symmetry(shift_to_apply) # shift the model
print(m) # now the origin shift is zero (original location)
Let's put this model into our boxed map_model_manager:
boxed_mmm.add_model_by_id(model_id = 'model', model = m) # load the model in
print(m) # automatically shifted to match the map_model_manager origin
We could have shifted the model without loading it into the map_model_manager if we wanted to. We use the function:shift_any_model_to_match:
First let's create a model with an origin shift of zero:
m.shift_model_and_set_crystal_symmetry(shift_to_apply) # shift the model again
print(m) # now the origin shift is zero (original location)
Now let's use shift_any_model_to_match to shift it to match the
boxed_mmm map_model_manager:
boxed_mmm.shift_any_model_to_match(m) # shift this model to match the map_model_manager
print(m) # automatically shifted to match the map_model_manager origin
There are some situations where we we want to keep the coordinates in a model exactly the same, but specify that there is an origin shift that matches a particular map_manager. This could happen if we use the map_data from a map_manager to create some coordinates (relative to the origin of that map) and create a model from these coordinates. Now that model has an origin of (0,0,0) but the map_manager may have some other origin. We just want to say that the model refers to the same origin as the map_manager. We can use the set_model_symmetries_and_shift_cart_to_match_map method to do this:
mm = boxed_mmm.map_manager() # get a map manager
mm.set_model_symmetries_and_shift_cart_to_match_map(m) # set the symmetry and origin
Keep in mind that the set_model_symmetries_and_shift_cart_to_match_map function should rarely be used. Most of the time you want to shift the origin and the coordinates at the same time with shift_any_model_to_match
unit cell. On the other hand, cryo-EM maps are only defined inside the region occupied by the original map. The parameter called wrapping is used to specify whether values at a grid point outside the boundaries of a map are to be calculated by wrapping the grid point back into the original map using unit cell translations. In general, crystallographic maps that are full size can be wrapped in this way, but cryo-EM maps and any parts of crystallographic maps cannot. Wrapping is normally guessed from the way a map is created, but it can be specified as part of the initialization of map_manager objects and map_model_manager objects.
The map_model_manager and tools that check to see whether the map and model are compatible by comparing their crystal_symmetry values, and if present, unit_cell_crystal_symmetry values. The crystal_symmetry values reflect the dimensions of the boxed map or model, and the unit_cell_crystal_symmetry values (if present) reflect the dimensions of the original map and model unit cells.
In most cases, as long as the map_model_manager has been used to keep track of origin shifts and unit_cell_crystal_symmetry, a map and model that have gone through the same boxing steps will be compatible and nothing further needs to be done.
However in some cases the crystal_symmetry of a model (for example) might be changed outside the map_model_manager, in which case an attempt to combine that model with a map that is otherwise compatible would fail with an error message that says that the symmetry of a map does not match that of a model and that they therefore cannot be combined successfully.
If the model coordinates really are compatible with the map and the crystal symmetry of the model is incorrect, the map_model_manager can be instructed to ignore the symmetry of the model with the ignore_symmetry_conflicts=True keyword.
If you have a map_manager and a model that is in the same frame of reference as the map_manager but for some reason has a different crystal_symmetry or unit_cell_crystal_symmetry, you can match the symmetry of the map_manager (without changing any coordinates in the model) with the map_manager member function set_model_symmetries_and_shift_cart_to_match_map().