Classic EM Tutorial, Part II
In part II of the classic tutorial we demonstrate more advanced features of the matchpt tool which supports the registration of point clouds even in cases where they are of unequal size. Here we dock a single ncd monomer simultaneously to both monomer densities. One of the problems encountered in this project is the weaker density of the detached ncd monomer compared to the attached one (see part I). Therefore, we demonstrate first the docking of two maps with vol2pdbmatchpt, and pdb2vol, followed by the advanced histogram matching with the volhist tool, and finally map averaging with volaver, to generate a suitable docking target.

The usage of these routines is demonstrated in the following, using the examples of part I. As a prerequisite at least the installation step of part I must be followed. The rest of part I can be skipped if the input files are copied from the "solutions" directory (or generated with the run_tutorial.bash script). More documentation is available in the user guide, on the methodology page, and in the published articles.
Content:

Registration of Two Maps

As a prerequisite to density histogram matching, two maps must be brought into register to maximize the density cancellation zone (Wriggers et al, 2011). Here we use the output files 4_s1.situs and 4_s2.situs from part I as target and probe maps. Situs docking tools require an atomic structure (probe) and a density map (target), so the probe map must be first made movable by creating a pseudo-PDB file with vol2pdb

Here is the shell window input and output for the entire vol2pdb transformation:

% ./vol2pdb 4_s2.situs 6_s2.pdb
lib_vio> Situs formatted map file 4_s2.situs - Header information:
lib_vio> Columns, rows, and sections: x=1-11, y=1-12, z=1-13
lib_vio> 3D coordinates of first voxel: (414.000000,222.000000,120.000000)
lib_vio> Voxel size in Angstrom: 6.000000
lib_vio> Reading density data...
lib_vio> Volumetric data read from file 4_s2.situs
lib_vwk> Map density info: max 30.000000, min 0.000000, ave 3.154429, sig 6.621312.
vol2pdb> Warning: Densities below cutoff 0.005000 will be ignored.
vol2pdb> The PDB occupancy field contains the densities above cutoff rounded to two decimals.

Map densities have now been written to the PDB occupancy field.

Next, we perform a matching of the probe to the target map using matchpt, as shown in part I, saving the best solution to a new file 6_s2.dock.pdb. At the command prompt enter:

./matchpt NONE 4_s1.situs NONE 6_s2.pdb
mv mpt_001.pdb 6_s2.dock.pdb
rm mpt_*

Finally, we map the rotated and translated probe back to the grid of the target map using pdb2vol. In the following session we select a voxel spacing of 6Å as before:

./pdb2vol 6_s2.dock.pdb 6_s2.dock.situs
lib_pio> 296 atoms read.
pdb2vol> Found 0 hydrogens, 0 water atoms, 0 codebook vectors, 296 density atoms
pdb2vol> Mass-weighting on.
pdb2vol> B-factor thresholding off.
pdb2vol> 296 out of 296 atoms selected for conversion.
pdb2vol>
pdb2vol> The input structure measures 46.709 x 53.183 x 66.460 Angstrom
pdb2vol>
pdb2vol> Please enter the desired voxel spacing for the output map (in Angstrom): 6
pdb2vol>
pdb2vol> Projecting atoms to cubic lattice by trilinear interpolation...
pdb2vol> ... done. Lattice smoothing (sigma = atom rmsd):  4.247 Angstrom
lib_vio> Writing density data...
lib_vio> Volumetric data written in Situs format to file 6_s2.dock.situs
lib_vio> Situs formatted map file 6_s2.dock.situs - Header information:
lib_vio> Columns, rows, and sections: x=1-14, y=1-15, z=1-17
lib_vio> 3D coordinates of first voxel: (384.000000,240.000000,126.000000)
lib_vio> Voxel size in Angstrom: 6.000000
pdb2vol>
pdb2vol> Map projected to lattice only, and written to file 6_s2.dock.situs.
pdb2vol> Effective spatial resolution (2 sigma) of lattice projection:  8.494A

Density Histogram Matching of Two Registered Maps

Now we can interactively match the map densities with
volhist using an affine (linear) model with gain and bias (Wriggers et al, 2011). This mode is activated when calling volhist with three arguments. After entering

