Deformed alignment of super-resolution images for semi-flexible structures

Due to low labeling efficiency and structural heterogeneity in fluorescence-based single-molecule localization microscopy (SMLM), image alignment and quantitative analysis is often required to make accurate conclusions on the spatial relationships between proteins. Cryo-electron microscopy (EM) image alignment procedures have been applied to average structures taken with super-resolution microscopy. However, unlike cryo-EM, the much larger cellular structures analyzed by super-resolution microscopy are often heterogeneous, resulting in misalignment. And the light-microscopy image library is much smaller, which makes classification challenging. To overcome these two challenges, we developed a method to deform semi-flexible ring-shaped structures and then align the 3D structures without classification. These algorithms can register semi-flexible structures with an accuracy of several nanometers in short computation time and with greatly reduced memory requirements. We demonstrated our methods by aligning experimental Stochastic Optical Reconstruction Microscopy (STORM) images of ciliary distal appendages and simulated structures. Symmetries, dimensions, and locations of protein complexes in 3D are revealed by the alignment and averaging for heterogeneous, tilted, and under-labeled structures.


Introduction
In the past decade, the development of localization-based super-resolution microscopy has brought the light microscopy to nanometer scales. Imaging beyond the diffraction limit addressed structural and functional biomedical questions of subcellular organelles that could not be resolved by conventional light microcopy. For example, the in situ dissection of macromolecular protein complexes includes ciliary transition zone [1][2][3], neuronal synapses [4], nuclear pore complex [5,6], focal adhesion complex [7], clathrin-coated pits [8], centrosome [9,10], and the escort complex at viral budding sites [11].
To obtain optical super-resolution images, the target proteins or DNA/RNA sequences are fluorescently labeled by organic dyes or fluorescent protein tags. However, the targets are usually under labeled. The cause of the low labeling efficiency includes low affinity antibodies, dye PLOS  quenching, immature fluorescent proteins, and less-than-one dye to antibody conjugation ratios, etc. In certain cases, it might not be possible to achieve full labeling experimentally. Consequently, image alignment and averaging are often required to make accurate conclusions about the biological system of interest. Template-free cryo-EM image alignment procedures have been applied and adapted to average structures taken with super-resolution microscopy in 2D [12][13][14], and in 3D [15]. Unlike cryo-EM, the much larger cellular structures analyzed by super-resolution microscopy are often not completely rigid and homogenous. For instance, the diameter of the rings of ciliary distal appendages varies from 369 to 494 nm [1]. And the shape for the rings are often elliptical due to the flexibility of cilia. In general, many cellular organelles larger than 200 nm are semi-flexible. These semi-flexible structures keep the common symmetry and angular arrangement, but have heterogeneous size and shape, either because of elastic deformation or because of variations in molecular composition. A wide range of cellular structures have such characteristics, including centrioles [9], ciliary transition zones, ciliary distal appendages [1][2][3], and pre-and post-synapse density [4]. In this case, direct alignment and averaging will lose structural details. Cryo-EM particle alignment deals with structural heterogeneity using classification and class-averaging [16][17][18]. However, super-resolution structure analysis often has much fewer images to start with for the following two reasons. First, there are only a few target organelles in a cell. For instance, usually one embryonic fibroblast cell only develops one primary cilium [1], and each cell only has two centrioles [9]. Second, different from EM, superresolution light microscopes detect fluorescence labels instead of electron density. Many structures seriously under-labeled cannot be used for the direct alignment and averaging. The small number of usable super-resolution images is often not enough to meet the minimal amount required for classification. Consequently, heterogeneity in semi-flexible structures will cause misalignment in direct rigid-registration used in existing algorithms for super-resolution image alignment [12][13][14][15], leading to substantial degradation of the resolution of the aligned images.
To address these problems, we take structural flexibility as one degree of freedom for image registration. We designed algorithms to first deform the semi-flexible structures to fairly uniform shape and size, and then to align all deformed structures based on cross correlations. As characterized by Fourier Ring Correlation (FRC) analysis, we have shown that our deformed alignment algorithm achieved better aligned image resolution for semi-flexible structures compared to the state-of-the-art rigid registration algorithm [12]. Although here we specifically developed our deformed algorithm for the ciliary transition zone and distal appendages [1], this procedure could be generalized for many other subcellular structures such as the nuclear pore complex and the centriole.

