A tessellation-based colocalization analysis approach for single-molecule localization microscopy.

Validation on simulated data and comparison with other methods

As already discussed, one major limitation of existing colocalization methods dedicated to λSMLM is, in addition to be time-consuming and restricted to 2D data, their sensitivity to the molecular organization, making them difficult to parametrize for routine analysis of different biological samples. This is an important limiting factor since molecular densities of each channel can strongly fluctuate experimentally depending on the labelling strategy, acquisition parameters and biological models of investigation. In order to validate the robustness of our method with respect to the density parameter, we simulated and analyzed 2-color SMLM data of well separated 50 nm radius clusters of circular and square shapes, with varying respective density ratios between 1:1 and 1:5 and varying colocalization ratios between 0 to 100% (Fig. 2a ). We computed the Manders coefficients ( M A and M B ) and the Spearman rank correlation coefficients ( S A and S B ) from the scatterplots of each channel for all the conditions using Coloc-Tesseler. Both quantitative analyses always successfully assessed the level of colocalization of the simulations regardless of the respective molecular densities, and without changing any parameters (Fig. 2b ). All the quantifications remained stable even in the case of strong density ratios between two channels (Supplementary Fig. 4 ). We also validated the efficiency of Coloc-Tesseler to quantify 3D localization data with the same formalism used to analyze 2D data. We more especially investigated the effect of the anisotropic localization precision by scrambling the localizations coordinates with different amplitudes ranging between 0 and 20 nm laterally and between 0 and 60 nm axially. For all the simulated conditions, Coloc–Tesseler systematically performed accurate colocalization analysis, illustrating its capability to analyze realistic 3D λSMLM data with a good robustness to the localization accuracy anisotropy (Supplementary Fig. 5 ). However, these simulations also point-out that Spearman coefficients are more sensitive to localization accuracy, which is expected since degrading the localizations accuracy homogenizes the density distribution of the two channels and increase their correlation. In order to further mimic the heterogeneity existing in biology, we simulated data sets of fully colocalized clusters with varying molecular density, cluster density (i.e. number of clusters per surface unit) and cluster size (Supplementary Fig. 6 ), from which no significant change is expected in the colocalization analysis. This variability is a common situation that can be found in many biological systems, where receptors are organizing and clustering in a highly dynamic manner upon activation to trigger molecular signaling 20 – 22 . We then compared the capability of Coloc-Tesseler (CT), Clus-Doc 10 (CD) and Getis & Franklin 8 (GF) methods to quantify these colocalization data. We first analyzed the impact of changing the density of clusters, ranging between 50 and 200 clusters per fixed field of view (Supplementary Fig. 6a ). Modulating the density of clusters in a fixed background impacts the ratio between clustered and non-clustered molecules, destabilizing the methods quantifying the colocalization by ratiometric computation between the colocalized molecules inside and outside the clusters. As an illustration, GF colocalization showed an important variability and dropped from 56 to 27% when increasing the density of clusters. Combined with DBSCAN to segment the clusters and restrict the colocalization analysis to the clustered molecules, CD allowed stabilizing the quantifications between 69 and 63%. Nevertheless, despite its stability, the computed colocalization value remains low for a 100% simulated colocalization. Since CT uses a similar concept to compute the Manders coefficients, it exhibited strong robustness with respect to the cluster density with a colocalization ranging between 87% and 80.5%. We then investigated the stability of the three methods with respect to the relative molecular densities between channels, ranging from 1:1 to 1:8 (Supplementary Fig. 6b ). This variation in density is common in multicolor SMLM, where molecular densities of clusters in each color may differ due to differences in protein aggregation, antibody specificity, and fluorophore photophysics. Since DBSCAN is known to be sensitive to the molecular density and background 12 , CD method was not stable when varying the respective molecular densities. It resulted in strong differences in the case of high relative density ratio between the channels, with 89% colocalization for channel A and 65% for channel B in the case of 1:8 density ratio. On the contrary, GF was very stable to this parameter thanks to its intrinsic normalization formalism, with a colocalization varying between 60 and 62%. However, GF being more sensitive to background, the colocalization values remained too low for a 100% colocalization condition. CT exhibited a high and stable colocalization ranging between 82 and 96% for all the conditions, thanks to its normalization mechanism and its robustness to noise. Finally, we analyzed the robustness of the 3 methods with respect to the cluster size heterogeneity. We simulated a mixed population of clusters of 100 and 200 nm diameter, with varying percentage of each population between 0 to 100% (Supplementary Fig. 6c ). In this case, both GF and CD exhibited a strong sensibility to this parameter, with colocalization values ranging between 56 and 29% for GF and between 69 and 34% for CD. On the contrary, CT remained highly robust with colocalization values ranging between 89 to 87% for all the conditions. The colocalization analysis performed on these synthetic data illustrates the limits of existing techniques and the versatility and robustness of Coloc–Tesseler to properly quantify the heterogeneity found in single-molecule localization microscopy and cell biology. This heterogeneity is a very important parameter to consider when aiming to decipher the molecular co-organization between different molecular species, under various biological conditions.