./volhist 4_s1.situs 6_s2.dock.situs tmp.situs

enter the number of histogram bins (30). Next, enter the surface value 20 for file 4_s1.situs and a value of 12 for file 6_s2.dock.situs (see part I of this tutorial). These surface isovalues will be automatically matched. Subsequently, the user has the uption of entering various bias values to interactively center the central peak of the difference histogram. We find that a bias value of -4.0 (corresponding to a gain factor of 2.0) provides for a well centered central peak. To save the transformation of 6_s2.dock.situs to the file tmp.situs, enter any value above 20.

Here is the appearance of the central peak for bias value -4.0:

volhist> Enter new bias value < 20.000000 to try again, or any value >= 20.000000 to exit and save current transformation: -4
volhist> Computing difference histogram of (isovalue thresholded) infile1 - (infile2 * gain + bias).
lib_vwk> Printing centered difference histogram (31 histogram bins used):
lib_vwk> Difference density range: -36.265322 to 35.000000, lower tail clipped, showing 31 histogram bins about center:
 -35.000 | 0   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
 -32.667 |= 1                                                                          
 -30.333 |= 1  .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
 -28.000 |===== 3                                                                      
 -25.667 |=== 2    .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
 -23.333 |======= 4                                                                    
 -21.000 |===== 3  .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
 -18.667 | 0                                                                           
 -16.333 |= 1  .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
 -14.000 |=== 2                                                                        
 -11.667 |======= 4    .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
  -9.333 |============== 8                                                             
  -7.000 |=================================== 19   .   .   .   .   .   .   .   .   .   .
  -4.667 |======================================== 22                                  
  -2.333 |============================================== 25    .   .   .   .   .   .   .
   0.000 |====================================================================== 38 ....
   2.333 |============================================ 24  .   .   .   .   .   .   .   .
   4.667 |========================================== 23                                
   7.000 |========================================== 23    .   .   .   .   .   .   .   .
   9.333 |===== 3                                                                      
  11.667 |================ 9   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
  14.000 | 0                                                                           
  16.333 |=== 2    .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
  18.667 | 0                                                                           
  21.000 |============================================================== 34    .   .   .
  23.333 |========================= 14                                                 
  25.667 |========================= 14 .   .   .   .   .   .   .   .   .   .   .   .   .
  28.000 |======================= 13                                                   
  30.333 |===== 3  .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
  32.667 |=== 2                                                                        
  35.000 |= 1  .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
volhist> Gain, bias used for affine transformation that matches infile2 to infile1: 2.000000, -4.000000 


The matching of the docked probe map was performed here to obtain the gain and bias values that will be applied to the original probe map in the following map averaging. In the run_tutorial.bash script this entire step is ignored because it requires interactive exploration that cannot be automated.  We assume the gain and bias values (2.0 and -4.0) are known in the following. Also, the output file tmp.situs is no longer needed.
Map Averaging

We wish to use the map averaging tool volaver to build up a combined map from the (original, undocked) 4_s1.situs and 4_s2.situs monomer maps, where the "weaker" density of 4_s2.situs is subject to affine transformation with gain 2.0 and bias -4.0. This "strengthening" of the weaker monomer is necessary to give both of them equal weight in the docking below. To apply the affine transformation to 4_s2.situs, call volhist with only two arguments:

./volhist 4_s2.situs 7_s2.situs

Enter the number of histogram bins (30). Then enter the gain factor (2.0) and offset value (-4). The output file is written to disk. Here is the complete volhist session:

