An advanced workflow for single-particle imaging with the limited data at an X-ray free-electron laser

The 3D structure of bacteriophage PR772 was determined from a single-particle imaging experiment at the Linac Coherent Light Source, using the limited number of data, with a spatial resolution of 6.9 nm.


Introduction
The recent worldwide outbreak of the COVID-19 pandemic has indicated urgency for the development of complementary imaging techniques for the study of virus structures at high resolution. Currently, cryogenic electron microscopy (cryo-EM) outperforms other methods of analysis of single virus particles, although to determine the structure one needs to freeze the samples at cryogenic temperatures. Understanding the nature of the same viruses at room temperature as well as in aerosols may provide information on the actual size distribution of virus particles in aerosols, a factor that may affect infection of humans by these viruses.
As suggested about two decades ago, intense femtosecond pulses from X-ray free-electron lasers (XFELs) may outrun conventional radiation damage, found at synchrotron X-ray sources, or the Coulomb explosion, found at the intensities of XFELs. They therefore allow, in principle, atomic resolution structure determination of isolated macromolecules in their native environment at room temperatures (Neutze et al., 2000). After the start of the first hard X-ray FELs (Emma et al., 2010;Ishikawa et al., 2012), it became clear that atomic resolution could be achieved with intense X-ray pulses from an XFEL on micro-and nano-crystals of proteins, leading to the serial femtosecond crystallography (SFX) technique Aquila et al., 2012;Boutet et al., 2012). However, progress toward high-resolution 3D electrondensity images of non-crystalline biological particles using the single-particle imaging (SPI) approach has been slow compared with SFX Seibert et al., 2011;Hantke et al., 2014;Ekeberg et al., 2015).
There are several reasons for the lower resolution achieved in SPI experiments in comparison with the SFX technique. The most important include: lack of crystalline periodicity to amplify the signal, weak single-particle signal from non-periodic nanoscale objects compared with instrumental background, limited number of usable data frames collected and heterogeneity of the samples. Radiation-damage processes initiated by X-ray photoionization may also play an important role at high resolution, but this limitation is minimized at XFEL sources by performing experiments in the 'diffract and destroy' mode (Chapman et al., 2006). One may expect that in order to further enhance resolution, one could increase the power of XFEL pulses to boost the signal-to-noise ratio for bio-particles, such as single proteins or virus particles that scatter very weakly. Unfortunately, increased XFEL fluence strips electrons from atoms more efficiently and the scattered power from the bound electrons (which contains structural information) does not increase in proportion to the X-ray fluence (Lorenz et al., 2012;Gorobtsov et al., 2015). Decreasing the pulse length below 1 fs may help to outrun Auger decay, but further decrease of the pulse duration may not lead to the desired suppression of all ionization channels .
Beside these fundamental limitations, there are a few other limiting factors. Background scattering originating from the beamline obscures the weak scattering signals from biological samples, and fluctuations of the beam position and intensity can cause additional variables that need to be accounted for in reconstruction algorithms. Moreover, single-particle sample delivery remains a challenging topic. In addition, despite strong development of detector technology [see, for example, Allahgholi et al. (2019)], the dynamic range of the present detectors is still not sufficient to capture a full diffraction pattern in a strong single shot. To tackle all these limitations of the SPI technique and to push the methodology further, a dedicated SPI consortium was formed at the Linac Coherent Light Source (LCLS) .
Several results were reported from the SPI consortium, covering both hard and soft X-ray experiments and focused primarily on well characterized viruses with sizes of a few tens of nanometres. Different analysis methods have been applied to virus structure determination [see, for example, Munke et al. (2016), Reddy et al. (2017), Kurta et al. (2017), Rose et al. (2018) and Shi et al. (2019)]. In the present work, we advance the previous studies performed on bacteriophage PR772 (Rose et al., 2018) by extending the general analysis pipeline of the SPI data as first proposed by Gaffney & Chapman (2007). The main upgraded features include: a classification method based on the expectation-maximization (EM) algorithm, developed in cryo-EM (Dempster et al., 1977), and mode decomposition for the final virus structure determination (Vartanyants & Singer, 2010;Thibault & Menzel, 2013). EMbased algorithms were first applied in SPI data analysis for 3D orientations recovery (Loh & Elser, 2009;Ayyer et al., 2016). Below we present a detailed description of all the steps that allowed us to obtain the electron density of the PR772 virus with an improved resolution compared with the previous SPI studies on the same virus (Rose et al., 2018).

