Hi Jan,
As Paul pointed out, you can use COOT to accomplish what you want. Particularly
you can take a look at the following functions in section 6 of the COOT manual
(these are in scheme):
(set-map-mask-atom-radius radius)
(mask-map-by-molecule imol-map imol-model invert-mask?)
(transform-map imol rotation-matrix trans point radius)
(transform-map-using-lsq-matrix imol-ref ref-chain ref-resno-start
ref-resno-end imol-mov mov-chain mov-resno-start mov-resno-end imol-map
about-pt radius)
Please take a look at the attached python example. It can be run with COOT
graphical interface by:
coot scr.py
You can run it without COOT graphical interface by un-commenting the exit()
line at the end, then:
coot --no-graphics -s scr.py
If you see any part of the resulting transformed map being cut in the protein
region, play with the rotation center (supply your coordinates [x,y,z] to
replace the rotation_center() function. In this script, the current rotation
center is set upon reading the last pdb.).
BTW, if you see jagged map after rotation, that's probably because the map
sampling frequency before the rotation is too low. When you rotate a map, there
is always a resolution loss, because of the rastering nature of the maps
(imagine rotating a mesh then resample it to become another mesh of the same
spacing, how can you keep all the information?). So you had better use a very
high sampling frequency when you are preparing a map that is going to be
rotated. Since we often start from a MTZ file, you may even start by computing
the starting map by specifying the sampling frequency to be more than 5x of the
highest resolution to be very safe. On the other hand, if you use COOT i guess
an interpolation is going to be done internally before the rotation (however
interpolation is, maybe to some degree, not as good as computing a very
over-sampled map from MTZ, I think.).
Zhijie
________________________________
From: CCP4 bulletin board <[email protected]> on behalf of Jan Abendroth
<[email protected]>
Sent: Tuesday, March 26, 2019 2:05 AM
To: [email protected]
Subject: Re: [ccp4bb] map rotation
Thanks, Eleanor!
what i really want to do with this script is to compare ligand binding sites
from many structures in several space groups. So, I want to move coordinates
and density onto one reference molecule. Since we ran into issues, I used this
dimer structure as a simple test case.
The script you outlined, only recreates the same density w/o any rotation.
mapmask \
mapin AB_full-cell.ccp4 \
xyzin B.pdb \
mskout B.msk \
<< eof
border 5
eof
maprot \
mapin AB_full-cell.ccp4 \
mskin B.msk \
mapout B_rot.map \
<< eof
MODE to
AVER
rota euler 152.440 110.243 28.112
TRANS -42.212 5.510 -57.243
end
eof
Btw, I had to remove the rota euler 0 0 0 / trans 0 0 0 lines. Maprot
complained about too many operators.
The odd thing -to me- is that when I shift a map around molecule B or the dimer
(mapmask with border 5) with a small amount ( euler 1 0 0 or trans 0.5 0.5
0.5) the map does not look jagged at the edges of the molecule, while it does
when I rotate the full amount to match B on A:
maprot \
wrkin AB-5.map \
mapout AB_rot.map \
<< eof
CELL xtal 61.0100 142.3600 68.2800 90.0000 97.1980 90.0000
GRID xtal 100 228 112
MODE to
AVER
rota euler 1 0 0
TRANS 0 0 0
end
eof
Cheers,
Jan
On Mon, Mar 25, 2019 at 3:49 AM Eleanor Dodson
<[email protected]<mailto:[email protected]>> wrote:
Hmm - I find maprot extremely confusing, but remember a wrkmap does not use any
symmetry so maybe that is why some is lost.
I would have done this, but I havent tested it. And the documentation is
SERIOUSLY confusing!!
Do I understand you want to ADD the density for mol B to that of Mol A
mapmask mapin whole-cell.map
xyzin A.pdb
mskout A-mol.msk
Then
maprot
mapin whole-cell.map
mskin A-mol.msk ! only density withing this mask is
interesting
mapout A-mol.map
MODE to
AVER
rota euler 0 0 0 ! to pick up mol A density
trans 0 0 0
rota euler 152.440 110.243 28.112 ! to rotate map by B to A rotation
TRANS -42.212 5.510 -57.243
end
On Mon, 25 Mar 2019 at 04:58, Jan Abendroth
<[email protected]<mailto:[email protected]>> wrote:
Hi all,
thanks for the feedback. Suggestions like coot or pymol won't work for us well,
since we will have to do this with dozens of structures/maps.So, I'd rather
have this scripted.
Still running into some issues that I think relate to maprot.
My understanding is that I first have to create a map covering molecule B that
I want to map on A. Checking the extend of the map in chimera confirms that
this worked:
mapmask \
mapin 2mol_2mFo-DFc.map \
xyzin 2mol_B.pdb \
mapout 2mol_2mFo-DFc_B.map \
<< eof
border 5
eof
Next, I need to rotate/translate the map in maprot. Since in maprot, mapin
requires a map that covers the unit cell, I use wrkin and 'mode to' as below.
In this script, the cell and grid values are the same mapdump provides me for
the map. The rotation and translation are from superpose, RMSD of that
superposition is 0.5Å.
maprot \
wrkin 2mol_2mFo-DFc_B.map \
mapout 2mol_2FoFc_rot.map \
<< eof
CELL xtal 61.0100 142.3600 68.2800 90.0000 97.1980 90.0000
GRID xtal 100 228 112
MODE to
AVER
rota euler 152.440 110.243 28.112
TRANS -42.212 5.510 -57.243
eof
The issue now is that the superposed map for the center of molecule A looks
great. Towards the edges of the molecule it gets weaker, does not match up with
the molecule or stops entirely. Again, molecule and maps between A and B, as
visualized in Coot by NCS hopping, are very similar.
I am still quite puzzled by what is happening. I guess I am missing something
in maprot. Any input would be appreciated. This is public data, so I would be
happy to share the data.
Cheers,
Jan
On Wed, Mar 20, 2019 at 9:26 PM Jan Abendroth
<[email protected]<mailto:[email protected]>> wrote:
Hi all,
this should be easy, scripting the rotation of a map.
Purpose for this is: Superimpose several structures of the same protein that
crystallized in different space groups, and then drag the maps along.
As a simple test, I took a dimeric protein and try to superimpose molecule B
along with the map on molecule A.
The execution should be straightforward:
a) take a map that covers the unit cell (fft),
b) generate a mask around molecule B (mapmask),
c) apply rotation/translation that I obtain from superimposing molecule B on
molecule A.
The issue is that the obtained map covers both molecule A and B (not a big
deal), more importantly, it cuts of certain areas on both molecules. Molecule A
and B have low RMSDs (0.5Å).
I must be missing something fairly obvious, have not been able to see what.
Feedback would be much appreciated. Scripts are below.
Thanks!
Jan
mapmask \
mapin 2mol_2mFo-DFc.map \
xyzin 2mol_B.pdb \
mskout 2mol_2mFo-DFc_2mol_B.msk \
<< eof
border 2
eof
maprot \
mapin 2mol_2mFo-DFc.map \
mskin 2mol_2mFo-DFc_2mol_B.msk \
wrkout 2mol_2mFo-DFc_rot.map \
<< eof
MODE from
AVER
ROTA euler 152.440 110.243 28.112
TRANS -42.212 5.510 -57.243
eof
--
Jan Abendroth
Seattle / Bainbridge Island, WA, USA
home: Jan.Abendroth_at_gmail.com<http://Jan.Abendroth_at_gmail.com>
--
Jan Abendroth
Emerald Biostructures
Seattle / Bainbridge Island, WA, USA
home: Jan.Abendroth_at_gmail.com<http://Jan.Abendroth_at_gmail.com>
work: JAbendroth_at_embios.com
http://www.emeraldbiostructures.com
________________________________
To unsubscribe from the CCP4BB list, click the following link:
https://www.jiscmail.ac.uk/cgi-bin/webadmin?SUBED1=CCP4BB&A=1
--
Jan Abendroth
Emerald Biostructures
Seattle / Bainbridge Island, WA, USA
home: Jan.Abendroth_at_gmail.com<http://Jan.Abendroth_at_gmail.com>
work: JAbendroth_at_embios.com
http://www.emeraldbiostructures.com
________________________________
To unsubscribe from the CCP4BB list, click the following link:
https://www.jiscmail.ac.uk/cgi-bin/webadmin?SUBED1=CCP4BB&A=1
########################################################################
To unsubscribe from the CCP4BB list, click the following link:
https://www.jiscmail.ac.uk/cgi-bin/webadmin?SUBED1=CCP4BB&A=1
handle_read_draw_molecule_with_recentre("a.pdb", 1) #the reference pdb, imol 0
handle_read_draw_molecule_with_recentre("H.pdb", 1) #the pdb used for masking and as the moving molecule, imol 1
handle_read_ccp4_map('1.map', 0) #input map, imol 2
#set_last_map_contour_level(0.15)
#set_last_map_sigma_step( 0.010)
#set_last_map_colour( 0.80, 0.00, 0.80)
set_map_mask_atom_radius(3)
mask_map_by_molecule(2,1,1) #mask imol2 around imol 1 to generate imol 3, the masked map around H.pdb
transform_map_using_lsq_matrix(0,'L',1,500,1,'H',1,500,3,rotation_center(),100)
#lsq match imol 1 H chain to imol 0 L chain. Then rotate&translate imol 3, generate a new map centered at the current rotation center
#The current rotation center is set on line 2 when reading 'H.pdb'. The transformed map has a maximum radius of 100A around this point. However the saved map is also limited by the unit cell.
#This generates map imol 4
#Alternatively one can set the center of the new map by specifying the coordinates for the center point:
#transform_map_using_lsq_matrix(0,'L',1,500,1,'H',1,500,3,[0,0,0],100) #if the point of interest is at [0,0,0]
#improperly set center will cause map to be cut inside protein
export_map(4,'transposed_masked.map') #save imol 4, the transformed, masked map.
#exit(0)
#uncomment the above line to run without graphics as: coot --no-graphic -s scr.py
#otherwise one has to type (exit 0) to quit coot