./volhist 4_s2.situs 7_s2.situs              
llib_vio> Situs formatted map file 4_s2.situs - Header information:
lib_vio> Columns, rows, and sections: x=1-11, y=1-12, z=1-13
lib_vio> 3D coordinates of first voxel: (414.000000,222.000000,120.000000)
lib_vio> Voxel size in Angstrom: 6.000000
lib_vio> Reading density data...
lib_vio> Volumetric data read from file 4_s2.situs
lib_vwk> Density information. min: 0.000000 max: 30.000000 ave: 3.154429 sig: 7.216357 norm: 7.875673
lib_vwk> Above zero density information. ave: 18.287162 sig: 5.016389 norm: 18.962712
lib_vwk> Please enter the of number histogram bins: 30
lib_vwk> Printing voxel histogram, 30 histogram bins
lib_vwk> (density value; voxel count; top-down cumulative volume fraction):
   0.000 |=======================================================================-> 1420 | 1.000e+00
   1.034 | 0                                                                             | 1.725e-01
   2.069 | 0 .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 1.725e-01
   3.103 | 0                                                                             | 1.725e-01
   4.138 | 0 .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 1.725e-01
   5.172 | 0                                                                             | 1.725e-01
   6.207 | 0 .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 1.725e-01
   7.241 | 0                                                                             | 1.725e-01
   8.276 | 0 .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 1.725e-01
   9.310 | 0                                                                             | 1.725e-01
  10.345 | 0 .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 1.725e-01
  11.379 | 0                                                                             | 1.725e-01
  12.414 |======================================================== 28    .   .   .   .   | 1.725e-01
  13.448 |====================================================================== 35      | 1.562e-01
  14.483 |==================================================== 26    .   .   .   .   .   | 1.358e-01
  15.517 |====================================================================== 35      | 1.206e-01
  16.552 |========================================================== 29  .   .   .   .   | 1.002e-01
  17.586 |============================================== 23                              | 8.333e-02
  18.621 |========================== 13  .   .   .   .   .   .   .   .   .   .   .   .   | 6.993e-02
  19.655 |==================================== 18                                        | 6.235e-02
  20.690 |======================== 12    .   .   .   .   .   .   .   .   .   .   .   .   | 5.186e-02
  21.724 |==================== 10                                                        | 4.487e-02
  22.759 |====================== 11  .   .   .   .   .   .   .   .   .   .   .   .   .   | 3.904e-02
  23.793 |====================== 11                                                      | 3.263e-02
  24.828 |========================== 13  .   .   .   .   .   .   .   .   .   .   .   .   | 2.622e-02
  25.862 |============== 7                                                               | 1.865e-02
  26.897 |============== 7   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 1.457e-02
  27.931 |======== 4                                                                     | 1.049e-02
  28.966 |================== 9   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 8.159e-03
  30.000 |========== 5                                                                   | 2.914e-03
lib_vwk> Maximum at density value 0.000
volhist> Enter a scaling factor (gain) by which map densities will be multiplied: 2
volhist> Enter offset density value (bias), this will be added after scaling by gain factor: -4
volhist> Calculating new voxel histogram
lib_vwk> Density information. min: -4.000000 max: 56.000000 ave: 2.308858 sig: 14.432713 norm: 14.616225
lib_vwk> Above zero density information. ave: 32.574324 sig: 10.032778 norm: 34.084355
lib_vwk> Printing voxel histogram, 30 histogram bins
lib_vwk> (density value; voxel count; top-down cumulative volume fraction):
  -4.000 |=======================================================================-> 1420 | 1.000e+00
  -1.931 | 0                                                                             | 1.725e-01
   0.138 | 0 .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 1.725e-01
   2.207 | 0                                                                             | 1.725e-01
   4.276 | 0 .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 1.725e-01
   6.345 | 0                                                                             | 1.725e-01
   8.414 | 0 .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 1.725e-01
  10.483 | 0                                                                             | 1.725e-01
  12.552 | 0 .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 1.725e-01
  14.621 | 0                                                                             | 1.725e-01
  16.690 | 0 .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 1.725e-01
  18.759 | 0                                                                             | 1.725e-01
  20.828 |======================================================== 28    .   .   .   .   | 1.725e-01
  22.897 |====================================================================== 35      | 1.562e-01
  24.966 |==================================================== 26    .   .   .   .   .   | 1.358e-01
  27.034 |====================================================================== 35      | 1.206e-01
  29.103 |========================================================== 29  .   .   .   .   | 1.002e-01
  31.172 |============================================== 23                              | 8.333e-02
  33.241 |========================== 13  .   .   .   .   .   .   .   .   .   .   .   .   | 6.993e-02
  35.310 |==================================== 18                                        | 6.235e-02
  37.379 |======================== 12    .   .   .   .   .   .   .   .   .   .   .   .   | 5.186e-02
  39.448 |==================== 10                                                        | 4.487e-02
  41.517 |====================== 11  .   .   .   .   .   .   .   .   .   .   .   .   .   | 3.904e-02
  43.586 |====================== 11                                                      | 3.263e-02
  45.655 |========================== 13  .   .   .   .   .   .   .   .   .   .   .   .   | 2.622e-02
  47.724 |============== 7                                                               | 1.865e-02
  49.793 |============== 7   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 1.457e-02
  51.862 |======== 4                                                                     | 1.049e-02
  53.931 |================== 9   .   .   .   .   .   .   .   .   .   .   .   .   .   .   | 8.159e-03
  56.000 |========== 5                                                                   | 2.914e-03