Validation on well-known experimental data

We validated our tessellation-based colocalization analysis on 2D and 3D experimental data using well-known biological structures, such as microtubules and nuclear pore complex, as well as actin cytoskeleton regulators. As a control of our method to efficiently analyze 2D and 3D experimental λSMLM, we performed several SMLM acquisitions of samples with varying labeling and molecular densities. For each data set, we quantified, with all the colocalization methods, the entire image as well as different ROI displaying different molecular densities. For the 2D analyses, we only considered the lateral coordinates of the localized molecules. As a positive control, we first acquired and analyzed astigmatism-based 3D single-color DNA-PAINT Tubulin data composed of 7,871,312 localizations, which we randomly split over time to mimic colocalized data with varying relative labelling density ratios ranging between 50%/50 and 10%/90% (Fig. 3a–c ). The normalized density scatter-plots illustrate both the perfect colocalization and correlation between the two channels (Fig. 3d ), as well as the robustness to molecular densities (Supplementary Fig. 7 ). As expected from the simulations, in 2D, CT performed very well for all the density conditions, with a colocalization ranging between 96 to 91% for Manders, and between 76 to 61% for Spearman (Fig. 3e ). CT also exhibited a very strong stability for the different zones on the image, with less than 17% (resp. 26%) fluctuation between ROIs for Manders (resp. Spearman) coefficients. GF also performed quite well, with a colocalization ranging between 83 and 82%, but a higher variability between the different ROI up to 28% (Fig. 3e ). However, CD’s sensitivity to molecular density, made it fastidious to adjust the parameters for the different relative densities. In addition, it could only achieve between 51 and 22% colocalization at the best, with up to 17% fluctuation between ROIs. As observed from the simulations, variabilities between ROIs on the same data illustrates the sensitivity of GF and CD methods, that rely on ratiometric computation between localizations inside the structures and background localizations. The 3D colocalization analysis performed with CT provided very similar results to 2D, with colocalization ranging between 91 to 80% (resp. 83 to 62%) for Manders (resp. Spearman), and an improved variability below 10% (resp. 8%) for all the ROI (Fig. 3f ). Fig. 3 Colocalization analysis on well-characterized biological structures. a – f Colocalization analysis of 3D single-color Tubulin data. ( a ) Random split of the localizations to obtain 2 channels with a density ratio of 50/50 (scale bar 5 µm). Left: 2D projection overlay of the 2 channels. Right: 3D color coding projection of all the localizations. b Magnified views of 3 ROI in a (scale bar 2 µm). c Point rendering for different density ratios within the magnified view in b (left 33 : 67%, middle 20 : 80%, right 10 : 90%, scale bar 250 nm). d Scatterplot of the normalized densities of the 2 colors. e 2D colocalization analyses performed with Coloc-Tesseler (Manders and Spearman), Clus-DoC and Getis & Franklin. f 3D colocalization analysis performed with Coloc-Tesseler (Manders and Spearman). g – k Colocalization analysis of a 3D two-color Tubulin/Lamin dataset (scale bar 5 µm). h Lateral (top) and axial (bottom) views of two magnified ROI of the dataset in f (scale bar 2 µm). i Scatterplot of the normalized densities of the 2 colors. j 2D colocalization analyses performed with Coloc-Tesseler (Manders and Spearman), Clus-DoC and Getis & Franklin. k 3D colocalization analysis performed with Coloc-Tesseler (Manders and Spearman). l – o Colocalization analysis of 2D and 3D two-color acquisitions of nuclear pore complexes (NPC). l 3D NPC dataset of genome-edited expressing Nup107-SNAP cells, labelled with BG-AF647 (magenta) and WGA-CF680 (green) (scale bar 500 nm). m Magnified views of a 2D (left, middle) and 3D (right) NPC (scale bar 50 nm). The same 2D NPC was used for the intensity-based (left) and point (middle) rendering. The point representation (middle) illustrates the result of the 5-class classification obtained using Color-Tesseler. n Scatterplot of the normalized densities of the 2 colors for a 2D (left) and 3D (right) NPC data sets. o 2D (left) and 3D (right) colocalization analyses performed with Coloc-Tesseler (Manders and Spearman), Clus-DoC and Getis & Franklin. In all box plots the center line is the median, the square is the mean and the bounds of the boxes are the 75 and 25% percentiles i.e., the interquartile range (IQR) As a negative control, we performed 3D two-color DNA-PAINT experiments of non-overlapping microtubule (4,086,102 localizations) and lamin (828,756 localizations) structures (Fig. 3g, h ). The normalized density scatterplot of the entire cell nicely illustrates the non-colocalization of the 2 channels. At the entire image level, CT performed very well, with a colocalization of 1% for Manders and −5% for Spearman (Fig. 3i, j ). GF didn’t perform very well in this condition, quantifying up to 20% colocalization and CD performed very well with only 3% colocalization. However, while both GF and CD remained quite stable for all the regions, Manders CT colocalization increased to 20% in the nucleus region, due to 2D projection artefacts (Fig. 3j ). This was corrected using the 3D colocalization analysis (Fig. 3k ), providing less than 0.3% colocalization for all the ROI, illustrating the importance of 3D colocalization analysis in the case of 3D protein distribution. As a more challenging negative control, we performed 2D and 3D dual-color dSTORM experiments of the nucleoporin Nup107 (34,362 localizations in 2D and 46,081 localizations in 3D) and wheat germ agglutinin (WGA) which binds to the disordered and glycosylated regions in the center of the nuclear pore (79,990 localizations in 2D and 35,321 localizations in 3D). Differences in the total number of localizations compared to the microtubule experiments are mainly due to the protein organization, labelling strategies and processing of the data, illustrating the broad range of densities that can be found in experimental SMLM. While these two structures are known for not being colocalized at the nanoscale level, they are very close to each other and organized in a concentric manner (Fig. 3l, m ). The normalized density scatterplots illustrate both the complexity of the data and the difference between 2D and 3D normalized density distributions (Fig. 3n ). While CT could still perform quite well in 2D, with 14% for Manders and 28% for Spearman, and classify efficiently the localizations (Fig. 3l ), GF failed to quantify properly the colocalizations, with more than 80% colocalization (Fig. 3o ). CD performed very well for this dataset with only 1.5% colocalization, even if overall CD always provided much lower values in our hands compared to CT and GF. Interestingly, the 3D quantifications provided by Coloc-Tesseler still performed well, improving the Mander’s colocalization from 14% to 7.5%, but degrading the Spearman rank coefficient from 28 to 38%. This loss of efficiency in 3D versus 2D is certainly due to the decrease of resolution, as illustrated on simulations (Supplementary Fig. 5 ).