Experiment
The experiment was performed in the Atomic, Molecular and Optical Science (AMO) instrument (Bozek, 2009;Osipov et al., 2018) at the LCLS at SLAC National Accelerator Laboratory using the LAMP end station [see Li et al. (2020) for experimental details]. PR772 bacteriophage growth and purification were carried out essentially as previously described (Reddy et al., 2017). Viruses in ammonium acetate volatile buffer were aerosolized in a helium environment using gas dynamic virtual nozzles (GDVN) that were 3D printed via 2 photon polymerization photo-lithography with a Nanoscribe Professional GT printer . Laboratory measurements showed reproducible jet diameters in the range of 0.5-2.0 mm (DePonte et al., 2008;Weierstall et al., 2012). The particles then passed through a differentially pumped skimmer for pressure reduction and were then injected/focused into the sample chamber using an aerodynamic lens injector Hantke et al., 2014). The focused particle stream intersected the focused and pulsed XFEL beam.
The LCLS had a repetition rate of 120 Hz, for this experiment the average pulse energy was $4 mJ, with a focal diameter of $1.5 mm and a photon energy of 1.7 keV (wavelength 0.729 nm). Diffraction patterns were recorded by a pnCCD detector (Strü der et al., 2010) mounted at a 0.130 m distance from the interaction region. In the experiment, a silver behenate salt was used as a calibration agent to determine sample detector distance and panel position. The scattering signal was only recorded by one of the two detector panels (one panel was not operational during the experiment because of an electronic fault). The size of the working panel was 512 by 1024 pixels with a pixel size of 75 Â 75 mm 2 , with research papers IUCrJ (2020). 7, 1102-1113 Assalauova et al. An advanced workflow for single-particle imaging the long edge closest to the interaction point. In the middle of the experiment (at run 205) the detector panel was moved 1 mm vertically up relative to the incoming X-ray beam to reduce background scattering. Our work is based on the data obtained from this panel covering a part of reciprocal space, as shown in Figs. 1(a)-1(d).

Initial classification steps
SPI data analysis involves many subsequent steps leading to the 3D reconstructed particle structure from a large set of 2D diffraction patterns (Gaffney & Chapman, 2007). Improvements in the data-analysis pipeline at early stages can result in a significant enhancement in reconstruction quality. Therefore, several pre-processing methods were developed to avoid experimental artifacts on the collected diffraction patterns. The pre-processing stages include: hit finding, background correction, beam-position refinement and particle-size filtering.
The initial experimental dataset, as collected at the AMO instrument, consists of $1.2 Â 10 7 diffraction patterns (9 Â 10 6 patterns before and 3 Â 10 6 patterns after moving the detector panel) (see Li et al., 2020). The hit finding was performed using the software psocake in the psana framework (Damiani et al., 2016). As a result, 1.9 Â 10 5 diffraction patterns were identified as hits from the initial set of diffraction patterns (see Table 1) and the signal from these hits was converted to photon counts. Examples of selected diffraction patterns are shown in Fig. 1.
Visual inspection of the selected patterns revealed the presence of additional instrumental scattering near the center of the diffraction patterns [see Fig. S1(a) in the Supporting information]. The scattering was caused by the interaction of the tails of the FEL beam with the upstream apertures or from the sample injector. To subtract this additional signal, histograms of intensity for each pixel were analyzed (Figs. S1 and S2). It was assumed that for pixels with zero photon counts, the histogram of intensity distribution should have a Gaussian form and a mean value equal to zero. The mean value of this distribution was subtracted from each diffraction pattern. Subtraction of this additional scattering contribution was an important step prior to finding the beam center position (see the Supporting information for details).
The particle-size filtering was performed in two steps. First, the power spectral density (PSD) function, i.e. angular averaged intensity, of each diffraction pattern was fitted with the form factor of a solid sphere in a wide range of sizes from 30 to 300 nm. Next, a fidelity criterion was introduced to distinguish research papers 1104 Assalauova et al. An advanced workflow for single-particle imaging IUCrJ (2020). 7, 1102-1113

