SMLocalizer, a GPU accelerated ImageJ plugin for single molecule localization microscopy

Abstract Summary SMLocalizer combines the availability of ImageJ with the power of GPU processing for fast and accurate analysis of single molecule localization microscopy data. Analysis of 2D and 3D data in multiple channels is supported. Availability and implementation Plugin freely available for Fiji and ImageJ2.0 through https://sourceforge.net/projects/smlocalizer/. Plugin also available for continuous updates through ImageJ update system, add http://sites.imagej.net/Cellular-Biophysics-KTH/ as update site in ImageJ. Java and CUDA source code freely available on the web at https://github.com/KristofferBernhem/SMlocalizer. Supplementary information Supplementary data are available at Bioinformatics online.


OS and hardware requirements
SMLocalizer will run on any system capable of running ImageJ or Fiji. At time of writing this includes (taken from http://imagej.net/Downloads) : ImageJ will run on any system that has a Java 8 (or later) runtime installed. This includes, but is not limited to: 1. Windows XP, Vista, 7 or 8 with Java installed from java.com 2. Mac OS X 10.8 "Mountain Lion" or later with Java installed from java.com 3. Ubuntu Linux 12.04 LTS or later with OpenJDK 8 installed SMLM images are typically large, as such, SMLocalizer requires that the java heap space memory be set high enough for the files the user intends to load. Typical values should be 3x the raw data file size. Java heap space memory is set in /Edit/Options/Memory & Threads.

GPU acceleration
SMLocalizer by default runs all processes on the system available CPU cores. If a compatible NVIDIA card is available, the user can select to transfer functions to the GPU instead. For GPU accelerated computation SMLocalizer 2.
x.x requires a NVIDIA GeForce gtx970 or later graphics card. Older version may work but none have been tested. The main limitation is memory available, with the current code expecting 4 GB of memory available. Ensure that the latest drivers for the graphics card are installed. For windows, timeout can be an issue for larger data files, this can be fixed through changing the graphic device timeout to 10+ seconds, see http://answers.microsoft.com for details on this.
The OS specific CUDA archives needs to be places in \plugins\jars directory of the main ImageJ directory (or into \jars for Fiji). Download the Fiji or ImageJ specific .rar and extract it into the main directory.

INSTALLING SMLOCALIZER
In ImageJ or Fiji, go to the help menu and select Update (see fig 1). This will after a checking of your current plugins take you to the updater window.
Select Manage update sites and Add update site. A new line will be added to the list of update sites (see fig 2). Replace the name New with SMLocalizer and add the URL: http://sites.imagej.net/Cellular-Biophysics-KTH/. ImageJ will now keep your SMLocalizer up to date.

PROCESSING SMLM DATA
The following section aim to explain the details required to properly use SMLocalizer. For a more detailed description of the algorithms, see section 5.
As ImageJ does not have a clear way to display processing progress, the keys used for the different algorithms will remain depressed whilst their respective algorithms are running.

PROCESSING
SMLocalizer can fit one-color 2D experiments without calibrations by:

LOADING IMAGES AND SETTINGS
Loading images for processing by SMLocalizer is done through the main ImageJ or Fiji interface. Some image formats lack compression and require additional space. A suggestion for performance optimization is to resave these as .tif stacks before proceeding.
Once installed and an image stack has been loaded, the plugin can be started. In the Plugins menu, select SMLocalizer to start the main plugin. Upon startup SMLocalizer will load the settings last stored (or default if no new settings have been stored). New settings can be stored or loaded using the two buttons in the bottom right corner (see fig 3).

CALIBRATION
SMLocalizer require bead images for calibration of 3D fits and can use multispectral beads for calibration of channel offset in x-y (2D) or x-y-z (3D). For 2D one-color data calibration is not required. For 2D multicolor data SMLocalizer will run but will not compensate for spectral shift between channels without calibration.
For 3D fit calibration image stacks with known z-step size with the same channel setup as for the intended experiment using sub 100 nm beads will need to be acquired. The beads need to emit light in wavelengths covering all filters that will be used in the experiments that will be fitted. So if an experiment calls for the use of the use of Atto-488 and AlexaFluor-647 as STORM dyes a bead emitting photons in the 510-525 and 650-670 nm range is required. Center focus on the beads and include in 10-20 nm steps 1 µm down and 1 µm up from this imaging plane in the image stack, generating a 100-200 frame image with two channels in the dye example above. Ensure that the order of filter imaging is the same as for the experiment, so if Atto-488 will be imaged first, place it first in the calibration experiment as well. 10-20 well-separated beads should be included in the image stack for best results. Also, ensure that you have good signal through the entire stack.
Once the stack has been loaded into ImageJ, select the correct 3D modality (or single slice multichannel image for 2D calibration) and press Calibrate (see fig 4). For 3D modalities, the software will  ask you for the z-step size in nm. Once you press OK, the calibration will commence. Once completed you will be presented with a list of fitted data points in a result table. What is common is that the very edge of the calibration file is a more uncertain then the more central parts, and subsequently any errors will be larger at the edges of the calibration. The calibration algorithms will generate a lookup table for the modality for z fitting and a series of offset values for x-y-z against the first channel for multichannel images. This offset will be applied to all fitted 3D data before being presented to the user.
As small shifts in system alignment will affect the performance, new calibration files should be generated often, preferably weekly for good performance.

BACKGROUND CORRECTION
SMLocalizer performs background correction based on the work of Hoogendoorn el al 3 . In short, static signal (for longer than half the filter width (see fig 5) from Basic input) is removed and transient increases remain (an example before and after image can be seen in fig 6). Background correction can be done through the Process button that will start with background correction or through Correct background which will background correct and replace the current open image stack (see fig 6). SMLocalizer will use a default 101 frame width for background corrections if the checkbox is not checked and the value for the channel altered.   SMLocalizer performs 2D Gaussian and 3D multi-modality fitting on regions extracted from the open image stack through minimizing least square errors of fit. Two settings (Pixel size and Total gain) needs to be looked over for the Gaussian fitting with an additional parameter that can be modified for experienced users, Minimal signal (see fig 7). This will override SMLocalizers standard region identification to only include areas with center pixels stronger then Minimal signal (after background correction).

Modality selection
SMLocalizer can fit both 2D and 3D SMLM data. In the dropdown menu at the bottom half of Basic input (see fig 7), the user can select between 2D and different 3D modalities. In order for fit 3D data, a calibration needs to be generated for that modality. See Calibration (3.2) for details concerning how to generate calibrations. For 3D, select the correct modality from the list: SMLocalizer will proceed to use the calibration for that modality and fit the 3D data, returning x-y-z coordinates for all events that fall within the calibration range. 2D and 3D data are then handled in the same way with regards to parameter filtering, image rendering, channel alignment and drift correction. For details concerning how 3D modalities are processed, see section 5.2. For Biplane data the input data needs to be stitched so that both planes are in a single image next to each other.

Image pixel size
User set parameter that is camera dependent. Change this value to the value that your system uses.

Total gain
This is the total conversion rate from a photon being read by the camera to the resulting pixel intensity. Camera manufacturers provide information on how this number is obtained.

Minimal signal
Minimal intensity of center pixel of a region. Pixels below this value will not be considered a possible particle. By default this is not used but rather particles will be located automatically if the center pixel is frame mean + 0.7σ or stronger.

Figure 9 Relevant buttons and fields for parameter range settings.
Parameter range (see fig 9) is used to select which fitted particles should be included in downstream algorithms. By selecting a parameter type, that type will be used and only particles with values between min and max will be included. Setting parameter range is non-destructive, that is, particles outside of the selected range are not removed but only excluded. An example list of what the parameter distribution can look like after fitting can be found in Figure 8.
The distribution of each parameter can be obtained from the result table. Select the Result menu and Distribution. Select the parameter of interest to plot the selected distribution (see Figure 10).

Clean table
By clicking Clean table SMLocalizer will remove particles that fall outside the selected parameters range. This is destructive.

Photon count
This is the number of photons SMLocalizer has calculated a particle has emitted. For PALM it is reasonable to require 100+ photons per particle whilst for STORM imaging 500+ photons should be required.

Sigma x y
This is the sigma in x and y for the Gaussian fit. Sigma of the Gaussian fit relates to the full width at half maximum (FWHM) as = 2√2 2 (wide field). The width of an in focus emitter will change depending on the optics and filters used but is typically in the range of 100-200 nm. Out of focus emitters will have a broadening of their PSF in 2D which will be reflected by larger values of sigma.

R^2
This is the goodness of fit for a given particle. A value of 1.0 is a perfect fit and 0.0 would indicate no correlation between fit and real data. With decreased R^2 values come an uncertainty in the localization. In our experience values above 0.85 is a good compromise.

Precision x y
This is the precision of the Gaussian fit in x and y. Dependent on sigma x y and photon count. Precision of the fit scales with the square root of the number of photons. This will be highly probe dependent but reasonable values typically fall between 5 -50 nm.

Precision z
This is the precision of the Gaussian fit in z. Precision of the fit scales with the square root of the number of photons for astigmatism and biplane modalities. For PRILM and Double Helix modalities, the error in x-y determination is used to calculate the precision in z. This will be highly probe dependent but reasonable values typically fall between 5 -75 nm. If the precision calculations for PRILM or Double Helix falls outside the calibration range a value of 1000 will be reported.

Frame
Which frames to include. Discarding early portions with to high density of events can improve upon results.

Z
Which z coordinates to include. SMLocalizer will during 3D fitting set z = 0 at the plane of focus with z decreasing towards the beginning of the calibration stack and increasing towards the end of the calibration stack. Any 3D modality will be more accurate in the region surrounding z = 0. SMLocalizer will perform drift correction on the fitted particles during processing if the checkbox is checked (see fig 11). By clicking the Drift correct button (see fig 11), a result table of fitted particles can be drift corrected outside of the main sequence of processing performed by the Process button (see fig 11). Drift correction is performed by maximizing the correlation between two adjacent bins of particles, separated based on frame number. Drift correction is done on each channel separately. Drift correction is prone to artefact introduction if too few particles are included.

DRIFT CORRECTION AND CHANNEL ALIGNMENT
Alignment of channels is performed in the same manner, accounting for drift between start of channel acquisition and chromatic shifts.

Max drift
The maximum drift from one bin to the next in X/Y and Z separately. Smaller allowed drifts will reduce computational cost at the risk of not finding the correct drift.

Number of bins
The number of bins (in time) the fitted particles should be divided into.

Figure 12 Buttons and fields relevant for image rendering.
SMLocalizer will, if the checkbox for Render image is selected (see fig 12), perform this action during Process. Render image uses the particles within the selected parameter range to render an image with a pixel size set by the user. For 3D data an image stack will be created with the voxel size chosen by the user (see fig 12). Default is that each particle adds a value of one to the image, allowing for counting particles in the image. All image editing ImageJ or Fiji is capable of performing can be performed on these result images. Gaussian smoothing multiplies the image values and applies a two pixel wide Gaussian filter. SMLocalizer can perform 2D cluster analysis on the localized particles through the DBSCAN 8 algorithm. A results image of the clusters found will be generated based on the pixel size set in Render image (see 3.8). Particles will be considered to be part of a cluster if they have Min connections neighbors within Epsilon (see fig 13). Figure 14 Multichannel selection.

MULTICHANNEL IMAGES
SMLocalizer can analyze multichannel (but not currently multi frame and slice) images. At the bottom part of the Basic input tab (see fig 14), a list of channels is available. Default is a single channel but more can be added in the list. Each channel has its own settings for all fields and checkboxes that needs to be set.

SUBSEQUENT ANALYSIS IN MATLAB
Once a results table has been obtained from SMLocalizer the data can be transferred to any software that accepts .tif or tab separated tables as input. For Matlab a short function is available online on https://sourceforge.net/projects/smlocalizer/ that will translate the results table into a struct in Matlab containing vectors with all data.
To transfer SMLM data processed in SMLocalizer to Matlab, click File/SaveAs in the results table and store the file, adding .txt at the end of the file name. Once this has been done, download the matlab function LoadSMLocalizer.m from Sourceforge, start Matlab and add LoadSMLocalizer.m to the path (through the Set Path button or by copying the file to the folder containing your other custom Matlab functions). To load the data, use the call:
Start by downloading the dataset and install SMLocalizer into ImageJ/Fiji (this tutorial uses Fiji).
All color maps have been changed to Red hot.
1. Load the dataset into Fiji. 2. Zoom in image to 800 %. 3. Go to Image/Adjust/Brightness / Contrast and press Reset in the B&C panel now displayed.
Next zoom in on the image to a 800% zoom (see fig 18).    7. For this set there are not enough events to do any meaningful drift correction. Pressing Drift correct will tell the user that no correction could be done under these settings. We do not suggest reducing the demands as artifacts are likely to appear. 8. Instead, render the result table using Gaussian smoothing and a R^2 range of 0.85-1.0. To do activate Gaussian smoothing in the Render tab, activate R^2 selection and set its range to 0.85 -1.0 in the parameter selection tab and finally press Render Image (see fig 22).    23). Proceed to select the correct modality in the dropdown menu (Double Helix) and press Calibrate (see fig 24).