Overview of the workflow
We designed two algorithms to align 2D and 3D semi-flexible ring structures, respectively. The 2D deformed alignment algorithm is suitable for heterogeneous ring-shaped structures in which the heterogeneity is mainly caused by the flexibility of the structure, for example, transition zone and centrioles. The algorithm first deforms individual structures by circulating the structures, and then aligns the images. The 3D alignment algorithm was designed for randomly oriented flat structures in 3D. The random orientation causes the heterogeneity in the x-y projection. This 3D algorithm rotates the structures in 3D, and aligns the projection on the x-y plane, without structure deformation. Both deformed alignment algorithm and 3D alignment algorithm share the same alignment metric, the cross correlation with the reference in the frequency domain. And the initial reference is just the average image of the structures to be aligned. After iterations of alignment, all the images in the frequency domain are inverse Fourier transformed to the space domain. The final alignment information is used to deform, translate and rotate the coordinates in the molecule lists of the super resolution images. Both deformed alignment algorithm and 3D alignment algorithm achieved sub-pixel precision and fast alignment for linear computing time. On top of the single-color alignment algorithms, we facilitated multicolor alignment by applying the 3D alignment parameters for the reference channel to other imaged channels.
We aligned experimental data, single-color STORM images, to validate and demonstrate the utility of our deformed 2D alignment algorithm, and 3D algorithm. We used simulated structures to further test these algorithms' capabilities to handle various image resolution, labeling efficiency, and structural complexity. In addition, we validated our two-color alignment algorithm using simulated two-color images. Because the quality of the alignment performance is reflected in the resolution of the average image of aligned super-resolution structures we employed the FRC resolution [19] to quantify the performance of different alignment methods. All software codes and sample data in this manuscript are available at https://github. com/Huanglab-ucsf/Deformed-Alignment.

STORM image acquisition
We have described the technical details of the STORM image acquisition and reconstruction methods in our previous work on ciliary transition zone [1]. The same instrument and software were used to acquire the STORM data for algorithm testing in this work. Briefly, the photoswitchable dye pair Alexa 405/Alexa 647 was used for STORM imaging with a ratio of 0.8 Alexa 647 molecules per antibody. During the imaging process, the 405 nm activation laser was used to activate a small fraction of the Alexa 647 at a time, and individual activated fluorophores were excited with a 642 nm laser. The typical power for the lasers at the back port of the microscope was~1 kW/cm 2 for the 642nm imaging laser and 0-20 W/cm 2 for the 405 nm activation laser.
Analysis of STORM raw data was performed in the Insight3 software, which identifies and fits single molecule spots in each camera frame to determine their x, y and z coordinates as well as photon numbers. Sample drift during data acquisition were corrected using imaging correlation analysis. The lateral and the axial localization precision was calculated to be 17 nm (standard deviation).

Coordinate-based image deformation
Taking structural flexibility as one degree of freedom for image registration, we deform the semi-flexible structures before aligning the images (Fig 1A). Compared to pixel-based images, it is much easier to apply complex deformation functions to coordinate-based images. Fortunately, the data format of localization microscopy is a list of coordinates. Each coordinate (x, y, z) for a 3D image or (x, y) for a 2D image locates a fluorophore. First, we used robust fitting to fit the coordinates of individual structures to ellipses described as where the algebraic parameters A to E are converted to geometric parameters center (x 0 , y 0 ), long axis a, short axis b, and the rotation angle f. Robust fitting is an alternative to least squares fitting when data are contaminated with outliers. It is more resistant to outliers in the data because it uses least absolute deviations rather than least squares. Because STORM images are composed of noisy localizations, robust fitting should work better than least square fitting to find a function which closely approximate the localizations. In Fig 1B-1D, we compared the ellipse fittings of a STORM image of ciliary distal appendages with robust fitting and with least square fitting, respectively. Fig 1B is a STORM image of the distal appendages of a cilium. The structure looks like an ellipse with the long axis titling towards northwest. The robust fitting results in an ellipse with the long axis titling towards northwest (red ellipse in Fig 1C), agreeing with the structure. However, the least square fitting gives an ellipse with the long axis slightly titling towards northeast (blue ellipse in Fig 1D), disagreeing with the orientation of the structure. This comparison demonstrated the advantage of robust fitting for STORM data analysis. In this algorithm, coordinates out of 1.5 standard deviations are removed to clean the images. With the geometric parameters obtained from the ellipse robust fitting, we center the structure to (x 0 , y 0 ) and rotate all the coordinates by-f around the center. Finally, we circularize the structure by simply scaling (x, y) with Eqs [2] and [3], where a is the long axis, b is the short axis, and R is the average radius of all structures to be aligned (Fig 1A).

