Virions

P22 virion structure

Test
renderings of the P22 virion EMD-1220 data set

"For
a physicist, it is indeed a great joy to learn how we can use beautiful
mathematics to understand the real world." In Stephen Weinberg's "Without God" (9/25/08, The New York Review of Books)

The EMD_1220 volume
data set was reconstructed from cryo-electron
microscopy images. It represents, roughly speaking, a three
dimensional estimate of electron densities of an average P22 virion
particle.

These notes outline post-processing that I performed in order to do renderings of the structure and sub-structures represented in the EMD_1220 volume.

These notes outline post-processing that I performed in order to do renderings of the structure and sub-structures represented in the EMD_1220 volume.

Data source

Processing notes

Other structural questions

Processing notes

Convert
from CCP4 map (.map/.xml) to Space
volume (.vol) file format

Allign the virion's raw data longitudinal axis with the volume's z-axis accurately

Coregister over-sampled original volume data to axially symmetric target

Rotational symmetry decomposition

Spectral density segmentation

Programs and codeAllign the virion's raw data longitudinal axis with the volume's z-axis accurately

Coregister over-sampled original volume data to axially symmetric target

Rotational symmetry decomposition

Spectral density segmentation

Other structural questions

The homogeneous (coordinate) transformation matrix for 3D points

2008-09-08 I just took a look at the EMD_1220 volume data set and it is very nice. My goal is to fully segment and color render, with cutaways, this volume data.

