Adaptive FDR thresholding of Fourier shell correlation for resolution estimation of cryo-EM maps

Fourier shell correlation (FSC) has become a standard quantity for resolution estimation in electron cryo-microscopy. However, the resolution determination step is still subjective and not fully automated as it involves a series of map interventions before FSC computation and includes the selection of a common threshold. Here, we apply the statistical methods of permutation sampling and false discovery rate (FDR) control to the resolution-dependent correlation measure. The approach allows fully automated and mask-free resolution determination based on adaptive thresholding of FSC curves. We demonstrate the applicability for global, local and directional resolution estimation and show that the developed criterion termed FDR-FSC gives realistic resolution estimates based on a statistical significance criterion while eliminating the need of any map manipulations. The algorithms are implemented in a user-friendly GUI based software tool termed SPoC (https://github.com/MaximilianBeckers/SPOC).


Abstract
Fourier shell correlation (FSC) has become a standard quantity for resolution estimation in electron cryo-microscopy. However, the resolution determination step is still subjective and not fully automated as it involves a series of map interventions before FSC computation and includes the selection of a common threshold. Here, we apply the statistical methods of permutation sampling and false discovery rate (FDR) control to the resolution-dependent correlation measure. The approach allows fully automated and mask-free resolution determination based on adaptive thresholding of FSC curves. We demonstrate the applicability for global, local and directional resolution estimation and show that the developed criterion termed FDR-FSC gives realistic resolution estimates based on a statistical significance criterion while eliminating the need of any map manipulations. The algorithms are implemented in a user-friendly GUI based software tool termed SPoC (https://github.com/MaximilianBeckers/SPOC).

Introduction
Electron cryo-microscopy (cryo-EM) is becoming established as one of the methods of choice for macromolecular structure determination. The technique has undergone major improvements in hardware and software, allowing to determine 3D structures routinely at close-to-atomic resolution (Kühlbrandt et al., 2014;Bartesaghi et al., 2015;Weis et al., 2019). These resolutions provide the landmark features for atomic model building and the resulting atomic coordinates rationalize biological mechanism and function. At the core of all these developments are improvements in the resolvability of molecular detail that can be achieved. Therefore, the resolution is a reported value that makes the data comparable with other structural biology methods, and more importantly, gives guidance how confidently we can interpret the density of a given feature in 3D space. Fourier shell correlation (FSC) curves are the standard metric for resolution estimation of cryo-EM maps. FSCs were originally introduced to the field of cryo-electron microscopy (Harauz & Van Heel, 1986) and more recently were also applied to related fields such as super-resolution microscopy (Nieuwenhuizen et al., 2013;Banterle et al., 2013).
The FSC measures the correlation between Fourier coefficients within every resolution shell of two independently determined half-maps. A typical curve shows high correlations at low resolutions until it drops to zero at higher resolutions when noise starts to dominate the signal. In order to report a resolution value of structures based on FSC curves, a threshold value needs to be selected. A fixed threshold of 0.143 has been proposed (Rosenthal & Henderson, 2003), which is widely used for resolutions better than 10 Å when map features can be used to validate the obtained resolution.
For lower resolutions, a more conservative 0.5 threshold has typically been used (Spahn et al., 2004). The 0.5 threshold is also favored when local FSCs within small windows are computed and local resolution is estimated (Cardone et al., 2013).
Nevertheless, fixed value thresholds ignore the effective number of Fourier coefficients in the respective resolution shell ( Van Heel & Schatz, 2005). This effect becomes particularly relevant at local resolution estimation with small window sizes (Cardone et al., 2013) and for 3D FSC with small number of Fourier pixels (Zi Tan et al., 2017). Alternatively, other criteria like the 2 and 3 (Saxton & Baumeister, 1982;Orlova et al., 1997) as well as the half-bit criterion have been proposed that compensate for these effects ( Van Heel & Schatz, 2005). More generally, applying thresholds such as the 0.143/ 0.5 or half-bit criterion are to be interpreted as identifying the highest resolution shell that still contains a minimum of interpretable signal.
In order to get correct estimates of resolution-dependent information measures of the molecular density within a 3D reconstruction, solvent noise needs to be removed from the volume, as big parts of the reconstructed volume correspond to noise only and thus bias the resolution towards lower values. Therefore, mask application is critical, yet poses the danger of introducing artificial correlations. This effect has been proposed to be corrected by mask deconvolution (Chen et al., 2013). In practice, this approach still involves testing of several masks until the user reaches an agreement between the displayed molecular features and the estimated resolution. Although the calculation of corrected FSC curves considering the molecular mass have been proposed to circumvent this problem (Sindelar & Grigorieff, 2012), such an approach is also used in combination with a defined threshold criterion.
Alternatively, -thresholds were proposed in order to provide cutoffs when the FSC exceeds the random correlations of pure noise (Saxton & Baumeister, 1982;Van Heel & Schatz, 2005). However, statistics based on simple -thresholds suffer from several drawbacks as linking -levels to probabilities requires strict assumptions about the underlying probability distributions, which are usually not known. In order to minimize subjectivity and increase robustness of FSC computations, we developed a new approach of adaptive thresholding for resolution determination using a combination of permutation sampling and p-value correction by false discovery rate control. We show that the procedure is able to detect signal with high precision over a wide range of noise levels that make masking of half-maps obsolete. We validate the approach by comparing a large batch of structures with reported resolution values of the EM databank and further extend the approach successfully to challenging cases of local and directional resolution determination.

Resolution estimation using statistical thresholding of the FSC
In order to circumvent principal and practical problems of thresholding FSC curves, we developed a procedure for identifying the highest resolution shell of interpretable signal based on parameter-free permutation sampling and subsequent statistical inference. Permutation sampling tests are statistical procedures for estimating the nullhypothesis from the data itself and thus do not require any assumptions about the underlying distributions (Lehmann & Romano, 2005). Permutation testing of FSC coefficients for each resolution shell is straightforward: under the null-hypothesis the two resolution shells from the half-maps are independent, we can thus generate new samples by permutation, i.e. scrambling the Fourier coefficients of the second resolution shell (Figure 1a). As a result, we compute a large series of noise FSCs.
Effectively, this permutation allows sampling of the null distribution and as a result FSC-values for each resolution shell can be conveniently transformed into p-values.
In order to account for this multiple testing problem, p-values are then corrected by means of false discovery rate (FDR) (Benjamini & Hochberg, 1995) control and thresholded at 1 %, as previously applied to cryo-EM map thresholding . We term this approach the FDR-FSC method of resolution determination.
In order to test whether the permutation approach captures the principal distribution of FSC coefficients, we simulated two pure noise half-maps with an average of 0 and a standard deviation of 1 and applied the permutation sampling as described. To assess the true distribution of FSC coefficients, we simulated 5,000 noise half-maps (again with an average of 0 and a standard deviation of 1) and calculated FSC curves for each pair of 5,000 half-maps. Comparison of the simulated with the permuted histograms generated from the FSC values of the half-Nyquist resolution shell (1/4 pixel) and the corresponding empirical cumulative distribution functions (ECDF) show that the distributions from the permutation approach and the simulation follow each other closely (Figure 1b), in particular at the tail of the distributions (Figure 1c).
However, the comparison within the half-Nyquist shell shows that simple standard deviations computed from non-permuted FSC distributions as they are used by the proposed 3 (Orlova et al., 1997) and modified 3 criteria ( Van Heel & Schatz, 2005) for resolution estimation can give rise to overestimates in comparison with permuted distributions. Next, we wanted to assess the performance over all resolution shells and plotted the right-sided 10, 5.0, 1.0 and 0.15 % cutoff values (percentiles) for each resolution shell. The percentiles of the distributions from the permutation and the simulation follow each other closely (Supplementary Figure 1a). Together, the permutation approach within each resolution shell is able to accurately represent the FSC distribution under the null hypothesis and the standard deviations computed from the permuted distributions are more robust than simple standard deviations as they are used for traditional 3 criteria.
When symmetry is applied during the image reconstruction, additional dependencies between half-maps are introduced. These effects on the FSC computation and cutoff determination have been discussed previously and can be compensated by the reduction of the sample size ( Van Heel & Schatz, 2005). Therefore, in the permutation approach we took this effect into account and reduced the sample size from each resolution shell by the factor of available symmetry-related asymmetric volume units, i.e. D4 symmetry has 8 asymmetric volume units and reduces the effective sampling to 1/8 or 12.5 % of the initial sample size. In order verify the approach, we tested permutation-based sampling of two symmetry-imposed half-maps in the presence and absence of the sample size correction factor. In analogy to the test on the permutation approach above, we compared the permuted distributions with the ones obtained 5,000 noise symmetrized half-maps by plotting the 10, 5.0, 1.0 and 0.15 percentiles as a function of spatial frequency. In the absence of any sampling correction, the permuted and simulated percentiles do not overlap whereas as in the presence of a D4 and D7 symmetry correction factor, the permuted and simulated percentile lines approach each other more closely (Supplementary Figure 1b and c). Similar effects of introducing additional dependencies between half-maps occur during masking, which most image reconstruction algorithms make use of already by reconstructing density within a sphere. This has been discussed in ( Van Heel & Schatz, 2005) as the filling degree of the reconstructed object in the volume. In order to compensate for such additional effects, we estimated the effective sample sizes from spherical masking. We minimized the deviation between the true distribution and the sample-size corrected distribution by testing systematically effective sample size factors between 0 and 1 and computing the Kolmogrow-Smirnow distance of the two distributions (Massey, 1951) (see Methods for more details). Using this approach, we found that a soft circular mask of volume size diameter reduces the effective sampling to 70% of the initial sample size (Supplementary Figure 2a) Figures 3 and 4). In turn, these tests also show that any map manipulations e.g. by user-defined masking will affect the FDR-FSC measurements and should be avoided. In conclusion, operations that introduce additional dependencies between half-maps such as volume symmetrization and masking require compensation by using a sample size correction factor.
Based on the proposed FDR-FSC approach, the resolution estimate will now be assigned to the highest spatial frequency that contains significant signal at 1% FDR (from now on referred to as FDR-FSC). For conventional FSC computation, the halfmaps are masked by the user, whereas for the here proposed FDR-FSC approach the untreated half-maps are used for FSC determination. The resolution value will be assigned to the resolution shell that crosses the significance threshold for the first time  (Table 1). In conclusion, the determined resolution estimates using 1 % FDR-FSC are in close agreement with previously reported resolution values from the EMDB and can be computed in a fully automated fashion without any user interference.

Application of FDR-FSC to local resolution estimation
Cryo-EM maps usually exhibit local variations of resolutions and estimating local resolutions has become critical in order to evaluate the structures (Cardone et al., 2013;Kucukelbir et al., 2014;Vilas et al., 2018). We reasoned that the adaptive FDR-FSC thresholding approach could provide benefits over fixed threshold approaches.
Therefore, we extended our approach by computing local FSC curves and tested it on cases with large resolution variations. For initial testing, we simulated a composite

Application of FDR-FSC to directional resolution estimation
In addition to local resolutions, directional resolutions have been recently evaluated to assess the effect of preferred particle orientations (Zi Tan et al., 2017). In practice, the FSC curve from a conical 3D Fourier transform is calculated including voxels of a specified angle. In analogy to local resolution windows, due to the limited number of Fourier coefficients inside a conical volume, directional FSCs suffer from poor statistics and therefore the FDR-FSC approach could provide similar benefits as demonstrated in the case of local FSCs. We tested the approach in more detail using three different cases: a map of the soluble portion of the small influenza hemagglutinin (HA) trimer with highly preferred orientations (EMD8731, Figure 4a) (Zi Tan et al., 2017), a highly symmetric apo-ferritin map (EMD0144, Figure 4b) (Zivanov et al., 2018) and an asymmetric map of -secretase (EMD3061, Figure 4c) (Bai et al., 2015).
Inspection of the directional resolution plots reveals that the 0.143 FSC-thresholds tend to give more optimistic resolution estimates compared with the 1% FDR-FSC approach. This effect is more pronounced at lower resolutions. As a result, the 1% FDR-FSC criterion displays stronger resolution differences with lower resolutions up to 8 Å in the vertical direction and higher resolutions up to 4.2 Å in the horizontal direction (Figure 4a and 4c). Due to highly preferred orientations, such a result can be expected in this sample (Zi Tan et al., 2017). In contrast, the directional values derived by the 0.143 FSC threshold are not able to resolve the differences in directional resolution apparent in the HA map. For the 1.6 Å resolution map of apoferritin, clear directional resolution differences are not displayed due to high symmetry and homogenous particle orientations (Figure 4b). Similarly, the asymmetric map of -secretase only exhibits minor directional resolution differences (Figure 4c). We attribute the larger detected resolution differences by the FDR-FSC approach to the more faithful resolution determination in cases of smaller number of Fourier coefficients in analogy to the benefits observed for local resolution. In conclusion, FDR-FSC is suitable for local resolution determination as well as for directional resolution measurements as it estimates resolution at different noise levels more robustly in comparison with fixed threshold approaches.

Discussion
Resolution estimation is one of the essential tasks to assess the experimental quality and confidence for the interpretation of cryo-EM maps. Therefore, an automated procedure with least user interference delivering robust results is desirable. Here, we propose a new approach of thresholding of FSC curves by non-parametric permutation sampling followed by FDR control and show that this approach is suitable for reproducible resolution determination. The approach does not have any free parameters. The only additional information that needs to be supplied is the volume symmetry applied. This way, the FDR-FSC approach enables resolution estimation without masking and thus eliminates the requirement of the mask including the optional deconvolution process (Chen et al., 2013). Although mask-free approaches have been proposed (Sindelar & Grigorieff, 2012), they require estimates of the expected molecular volume, which corresponds to information that may not be routinely available from the cryo-EM experiment in particular when heterogeneity is involved. Further application to local resolution estimation showed that judged by visual map features the adaptive FDR-FSC method captures locally high-resolution features in the core of protein structures as well as lower resolution in the periphery of a protein complex. Complementary to local resolutions, directional FSCs can also be analyzed by FDR-FSC and resulted in well-estimated resolutions in cases of anisotropic reconstructions. The algorithm is implemented in a Python program that takes a few seconds for small maps (<200 pixels volume size) up to a few minutes for big maps (>400 pixels volume size).
One potential conceptual advantage of the proposed FDR-FSC approach is that the method only assesses whether any significant correlations beyond random noise can be detected in a given resolution shell. This is different to the commonly used approach of determining the absolute spectral information content of the structure of interest. In this way, the FDR-FSC approach tolerates even high noise levels (see  (Benjamini & Yekutieli, 2001). Both of these issues occur in FSCs computed from cryo-EM maps. First, correlations and dependencies between neighboring resolution shells in cryo-EM are known to exist due to uncertainties in alignment and are propagated to the reconstruction volume.
Second, limited number of Fourier coefficients are present in smaller windows in cases of FSC-based local and directional resolution estimation. The FDR control procedure automatically adjusts the effective p-value threshold to the number of available resolution shells. For such applications, however, we also showed that the window size remains a free parameter that affects the computation of local resolutions regardless of which thresholding procedure is being used (see Supplementary   Figure 6). This is a general property when local FSC calculations are computed and is independent from the respective threshold criterion as it is already elaborated in detail in (Cardone et al., 2013). Irrespective of the window size, the here presented FDR-FSC method has advantages over 0.5 FSC, as it is less sensitive to differences in noise levels in the center and periphery of the protein that can affect the resolution determination.
Given its robustness and minimal input requirement, it is conceivable that the proposed FDR-FDR approach can be used to re-evaluate many deposited EMDB structures when raw half-maps were available. In our evaluation, we found very good agreement between the reported resolution values, mostly based on the 0.143 FSC criterion, and the FDR-FSC resolution values. This correlation shows that the overwhelming majority of depositors submits their resolution value with care in correspondence to the presented visual features such as beta-strand separation and side-chain densities. For future EMDB depositions, the FDR-FSC method could be used as an additional resolution assessment in the context of other standard validation tools. Due to the increased statistical robustness, we see particular use for extending resolution determination into more challenging applications such as local and directional resolution where resolution reporting is often critical for interpretation but technically less standardized. In the current manuscript, we presented the FDR-FSC method with a focus on cryo-EM map evaluation. It has to be noted that the proposed framework is applicable to 2D Fourier ring correlations as used in super-resolution microscopy (Nieuwenhuizen et al., 2013). Due to the demonstrated ease of use and robustness, we anticipate a common applicability of the FDR-FSC method to image analysis and processing tools in and beyond the cryo-EM field.

Permutation testing for Fourier shell correlation coefficients
Let ' ( and ' ( complex random variables and we denote the corresponding Fourier coefficients at the specific location + , = 1, … , , in resolution shell of half-map 1 and 2, where ∈ ℕ is the number of Fourier coefficients in the respective shell. Due to the Friedel symmetry inherent to Fourier transforms, we restrict ourselves to one half of all locations + . The pair is simply represented as the complex conjugate. We denote the Friedel symmetry corrected sample size with . In the case of two independent reconstructions, dependence between ' ( and ' ( is introduced through an effect of the signal at position + , which we denote with ( + ) ∈ ℂ and which is the same in both volumes. Thus, Fourier coefficients are usually modelled as ' ( = + + : ( + ) and ' ( = + + ; + (1), where : ( + ) and ; ( + ) are complex valued noise variables ( Van Heel & Schatz, 2005). No assumptions about the distribution of the noise in resolution shell are made. Furthermore, we denote by < = ( ' = , … , ' > ) and similarly < = ( ' = , … , ' > ) the complete set of resolution shells.
Initially, we have to clarify whether the Fourier coefficients in the respective resolution shell are dependent with respect to the location + , i.e. that groupings ( ' ( , ' ( ) show statistical dependencies. Expressed in statistical terms: the null hypothesis @ is that ' ( and ' ( are independent of each other, and the alternative hypothesis A is that there are dependencies. In order to test such a hypothesis, correlation coefficients can be used (Lehmann & Romano, 2005 and it becomes clear that < , < is a real valued random variable. Statistically, < , < is the estimator of the true value of the Fourier shell correlation, which we denote from now on with FSC. It is important to emphasize here that we are testing for dependencies by using a correlation coefficient; i.e. the null hypothesis is not FSC = 0.
In order to test @ , we calculate a p-value as follows. Using we denote the estimated value for the FSC in resolution shell . Moreover, using we denote the random variable describing the FSC as above. The p-value of this observation is then given as the probability that is bigger than under the null hypothesis, i.e.
where ℙ is the true probability distribution under the null hypothesis.
A permutation test is performed as follows. Under @ , the Fourier coefficients ' ( and ' ( are independent for = 1, … , . Newly paired samples of Fourier coefficients can be generated by permutation of < . Thus, the sampled null distribution is the distribution of from independent half-maps. For a detailed discussion regarding permutation of correlation coefficients we refer to (DiCiccio & Romano, 2017 , where is the cardinality of the set , i.e. the number of permutations selected for the Monte-Carlo estimate. It is important to note that the null distribution of the estimated by the permutation approach is not necessarily the distribution of the in the absence of any signal. As we are permuting in the presence of possible signal, which adds additional variation, the permutation distribution of the will have larger tail probabilities than the distribution without any signal. Consequently, this could give rise to more conservative resolution estimates. However, simulations in the presence of signal at high and low signal-to-noise ratios showed that this effect seems to have less practical relevance (Supplementary Figure 7). Moreover, the amount of signal in the most critical resolution shells close to the true resolution of the structure is low and will thus have limited influence on the actual distribution. As a compromise between computational efficiency and statistical accuracy, we restrict the number of permutations to a maximum of 1,000. Permutations are only performed for resolution shells with an effective sample size >10, which allows for more than 1,000 permutations. In cases of insufficient sampling below 10, the program will generate a warning message.

Multiple testing correction
Testing the complete set of resolution shells of a volume results in a multiple testing problem. In order to correct for the multiple testing problem, estimated p-values are subsequently adjusted for the false discovery rate (FDR) of detected resolution shells, i.e. the maximum amount of false positive resolution shells. The FDR is given as with ∈ ℕ @ the number of false rejections, ∈ ℕ @ the number of true rejections and () denotes the expectation value. For a more detailed description of the exact adjustment in the context of cryo-EM we refer to . The false discovery rate (FDR) is controlled with the Benjamini-Yekutieli procedure (Benjamini & Yekutieli, 2001), which controls the FDR under arbitrary dependencies.

Effective sample size corrections
Correction of the sample size has been found to be an important factor in presence of symmetry and masking ( Van Heel & Schatz, 2005) (Results section), as this leads to additional dependencies between Fourier coefficients and thereby reduces the true sample size to an effective sample size jkk . The effect of imposed symmetry during reconstruction can be conservatively corrected by where is the number of Fourier coefficients in the respective resolution shell and lm is the number asymmetric volume units at the given symmetry ( Supplementary   Figure 1b). Effective sample sizes are incorporated in the permutation framework by sub-sampling of Fourier coefficients. Effects of masking on the effective sample size are more complicated and depend on the specific shape and volume of the mask.
Here, we estimate the effective sample size jkk ≈ * by finding the factor ∈ (0,1] through minimization of the mean Kolmogorov-Smirnov distance q over the resolution shells. The two-sample Kolmogorov-Smirnov statistics (Massey, 1951) q is a measure of similarity of two empirical cumulative distribution functions (ECDF) and is defined as where in our setting A,q ( ) denotes the ECDF of the permutation approach. This approach is estimated using an effective sample size jkk ≈ * , and I,m+z ( ) an estimate of the true cumulative distribution function in presence of the respective mask. I,m+z ( ) can be obtained by repetitive simulation of two masked noise maps and subsequent FSC calculation. q is then calculated as the mean of the Kolmogorov-Smirnov statistics q,' of the resolution shells , i.e.
where is the number of resolution shells. An estimate for is then given by = arg min q∈(@,A] q 10 . In the presented algorithm, we apply a soft circular mask for global resolution estimation by default, as the resulting effect can be corrected by an effective sample size of jkk ≈ 0.7 , i.e. q was minimized for ≈ 0.7 (Supplementary Figure 2a), which allows accurate calculation of tail probabilities for various box sizes (Supplementary Figure 3). Application of windowing functions, as used for local resolution estimation (see below), also leads to reduced effective sample sizes. In a similar way as for the soft circular mask, we found that a Hann-window leads to an effective sample size jkk ≈ 0.23 (Supplementary Figure 2b, Supplementary   Figure 4).

Local Resolution estimation
Local resolutions are estimated by a moving window across both half maps and subsequent calculation of local FSC thresholds (Cardone et al., 2013). In order to account for high-resolution artifacts introduced through spectral leakage, a Hannfunction is used as a windowing function. The masking effects, which are introduced from the Hann window, are corrected as described above with a correction factor ≈ 0.23. Moreover, effects of symmetry do not have any influence on the effective sample sizes when the size of the sliding windows is smaller than the size of the asymmetric volume unit in the map. In order to speed up the calculations, a step-size option is implemented allowing movement of the sliding window of more than a single voxel.
Moreover, in order to avoid repeating permutations on the same map, the permutations are only done on 10 random locations of the sliding window. The resulting samples of the null distribution are merged and subsequently used for pvalue calculation at all locations. We found that window sizes around seven times the estimated global resolution provide a reasonable compromise between locality and resolution sampling in Fourier space for most cases. Presented experiments were performed using a window size of 25 pixels.

Directional resolution estimation
The implementation of the directional resolution estimation follows Lyumkis and colleagues (Zi Tan et al., 2017). For each direction, the FSC is calculated by taking those voxels from each half-map that are included by a specified angle at the respective direction. This results in rotating an inverse cone over one half of the 3D

Software
The procedures are implemented together with the previously published Confidence Map tools in a GUI based software named SPoC -Statistical Processing of cryo-EM maps (Supplementary Figure 8). Code is written in Python3 based on NumPy (Van Der Walt et al., 2011), matplotlib (Hunter, 2007), SciPy (Oliphant, 2007), mrcfile (Burnley et al., 2017) and parallelized with multiprocessing. The software is available at https://github.com/MaximilianBeckers/SPOC. Figures were prepared with UCSF Chimera (Pettersen et al., 2004).