Figure 1
Examples of diffraction patterns from the SPI experiment. (a), (b) Diffraction patterns corresponding to a single PR772 virus hit by an XFEL beam. (c) A diffraction pattern corresponding to a non-single hit which was sorted out from the analysis at the classification step. (d) A sum of 1.9 Â 10 5 diffraction patterns identified as hits. White regions in the center of the diffraction patterns as well as white stripes correspond to a mask introduced to reduce artefacts owing to the signal exceeding the detector capabilities. The mask was the same before and after the move of the detector panel. (e) The PSD function of the scattered intensity for the sum of all diffraction patterns identified as hits collected in the experiment. The signal until the corner and edge of the detector is indicated by the green and black vertical dashed lines, respectively. For orientation determination the data until momentum-transfer values of 1 nm À1 (shown by the red dashed line) were used. between successful and unsuccessful size findings (see the Supporting information for details). As an outcome of this analysis, the number of selected diffraction patterns was reduced from 1.9 Â 10 5 to 1.8 Â 10 5 , which were used for the final particle-size filtering. A histogram for different particle sizes in the selected size range is shown in Fig. 2. An extended range of sizes observed in this figure corresponds to clusters of particles stuck together and the varying thickness of a hydration layer. One can also identify a peak in this histogram in the range from 55-84 nm. This size range agrees well with the expected virus particle size of $70 nm (Rose et al., 2018;Reddy et al., 2019). The considered range contains 1.8 Â 10 4 diffraction patterns, which were selected for further single-hit classification (see Table 1).