2008-09-18 The post-processing has morphed into a somewhat separate project: make an average over axial n-fold rotations, where n is point dependent and point-wise derived. I considered this approach to "cleaning" the data (averaging out noise that doesn't have a rotational symmetry) while rendering EMD-1222 [To Do: link]. I've tested the idea (see rotational averages) and now can write the MATLAB code that implements it with precision.

2009-10-08 These notes are reconstructed, fixed up and reformatted from a mangled version that was dumped from the previous site, Enzymind.com. Currently I'm reinstalling FSL so I can make more progress on this.

"for a higher resolution (better) reconstruction check out EMD-1220"

This EMD-1220 data set (a cryo-EM assymetric reconstruction, 3-D grayscale) is available (accessed 2008-09-08) at http://www.ebi.ac.uk/msd-srv/emsearch/atlas/1220_downloads.html, in the European Bioinformatics Institute's Protein Data Bank Europe (PDBe). Also see the EM Data Bank, which provides three-dimensional density maps for non-viruses.

From the data set web page description:

The structure of an infectious P22 virion shows the signal for headful DNA packaging

Gabriel C. Lander, Liang Tang, Sherwood R. Casjens, Eddie B. Gilcrease, Peter Prevelige, Anton Poliakov, Clinton S. Potter, Bridget Carragher, and John E. Johnson:

The structure of an infectious p22 virion shows the signal for headful DNA packaging J Mol Biol (2006) 312, pp. 1791-1795 [PubMed entry 16709746]

3-Jun-2006

singleParticle, 17 Å (resolution determined by FSC at 0.5 cut-off

** **

2008-09-08

MATLAB program convert_map_volume.m does this (see attached). First read the .xml file (as text) to find the volume dimensions (the EMD_1220.map is 256 x 256 x 256). A hardcoded black and white point must be chosen for each data set conversion. A histogram of final values is displayed that can help with this.

I'm not sure I captured the full dynamic range, particularly at the lowest values. Add gamma conversion control? Or redo with 32-bit NIFTI format?

MATLAB program convert_map_volume.m does this (see attached). First read the .xml file (as text) to find the volume dimensions (the EMD_1220.map is 256 x 256 x 256). A hardcoded black and white point must be chosen for each data set conversion. A histogram of final values is displayed that can help with this.

I'm not sure I captured the full dynamic range, particularly at the lowest values. Add gamma conversion control? Or redo with 32-bit NIFTI format?

MATLAB command:

>> volData = convert_map_volume( 'emd_1220.map', 'emd_1220_8-bit', [ 256 256 256 ] );

2008-09-14

Rewrote convert_map_volume.m for NIFTI 32-bit output.

Rewrote convert_map_volume.m for NIFTI 32-bit output.

MATLAB command:

>> volData = convert_map_volume( 'emd_1220.map', 'emd_1220', [ 256 256 256 ] ); [To Do: link]

>> volData = convert_map_volume( 'emd_1220.map', 'emd_1220', [ 256 256 256 ] ); [To Do: link]

2008-09-16

The approach I chose was to find a best fit longitudinal axis took advantage of the high radial 5-fold symmetry of the capsid. The capsid was first roughly alligned with the volume's z-axis (within about an angle of 1 px./80 px., just an eyeball resampling at this orientation in Space), and roughly segmented (a binary mask that included almost all the capsid, but not the interior or tail structures). The masked volume (just the capsid) was then averaged with its four roughly equivalent 2*pi/5 rotations. This average volume is axially symmetric, and served as a target for coregistration of the raw data.

The approach I chose was to find a best fit longitudinal axis took advantage of the high radial 5-fold symmetry of the capsid. The capsid was first roughly alligned with the volume's z-axis (within about an angle of 1 px./80 px., just an eyeball resampling at this orientation in Space), and roughly segmented (a binary mask that included almost all the capsid, but not the interior or tail structures). The masked volume (just the capsid) was then averaged with its four roughly equivalent 2*pi/5 rotations. This average volume is axially symmetric, and served as a target for coregistration of the raw data.

2008-09-08

Cubic resampled (x2) and linearly resampled in an approximately axial orientation. Binary mask for substructures (capsid, inner head, outer head, DNA shells 1, 2, 3, and 4) using 2 and 3-D fill on the full binay image. The capsid has the primary orientation information, at least for any 5-fold symmetry structure. All using Space Software.

Cubic resampled (x2) and linearly resampled in an approximately axial orientation. Binary mask for substructures (capsid, inner head, outer head, DNA shells 1, 2, 3, and 4) using 2 and 3-D fill on the full binay image. The capsid has the primary orientation information, at least for any 5-fold symmetry structure. All using Space Software.

After masking the high-value 8-bit volume with the rough capsid binary segmentation (using Space), FSL was used to make 5 copies, at n*2*pi/5 rotations.

- 2*pi/5 rotation transform matrix:

MATLAB calculation for xt and yt
coordinate transformation matrix:

>> xt = xc + ( sqrt( xc*xc + yc*yc ) )*sin(
(2*pi)/5 - atan( xc/yc ) )

>> yt = yc - ( sqrt( xc*xc + yc*yc ) )*cos(
(2*pi)/5 - atan( xc/yc ) ) ,

where xc = 165 and yc = 163 are the
center coordinates with respect to the corner origin.

The full transformation matrix, in homogenous
coordinates:

```
0.309017
-0.951057 0 269.0934
```

0.951057 0.309017 0 -44.2941

0 0
1 0

0 0
0 1

- Apply rotation transform matrix four times to the masked capsid volume:

FLIRT
| Utils | Apply FLIRT transforms
(ApplyXFM)

or

flirt
-applyxfm -init -out -ref

- These 5 volumes were averaged to create capsid_mean.nii:

FSL merge of the five capsid
rotations to 4-D volume and collapse across the "time" (n=5) dimension:

fslmerge -t
capsid_avg.nii.gz emd_1220_rot_0 emd_1220_rot_72 emd_1220_rot_144
emd_1220_rot_216 emd_1220_rot_288

fslmaths
capsid_avg -Tmean capsid_mean

Note that this took a long time,
apparently due to memory limitations. I think this time could be
reduced using fslmaths -add, but would require a conversion to a
floating point format, and adding each one sequentially. Adding in the
8-bit format caused value wrap-around (8-bit overflow).

I think that adding sequentially will
be required for the the alligned
image rotation averages, to maintain the bit depth without memory
limitations.

[To Do: Screenshot of capsid and capsid_mean section.]

An assymetric capsid slice (left, high
contrast image of an angled slice, resampled using linear
interpolation) compared with the coregisitration target (right), a
5-fold mean image, on orthogonal slice planes:

Comparison of the nearly alligned
roughly segmented capsid (left) and the 5-fold rotational average
(right). At top is an axial slice through the vertex at top of capsid
(magnified x2 relative to bottom images), where the mis-allignment is
most pronounced.

2008-09-15

Raw data was cropped
(to limit memory footprint) and over-sampled (tri-cubic) by a factor of
x2, in Space Software.

The raw
256x256x256x4 byte = 64 MB volume is too large to resample x23. So need
to crop (in Space after Volume |
Volume Info | Center, exact), both this raw data and the capsid
weighting volume, emd_1220_8-bit_high-125-230_x2_capsid_weight.nii.gz.

The raw data crop at native resolution:

A binary coregistration weighting volume was created (emd_1220_8-bit_high-125-230_x2_capsid_weight_cropped.nii) by segmenting the capsid at the original orientation, in a procedure identical to that above (rough segmentation of capsule head).

This was using a high values (125-230) 8-bit volume, but at the original orientation. This tilted orientation made it more difficult to segment the tail/capsid boundary, but it doesn't need to be perfect -- even a majority portion of the capsid would coregister well.

The raw data crop at native resolution:

x: [ 41 : 212 ] (even width s.t.
the resampled dimensions will be odd)

y: [ 35 : 218 ]

z: [ 39 : 256 ]

y: [ 35 : 218 ]

z: [ 39 : 256 ]

A binary coregistration weighting volume was created (emd_1220_8-bit_high-125-230_x2_capsid_weight_cropped.nii) by segmenting the capsid at the original orientation, in a procedure identical to that above (rough segmentation of capsule head).

This was using a high values (125-230) 8-bit volume, but at the original orientation. This tilted orientation made it more difficult to segment the tail/capsid boundary, but it doesn't need to be perfect -- even a majority portion of the capsid would coregister well.

Weighting origin
tweaked to exactly overlap raw data:

Note the value wrap flaw on the high luminance "needle" center at bottom-left. This might be due to a flaw in the Space tri-cubic resampling algorithm, or in the MATLAB conversion code. [ 2008-09-15 I can minimize the effect of this, but it seems to be a problem with Space's storage or interpretation of NIFTI range information.]

FSL coregistration:

Target: capsid_mean.nii

Object: emd_1220_cropped_x2.nii

Weighting of object: emd_1220_8-bit_high-125-230_x2_capsid_weight_cropped.nii

Coregistered output: emd_1220_cropped_x2_alligned.nii

Normalized correlation method. Straight correlation method doesn't work, because the raw values (32-bit) are a very different range than the mean capsid (8-bit).

Tri-linear resampling.

flirt -in /emd_1220_cropped_x2.nii.gz -ref capsid_mean.nii.gz -out /emd_1220_cropped_x2_alligned.nii.gz -omat /emd_1220_cropped_x2_alligned.mat -bins 256 -cost normcorr -searchrx 25 35 -searchry -5 5 -searchrz -5 5 -dof 9 -inweight /emd_1220_8-bit_high-125-230_x2_capsid_weight_cropped.nii.gz -interp trilinear

Object: emd_1220_cropped_x2.nii

Weighting of object: emd_1220_8-bit_high-125-230_x2_capsid_weight_cropped.nii

Coregistered output: emd_1220_cropped_x2_alligned.nii

Normalized correlation method. Straight correlation method doesn't work, because the raw values (32-bit) are a very different range than the mean capsid (8-bit).

Tri-linear resampling.

flirt -in /emd_1220_cropped_x2.nii.gz -ref capsid_mean.nii.gz -out /emd_1220_cropped_x2_alligned.nii.gz -omat /emd_1220_cropped_x2_alligned.mat -bins 256 -cost normcorr -searchrx 25 35 -searchry -5 5 -searchrz -5 5 -dof 9 -inweight /emd_1220_8-bit_high-125-230_x2_capsid_weight_cropped.nii.gz -interp trilinear

The final alligned
volume, emd_1220_cropped_x2_alligned.nii,
is 32-bit, 333 x 339 x 435 vx. (after post-cropping), origin at (164,
169, 217), 188 MB:

EMD-1220, orthogonal slices of alligned volume, high contrast. Averaged 5-fold rotation of capsid/DNA and 6-fold rotation of tail structures.## Rotational averages

2008-09-16

2008-09-17

Write a specialized MATLAB program for radial symmetries analysis and n-fold rotational averaging. Started rotation_average_volume.m.

Here's a proof of concept image using a test version of rotation_average_image.m. The concentric DNA rings near the axial capsid center are nearly circular. By rotational averaging the central radial regularities are highlighted:

The original image (left) and an average with its pi/2 radian rotation about the center (right). The ring definition near the center is enhanced. Note that the capsid (pink at left) does not have a 2-fold symmetry on this plane, so its rotation has only a small overlap (pink at right).

DNA rings averaged over 5-fold (left) and 90-fold (right) rotations. The 5-fold rotation matches the capsid symmetry. Nine DNA rings are clear, but the central features are likely artifacts (noise, or different rotational symmetry).

2008-09-21

Tentative procedure to find center accurately and average over 2*pi/5 and 2*pi/6 rotations. I've eperimented with this procedure on 2-D slices, but only minor modification need to be made to do 3-D serial slices:

2008-10-01

- Remove background (set background to marker symbol, a low value). See remove_background_peak_image for 2-D demo. [To Do: Don't remove background values from central features.]
- For 3-D volume, use Space 3-D fill tool to accomplish this. Make an 8-bit mask and apply to 32-bit volume.

Did this process on emd_1220_cropped_x2_alligned.nii. Note that it is not needed for the rotational averaging algorithms (below). But for finding the center in single slices, it worked well to remove the contribution from noise, only taking into account the standard deviation of histogram values that were not on the peripheral noise.

- Exterior mask Overlay | New Blank Overlay | 8-bit with Edit | 3-D Fill (range = 1). Highest contrast on original, at a black point that best segregated interior an exterior. Fill in all exterior to the capsid (except for discontinuous blobs) without much attention to the exact background level. At lower tail planes the background values were contiguous, so used a wide pen tool to manually segment (circle) ROIs.
- Interior mask: Close original, make a new blank overlay to the exterior mask. Fill interior. There is some bleed into random noise contiguous with the interior, but it is not important for rotational averaging as at least n-1 points must be non-zero for an n-fold rotation.
- Overlay | Mask original (16-bit) with interior mask (8-bit, on top), and Overlay | Merge. Save as NIFTI because there is some BP/WP/Min/Max (Colors | Transforms and Histograms) interpretation bug with Space after re-opening.
- Sub-pixel approximation of center by applying rotations (5 and/or 6-fold) an minimizing standard deviation. Iterative hill climbing.
>> rotation_average_image( 'test_remove_background-zeroback.png', [169 169], 5, true, 'test_' );

2008-09-22

To Do: Ignore zero values (background), and require
a
minimum number of low deviation points. Ignore background point in net
deviation calculation. Set map points to zero if not enough forground
points on rotation.

2008-09-22

To Do: Accept an n-fold rotation map and only
calculate
deviation of rotation points for each corresponding n-fold rotation.

- Use center estimate to find best n-fold rotation by finding minimum standard deviation of each rotated point with each n-fold rotation.

>>rotation_nfold_image( 'test_remove_background-zeroback.png', [169.2383 170.0625], [5 6 10 12 15 18 24 25 30], true, '' );

2008-09-22

To Do: Jitter the center point for each point
rotation.
Save the maximal point for the n-fold rotation that was most stable
with jitter. Also save the jitter maps.

2008-09-24

Experimented with this a lot. This point by point
method
(evaluating best n-fold fit at each point) is not the way to go -- to
noisy and not easy to integrate neighborhood and global n-fold
information. See more simple approach below.

Repeat center
finding,
using the new n-fold rotation map. Different centers for different n?

Rough 2 part segmentation of alligned volume by symmetry (5-fold for capsid and DNA, 6 or 12-fold for tail). Include some overlap at intersections. Find best fit n-fold rotation and center for all of one segmented slice, checking for other possible symmetries (circular symmetry, nearby odd symmetries). Repeat on all slices. Trim and blend two (or more) symmetry averaged volumes.

Rough 2 part segmentation of alligned volume by symmetry (5-fold for capsid and DNA, 6 or 12-fold for tail). Include some overlap at intersections. Find best fit n-fold rotation and center for all of one segmented slice, checking for other possible symmetries (circular symmetry, nearby odd symmetries). Repeat on all slices. Trim and blend two (or more) symmetry averaged volumes.

5-fold
rotational averaging of capsid (white on right) and 6-fold averaging of
tail (pink on right). Some (unknown) parts of the 12-fold symmetric
"bumps" at the tail periphery are artifacts of the rotational
averaging. Ameliorate by excluding oulliers on annular rings before
averaging.

Used rotational center finding and averaging on several planes of the
16-bit alligned data with exterior noise zeroed (see Exterior mask
above). It is beautiful, both 5 and 6-fold components on the critical
boundary planes. The calculated centers were dead on for all capsid
only planes, and only a possible +-.25 pixel shift of tail only planes,
with little or no axis skewing.

Currently used this program, with only one plane of interest and one symmetry hardcoded, to center find and rotationally average:

Currently used this program, with only one plane of interest and one symmetry hardcoded, to center find and rotationally average:

>> rotation_average_volume( 'emd_1220_cropped_x2_alligned_background_masked.nii', [170 165] );

[To Do: Example average images, derived center estimates.]

2008-10-01

- Extract an annulus that has 5 and 6-fold components from a single plane (radial components_image.m). Do component analysis manually and code annulus_components.m
- Average over m, n, ... -fold increments, excluding outliers. Hold in array with dimensions radius, angle, [1, m, n, ...].
- Filter high frequencies ( > m*n cycles/annulus) and normalize.
- Least squares fit of these models to the original. Find residual.

2008-10-06

Testing radial_components_image.m

5-fold symmetric capsid structures (blue) overlaid with 12-fold injectosome structure. This axial plane is near where these two symmetries come in contact, or overlap.2008-10-15

>> fPwr = radial_components_volume( 'emd_1220_cropped_x2_alligned.nii', [170 165] );

--> calls: fpwower = radial_components_image( imPlane, xyCenter, false, num2str( iP ) );

`% fPwr = radial_components_volume( strInputFilename, xyCenter )`

% ------------------------------------------------------------------------------

%

% Fourier decomposition of annuli about the center of a selected serial

% grayscale z-planes of a NIFTI volume.

`%%%%%%%%%%%%%%%%%%%%%%%%`

`% Hardcoded information:`

bShow = false; % [true, false] If true, displays each raw slice.

bClear = false; % [true, false] Unload volume from memory while analyzing each plane.

iPlanes = 1:435 % [ n : m ] z-coordinates of planes to be analyzed.

rRes = 2;

`%%%%%%%%%%%%%%%%%%%%%%%%`

% fPwr = radial_components_image( strInputFilename, xyCenter, bShow, strOutputFilePrefix )

% ------------------------------------------------------------------------------

%

% Fourier decomposition of annuli about the center of a grayscale image.

%

% This version is specialized for analysis of the 5 and 6-fold symmetry

% of the P22 virion (particularly the EMD-1220 data set), but it could be

% modified for other annular symmetries.

%%%%%%%%%%%%%%%%%%%%%%%%

% Hardcoded information:

bWrite = true;

bDoReconstruct = true;

rRes = 2;

% To Do: Calculate rmax from the given center and image size, use this

% parameter as a second hard maximum.

rmax = 163;

% Annular sampling rate; the number of sample points on each annulus.

% Lower integer factors have some computational speed advantage.

% This is longer than needed, and takes up considerable memory.

np = 2*2*3*3*5*7*11 % = 13860, excludes 8 as a factor

%

`%%%%%%%%%%%%%%%%%%%%%%%%`

% fPwr - total Fourier power of annuli at each radius.

% 1 - 1 cycle/annulus

% 2 - 2*n, excluding mod 5 and mod 6

% 3 - 3*n, excluding mod 5 and mod 6

% 4 - 4*n, excluding mod 5 and mod 6

% 5 - 5*n cycles/annulus, excluding 30, 60, ...

% 6 - 6*n cycles/annulus, excluding 30, 60, ...

% 7 - 7*n, excluding mod 5 and mod 6

% 8 - all nonzero and not in 5 or 6-fold

% 7 - [not used, empty]

% 10 - zeroth (DC, or constant) component

% 11 - 11*n, excluding mod 5 and mod 6

% 12 - 12

% 13 - 13*n, excluding mod 5 and mod 6

% 14 - 2

% 15 - 3

% 16 - 4

% 17 - 6

% 18 - 8

% 19 - 30*n

% 20 - 60*n

>> serial_mat_to_nii( 'test_', 'Test_out.nii' )

% serial_mat_to_nii( strInputFilestem, strOuputFilename )

% ------------------------------------------------------------------------------

%

% Read a series of 2-D matrices from .mat format into a 3-D volume and

% save as a NIFTI volume.

%

% Currently it is used for assembling the planes output by

% radial_components_image.m.

% This version is specific to a set of .mat files that

% contained an image in a structure field .im56fold.

`%%%%%%%%%%%%%%%%%%%%%%%%`

`% Hardcoded information:`

bShow = false; % [true, false] If true, displays each raw slice.

iPlanes = 1:435 % [ n : m ] z-coordinates of planes to be analyzed.

`%%%%%%%%%%%%%%%%%%%%%%%%`

Before (left) and after (right) rotational symmetry filtering. High contrast to emphasize DNA packing structures. The 5 and 6-fold filtering parameters are just a rough guess here. There are some flaws, for example the blobs at far right are artifacts because I neglected to include overtones of 5-fold symmetries that are also overtones of 6-fold symetries, even though there are negligable 6-fold components at that location.

Spectral powers above a threshold, for several sets of frequencies using spectral_segmentation_EMD_1220.m:

MATLAB commands:

>> load emd-1220_RPwr

>> spectral_segmentation_EMD_1220( 'Cross_section_half.png', fPwr )

where emd-1220_RPwr.mat was generated by radial_components_volume.m (see above), containing the fPwr 3-D array.

`%%%%%%%%%%%%%%%%%%%`

`% Hardcoded information:`

`bShow = false; % [true, false] If true, displays each raw slice.`

`bWrite = true; % [true, false] If true, write segmentation image.`

`strOutFilestem = 'test_seg_fundamental_6'`

`iChannels = [ 21:33 ]`

`% iChannels = [ 32 ]`

`%ndev = 5.0; % Segmentation threshold, n standard deviations from median cut.`

`ndev = 6.0; % Segmentation threshold, n standard deviations from median cut.`

`%%%%%%%%%%%%%%%%%%%`

Reperesentative cross section (left) compared with spectral segmentation (right).

Red: f = 5 or 10 cycles/annulus

Green and cyan: f = 6 cycles/annulus

Blue: f = 12 cycles/annulus

Magenta: f = 5 and 12 cycles/annulus

Yellow: f = 6 and 12 cycles/annulus

White: Other frequencies, primarily 7 and 13 cycles/annulus, and overlap of 5, 6 and 12 cycles/annulus.

The white in the bottom right corner is an artifact due to constant (zero) filling where there was no image data.

## Spectral density segmentation

2008-10-20

For second itteration of spectral density segmentation:

2008-10-24

- Use radial_components_volume.m to generate the f = 0 (DC, or constant) component volume.
- Output in 16-bit, format instead of 32-bit, for convenience.
- Add f = 14, 17, 18 and 19 components to spectral powers array.
- Pre-filter super-Nyquist frequencies?
- Use a map of significant spectral power for generating component volumes: imPwrSigMap( r, z, fa:fz ), where ( rn, zn, fa:fz )n <= 1

Done and tested. Need to add logic for overlapping regions -- what ratio of overtones to use.

2008-10-21

- Test generation of another frequency and harmonics volume. fa = n*7, fa = n*13, or fa = 1 (fundamental frequency only).

Made imPwrSigMap. Find clusters with at least one value == 1 and remove all other values (set to zero).

2008-10-24

- Fill center point of reconstructions properly

It is apparent that a large fraction of the non-5-6-12-fold power spikes (1-7-13-fold, larger white blobs in the spectral segmentation map above) are due to mis-estimation of the center(s) of rotation. Remember that the original (single) estimate was based on the capsid segmentation, all 5-fold.

Re-estimate the centers for the 5-fold, 6-12-fold and central (about circularly symmetric) regions independently. Use hill climbing and spectral decomposition. For which center is there the least residual power?

2008-11-03

Currently testing rotation_center_volume/image.m for center estimation by n-fold spectral power peak , about the longitudinal axis.

rotation_center_image.m calls radial_components_image.m

with the hardcoded parameter:

nCriterion = 2

2008-09-10

- Full final segmentation and mask.
- Separate thresholds for some or all structures, as separate 8-bit images.
- Wedge mask. Apply to elements.

Tested a 90 degree (volume corner) wedge, it works well. In Space Software, simpy create an all white volume and adjust its x/y origin (Volume|Volume Info) until its corner is near the origin of the virion volume object. Apply Overlays|Mask with this overlay and Overlays|Merge All.

- Render combinations.

My free Windows GUI
for display, navigation,
rendering, editing, and measurement of 3-D data sets.

Space was used for
some of
the processing (resampling, value scaling, etc), most of the image
slices, and all of the rendering.

Free linux based analysis tools intended for for FMRI, MRI and DTI brain imaging data, but the coregistration and volume math tools are suitable for use with other volume data. GUI and command line tools.

Used for
coregistration, application of some rotation transforms, and volume
averaging.

2008-09-12

To avoid writing code that allows finding the axis of rotation accurately (see above, iterated least squares self-coregistration fit about capsule's 5-fold symmetric axis, and head's six-fold rotational axis), I decided to try a well tested coregistration program. The FLIRT tool in the Linux FSL toolbox does volume coregistration with a variety of options, and also has functions that can resample, transform and average volumes.

To avoid writing code that allows finding the axis of rotation accurately (see above, iterated least squares self-coregistration fit about capsule's 5-fold symmetric axis, and head's six-fold rotational axis), I decided to try a well tested coregistration program. The FLIRT tool in the Linux FSL toolbox does volume coregistration with a variety of options, and also has functions that can resample, transform and average volumes.

As I'm using a
Windows
operating system, I needed to install VMware Player to
run the FSL
virtual machine. This was a huge hassle, taking a lot of time to
get everything working about right. [To Do: link to FSL
instructions, outline installation.]

Space Software still has some compatibility issues (and a bug or two) with respect to the NIFTI volume format. In particular, the first value of the dim[8] header item gets written as 1 (short integer), while it should be 3 (the number of dimensions of the data set. The work-around is to manually set this value in the Space Volume | Volume Information | View/Edit Header dialog box. Also, the header length appears to by off by 4 bytes, but this didn't crash FSLView (v3.0) or FLIRT

4x4 transforms can by applied to volumes using the FLIRT utility (flirt -applyxfm, -init and -out). This transform program can be accessed through the FLIRT GUI using the Util | Apply FLIRT Transform (ApplyXFM) button (lower right).

Space Software still has some compatibility issues (and a bug or two) with respect to the NIFTI volume format. In particular, the first value of the dim[8] header item gets written as 1 (short integer), while it should be 3 (the number of dimensions of the data set. The work-around is to manually set this value in the Space Volume | Volume Information | View/Edit Header dialog box. Also, the header length appears to by off by 4 bytes, but this didn't crash FSLView (v3.0) or FLIRT

4x4 transforms can by applied to volumes using the FLIRT utility (flirt -applyxfm, -init and -out). This transform program can be accessed through the FLIRT GUI using the Util | Apply FLIRT Transform (ApplyXFM) button (lower right).

MATLAB is an
expensive high-level language and interactive
environment that enables fast coding with matrix based data. I wrote
these MATLAB programs for several specific purposes.

convert_map_volume.m

Reads 32-bit
floating point volume data from a CCP4 format file, and
writes the 3-D array into an 8-bit unsigned Space Volume Type 1 or a
32-bit single file format NIFTI(.nii) file. Both formats are readable
by Space Software, and FSL and other programs read the NIFTI format.

[To Do: Add argument to convert_map_volume.m for desired output volume format. Document and include the NIFTI MATLAB tools I used: make_nii.m and save_nii.m

"Part of this file is copied and modified under GNU license from MRI_TOOLBOX developed by CNSP in Flinders University, Australia NIFTI data format can be found on: http://nifti.nimh.nih.gov

- Jimmy Shen (pls@rotman-baycrest.on.ca) "

2008-09-14

Rewrote convert_map_volume.m to include NIFTI 32-bit output option.

[To Do: Document and attach this new code.]

[To Do: Add argument to convert_map_volume.m for desired output volume format. Document and include the NIFTI MATLAB tools I used: make_nii.m and save_nii.m

"Part of this file is copied and modified under GNU license from MRI_TOOLBOX developed by CNSP in Flinders University, Australia NIFTI data format can be found on: http://nifti.nimh.nih.gov

- Jimmy Shen (pls@rotman-baycrest.on.ca) "

2008-09-14

Rewrote convert_map_volume.m to include NIFTI 32-bit output option.

[To Do: Document and attach this new code.]

Average n-fold
rotations of an image about a near-central point. The
location of the center is refined using a hill climbing algorithm that
maximizes the standard deviation of values in the average image.

[To
Do: Example usage.]

rotation_average_volume.m

rotation_nfold_image.m

See MATLAB code file, remove_background_peak_image.m.

[ vPeak fFWHM_range
imZeroBackground ] =
remove_background_peak_image( strInputFilename, fFWHM, bShow )

USAGE: [ vPeak fFWHM_range imZeroBackground ] = remove_background_peak_image( 'test_remove_background.png', 1.5, true );

USAGE: [ vPeak fFWHM_range imZeroBackground ] = remove_background_peak_image( 'test_remove_background.png', 1.5, true );

For an image with a relatively uniform noisy background (resulting in a luminance histogram peak), find the peak luminance histogram hump and the range of values above a given fraction of the full-width at half-maximum FWHM. For example, with fFWHM = 1.5 the width of the range of values that is set to zero is 1.5 times the derived FWHM.

For an image with a relatively uniform noisy background (composing a luminance histogram peak), find the peak luminance histogram hump and the range of values above a given fraction of the full-width at half-maximum FWHM. For example, with fFWHM = 1.5 the width of the range of values that is set to zero is 1.5 times the derived FWHM.

2008-09-24

I'm not going to use this for this project, although this is why wrote it. I think it will be easier to mask the backgound on the whole volume using the Space Software 3-D Fill tool. This avoids dealing with interior pixels/voxels that this program removes from the interior.

I'm not going to use this for this project, although this is why wrote it. I think it will be easier to mask the backgound on the whole volume using the Space Software 3-D Fill tool. This avoids dealing with interior pixels/voxels that this program removes from the interior.

Example usage:

This example uses a hardcoded parameter used to remove the uniform image area (at bottom of original image below) before finding the histogram peak:

valueIgnore = 125; % This values will be set to zero before processing and in output.

This example uses a hardcoded parameter used to remove the uniform image area (at bottom of original image below) before finding the histogram peak:

valueIgnore = 125; % This values will be set to zero before processing and in output.

MATLAB command and command line output:

>> remove_background_peak_image( 'test_remove_background.png', 1.5, true );

>> remove_background_peak_image( 'test_remove_background.png', 1.5, true );

szInput =

339 333

Peak value: 128.5

Low range value: 118.3749

High range value: 137.8805

Original and final images, and intermediate (image and histogram) outputs:

original (input) image

final (output) image

histogram of original after ignored value removed

histogram of final image

image after ignore value set to zero

2008-10-15

Central and below axial center
rings don't fit the rest of the DNA spool pattern. What are they and
how do they just float, apparently disconnected?

2008-10-15

Some of the peripheral "stacked rings" of
DNA are very robustly detected, indicative of "spooling". But the
non-uniform breaks and bridges between the rings also appear to be
real. How could this be true for a rotationally averaged image? There
doesn't seem to much rhyme or reason - not correlated across concentric
rings or capsid structures.

2008-10-15

Central axial structures appear to
be real but not "lead DNA". There's a hyper-density plug, and a
hypo-density cap.

2008-10-08

What do small Fourier components
that are not multiples of 5 or 6 mean? Could there be admixtures of,
say, 11-fold tail structures, and would these show up as 11-fold
components? Would their phase be coherent if there were an admixture?
Can non-harmonic distortions be distinguished from real n-fold
regularities?

2008-10-07

The 12-fold portion of the tail is
fit into a 5-fold vertex of the capsid, like a square peg in a round
hole. The Fourier decomposition shows the boundary has a narrow region
where both annular frequencies are significant. If the 5-fold and
6-fold symmetries are independently reconstructed, are the overlapping
features due to distortion of one or both structures, or are the
"electron densities" (or whatever is measured by EM) of both structures
are interpenetrating? Can these two mechanisms be distinguished in some
ideal case?

2008-11-25

I got
distracted fom this project in making an animation that I naimed "Mixed
symmetries", because I combined some 3-D periodic, 3-D perceptual
rotation, and 2-D periodic symmetries. After I posted it I googled
"mixed symmetries", yielding some physics energy level papers. But the
top hit on the "images" search with that phrase was a rendering of a Caudovirus, Epsilon 15, a virus
that infects the bacterium Salmonella based on Cryo-EM reconstruction!
There's a good 2006 paper [To Do: Refind this link in
"Virus_structure_and_mechanics_information"] with very nice
illustrations that explores just this question. Quite a coincidence as
I had not looked into this close relative of P22, and it has obvious
distortions of molecular configuration at the 6/5-fold boundary. I
though that there would be a long list of topics that used the phrase,
but the top two happen to be obscure side topics I've worked on.

2008-09-24

Is the axis of the capsid alligned with
that of the tail?

The allignment was based
on the capsid symmetry. I expect the capsid center derived from each
axial slice will be close to parallel with the volume's z-axis. But the
tail structures might be skewed due to the 5 to 6/12-fold interlocking.

2008-09-24

Is there a distortion of the icosohedral
capsid tail vertex due to the tail?

After the full n-fold
segmentation is done, a difference with one (or an average of several)
vertex will demonstrate any distortion.

2008-10-01

Rotional symmetry program, for EMD-1220:

I've looked closely at how
to symmetrically average planes normal to the axis, particularly near
the boundary of 5 and 6-fold symmetry. This tail-capsid boundary has
components of both symmetries, and the image can be partitioned with a
simple linear model.

The tail-capsid interface is particularly interesting with respect to molecular machines. [2008-10-15 1950's molecular machine book by ???von Hippel] The dsDNA flows through the axis of the tail in both directions, while being packaged and while being inserted in a new host cell. The mechanism of both processes is not known. It is likely that the central tail structure is the basis of a molecular motor. It has been suggested that its incommensurate 5 and 6-fold boundary at the capsid vertex might form a kind of ratchet.

"Viral DNA packaging motors inject viral genomic DNA into capsids as part of their replication cycle, packing it very tightly. [5]"

The tail-capsid interface is particularly interesting with respect to molecular machines. [2008-10-15 1950's molecular machine book by ???von Hippel] The dsDNA flows through the axis of the tail in both directions, while being packaged and while being inserted in a new host cell. The mechanism of both processes is not known. It is likely that the central tail structure is the basis of a molecular motor. It has been suggested that its incommensurate 5 and 6-fold boundary at the capsid vertex might form a kind of ratchet.

"Viral DNA packaging motors inject viral genomic DNA into capsids as part of their replication cycle, packing it very tightly. [5]"

1. Roll by | |

2. Pitch by | |

3. Yaw by | |

4. Translate by |

------------------------------------

There are no restrictions on use of the images or code on this page. Claiming to be the originator of the material, explicitly or implicitly, is bad karma. A link (if appropriate), a note to dow[at]uoregon.edu, and credit are appreciated but not required.

Comments are welcome (dow[at]uoregon.edu).