Methods 2D simulations

We simulated several single-molecule colocalization data sets organized in clusters of various stoichiometry ratios, with R defined as the enrichment ratio between the densities of molecules inside and outside the clusters. Simulations consist of randomly placed, non-overlapping circular clusters of 100 nm-diameter (channel A) and 100 nm square clusters (channel B). A reference condition was defined with a cluster density of 0.013 mol nm −2 , an enrichment factor R = 10 and an image dimension of 2.5 × 2.5 µm. Channel A was fixed with the reference condition while channel B was defined with cluster densities varying linearly between 0.013 mol nm −2 and 0.065 mol nm −2 with 0.0026 mol nm −2 steps, corresponding to ratios between 1:1 and 1:5 with 20 steps. We simulated three colocalization conditions by varying the inter-distance between clusters of the two-colors of 0 nm, 50 nm and 125 nm. Each condition, corresponding to a given enrichment factor and colocalization, was simulated 10 times, leading to a total of 200 simulations. 3D simulations 3D simulations were analogous to 2D simulations conditions, extended with the axial (z) coordinate. They consisted of 100 nm-diameter spherical clusters of density 1.9 × 10 –4 mol nm −3 randomly distributed in a 2.5 × 2.5 × 1 µm volume with an enrichment factor R = 146 . In order to test the efficiency of our method with respect to the localization accuracy, we degraded our simulation data by scrambling the molecule positions with 5 different localization accuracy couples of (Δ xy = 0 nm, Δ z = 0 nm), (Δ xy = 20 nm, Δ z = 20 nm), (Δ xy = 20 nm, Δ z = 40 nm) and (Δ xy = 20 nm, Δ z = 60 nm). We then simulated the three colocalization conditions defined in the 2D case. Each condition, corresponding to a given scrambling mode and colocalization, was simulated 10 times, leading to a total of 120 simulations. Voronoï-based colocalization analysis The Voronoï diagram V is a space-partitioning technique subdividing a space documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$S in {Bbb R}^n$$end{document} S ∈ R n containing an ensemble of localizations documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$L = left{ {s_i,2 le i le n} right} in S$$end{document} L = s i , 2 ≤ i ≤ n ∈ S into polytopes (i.e. polygons in 2D and polyhedrons in 3D). Each localization s i is described by a unique polytope P i ∈ S centered on s i , with s i ∈ P i , from which several parameters can be computed such as its area A i (or volume V t in 3D), n i the number of direct neighbors (i.e. polytopes sharing a common edge with P i ) of s i , d ( s i , s i,j ) the Euclidian distance from s i to one of its direct neighbor s i,j or its 1 st rank local localization density δ i (Supplementary Fig. 1a ), other parameters and their definitions can be found in Levet et al. 12 ). The normalized 1 st rank density is defined as documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$widehat {delta _i} = delta _i/delta$$end{document} δ i ^ = δ i ∕ δ , with δ the average density of a spatially random reference distribution. In the case of the colocalization analysis between two channels, A and B , there are two ensembles of localizations documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$L^{rm{A}} = left{ {s_i^{rm{A}},2 le i le n^{rm{A}}} right} in S$$end{document} L A = s i A , 2 ≤ i ≤ n A ∈ S and documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$L^{rm{B}} = left{ {s_i^{rm{B}},2 le i le n^{rm{B}}} right} in S$$end{document} L B = s i B , 2 ≤ i ≤ n B ∈ S with their corresponding Voronoï diagrams V A and V B (Fig. 1a ). Each Voronoï diagram can be analyzed independently using a unique threshold T ≥0, allowing classifying the localizations in 3 orthogonal classes: two high-density classes C A and C B , and one background class documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{bar {rm{A}}bar {rm{B}}}$$end{document} C A ¯ B ¯ , defined by: 1 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{rm{A}} = left{ {s_i^{rm{A}} in L^{rm{A}}{mathrm{|}}widehat {delta _i}^{rm{A}} ge T} right}$$end{document} C A = s i A ∈ L A ∣ δ i ^ A ≥ T 2 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{rm{B}} = left{ {s_i^{rm{B}} in L^{rm{B}}{mathrm{|}}widehat {delta _i}^{rm{B}} ge T} right}$$end{document} C B = s i B ∈ L B ∣ δ i ^ B ≥ T 3 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{bar {rm{A}}bar {rm{B}}} = left{ {s_i^{rm{A}} in L^{rm{A}}{mathrm{|}}widehat {delta _i}^{rm{A}} < T} right} cup left{ {s_i^{rm{B}} in L^{rm{B}}{mathrm{|}}widehat {delta _i}^{rm{B}} < T} right}$$end{document} C A ¯ B ¯ = s i A ∈ L A ∣ δ i ^ A < T ∪ s i B ∈ L B ∣ δ i ^ B < T with documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{rm{A}} + C^{rm{B}} + C^{bar {rm{A}}bar {rm{B}}} = L^{rm{A}} + L^{rm{B}}$$end{document} C A + C B + C A ¯ B ¯ = L A + L B (Fig. 1b ). The surfaces (or volumes in 3D) occupied by each high-density class O A and O B , can be computed by segmentation from the neighboring polygons of C A and C B as described in 12 (Supplementary Fig. 2b ). A simplified version of the Manders’ fractional overlapping coefficients, only accounting for overlapping surfaces can then be defined as: 4 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$M^{rm{A}} = frac{{O^{rm{A}} cap O^{rm{B}}}}{{O^{rm{A}}}}$$end{document} M A = O A ∩ O B O A 5 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$M^{rm{B}} = frac{{O^{rm{A}} cap O^{rm{B}}}}{{O^{rm{B}}}}.$$end{document} M B = O A ∩ O B O B . This formulation requires segmenting both channels in order to extract the surface (or volume in 3D) of each channel. However, the subdividing space Voronoï tessellation architecture makes it possible to avoid this step. Indeed, since V A and V B are defined on the same spatial domain S , it is possible to pair the molecules in both channels using overlapping polytopes. Each localization documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_i^{{{rm{A}}}}$$end{document} s i A (resp. documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_i^{rm{B}}$$end{document} s i B ) can therefore be associated with its corresponding molecule documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_j^{rm{B}}$$end{document} s j B (resp. documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_j^{rm{A}}$$end{document} s j A ) in the other channel, with documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_j^{rm{B}}$$end{document} s j B (resp. documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_j^{rm{A}}$$end{document} s j A ) determined such as documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_i^{rm{A}} in P_j^{rm{B}}$$end{document} s i A ∈ P j B (resp. documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_i^{rm{B}} in P_j^{rm{A}}$$end{document} s i B ∈ P j A ). We introduce a pair-density descriptor ( documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$widehat {delta _i}^{rm{A}},widehat {delta _j}^{rm{B}}$$end{document} δ i ^ A , δ j ^ B ) (resp. ( documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$widehat {delta _j}^{rm{A}},widehat {delta _i}^{rm{B}}$$end{document} δ j ^ A , δ i ^ B )) to quantify the spatial co-organization of the localization documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_i^{rm{A}}$$end{document} s i A (resp. documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_i^{rm{B}}$$end{document} s i B ). The high-density classes C A and C B can then be divided into 4 classes defined by: 6 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{{rm{AB}}} = left{ {s_i^{rm{A}} in C^{rm{A}}{mathrm{|}}widehat {delta _j}^{rm{B}} ge T} right}$$end{document} C AB = s i A ∈ C A ∣ δ j ^ B ≥ T 7 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{{rm{A}}bar {rm{B}}} = left{ {s_i^{rm{A}} in C^{rm{A}}{mathrm{|}}widehat {delta _j}^{rm{B}} < T} right}$$end{document} C A B ¯ = s i A ∈ C A ∣ δ j ^ B < T 8 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{{rm{BA}}} = left{ {s_i^{rm{B}} in C^{rm{B}}{mathrm{|}}widehat {delta _j}^{rm{A}} ge T} right}$$end{document} C BA = s i B ∈ C B ∣ δ j ^ A ≥ T 9 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{{rm{B}}bar A} = left{ {s_i^{rm{B}} in C^{rm{B}}{mathrm{|}}widehat {delta _j}^{rm{A}} < T} right}$$end{document} C B Ā = s i B ∈ C B ∣ δ j ^ A < T with documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{{rm{AB}}} + C^{{rm{A}}bar {rm{B}}} = C^{rm{A}}$$end{document} C AB + C A B ¯ = C A and documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{{rm{BA}}} + C^{{rm{B}}bar {rm{A}}} = C^{rm{B}}$$end{document} C BA + C B A ¯ = C B . These four classes describe whether a high-density localization in channel A (resp. B) lies ( C AB , resp. C BA ) or not ( documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{{rm{A}}bar {rm{B}}}$$end{document} C A B ¯ , resp. documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{{rm{B}}bar {rm{A}}}$$end{document} C B A ¯ ) into a high-density polytope on the other channel (Fig. 1c ).