Single-hit diffraction patterns classification
The key step of data selection in this work was single-hit classification. The angular X-ray cross-correlation analysis classification that was used in the previous work (Bobkov et al., 2015;Rose et al., 2018) was not as effective with the present experimental data owing to the absence of the scattering signal on one half of the detector and the low single-hit rate.
The single-hit classification approach was based on the EM algorithm (Dempster et al., 1977;Loh & Elser, 2009;. This algorithm allows for unsupervised clustering of data when neither initial data assignments to clusters nor cluster parameters are known. The dataset is distributed into a pre-defined number of clusters and cluster parameters are retrieved automatically at the same time. Later, manual input is used to classify each of the cluster models, based on symmetry considerations or expected fringe patterns, in order to perform data selection. This algorithm is commonly used in cryo-EM for unsupervised single-particle clustering [see, for example, the software RELION (Scheres et al., 2005)]. Low contrast and low signal-to-noise level are common problems for cryo-EM and SPI so we implemented in this work the original EM algorithm developed in cryo-EM for clustering of SPI data.
Unsupervised EM clustering is an iterative process; the algorithm starts from a random model for each class. At each iteration the probability for each diffraction pattern to belong to a certain class is calculated. To accommodate for random particle orientations in the SPI experiment, the cluster model is compared with the 2D diffraction pattern rotated in-plane. After probabilities are evaluated, a new cluster model is calculated by weighted averaging of patterns that belong to each cluster in each orientation. The weights are defined by computed probabilities. In fact, the algorithm imitates the expand-maximize-compress (EMC) algorithm (Loh & Elser, 2009) but, instead of rotation in 3D space, an in-plane rotation and cluster distribution are used (Loh, 2014). When the EM algorithm converges, the supervised class selection takes place. The cluster models that correspond to single hits of an investigated particle are selected manually by an expert. As an example, in one of the EM clustering we selected classes 1 and 2 as classes containing six high-contrast fringes that we attribute to scattering from a PR772 virus particle [see Fig. 3

(a)].
To make classification more accurate we performed five independent EM clusterings starting from random cluster models. Each EM clustering produced slightly different results, which are summarized in Table 1. The intersection of all the results was considered as a stable single-hit selection [see Fig. 3(b)]. Finally, we ended with the dataset containing 1085 single-particle patterns (see Table 1). The PSD function of an average of all the selected diffraction patterns is shown in Fig. 3(c). The scattering signal from the virus particle extends up to the momentum-transfer value q max = 1 nm À1 in reciprocal space, which corresponds to a resolution (2/q max ) of 6.3 nm in real space.
A manual search and selection of single hits from the dataset of 1.9 Â 10 5 diffraction patterns produced a new dataset containing 1393 diffraction patterns (Li et al., 2020). From this selection, 574 diffraction patterns are also present in our EM-based selection. The PSD function for the manual selection shown in Fig. 3(c) has a lower contrast without visible fringes. We attribute this mostly to the absence of the research papers IUCrJ (2020). 7, 1102-1113 Assalauova et al. An advanced workflow for single-particle imaging 1105 Figure 2 A particle-size histogram after the PSD-function filtering. The blue area corresponds to 1.8 Â 10 5 diffraction patterns, satisfying the fidelity score criterion FS > 1.05 (see the Supporting information for details). The orange area corresponds to manually selected single-hit diffraction patterns. A range of particle sizes from 55 to 84 nm (1.8 Â 10 4 diffraction patterns) was selected for the further SPI analysis (green dashed area).

Table 1
Datasets selected at different stages of the analysis: hit finding selection, PSD-fitting quality filtering, particle-size filtering and single-hit diffraction pattern selection.
The percentage of the chosen dataset to the initial one S 0 is given in parentheses. size-filtering step in manual selection. The virus size fluctuations could be caused by a slight change in hydration-layer thickness, variation of the position where the particles interact with the X-ray beam, slightly upstream or downstream of the focus, or a real sample size distribution of the PR772 virus. We note that in our previous work (Rose et al., 2018) the number of single-hit diffraction patterns was $7.3 Â 10 3 , which is about one order of magnitude higher than in this work. The reason for the smaller number of patterns was a combination of downtime owing to detector troubleshooting and time needed for sample-injection optimizations. The GDVN flow and pressures needed to be optimized because of the relatively large initial droplet diameters of $1-2 mm.

Orientation determination and background subtraction
For orientation determination we used the EMC algorithm implemented in the Dragonfly software package . This iterative algorithm successfully combines 2D diffraction patterns into the 3D intensity distribution of the PR772 virus shown in Fig. 4(a) (see the Supporting information for further details). The results presented in Fig. 4(a) clearly show that the recovered 3D diffraction intensity contains a substantial highmomentum background that may be caused by scattering on helium from the carrier gas. To reduce the influence of this background contribution, a background subtraction was applied. Several approaches for the background correction in SPI experiment analysis were developed earlier (Rose et al., 2018;Lundholm et al., 2018;Ayyer et al., 2019) and it was understood that the background-correction method may affect the reconstructed structure. Here we used a combined approach for the background determination, first a constant background was subtracted and then, on the reconstruction step, the contrast of diffraction patterns was additionally enhanced by applying deconvolution algorithms (see the next section for details). In the first step, the background level was   defined as the mean signal value in the high-momentum region of the 3D intensity distribution, free of particle-scattering contributions (see the Supporting information for details). Reciprocal-space data with and without background subtraction are shown in Figs. 4(a)-4(c). As seen in Fig. 4(c), the power-law dependence of data after the background subtraction is the same up to a high q range of about one inverse nanometre. This is in contrast to the data after EMC orientation determination [the blue curve in Fig. 4(c)], which are clearly saturated at high q range.
As a result of the background subtraction, more features and higher contrast were revealed in the high momentumtransfer region. The fringe visibility or averaged contrast hi was defined as Here i is the local contrast, I i max and I i min are the PSD-function values in the local maxima and following minima, respectively, and N is a number of pairs of maxima and minima considered in this analysis. In our case, the average contrast values hi were calculated for the first N = 6 pairs. For the experimental data before the background subtraction, we obtained the value hi = 0.41. The fringe contrast after the background subtraction showed significant improvement with hi = 0.58.
To proceed further with the analysis, low (q < 0.12 nm À1 ) and high (q > 0.93 nm À1 ) momentum-transfer regions were eliminated from the final data selection owing to scattering artefacts present in these regions. After the background subtraction, negative values of intensity were set to zero. A 3D intensity distribution in reciprocal space after the background subtraction is shown in Fig. 4(d).

Phase retrieval and PR772 virus structure
To determine the electron-density distribution of the PR772 virus, the phase retrieval was performed using the 3D intensity distribution of the virus in reciprocal space described in the previous section [see Fig. 4(d)]. Iterative phase-retrieval algorithms are based on the Fourier transform between real and reciprocal/diffraction space using two constraints: in reciprocal space the amplitude of the signal is set equal to the experimentally measured values and in real space a finite support of the particle is used (Fienup, 1982;Marchesini, 2007).
The electron-density distribution of the PR772 virus was obtained with the following steps. First, the central gap, originating from the initial data masking of the diffraction patterns, was filled [see Figs. 4(a) and 4(b)]. This was accomplished by running multiple 3D reconstructions of the virus and leaving the central gap of the diffraction volume to freely evolve. All reconstructions produced practically identical continuous 3D functions in reciprocal space that smoothly covered the central gap (see Fig. S10). This reciprocal-space distribution with the filled central part was used for the next reconstruction step. On that step, 50 successive reconstructions using a combination of continuous hybrid input-output (Fienup, 2013), error reduction (Fienup, 1982) and shrinkwrap algorithms (Marchesini et al., 2003) were performed. In addition to these algorithms, we used the Richardson-Lucy deconvolution technique (Clark et al., 2012). Application of this method allowed us to take into consideration additional background and substantially enhance the contrast of the reconstructed diffraction pattern to the value of hi = 0.87. As a result, we obtained complex-valued real-space images for each of the 50 reconstructions (see the Supporting information for details).
On the final step, to identify the electron density of the virus we used mode decomposition. We first introduced a density matrix in the form where (r) relates to complex-valued real-space images and the brackets h . . . i indicate ensemble averaging over different reconstructions. We further decomposed the density matrix into orthogonal modes n (r) using principal component analysis (Khubbutdinov et al., 2019): where n refers to eigenvalues of this decomposition. This approach is advantageous in comparison with averaging that would make important object features blurry and possibly affect the final resolution. By considering the zero mode only, unique features present in all other reconstructions will be represented. In practice, we used all 50 reconstruction results and performed mode decomposition. The fundamental mode (with the weight factor 0 of 99%) was considered as the final result of reconstruction. Such a high weight-factor value indicates that all reconstructions converged essentially to the same result with an uncertainty level of only 1%. To determine the electron density, we took an absolute value of this fundamental mode complex-valued image. A three-times upsampled version of the real-space electron density is shown in Fig. 5 (see the Supporting information for details). As seen in Figs. 5(a)-5(c), the retrieved electron density of bacteriophage PR772 shows the expected icosahedral structure. A 2D cut of the virus structure is presented in Fig. 5(d), where a higher electron density in a thin outer shell is well resolved and is attributed to the capsid proteins arrangement. As can be seen in Fig. 5(d), the electron density in the center of the reconstructed virus particle is reduced. The reason for this may be the heterogeneity of virus particles present in solution and injected in the X-ray beam.
Bacteriophage PR772, like other members of the Tectiviridae family, contains an inner proteolipid membrane that facilitates delivery of the viral genomic DNA during infection (Mä ntynen et al., 2019). When the virus binds to a bacterial host cell the inner membrane is extruded from one of the viral vertices to form a nanotube that facilitates genome delivery. Spontaneous release of the viral genome has been reported (Peralta et al., 2013;Santos-Pé rez et al., 2017). Our bacteriophage PR772 preparations were also analyzed by cryo-EM (Fig. S6). Under the conditions used for plunge freezing, we observed particles that were more rounded than icosahedral. Preliminary volume analysis suggests that some particles had released their DNA. These particles appeared to be 'triggered' during the freezing process since transmission electron microscopy imaging did not reveal this. It may be possible from minor differences in virus preparation or upon XFEL sample delivery that some of the virus particles were similarly 'triggered' to release their genomes during cryo-freezing, thus resulting in a decrease in inner electron density. Previous analyses of XFEL snapshots on the same virus also suggested that some particles exhibit decreased inner density Hosseinizadeh et al., 2017).
To identify the particle size, the electron-density profiles in different directions were investigated. In this work, we used the same approach as developed earlier (Rose et al., 2018) and analyzed the electron-density profiles in the directions from facet to facet and from vertex to vertex of the reconstructed virus particle [see Figs. 6(a) and 6(b)]. For the particle-size estimate, we selected the electron-density threshold value of 0.2 as it was considered in the shrink-wrap algorithm during the phase retrieval. From this criterion we determined the particle size in different directions (see Table 2 and the Supporting information). Thus, the obtained mean particle size was 61 AE 2 nm between facets and 63 AE 2 nm between vertices. These sizes correspond well to the initial range of particle sizes (from 55 to 84 nm) considered at the initial research papers 1108 Assalauova et al. An advanced workflow for single-particle imaging IUCrJ (2020). 7, 1102-1113  classification step (see Section 4) and other XFEL data performed on the same virus (Kurta et al., 2017;Rose et al., 2018).
Similar to previous work (Kurta et al., 2017;Rose et al., 2018), we also observed a certain elongation of the particle shape, which might be inherently present in the viruses in solution or appear owing to the aerosol injection conditions. We defined the elongation of the particle by the following measure, where D max , D min and D mean are maximum, minimum and mean particle sizes, respectively. The virus structure obtained in the current work showed an elongation value of = 11% for sizes taken between vertices, which is similar to the result of the previous SPI experiments (Rose et al., 2018) with = 9%. The mean capsid thickness was obtained from the Gaussian fitting of the electron-density profiles with well pronounced features (Fig. S14) and was determined to be 7.6 AE 0.3 nm. The thickness of the capsid in recent cryo-EM studies (Reddy et al., 2019) was identified to be below 10 nm, which is in good agreement with the thickness determined in this experiment.
Finally, we determined the resolution of the reconstructed electron density of the PR772 virus. We used the Fourier shell correlation (FSC) approach (Harauz & van Heel, 1986) to determine the resolution of the reconstructed virus. For this method two independent sets of reconstruction were required, usually each is based on half of the available dataset. FSC measures the normalized cross-correlation coefficient between both reconstructions over corresponding shells in Fourier space. The results of the FSC analysis are shown in Fig. 7. To estimate the achieved resolution, we used half-bit threshold criteria ( van Heel & Schatz, 2005). Its intersection with the FSC curve gave a resolution value of 6.9 nm (see Fig. 7). The obtained result is better than the previously reported value of $8 nm for the same virus (Rose et al., 2018), though the number of diffraction patterns used for the final analysis was much lower (15% of the previous dataset). In this case, the resolution was limited by the number of diffraction patterns and moderate scattered intensity, whereas previous reconstructions were limited by the extent of the detector.