Fast image alignment using correlation in the Fourier domain
Different from our previously reported alignment framework using coordinate-based correlation, this work employs pixel-based correlation and DFT to accelerate computation speed. The DFT algorithm is a modification of the efficient subpixel image registration algorithm written by Guizar-Sicairos et al [20]. that achieved efficient subpixel image registration by up-sampled DFT cross-correlation. This algorithm registers images using 2D rigid translation. We considered sample rotation and implemented rotational registration ( Fig 1M). Briefly, the original algorithm first obtains an initial 2D shift estimate of the cross-correlation peak by fast Fourier transforming (FFT) an image to register to a reference image. Second, it refines the shift estimation by up-sampling the DFT only in a small neighborhood of that estimate by means of a matrix-multiply DFT. We added a loop to rotate the images to align from 0 to 359˚by 1˚at a step, and then performed the original first step after each rotation step to obtain 360 cross-correlation peaks. The biggest peak provides not only the initial 2D shift estimate but also the optimal angle of rotation. The images of single structure were aligned for ten iterations. To eliminate the bias that can be introduced by a template, we used the average of all the images to be aligned as the reference image for the first iteration. Each image was translated and rotated to the reference image to maximize the cross-correlation between the image to be aligned and the reference image in the Fourier domain. The average of the images aligned in one iteration was then served as the reference image for the next iteration. Our algorithm aligned the images with a precision of 1/100 of a pixel. The up-sampling rate can be adjusted higher without increasing computation time. We quantify the alignment using normalized root-mean-square error (NRMSE) E between f(x, y) and rotated g(x, y), defined by [4] E 2 ¼ min a;x 0 ;y 0 ;y P x;y ja � g y ðx À x 0 ; y À y 0 Þ À f ðx; yÞj where summations are taken over all image points (x,y), α is an arbitrary constant, g θ (x, y) is the image g(x, y) rotated by angle θ 0 as in Eq [5] x 0 Finding the x 0 , y 0 , and θ 0 for the minimum NRMSE is equivalent for the maximum crosscorrelation r fg , defined by [6] r fg x 0 ; y 0 ; y ð Þ ¼ P x;y f ðx; yÞg � y ðx À x 0 ; y À y 0 Þ ¼ where N and M are the image dimensions, ( � ) denotes complex conjugation, F(μ, ν) and G θ (μ,

3D alignment
Based on the 2D alignment algorithm, we developed a 3D alignment algorithm finding the minimum NRMSE (maximum cross-correlation) while rotating and translating the structures in 3D. According to Euler's rotation theorem, combinations of rotations around any two axes can reproduce the rotation around the third axis. For example, any rotation around y axis is equal to the combination of sequential rotations around z by 90˚, around x axis, and around z by -90˚. So, the 3D alignment algorithm adds only one for-loop out of the 2D alignment algorithm described above and excludes the deformation step. This for-loop rotates each individual structure around its x axis. The localization coordinates (x, y, z) are rotated by angle f in a preset range. The coordinates are transformed by operation

5: ½8�
The structures are aligned finding the x 0 , y 0 , θ 0 , and f 0 for minimum NRMSE or maximum cross-correlation r fg .

Two-color alignment
The two-color alignment includes two steps. We first picked the channel with higher resolution as the reference channel and aligned this channel with the algorithm described above. Second, the alignment parameters for the reference channel were applied to the other channel.

Coordinate-based image deformation and alignment
By deforming the super-resolution images, we correct the heterogeneity in semi-flexible structures which cause misalignment in rigid registration. The deformation improves the resolution and symmetry in the aligning result. With experimental data, we demonstrated this improvement, by comparing the average image of structures aligned by the deformed alignment algorithm (Fig 1H) to the aligning result of the rigid registration algorithm recently developed in our group [12] (Fig 1J). The FRC analysis quantified the resolution improvement made by the coordinate-based image deformation (Fig 1K & 1L). The FRC resolution of the average image of structures aligned with deformation is 41.4 ±1.2 nm, while the FRC resolution for the rigid registration is 45 ± 2.2 nm.
The algorithm was demonstrated with the experimental 2D STORM data of ciliary distal appendages of mouse tracheal epithelial cells (MTECs) with CEP164 labeled by Alexa Fluor 647. The deformed 2D alignment result (Fig 1H) of 31 under-labeled structures with semi-flexible shape and size (Fig 1A) showed clear 9-fold symmetry. The elongated distribution of CEP164 in each distal appendage was also represented in individual subunit of the average image. The elongated distribution of CEP164 in each subunit agrees with the fiber shape of the distal appendage imaged by EM [21]. 30 images can be aligned in a few minutes with a CPU at 2.8 GHz and 24 GB memory (our lab computer). And the computation time is linear to the number of the images to be aligned.
To test the capability of our algorithm to align images with various localization precision, labeling efficiency, and high structural complexity, we simulated three sets of localization images. (1) With the localization precision of 60 nm, we simulated twenty rings, each of which consists of 9 evenly distributed clusters with the diameter of 300 nm (Fig 2A average structure,  B single structure). The algorithm was able to resolve 9 clusters clearly by aligning simulated structures with 10 integrations in 363 seconds on a server with 2.8 GHz CPU and 24GB memory (Fig 2C). We used autocorrelation of aligned image as a function of the angle of rotation to quantify how well the symmetry is retrieved by the alignment. The peak at 38˚in the autocorrelation function (Fig 2J) agreed well with the 9-fold symmetry in the simulated images. (2) To test the alignment of structures with low labeling efficiency, we simulated 20 9-cluster rings with 5 random clusters labeled at the localization precision of 15 nm (Fig 2D average structure, E single structure). The average image of aligned under-labeled structures recovered clear 9 clusters (Fig 2F). The peak at 40˚in the autocorrelation function ( Fig 2K) retrieved with the 9-fold symmetry in the simulated images. (3) We also evaluated the alignment accuracy with simulated structures with higher structural complexity. We simulated 20 rings with 18 clusters which have alternative 15˚and 25˚angular spacing between the neighboring clusters ( Fig 2G  average structure, H single structure). The simulated labeling efficiency is 67%. The average image of aligned structures successfully represented the complex symmetry (Fig 2I). The peaks at 15˚, 25˚and 40˚in the autocorrelation function (Fig 2L) again perfectly retrieved the complex angular distribution of simulated structures.

3D alignment
We used the algorithm to align 29 experimental 3D STORM images of ciliary distal appendages with long axes varying from 376 to 470 nm (Fig 3B). These in situ structures were randomly oriented in 3D, and their labeling efficiency was about 60-85%. Before alignment, the average structure showed a ring shape but with no information on the angular distribution of the clusters. After 5 iterations of 3D alignment, the average image showed 9 clusters with approximately even angular distribution around the ring.
To test the capability of alignment of largely tilted structures, we simulated ten rings with 9 clusters. These simulated structures are randomly rotated around the x axis within 60˚and around the z axis within 360˚, at a localization precision of 15 nm (Fig 3F). The average of simulated ring images before alignment shows no information about the symmetry of the structures (Fig 3G). The 3D alignment result shows 9-fold symmetry (Fig 3H). The algorithm also provides good alignment results for simulated structures with a large tilting angle (90˚) around x-axis (Fig 3I and 3J), and for large localization precision (60 nm) in 3D (Fig 3K and 3L). These results indicate that our 3D alignment algorithm can accurately register largely titled and noisy images, and can efficiently extract the common structural pattern among these heterogeneous images.

Two-color alignment
To test the algorithm of two-color alignment, we simulated twenty two-color super-resolution images with 15 nm lateral localization precision and 30 nm axial localization precision ( Fig  4A). Each simulated structure is composed of two parallel rings with diameters of 200 and 300 nm, respectively. Each ring has 9 clusters that are evenly distributed. The average number of localizations in each cluster is 50. The angular distributions of the clusters of the smaller rings and the bigger ring are offset by 20˚. The two rings are 100 nm distal to each other. Our alignment algorithm successfully aligned the images. The aligned average image shows 9-fold symmetry and the 20˚angular offset between the two rings in the top view (Fig 4C), and 100 nm distal distance between the two rings in the side view (Fig 4E).
To validate the algorithm's capability of aligning two-color images that are tilted in 3D, we randomly rotated the structures used in the case above around x axis in a range of 30˚, and around z axis in a range of 360˚ (Fig 4F). The 3D alignment algorithm again efficiently aligned the titled structures and recovered the 3D arrangement in the average structure (Fig 4H  and 4J).

Conclusions
We have demonstrated a deformed 2D alignment algorithm that can accurately align semiflexible ring structures from both experimental STORM images of ciliary distal appendages and simulated images with high noise and complexity. Information on symmetry and common structural features were efficiently extracted from a few tens of heterogeneous structures. The cross-correlation based alignment algorithm is largely accelerated by DFT and registration in the frequency domain. We also demonstrated that our 3D alignment algorithm can accurately align images of structures tilted in space, using 3D STORM images of ciliary distal appendages, and simulated structures randomly rotated around x axis. For two-color alignment, we simply applied the alignment parameters for the reference channel to the other channel, and achieved sub-pixel alignment of two-color super-resolution images. The two-color algorithm can be easily expanded to multiple color alignment. Use of these algorithms makes accurate registration of super-resolution images within a hundredth of a pixel. Information on the general geometric feature of heterogeneous structures can be extracted with tens of images. All the algorithms including deformed 2D, 3D, and two-color alignments are computationally manageable on a regular desktop computer. The deformed 2D algorithm can be applied to any semi-flexible ring-shaped structures. The 3D and multicolor algorithms will provide substantial advantage for any applications that require aligning and averaging super-resolution images in 3D.