Edge effects correction

The localizations at the frontier between high- and low-density classes can lead to artifacts, an inherent limitation of Voronoï diagrams. Indeed, the region of influence of these border localizations, represented by their polytopes, is usually larger compared to the ones of the inside localizations, influencing the colocalization classification in a negative manner (Supplementary Fig. 3a-c ). In particular, it results in classifying localizations in C AB and C BA , even when high density regions of the 2 colors are only partially colocalized (Supplementary Fig. 3c ). The Delaunay triangulation, dual of the Voronoï diagram, is not prone to edge effect and can be used to correct these localizations by transferring them from C AB (resp. C BA ) to documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{{rm{A}}bar {rm{B}}}$$end{document} C A B ¯ (resp. documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C^{{rm{B}}bar {rm{A}}}$$end{document} C B A ¯ ). First, we determine Tr A and Tr B , the triangle sets that describe C A and C B . A triangle is part of Tr A (resp. Tr B ) if its three vertices belong to C A (resp. C B ) (Supplementary Fig. 3d ). To ensure a better stability of the correction, we remove outliers from Tr A and Tr B by using the interquartile range (IQR) method defined as: IQR = q75 - q25 (10) With q25 and q75 being the 25th and 75th percentiles of the triangle area distribution. Then outliers are defined as triangles with areas bigger than q75 + (IQR*1.5) (Supplementary Fig. 3e ). All outliers are then removed from Tr A and Tr B (Supplementary Fig. 3f ) and the final classes are defined as: 11 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${rm{corr}}left( {C^{{rm{AB}}}} right) = left{ {s_i^{rm{A}} in C^{{rm{AB}}}{mathrm{|}}s_i^{rm{A}} in Tr^{rm{B}}} right},$$end{document} corr C AB = s i A ∈ C AB ∣ s i A ∈ T r B , 12 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${rm{corr}}left( {C^{{rm{A}}bar {rm{B}}}} right) = C^{{rm{A}}bar {rm{B}}} + left{ {s_i^{rm{A}} in C^{{rm{AB}}}{mathrm{|}}s_i^{rm{A}} notin Tr^{rm{B}}} right}$$end{document} corr C A B ¯ = C A B ¯ + s i A ∈ C AB ∣ s i A ∉ T r B 13 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${rm{corr}}left( {C^{{rm{BA}}}} right) = left{ {s_i^B in C^{{rm{BA}}}{mathrm{|}}s_i^{rm{B}} in Tr^{rm{A}}} right}$$end{document} corr C BA = s i B ∈ C BA ∣ s i B ∈ T r A 14 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${rm{corr}}left( {C^{{rm{B}}bar {rm{A}}}} right) = C^{{rm{B}}bar {rm{A}}} + left{ {s_i^{rm{B}} in C^{{rm{BA}}}{mathrm{|}}s_i^{rm{B}} notin Tr^{rm{A}}} right},$$end{document} corr C B A ¯ = C B A ¯ + s i B ∈ C BA ∣ s i B ∉ T r A , where documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_i^{rm{A}} in Tr^{rm{B}}$$end{document} s i A ∈ T r B means that the localization documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_i^{rm{A}}$$end{document} s i A falls inside one triangle of Tr B (Supplementary Fig. 3g ). Scatterplot representation A scatterplot is a 2D histogram, commonly used in image-based colocalization analysis where the intensity of one channel is plotted against the intensity of the other channel for each pixel. It allows to visually investigate the degree of colocalization between two channels. We adapted the scatterplot representation to λSMLM data using the Tessellation-based architecture provided by the Voronoï diagrams, plotting each localization documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_i^{{{rm{A}}}}$$end{document} s i A (resp. documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s_i^{rm{B}}$$end{document} s i B ) with respect to its pair-density coordinates ( documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$widehat {delta _i}^{rm{A}},widehat {delta _j}^{rm{B}}$$end{document} δ i ^ A , δ j ^ B ) (resp. ( documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$widehat {delta _j}^{rm{A}},widehat {delta _i}^{rm{B}}$$end{document} δ j ^ A , δ i ^ B )). We defined one scatterplot for each channel, always keeping the density of channel A in abscises and the density of channel B in ordinates. The scatterplot of channel A (resp. B) is defined by n A (resp. n B ) points of coordinates ( documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$x_i^{rm{A}} = widehat {delta _i}^{rm{A}}$$end{document} x i A = δ i ^ A , documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$y_i^{rm{A}} = widehat {delta _j}^{rm{B}}$$end{document} y i A = δ j ^ B ) (resp. ( documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$x_i^{rm{B}} = widehat {delta _j}^{rm{A}},y_i^{rm{B}} = widehat {delta _i}^{rm{B}}$$end{document} x i B = δ j ^ A , y i B = δ i ^ B )). Voronoï Manders’ overlapping coefficients They can be derived from the thresholded scatterplot representations to quantify the spatial overlapping between the two channels. They are defined by: 15 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$forall ;s_i^{rm{A}} in C^{rm{A}},M^{rm{A}} = frac{{mathop {sum}nolimits_i {alpha _ix_i^{rm{A}}} }}{{mathop {sum}nolimits_i {x_i^{rm{A}}} }};{mathrm{with}}left{ {begin{array}{*{20}{c}} {alpha _i = 1,;if;s_i^{rm{A}} in {rm{corr}}left( {C^{{rm{AB}}}} right)} \ {alpha _i = 0,;if;s_i^{rm{A}} in {rm{corr}}left( {C^{{rm{A}}bar {rm{B}}}} right)} end{array}} right.$$end{document} ∀ s i A ∈ C A , M A = ∑ i α i x i A ∑ i x i A with α i = 1 , i f s i A ∈ corr C AB α i = 0 , i f s i A ∈ corr C A B ¯ 16 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$forall ;s_i^{rm{B}} in C^{rm{B}},M^{rm{B}} = frac{{mathop {sum}nolimits_i {alpha _iy_i^{rm{B}}} }}{{mathop {sum}nolimits_i {y_i^{rm{B}}} }};{mathrm{with}}left{ begin{array}{l}alpha _i = 1,;if;s_i^{rm{B}} in corrleft( {C^{{rm{BA}}}} right)\ alpha _i = 0,;if;s_i^{rm{B}} in {rm{corr}}left( {C^{{rm{B}}bar {rm{A}}}} right)end{array} right.$$end{document} ∀ s i B ∈ C B , M B = ∑ i α i y i B ∑ i y i B with α i = 1 , i f s i B ∈ c o r r C BA α i = 0 , i f s i B ∈ corr C B A ¯ Spearman’s rank correlation coefficients The normalized densities documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$widehat {delta _i}^{rm{A}}$$end{document} δ i ^ A (resp. documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$widehat {delta _i}^{rm{B}}$$end{document} δ i ^ B ) sorted in ascending order can be substituted by their rank documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$r(widehat {delta _i}^{rm{A}})$$end{document} r ( δ i ^ A ) (resp. documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$r(widehat {delta _i}^{rm{B}})$$end{document} r ( δ i ^ B ) ) in the ordered density distribution. The Spearman’s rank correlation coefficients quantify the similarities between documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$r(widehat {delta _i}^{rm{A}})$$end{document} r ( δ i ^ A ) and documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$r(widehat {delta _i}^{rm{B}})$$end{document} r ( δ i ^ B ) . Compared to the popular Pearson coefficients, they don’t impose any linear relationship between the densities of the two channels, making it more suitable for λSMLM data colocalization analysis. They can be derived from the scatterplot representations to quantify the spatial co-organization between the two channels. They are defined by: 17 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$forall s_i^{rm{A}} in L^{rm{A}},;S^{rm{A}} = 1 - frac{{6mathop {sum}nolimits_i {left( {rleft( {x_i^{rm{A}}} right) - rleft( {y_i^{rm{A}}} right)} right)^2} }}{{n^{rm{A}}left( {n^{{rm{A}}2} - 1} right)}}$$end{document} ∀ s i A ∈ L A , S A = 1 - 6 ∑ i r x i A - r y i A 2 n A n A 2 - 1 18 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$forall s_i^{rm{B}} in L^{rm{B}},;S^{rm{B}} = 1 - frac{{6mathop {sum}nolimits_i {left( {rleft( {x_i^{rm{B}}} right) - rleft( {y_i^{rm{B}}} right)} right)^2} }}{{n^{rm{B}}left( {n^{{rm{B}}2} - 1} right)}}$$end{document} ∀ s i B ∈ L B , S B = 1 - 6 ∑ i r x i B - r y i B 2 n B n B 2 - 1 with r ( x i ) and r ( y i ) the ranks of the scatterplot coordinates of the i th localization. We used the implementation provided by the Alglib library ( http://www.alglib.net/ ). DNA-PAINT Tubulin experiment MEF (Mouse Embryonic Fibroblast) cells were allowed to spread onto 1.5 H clean coverslips for 4 h prior fixation with 4% paraformaldehyde (Sigma) + 0.2% Glutaraldehyde (Sigma) + 0.3% Triton X-100 (Sigma) for 10 min, followed by the quenching of autofluorescence using 150 mM Glycin (Sigma) in PBS for 10 min and an additional permeabilization step with 0.3 % Triton X-100 in PBS for 10 min. Unspecific sites were then blocked using 5% BSA in PBS for at least 2 h prior incubation of primary antibodies for 2 h at room temperature. Anti-tubulin (rat, alpha-Tubulin #MA1–80017, Thermo Fisher) and anti-lamin (Goat, #sc-6217, Santa Cruz) were diluted at 1:500 in the blocking solution. After three washes, secondary antibodies were incubated 2 h at room temperature. Anti-rat-P1 and anti-Goat-P5 antibodies functionalized with two orthogonal oligonucleotides sequences, P1 and P5, were kindly given by Ralf Jungmann’s lab and were each used at 1:100 in blocking solution. Finally, coverslips were washed three times in blocking solution and then in PBS and were kept at 4 °C before use. Prior DNA-PAINT experiments, coverslips were incubated in a solution containing 100 nm fluorescent nanodiamonds (NDNV100 nmMd10 ml from Adamas Nanotechnologies, Inc) for 30 min and then washed three times in PBS. DNA-PAINT acquisitions were performed on a TiE Nikon microscope equipped with a TIRF illumination module (iLAS2, Roper Scientific) fiber coupled to a 635 nm and 561 nm lasers (Errol), a 100 × 1.49NA objective lens (CFI SR HP Apochromat TIRF 100XC Oil, Nikon), a quad-band filter set (ZET 405/488/561/640, Chroma), a N-STORM astigmatism lens (Nikon, France), and an EMCD camera (Evolve512, Photometrics). P1-Cy3b and P5-Cy3b imagers (kindly provided by Ralf Jungmann’s Lab) were used at a concentration of 0.2 nM and 0.5 nM respectively in imaging media (PBS supplemented with 500 mM of NaCl (Sigma)). For two-color experiments, sequential exchange PAINT approach was used: after acquisition using a first imager, coverslips were carefully washed several times with imaging media to remove all imagers before to add the second orthogonal imagers and perform a second acquisition sets. It enables to image both structures with the same Cy3b fluorophores avoiding possible chromatic aberrations. Acquisition sequences of 40,000 frames per channel were steered using MetaMorph software at 5 Hz in streaming mode. 3D single-molecule localization and super-resolution image reconstruction were achieved using the WaveTracer module (Molecular Devices) which uses a combination of wavelet-based localization and anisotropic Gaussian fitting methods 29 , 30 . dSTORM nuclear pore experiment Nuclear pores were stained and imaged as described previously 31 . Genome-edited U-2 OS cells that expressed Nup107–SNAP were cultured under adherent conditions in DMEM (high-glucose, without phenol red) supplemented with 10% (v/v) FBS, 2 mM l -glutamine, nonessential amino acids, and ZellShield at 37 °C, 5% CO 2 and 100% humidity. All incubations were carried out at room temperature. For nuclear pore staining, the coverslips were prefixed with 2.4% (w/v) formaldehyde (FA) in PBS for 30 s. Cells were permeabilized with 0.4% (v/v) Triton X-100 in PBS for 3 min and then fixed with 2.4% (w/v) FA in PBS for 30 min. Subsequently, the fixation reaction was quenched by incubation in 100 mM NH 4 Cl in PBS for 5 min. After being washed twice with PBS, the samples were blocked with Image-iT FX signal enhancer (Thermo Fisher Scientific, Waltham, MA, USA) for 30 min. The coverslips were incubated in staining solution (1 μM benzylguanine Alexa Fluor 647 (S9136S; NEB, Ipswich, MA, USA), 1 mM DTT, 1% (w/v) BSA in PBS) for 50 min in the dark. After being rinsed three times with PBS and washed three times with PBS for 5 min, the samples were stained with wheat germ agglutinin coupled to CF680 (29029, Biotium, Fremont, CA, USA). The coverslips were incubated in the staining solution (0.2 µg.mL −1 WGA-CF680 in 1% (w/v) BSA in PBS) for 5 min in the dark. After washing three times with PBS for 5 min, the sample was mounted for imaging. The sample was imaged in blinking buffer (50 mM Tris, pH 8, 10 mM NaCl, 10% (w/v) D-glucose, 35 mM 2-mercaptoethylamine, 500 μg.mL −1 GLOX, 40 μg.mL −1 catalase). SMLM image acquisition was performed at room temperature (24 °C) on a customized microscope equipped with a high-numerical-aperture (NA) oil-immersion objective (×160, 1.43-NA; Leica, Wetzlar, Germany) with homogenous multi-mode fiber illumination 32 . A closed-loop focus lock system was implemented, using the signal of a near-infrared laser reflected by the coverslip and its detection by a quadrant photodiode. The fluorescence emission was split by a 665 nm long-pass filter (AHF, Tübingen, Germany). Both color channels were imaged side-by-side on an EMCCD camera (Evolve512D; Photometrics, Tucson, AZ, USA) after filtering by a 685/70 and 676/37 filters (AHF), respectively. A cylindrical lens (f = 1000 mm, Thorlabs, Newton, NJ, US) introduced astigmatism for 3D localization. The pulse length of the 405 nm laser was automatically adjusted to retain a constant number of localizations per frame. The 3D positions of the fluorophores were determined with a MLE fit using an experimentally derived PSF model 31 . The color was assigned based on the relative intensity of the fluorophores in both spectral channels. A 3D redundant cross-correlation based drift correction was employed 31 and localizations persistent in consecutive frames were grouped into one localization. PALM-dSTORM experiments for cytoskeleton regulators We used dissociated rat hippocampal neurons transfected using Effectene (Qiagen) at 7 days in vitro (DIV) with F-actin regulators fused to mEos2 (mEos2::Abi1, mEos2::ArpC5A, mEos2::VASP). Neurons were co-transfected with constitutively active Rac1-Q61L and mEos2::Abi1 as a negative control for our single-molecule-based colocalization analysis. Data were acquired as in Chazeau et al. 22 by sequential dual-color SMLM using PALM for F-actin regulators fused to mEos2 and dSTORM for endogenous PSD95 immunostained with a mouse primary anti-PSD95 antibody revealed with an Alexa647-coupled anti-mouse secondary antibody.

Chromatic aberration correction

Localization errors induced from field-dependent chromatic aberrations were characterized and corrected using a bi-dimensional 3 rd order polynomial field correction. A single image containing ≥ 10 fiducial markers covering the entire field of view was acquired for each excitation wavelength used, and the localizations of the individual fiducial markers in each channel were paired and used to calculate the field of view transformation. After the acquisition, localizations from the second emission channel (typically 561 nm excitation, Channel B) were transformed into the space of the reference, far-red emission channel (640 nm excitation, Channel A), with an error smaller than the localization precision of the individual localizations. Chromatic and drift corrections We used 100 nm multicolor fluorescent microbeads (Tetraspeck, Invitrogen) or 100 nm fluorescent nanodiamonds (NDNV100nmMd10ml, Adamas Nanotechnologies) as fiducial markers to register multicolor experimental data and correct for lateral drifts. After drift correction, chromatic shift was automatically corrected using a two-stage process. First, the bead positions of the 2 channels were automatically computed as the barycenter of all the registration beads’ localizations. Second, all the localizations of the second channel (Channel B) were translated by the displacement vector computed between the 2 bead positions.

Bead filtering

Fluorescent beads were automatically excluded from the analysis after chromatic and drift correction. For each localization, we computed the number of neighboring localizations within a radius of 100 nm using kd-tree implementation of Jose Luis Blanco-Claraco ( https://github.com/jlblancoc/nanoflann ). Then, localizations with a number of neighbors greater or equal to 90% of the total number of frames were removed. Determination of the distance between clusters of F-actin regulatory proteins and PSD95 Clusters of F-actin regulator proteins were identified by a two-level segmentation process using SR-Tesseler 12 . First, we segmented the neuron contours by thresholding the localizations with documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$delta _i^1$$end{document} δ i 1 > 2 δ Ι , δ Ι being the average localization density of the whole image. Then we computed potential embedded clusters into the selected synapses using a threshold of documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$delta _i^1$$end{document} δ i 1 > 2 δ Ν , where δ Ν is the average localization density inside the neuron contour. PSD95 clusters were segmented inside the selected synapses using a single-level threshold of documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$delta _i^1$$end{document} δ i 1 > 2 δ Ν . For both F-actin regulators and PSD95, the cluster barycenter was defined as the centroid of all the localizations composing each cluster. The cluster distance between the 2 channels was computed as the shortest distance between the clusters’ barycenter on identified synapses having F-actin regulators and PSD95 clusters.

Implementation and benchmarking

Coloc-Tesseler software uses a combination of a multi-view OpenGL-based visualization with a C + + code optimization for the colocalization analysis. It is fast both for the colocalization analysis and the visual rendering and includes an efficient batch analysis mode. As an illustration, it took only 7:17 min, including 1:50 min to load the data, to analyze and display the complete experimental data sets, composed of 18 two-colors cells labelled with F-actin regulators fused to mEos2 and PSD95 immunostained with A647, corresponding to 8,276,861 localizations. The 3D colocalization analysis of the 3D DNA-PAINT Tubulin biggest dataset composed of 7,871,312 localizations took 19:01 min, including 1:52 min to load the data. Computation were done using a standard computer equipped with an Intel Xeon 2.40 GHz processor. Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Supplementary information Supplementary Information Peer Review File Reporting Summary