3D DOUBLE HELIX
The software will proceed to ask the user for the z-step size for the stack provided (see fig 25, here 10 nm).
The calibration algorithm will now find the optimal parameters for generating the calibration lookup table. Once complete a result calibration curve will be and an empty result table displayed (see fig 26).   The calibration curve will be stored in ImageJ for future use. The next step is to close all open windows and load the experimental file (sequence-as-stack-MT0.N1.LD-DH-Exp). Helix. Once all changes have been made, press Process. Alternatively press Correct background, once complete press Localize and once that is complete finally press Render image.

MITOFILIN PRILM
Provided on https://sourceforge.net/projects/smlocalizer/ are three datasets that are used in this tutorial. A bead stack for PRILM calibration, an overview widefield image and a PRILM modality dataset (see fig 31).
The sample is a U2OS cell, fixed using paraformaldehyde and labeled with a primary antibody against Mitofilin 9 (proteintech, id: 10179-1-AP) and a secondary antibody labeled with Alexa Fluor 647. The images are acquired on a Carl Zeiss Elyra PS.1 using 642 nm activation and increasing back-pumping with 405 nm. A Plan-Apochromate 100x/1.46 Oil (Zeiss) objective was used and emission was collected through a 655 nm long pass filter. Pixel size was 100 x 100 nm, integration time was 25 ms and the gain on the EMCCD camera was set to 100. Load the PRILM calibration stack into ImageJ and start SMLocalizer. Select PRILM as image modality, set Image pixel size to 100 and press Calibrate. When asked, set z-step to 10 nm and press ok. Once completed a plot with the resulting angle-z graph will be displayed along with an empty results table. Close both and proceed to load the Mitofilin PRILM dataset into ImageJ (see fig 32).
With this information, the basic settings in SMLocalizer can be set (Image pixel size: 100, Total gain: 100). Input these settings into SMLocalizer and press Process. Once a results table has been displayed, it is time to change the parameter ranges. Select Photon count and set it to 300-3000 and change the R^2 range to 0.7 -1.0. PRILM distorts the PSF and a great fit against a Gaussian will not work against the more distorted (far from focus) events. There are too few particles to do any reasonable drift corrections with so proceed to change Pixel size to 5 and 5 nm (in Render image section) and proceed to press Render Image (see fig 33a). Do not have the Gaussian smoothing checkbox ticked. Increase the pixel values for the image by the command Process/Math/Multiply and set the value to 5000 and press ok (see fig 33b&c).
For final image rendering, four new images will be generated, all using the just generated 4000x2000x725nm tif stack:  33e). 3. Start with the command Image/Stack/Z project and select Sum Slices and press ok. Use the command Process/Filters/Gaussian Blur and set the radius to 2.0 and press OK. Use the command Image/Lookup tables/mpl-inferno and finish by setting the histogram range (Image/Adjust/Brightness and Contrast) using auto (see fig 33f). 4. Start with the command Process/Filters/Gaussian Blur 3D and set all values to two. Press ok and proceed with finding the highest intensity frame (# 59) and press Auto in the B&C window. Proceed with the command Image/Hyperstacks/Temporal-Color code, select Thermal in the dropdown menu, and press ok. The resulting image will be color coded for depth according to the generated color time scale window. White is in focus, blue shifted colder colors are below the focus and red shifted warmer colors are above the focus. (see fig  33g&h)

2D GOLD BEAD TUTORIAL -DRIFT CORRECTION
Provided on https://sourceforge.net/projects/smlocalizer/ is a bead sample dataset, 50nm_Gold_bead.tif. Download the file for use in this tutorial, load it into ImageJ and start SMLocalizer (see fig 34a&b).

Figure 34 a) SMLocalizer GUI with all settings correct for this tutorial. b) Raw first frame of the input data zoomed in to 2400%. c) Initial, non corrected, output showing the introduced drift in of the bead with mpl-inferno lookup table. d) Drift correction in x and y calculated by SMLocalizer for 5 bins. e) Drift correction in x and y calculated by SMLocalizer for 10 bins. f) Drift correction in x and y calculated by SMLocalizer for 25 bins. g) Drift correction in x and y calculated by SMLocalizer for 50 bins.
As this is a bead sample, we can not background correct this image stack as that would remove the non-blinking bead from our dataset. Instead we directly localize the particles by hitting the Localize key. Wait for a results table to be displayed, change the Photon range to 2000 -4000, tick the checkbox for Gaussian smoothing and press Render image to display the non drift corrected results (see fig 34c). Next we calculate and apply the drift correction for this dataset. Press Drift correct and wait for the Drift corrections plot to appear. Conclude with pressing Render image again. In the figure 34 d-g above we used 4 different values for number of bins, both under and oversampling.