lib_vwk> Maximum at density value -4.000
lib_vio> Writing density data...
lib_vio> Volumetric data written in Situs format to file 7_s2.situs
lib_vio> Situs formatted map file 7_s2.situs - Header information:
lib_vio> Columns, rows, and sections: x=1-11, y=1-12, z=1-13
lib_vio> 3D coordinates of first voxel: (414.000000,222.000000,120.000000)
lib_vio> Voxel size in Angstrom: 6.000000

Note how the cutoff density in the histogram was moved from 12 to 20 (as in 4_s1.situs, which you can check with "./volhist 4_s1.situs" using only one argument).

Now we are finally ready to combine the monomer densities from both maps into a single "dimer" target map. The maps can be added using the volaver tool (subject to normalization by factor 2). At the prompt enter:

./volaver 4_s1.situs 7_s2.situs 7_combined.situs

The file 7_combined.situs now contains the dimer density.
Unequal Size Point-Cloud Registration

We now dock the high-resolution ncd structure into the larger dimer map with the matchpt utility. The dimer map fits exactly two input structures so we set the -units option to 2.0 (in general, this parameter may be non-integer). We also sort here the solutions by cross-correlation (CC) value using the "-ranking 1" option:

./matchpt NONE 7_combined.situs NONE 0_ncd.pdb -units 2.0 -ranking 1

The program selects again 7 codebook vectors per monomer using the RMSD criterion, and returns 10 best fits, ranked by the CC value. The program also shows the permutation of the vectors that determine the superposition.
The solutions and corresponding codebook vectors are written to disk. Here is the complete matchpt session:

./matchpt NONE 7_combined.situs NONE 0_ncd.pdb -units 2.0 -mincv 7 -maxcv 7
matchpt> File1 == NONE
lib_vio> Situs formatted map file 7_combined.situs - Header information:
lib_vio> Columns, rows, and sections: x=1-14, y=1-15, z=1-16
lib_vio> 3D coordinates of first voxel: (396.000000,222.000000,120.000000)
lib_vio> Voxel size in Angstrom: 6.000000
lib_vio> Reading density data...
lib_vio> Volumetric data read from file 7_combined.situs
matchpt> File3 == NONE
matchpt> Loading high-resolution structure PDB file: 0_ncd.pdb
matchpt> No codebook vectors available, will compute a series of vector
matchpt> sets of different sizes and rank them.
matchpt>
matchpt> Number of structure units contained in target volume set to 2
matchpt> Ranking option set to 1 (min RMSD / cross correlation)
matchpt>
matchpt> Computation and clustering of sets of (high res./low res.) codebook vectors.
lib_mpt> Variability of  8 codebook vectors in 8 statistically independent runs: 0.358 A
lib_mpt> Variability of 10 codebook vectors in 8 statistically independent runs: 1.269 A
lib_mpt> Variability of 12 codebook vectors in 8 statistically independent runs: 0.725 A
lib_mpt> Variability of 14 codebook vectors in 8 statistically independent runs: 1.889 A
lib_mpt> Variability of  4 codebook vectors in 8 statistically independent runs: 0.938 A
lib_mpt> Variability of  5 codebook vectors in 8 statistically independent runs: 1.152 A
lib_mpt> Variability of  6 codebook vectors in 8 statistically independent runs: 1.937 A
lib_mpt> Variability of  7 codebook vectors in 8 statistically independent runs: 0.301 A
lib_mpt> Variability of 16 codebook vectors in 8 statistically independent runs: 1.521 A
lib_mpt> Variability of 18 codebook vectors in 8 statistically independent runs: 2.673 A
lib_mpt> Variability of  8 codebook vectors in 8 statistically independent runs: 3.490 A
lib_mpt> Variability of  9 codebook vectors in 8 statistically independent runs: 2.141 A