Summary and outlook
We have presented an implementation of the SPI data-analysis workflow from diffraction patterns measured at the AMO instrument at the LCLS to reconstruct the electron density of bacteriophage PR772. This work is an important methodological step for development of SPI data analysis. We implemented several methods into the workflow including EMbased classification and mode decomposition that were crucial for the high-resolution final reconstruction. Although only half of the detector was operational, implementation of all of these steps allowed us to determine the PR772 virus structure with a higher resolution as compared with previous SPI studies.
From the initial set of 1.2 Â 10 7 experimentally measured diffraction patterns, the final single-hit selection set contained 1085 diffraction patterns. About 53% of this final set was also present in the single-hit selection made manually. The number of diffraction patterns classified for further analysis is substantially lower than in the previous experiment (Rose et al., 2018), improvements in sample delivery to increase the number of diffraction patterns are, therefore, crucial for advancements in reconstruction resolution.
The combination of all methods implemented in the workflow allowed us to obtain the 3D electron density of the virus with a resolution of 6.9 nm based on the FSC analysis, which is better than that obtained in the previous studies (Rose et al., 2018). The obtained mean PR772 virus size in this experiment was 61 nm (between facets) and 63 nm (between vertices), like in the previous SPI experiments performed on the same PR772 virus (Kurta et al., 2017;Rose et al., 2018). We also observed a similar elongation of $11% of the virus structure as was determined in the previous experiments. Electron-density profiles of the reconstructed virus PR772 normalized to the maximum value for the cut between vertices (a) and facets (b). The horizontal black dashed lines denote a particle-size threshold of 0.2. The mean virus size is 63 and 61 nm for the distance between vertices and between facets, respectively. Table 2 Particle-size analysis from facet to facet and from vertex to vertex for the reconstructed PR772 virus.
All distances between facets and vertices are given in the Supporting information.
Size ( The presented research is another step forward in SPI data analysis. Implemented steps may become especially important when SPI experiments are performed at high repetition rate XFELs, such as the European XFEL (Decking et al., 2020;Mancuso et al., 2019) and LCLS-II (Schoenlein et al., 2015) facilities. The first experiments performed at the European XFEL demonstrated the possibility of collecting the SPI data at a megahertz rate (Sobolev et al., 2020), which might be crucial for future progress in SPI experiments performed at XFELs.