Static events removal
SMLocalizer performs background correction based on the work of Hoogendoorn el al 3 . In short, each pixel is median filtered in time through the following steps: 1. The mean intensity in each frame is calculated. 2. Pixels are normalized to the frame mean.
3. The filter width number of pixels is used to calculate the running median for a specific pixel. 4. The median is subtracted from the original normalized pixel value and is subsequently rescaled using the frame mean. 5. The method removes static (over at least half the filter width) objects and retains transient increases (blinking events).
Pseudo code for performing the background correction (see Hoogendoorn el al 3 for details): Best results are obtained when the filter window width is > 10x the duration of blinking events.

Shot noise
Following median filtering described above each frame is run through a bicubic B-spline filter to further clean out the signal 10 . Each pixel is processed with the following 5x5 kernel: Where b is the background and amp the amplitude of the function. The residual error of fitting is calculated as: Where y is the measured pixel values and ̅ is the mean value.
Optimization is performed by minimizing SSE/SST after initial guesses of input parameters has been done.
 x0 and y0 are estimated as the weighted centroid.  σx and σy are evaluated between 80 nm and 200 nm and the combination yielding the best fit is used for initial estimate.  amp is estimated as the central pixel intensity.  b is estimated to 0.  ϑ is estimated as 0.

3D fitting
All 3D modalities start with 2D fitting using parameter settings obtained from the calibration file. The fitted particles are then translated using a lookup table and raw data to obtain 3D information. For all 3D modalities as part of the calibration, the x-y offset from the frame in focus is determined as function of z. This correction is applied to the final x-y coordinates.

