Mark Dow

Projects

Virions

P22 virion structure

Test renderings of the P22 virion EMD-1220 data set

EMD-1220 rendering notes

Data source

Processing notes

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.]

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
    Write a specialized MATLAB program for radial symmetries analysis and n-fold rotational averaging. Started rotation_average_volume.m.

2008-09-17
    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:

• 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.

2008-10-01
    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.

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.

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

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

Rotational symmetry decomposition

2008-10-01

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:

• 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

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

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

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

• Fill center point of reconstructions properly

2008-10-24
    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

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

2008-09-10
    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.

Programs and code

    Space Software

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

    FMRIB Software Library (FSL)

    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.

    MATLAB code

remove_background_peak_image.m

original (input) image	final (output) image
histogram of original after ignored value removed	histogram of final image
image after ignore value set to zero

Other structural questions

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:

The homogeneous (coordinate) transformation matrix for 3D points

1. Roll by
2. Pitch by
3. Yaw by
4. Translate by

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

 This work is dedicated to the Public Domain.
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