View source for Walkthrough on hexagonal lattice coherence through neighborhood analysis (I2PC – Instruct Course on Electron Tomography and Subtomogram Averaging 2025)

[[Category:WalkthroughsSpecific]]


==Neighborhood analysis==

Neighbourhood analysis stands at the forefront of methods used to investigate geometrical lattice parameters and offers a wide spectrum of visualization tools. This technique is precise, targeting specific elements, yet broad enough to facilitate overarching interpretive analysis. At its core, the method constructs a latent space of connections by mapping the relative positions of neighbouring particles onto a reference framework. This mapping involves a reorientation process, wherein the 'N' closest neighbours are repositioned within the reference particle's coordinate system, adjusted for its orientation. The richness of the latent space—characterized by the number and nature of these repositioned features—is dictated by a defined distance threshold. This threshold ranges from a minimal couple, representing the closest inter-particle distance, to the maximal extent, linking every particle within the neighbourhood to the reference point.

For additional insights into neighborhood analysis, visit [https://www.dynamo-em.org//w/index.php?title=Neighborhood_analysis this page.]  

==Data==

We will use data from the initial gold standard alignment derived from five tomograms of immature HIV-1 virus-like particles.

In this course, the data for this walkthrough is in:
 <tt >/home/student/data/HexagonalLattice/ </tt>
the tomogram with the data is in the link:
 <tt > /home/student/data/HexagonalLattice/template.em </tt>
the table linked with the data is in the link:
 <tt > /home/student/data/HexagonalLattice/table.tbl </tt>

===Inspecting data===

First, open a new local terminal and navigate to the directory containing the tutorial dataset and project.

Let's inspect the contents:

 >> dfile template.em
    
    filetype: volume
    
    size: 192 x 192 x 192 

 >> dfile table.tbl
  
    filetype: table
    
    size: 7068  x 35, aligns: 7068, averages: 7068 

These files represent an average as a template and its associated metadata (i.e., table). We will focus on the table data.

To visualize the particles along with their normals:

 tbl = dread('table.tbl');
 figure;
 h = dpktbl.plots.sketch(tbl,'haxis',gca());
 h.zlength.value=100;


[[File:tbl-positions-normals-vlps.png|centre|thumb|450px| VLPs particle positions and normals.]]

== Prepare (meta)data ==

Let's select a couple of Virus-Like Particles (VLPs):

 tblT2V2 = tbl(tbl(:,20) == 2 & tbl(:,22) == 2, :);
 tblT2V5 = tbl(tbl(:,20) == 2 & tbl(:,22) == 5, :);
 tblT2V3 = tbl(tbl(:,20) == 2 & tbl(:,22) == 3, :);

To visualize:

 figure;
 
 subplot(1, 3, 1);
 dpktbl.plots.sketch(tblT2V5);
 title('Ideal');
 pbaspect([1 1 1]);
 
 subplot(1, 3, 2);
 dpktbl.plots.sketch(tblT2V2); 
 title('Slightly flat surface');
 pbaspect([1 1 1]); 
 
 subplot(1, 3, 3);
 dpktbl.plots.sketch(tblT2V3); 
 title('1/4 missing, flat top');
 pbaspect([1 1 1]); 

[[File:selected-vlps.png|centre|thumb|800px|  Selected VLPs. From left to right: ideal, medium, and low quality instances.]]

Indeed, we will analyze three different VLP quality instances.

== Get mean distance to closest neighbor ==

To find the optimal distance threshold, we will use the following function to calculate the mean distance to the closest neighbor:

 function [meanDst, stdDsts] = mean_distance_to_closest_neighbor(pts);
 numPts = size(pts, 1);
 minDsts = zeros(numPts, 1);
 
 for i = 1:numPts
    dsts = pdist2(pts(i,:), pts);
    dsts(i) = Inf;
    minDsts(i) = min(dsts);
 end
 
 meanDst = mean(minDsts);
 stdDsts = std(minDsts);
 fprintf('Mean distance: %.4f\nStandard deviation: %.4f \n', meanDst, stdDsts);
 end

Invoke the function like so: 

 [meanDst, stdDsts] =   mean_distance_to_closest_neighbor(tbl(:,24:26)+tbl(:,4:6));

The mean distance to the closest neighbor is around 50-60 pixels.
 Mean distance: 50.42
 Standard deviation: 12.23

For a sanity check, you might also measure the lattice directly on the average volume:

 dshow template.em

[[File:average-overplot-hexagonal-lattice-measurements.png|centre|thumb|400px|  Overplot average and hexagonal lattice measurements.]]

== Neighborhood analysis computation ==

Here, we conduct the neighborhood analysis with a distance threshold set as the mean distance plus the standard deviation:

 tblT2V3N = dpktbl.neighborhood.analize(tblT2V3,'distance',meanDst+stdDsts,'modelColumn',22);
 tblT2V2N = dpktbl.neighborhood.analize(tblT2V2,'distance',meanDst+stdDsts,'modelColumn',22);
 tblT2V5N = dpktbl.neighborhood.analize(tblT2V5,'distance',meanDst+stdDsts,'modelColumn',22);

Then, we will need the couple table since it is our new representation (i.e., latent space):

 tblT2V3NC = real(tblT2V3N.coupleTable);
 tblT2V2NC = real(tblT2V2N.coupleTable);
 tblT2V5NC = real(tblT2V5N.coupleTable);

Now, let's visualize first-order distance clusters with synchronized viewing:

 datasets = {tblT2V3NC, tblT2V2NC, tblT2V5NC};
 titles = {"Incomplete + flat surface", "Slightly flat surface", "Ideal"};
 
 fig = figure;
 ax = gobjects(1, length(datasets));
 
 for i = 1:length(datasets)
    ax(i) = subplot(1, length(datasets), i, 'Parent', fig);
    pcshow(datasets{i}(:, 24:26), 'Parent', ax(i));
    pbaspect(ax(i), [1 1 1]);
    % Set titles using text, since the title function was not working
    xl = xlim(ax(i));
    yl = ylim(ax(i));
    zl = zlim(ax(i));
    titlePos = [mean(xl), yl(2), zl(2)];
    text(ax(i), titlePos(1), titlePos(2), titlePos(3), titles{i}, ...
        'HorizontalAlignment', 'center', 'Color', 'white', ...
        'FontSize', 14, 'FontWeight', 'bold', 'BackgroundColor', 'black');
 end
 
 link = linkprop(ax, {'CameraUpVector', 'CameraPosition', 'CameraTarget', 'CameraViewAngle', 'XLim', 'YLim', 'ZLim'});
 setappdata(gcf, 'StoreTheLink', link);
 
 for a = ax
    rotate3d(a, 'on');
 end

[[File:first-order-distance-clusters-three-different-vlps.png|centre|thumb|800px|   First-order distance clusters obtained vía neighborhood analysis. From left to right: low, medium, and high quality instances.]]


The color, by default, represents the Z-axis. Notice that the better the quality of the VLP, the less sparse each cluster is.

== Local analysis - lattice coherence ==

To quantify lattice coherence, we apply [https://www.youtube.com/watch?v=R2e3Ls9H_fc k-means clustering] and calculate the Euclidean norms for each centroid. The standard deviation of these norms provides a measure of cluster sparsity and hence, lattice coherence.

We encapsulate this process in a function:

 function [idx, cc_final, results, results_std] = perform_lattice_analysis(data, symm)
    [idx, cc_final] = kmeans(data, symm);
    results = vecnorm(cc_final, 2, 2); 
    results_std = std(results); 
 end

Then, execute the analysis:

 symm = 6; 
 [idx_ideal, cc_ideal_final, results_ideal, std_ideal] = perform_lattice_analysis(tblT2V5NC(:,24:26),symm);
 [idx_slightly_flat, cc_sf_final, results_sf, std_sf] = perform_lattice_analysis(tblT2V2NC(:,24:26),symm);
 [idx_incomplete_and_flat, cc_implt_final, results_implt, std_implt] = perform_lattice_analysis(tblT2V3NC(:,24:26),symm);

Display results neatly in a table

 resultsTable = table(std_ideal, std_sf, std_implt, 'VariableNames', {'Ideal', 'SlightlyFlat', 'IncompleteFlat'}, 'RowNames', {'StandardDeviation'});
 
 disp(resultsTable);

This will output:
        
                               Ideal             SlightlyFlat        IncompleteFlat                       
                         __________________    _________________    _________________
                           
    Standard deviation    0.0933777687916592    0.321031652546273    0.574732108201886

As demonstrated, the qualitative degree of imperfection correlates with our lattice coherence measurements.