PRILM PRILM fitting starts with generation of a calibration lookup table from an image stack of beads.
Once a calibration lookup table has been generated 3D fitting can be performed. The window used for fitting and the maximum distance allowed between lobes are taken from the calibration file, as is the chromatic corrections.

Z precision estimate for PRILM and double helix
The error in x-y localization propagates to an error in z determination. By shifting the centers for the two lobes by their precision in x and y and calculating a new angle we obtain the uncertainty in angle determination for a given double psf localization. By calculating the resulting z value we get boundaries for the uncertainty of z determination for the given localization. In figure 35a this shift and subsequent angle calculation is demonstrated. In figure 35b the angle is interpreted for a double helix calibration curve. If either the upper or lower border for the precision determination falls outside the calibration curve a value of 1000 will be reported. If both upper and lower bound falls within the calibration curve the mean z offset will be reported as z precision.

DRIFT CORRECTION AND CHANNEL ALIGNMENT
Drift correction is performed by calculating the cross correlation between different bins of particles (in time). For each bin an image (10x10x10 nm) (2D or 3D) is generated. The maximum correlation between two adjacent bins is obtained and the shift applied to the second bins data points.
For each bin the following pseudocode is executed:

IMAGE RENDERING
For single channel and multichannel a single image (with multiple slices for multichannel data) is generated. Particle coordinates are rounded to nearest multiple of pixel size and the resulting pixel coordinate value is increased by 1.
If Gaussian smoothing is selected the image values will be multiplied by 1000 and a Gaussian filter with 2-pixel filter radius will be applied to the image.
All images are calibrated to get accurate dimensions.

CLUSTER ANALYSIS
Cluster analysis will perform cluster analysis on the particles within a given channel that has parameters within a user set ok range. The current version implements DBSCAN for 2D clusters using the apache commons framework. See http://commons.apache.org/proper/commonsmath/apidocs/org/apache/commons/math4/ml/clustering/DBSCANClusterer.html for details on the functions called where eps = Epsilon and minPts = Minimum connections. The results are given as a new rendered image with only the clusters found present and an update to the result list with a column for cluster id. where x, y are 3 dimensional arrays, i is a 3 dimensional index, d is a 3 dimensional shift. and represent array means. Non overlapping voxels due to shift are discarded. output = find shift →max (xCorr)