matchpt>
matchpt> Point cloud matching for (high res./low res.) vectors:
matchpt> 4/8 vectors: 1 match found with RMSD 5.735 A, variabilities 0.938/0.358 A, CC: 0.410
matchpt> 5/10 vectors: 7 matches found with min RMSD 3.725 A, variabilities 1.152/1.269 A, max CC 0.716
matchpt> 6/12 vectors: 8 matches found with min RMSD 4.262 A, variabilities 1.937/0.725 A, max CC 0.702
matchpt> 8/16 vectors: 5 matches found with min RMSD 4.331 A, variabilities 3.490/1.521 A, max CC 0.737
matchpt> 9/18 vectors: 2 matches found with min RMSD 5.525 A, variabilities 2.141/2.673 A, max CC 0.679
matchpt> 7/14 vectors: 17 matches found with min RMSD 2.760 A, variabilities 0.301/1.889 A, max CC 0.742
matchpt>
matchpt> Based on min RMSD, selecting 17 matches from 7/14 vectors. Exploring top 10 CC solutions.
matchpt>
matchpt> Solution filename, codebook vector RMSD in Angstrom,
matchpt> cross-correlation coefficient, and permutation
matchpt> (order of low res fitted to high res vectors):
matchpt>
matchpt> mpt_001.pdb - RMSD:  5.404  CC:  0.742 - ( 6,11,12, 5, 7, 3, 9)
matchpt> mpt_002.pdb - RMSD:  2.760  CC:  0.696 - ( 4, 8, 1,10, 2,13,14)
matchpt> mpt_003.pdb - RMSD:  5.319  CC:  0.634 - (13, 1,10, 8, 2,14, 4)
matchpt> mpt_004.pdb - RMSD:  4.595  CC:  0.618 - (14,10, 8, 1, 2, 4,13)
matchpt> mpt_005.pdb - RMSD:  8.291  CC:  0.612 - (14, 2,10,13, 8, 1, 4)
matchpt> mpt_006.pdb - RMSD:  8.390  CC:  0.598 - ( 9, 8, 3, 6, 5,12,11)
matchpt> mpt_007.pdb - RMSD:  9.101  CC:  0.565 - ( 6, 7,11, 9,12, 5, 3)
matchpt> mpt_008.pdb - RMSD:  9.848  CC:  0.539 - ( 9, 4, 3, 6, 5,12,11)
matchpt> mpt_009.pdb - RMSD:  8.783  CC:  0.525 - ( 4, 6, 1, 8,13,10,14)
matchpt> mpt_010.pdb - RMSD:  7.775  CC:  0.509 - ( 8,10, 1,13,11, 4,14)
matchpt>
matchpt> All done.

The best solutions in this list (by CC value) are mpt_001.pdb and mpt_002.pdb. These are nearly identical to the solutions found earlier (files 5_s2_1.dock.pdb and 5_s1_1.dock.pdb). 

The following sequence of commands in the VMD text console (cf. VMD user guide) will load the two best solutions (red and orange), the microtubule surface (green), and the histogram matched dimer map (blue), and render them similarly to the final result in part I:

mol load pdb mpt_001.pdb
mol load pdb 
mpt_002.pdb
mol load situs 3_mt.situs
mol load situs 7_combined.situs
mol top 0
rotate stop
display resetview
display projection orthographic
mol modstyle 0 0 Cartoon 2.1 11 5
mol modstyle 0 1 Cartoon 2.1 11 5
mol modstyle 0 2 Isosurface 20 0 0 0 1 1
mol modstyle 0 3 Isosurface 10 0 0 1 1 1
mol modcolor 0 0 ColorID 3
mol modcolor 0 1 ColorID 1
mol modcolor 0 2 ColorID 7
mol modcolor 0 3 ColorID 0

Don't forget to hit "enter" after the last line! 

The result should look like this:

vmd result
(Click image to enlarge)

Note: As described in (Birmanns & Wriggers, 2007), the fast point cloud matching can be improved by a post-processing with the fast correlation-based tool collage. We have explained this workflow in the separate multi-dock tutorial.

Return to the front page .