The Artiatomi package

The following chapter aims in presenting each individual application in the ARTIATOMI package. This description is kept basic and focuses on giving an overview of the algorithms used. The order in which to use the tools in a processing pipeline does not necessarily correspond to the order here. Especially after the subtomogram-averaging step the order and iterative refinement depends largely on the dataset and its quality and thus can’t be described in a generic way. When best to perform for example 3D classification depends on the dataset and is not in the scope of this description.

1. Image viewer and editor: aiClicker

aiClicker can visualize most common image formats that are used in the field of cryo-electron microscopy. It offers a specialized user interface for four different kinds of images: dose-fractionation movie-stacks, single frame images, tilt-series and volumes. The dose fractionation viewer can show the images individually or performs the accumulation on-the-fly with or without image shifts. Single frame images are usually black-and-white for cryo-em, but aiClicker also supports color images in this mode. The tilt-series viewer and editor allows for selecting and tracing gold beads, perform marker alignment, alignment refinement based on shift information of individual sub-volumes, visualizes the tilt-series alignment in 3D with and without a deformation model. Finally, the volume volumes shows either the entire reconstructed tomogram and allows to select and visualize sub-volume positions. Or it shows the results of sub-tomogram averaging with FSC curves and achieved resolution.

2. Align dose fractionation movie stacks to a tilt series: aiImageStackAlignator

aiImageStackAlignator takes as input a single movie stack (single particle use-case) or a list of multiple stacks that form a tilt-series. As the first chain link, aiImageStackAlignator not only is responsible for aligning the image frames but also performs the data import for ARTIATOMI and creates the first .info metadata files.

For tomography, two import options are available: MDOC files as provided by SerialEM or a CSV file providing all needed parameters.

MDOC-import:

For the MDOC-import, aiImageStackAlignator relies on two input sources: the actual MDOC file and an additional microscope settings file. aiImageStackAlignator looks for a file named ‘artiatomi.settings’ in the ARTIATOMI_PATH or uses the settings file provided as a parameter. These files contain microscope specific information such as Cs value, acceleration voltage, the defocus field to use in the MDOC file, scaling factors, etc. in order to correctly build all the needed metadata for a tilt-series. The movie stacks must be in the same folder as the MDOC file.

CSV-file import:

A text file that contains all necessary information to create a tilt-series from movie stacks:

The first line must contain the number of image files to process. The second line gives a list of provided parameters, separated by ‘;’. The following values are possible:

  • AccumulatedDose
  • AccumulatedExposureTime
  • AcquisitionIndex
  • Binning
  • Cs
  • Defocus
  • ExposureTime
  • Filename
  • FlipOnYAxis
  • ImageRotation
  • MagnificationAnisotropy
  • MagnificationAnisotropyAngle
  • PixelSize
  • TiltAngle
  • TomogramNr
  • Voltage

The following lines then contain the actual data as described in line 2 of the file, all values are separated by ‘;’.

Usually, the accumulated dose and the accumulated exposure time are linearly dependent but are both used as separated metadata. The reason is that in ARTIATOMI, the accumulated dose is defined as the dose accumulated before the current acquisition, the accumulated exposure time is the time including the current acquisition. We have thus a zero-based and a one-based parameter for dose dependent processing. For example, the deformation model is based on accumulated dose: as the accumulated dose before the first frame is usually zero, the deformation model is thus the identity for the first frame, which means no deformation at all.

The imported frame parameter values can be stored directly in the tilt-series metadata info-file. But we encourage the use of the second option, i.e. to provide a filename for these frame parameters. Doing so, it is easy to share the settings across different versions of the same tilt-series, for example a binned or a denoised version.

Each input movie stack is processed individually and independently of the other stacks. For frame alignment, three different algorithms are implemented: Global alignment, Refine (One-To-All approach), or both where first the global alignment is computed and then the Refine algorithm is applied.

The resulting tilt-series is stored in a user-defined pixel type. Whereas floating point values are most accurate, it is in most of the cases sufficient to store the images as unsigned short (USHORT) with a scaling factor of 100 or 1000.

Global alignment:

The global alignment is based on the algorithm as described in doi:10.1038/nmeth.2472 and computes a matrix of CC-coefficients and a set of frame shifts is determined that minimizes the global shifts. As for large frame counts the matrix can get very large and many correlation operations have to performed, a block-wise strategy can be applied so that the CC-matrix is filled block wise, i.e. with correlation from only neighboring frames.

Refine:

The refine algorithms works in the following iterative way: first all frames except one are summed up with an initial shift, either zero or given from a previous global alignment. The left-out frame is then correlated with this sum and its shift is determined. This is repeated for all frames in the stack in multiple iterations. For extreme-low dose stack-frames as they occur in tomography, this algorithm usually performs better than the global alignment.

Skip:

For completeness also another option is mentioned: using no alignment algorithm at all, only the shifts from metadata are used to sum up the individual frames.

For correlation, aiImageStackAlignator uses phase correlation with a band-pass filter provided by the user. The correlation peak can be interpolated which results in sub-pixel shift determination.

The user can define a maximum allowed shift that discards individual frames that have a larger frame-to-frame shift than the limit. A large frame-to-frame shift implies that the electron beam or the sample were also moving during the acquisition of that frame and thus cancel out all high-frequency information. Further, one can define a fix number of frames to discard from the beginning of each movie stack, for example if one knows that the first frame always contains some artifacts or has too much shift. Removing frames does not alter the determined shifts, the frames are still used in the correlation routines. It is only in the last summation step where these frames get discarded, which ensures that a tilt-series with discarded frames has exactly the same geometry as the non-discarded one. The tilt-series’ metadata keep track of the number of frames used and dose-weighting can be adjusted accordingly.

3. Coarse pre-alignment: aiCoarseTiltSeriesAlignator

Having a tilt series, further alignment processing mainly depends on two parameters: image shifts and image rotation. Image rotation is assumed to be known with a precision +/- 1 degree for a given microscope/magnification and is usually stored for example in the MDOC file for each tilt series.

aiCoarseTiltSeriesAlignator then determines the shifts from tilt to tilt using phase correlation and tilt-angle dependent image stretching. Starting from the zero tilt once in the negative tilt range and once in the positive tilt range, each tilt is correlated with all frames with a lower tilt in the configurable angular range (usually around 7 degrees). If multiple frames fall into the angular range, the shift determined by each frame is averaged and weighted by the correlation coefficient. In case aiCoarseTiltSeriesAlignator is not able to find any meaningful correlation peak inside the allowed max-shift radius in the correlation map, an image can be automatically marked as “invalid” and is thus ignored in any further processing. This can happen for example if a grid bar enters the field of view in higher tilts.

In theory, a first reconstruction could be performed only with this coarse alignment.

We advise to check the quality of the tilt series and the processing done so far in aiClicker at this stage.

4. Marker tracking and tilt series alignment: aiBeadTracker or aiPatchTracker

For a precise reconstruction a tilt series must be aligned as good as possible and the coarse pre-alignment is not sufficient. In case that the tilt series contains enough visible gold bead markers, this fine alignment is done with aiBeadTracker. If no markers or not sufficient markers are available, aiPatchTracker is used.

aiBeadTracker:

aiBeadTracker works in the following fashion: On each tilt, marker candidates are determined. The user provides a range of gold bead sizes to scan and for all sizes, a marker reference is created, i.e., a black circular surface on white background with varying size. Each frame is then correlated with that reference, additionally with defocus dependent CTF distortions. As a strong low-pass filter is applied, the rough initial defocus values are accurate enough for this step. Further a gradient image is correlated with a ring-model of a gold bead, and if all three images result in a peak, we have strong confidence in the presence of a gold bead at the given position.

Finally, all marker candidates are checked for plausibility, which means that given the coarse pre-alignment, a track can be formed from one tilt to the other for a given marker position, whereas it is not necessary that a marker is visible on all tilts. This filters out any wrong candidates and usually only gold beads remain. Some other high contrast image features could remain as wrong markers, though.

aiPatchTracker:

Similar to aiCoarseTiltSeriesAlignator, aiPatchTracker uses phase correlation and tilt angle dependent-stretching. The difference is, that it uses image patches instead of the whole image so that individual image features are tracked from one tilt to the other. The user determines the patch size and the patch density.

Both tools provide the option to perform a marker alignment after marker determination / patch tracking, together with an additional filtering step where markers with large error score are discarded.

Marker alignment optimizes the following parameters:

  • Image shift in X and Y
  • In-plane image rotation (individually per tilt or constant for the entire tilt series)
  • Tilt angle, for each tilt with the option to penalize deviation for lower tilts
  • Beam declination angle
  • Magnification, with the option to fix the magnification of the first recorded frame to 1

So that the projection of a marker at global coordinate

\left(
\begin{array}{c}
X\\
Y\\
Z\\
\end{array}
\right)

and the matrices:

flipZ = \left(
\begin{array}{ccc}
1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & -1\\
\end{array}
\right)
\\
imgRot= \left(
\begin{array}{ccc}
cos(\psi) & -sin(\psi) & 0\\
sin(\psi) & cos(\psi) & 0\\
0 & 0 & 1\\
\end{array}
\right)
\\
tilt= \left(
\begin{array}{ccc}
cos(\theta) & 0 & sin(\theta)\\
0 & 1 & 0\\
-sin(\theta) & 0 & cos(\theta)\\
\end{array}
\right)
\\
beamDecl= \left(
\begin{array}{ccc}
1 & 0 & 0\\
0 & cos(\varphi) & -sin(\varphi)\\
0 & sin(\varphi) & cos(\varphi)\\
\end{array}
\right)
\\
adj= \left(
\begin{array}{ccc}
0 & 1 & 0\\
-1 & 0 & 0\\
0 & 0 & -1\\
\end{array}
\right)
\\
proj= \left(
\begin{array}{ccc}
1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 0\\
\end{array}
\right)
\\

gives us the pixel

\left(
\begin{array}{c}
px\\
py\\
\end{array}
\right)
=
\left(
\begin{array}{c}
imgSizeX/2+0.5\\
imgSizeY/2+0.5\\
\end{array}
\right)
+
\left(
\begin{array}{c}
shiftX\\
shiftY\\
\end{array}
\right)
+{1\over mag}
*
imgRot*proj*adj*beamDecl^T*tilt*beamDecl*flipZ*
def(
\left(
\begin{array}{c}
X\\
Y\\
Z\\
\end{array}
\right), accumulatedDose
)

With tilt angle θ, beam declination angle φ, image rotation angle ψ, magnification ratio mag, a deformation model def [https://doi.org/10.1016/j.jsb.2018.02.001] that displaces (X,Y,Z) depending on the position and the accumulated dose and estimated pixel value (px,py) is closest to the actually selected marker position for each individual tilt and each marker. Also note that imSizeX and -Y are the un-binned image sizes, the additional 0.5 sets the rotation center to the center of the central (un-binned) pixel. All positions given in pixel coordinates are assumed to be compensated for magnification anisotropy.

The flipZ and adjust matrix is necessary as the marker model is defined with a tilt axis on the Y-axis, the images in a tilt series on the other hand are usually aligned in a way that the X-axis is the tilt axis. These matrices simply switch from our global coordinate system to the different coordinate systems used in papers and publications for alignment.

Note that marker alignment with a deformation model is not available at this stage but only through aiClicker and aiToolbox. The reason for this is, that a deformation model allows to easily overfit a marker model even to erroneous markers. The user must thus first evaluate the quality and precision of the tracked markers before aligning with deformation models.

The fine alignment information can be stored directly in the tilt-series metadata info-file. But we encourage the use of the second option, i.e. to provide a filename for the alignment parameters. Doing so, it is easy to share the alignment across different versions of the same tilt-series, for example a binned or a denoised version.

The gold bead tracking gives good results for a large variation of different doses and gold beads, it is still worthwhile to verify the result in aiClicker.

5. CTF parameter estimation: aiCtfDetector or aiCtfDetectorGui

In order to reconstruct with CTF correction, the defocus and astigmatism parameters must be determined for each individual tilt. ARTIATOMI does three dimensional CTF-correction and thus the CTF parameters are not to be given only for 2D tilt-image of the tilt series but as a defocus value at a specific height in the 3D reconstruction volume. Per definition, CTF parameters are given for Z=0 in the global coordinate system.

In case that CTF is determined on the entire frame, no alignment information is used to estimate the deviation in Z as usually the difference is low especially for lower tilts. On the other hand, CTF can also be detected on an image patch, for example a carbon strip that has a well-defined position in the volume. In this case, the user only has to provide the position of that strip as a 3D coordinate, and both CTF detectors, GUI and CLI, track the point given the marker alignment and compensate for the Z height difference.

CTF is measured in the following way: The lowest tilt angle of the tilt series is used to determine astigmatism – the amount of astigmatism and its orientation. The angle convention is the same as used by GCTF. Once the defocus and astigmatism are determined for the minimum tilt image, the entire tilt series is scanned, but only for defocus assuming that astigmatism did not change from tilt to tilt.

Note that this is only the first step of CTF determination. If sub-tomogram averaging and sub-tomogram shift refinement is performed later on, the defocus and astigmatism values for each individual tilt can be further refined.

To estimate the defocus values used for recording, aiCtfDetector computes the power spectrum of the image and band pass filters it given the user provided parameters. Then the background is subtracted, using either the LP-approach [https://doi.org/10.1016/j.ultramic.2020.112950] or our quadrature-based method based on the quadrature operator introduced in [https://doi.org/10.1016/j.jsb.2012.12.006]. Then for a user defined range, the power spectrum is correlated with a modelled CTF for each value defined in that range. The defocus (and astigmatism) with the largest correlation peak is stored.

6. Tomographic reconstruction – full volume: aiRec

The tilt series now well aligned and with CTF properly determined can now be reconstructed. The reconstruction uses two configuration files: one like all the other application to define algorithmic parameters and another one that provides tomogram dependent settings, such as additional shifts, additional rotations or volume dimensions. Note that the dimensions are to be given in un-binned pixels. Doing so, one can define one common configuration for all tomograms in a dataset and include variations for each individual tomogram.

To position the reconstruction in the volume, one has multiple options: If a Z-shift is needed, it is best to incorporate this shift in the marker model and realign the tilt-series incorporating this shift. Otherwise, the volume can be shifted in X, Y and Z by any amount. This shift is taken into account later on when reconstructing or projecting sub-volumes. Further the reconstruction can also be rotated around the tilt-axis (adding a rotation to the tilt-angle) and around the X-axis to compensate for oblique orientations of the real sample. Note that the rotations effectively do not rotate the volume as the volume always must remain axis aligned. Instead, the projection geometry is rotated accordingly.

The user can choose in-between two fundamentally different reconstruction approaches: SART or weighted back-projection. Both support the same three-dimensional CTF correction and an alignment with a deformation model. The only restriction for SART is, that it is not possible to reconstruct with a deformation model when a volume is split over multiple GPUs. This is only possible using a single GPU as the forward projection would otherwise have to read from a memory portion allocated on another GPU while following the deformation.

7. Particle / sub-volume selection: aiTemplateMatching or manual selection

aiTemplateMatching correlates a small reference volume with the full reconstruction and performs an angular scan over a user-provided angular range. For each voxel in the full reconstruction the best correlation coefficient and the corresponding rotation angles are stored. Finally, the volume is filtered to eliminate local maxima nearby another maxima where the one with larger value is maintained. Maxima below a given threshold are removed. From this distance filtered three-dimensional CC-map a motivelist extracted. Finally, the user can provide a maximum number of particles to create, so that only the best N entries are stored in the motivelist. (A motivelist is a binary or JSON file that contains the position and orientation of individual sub-volumes of one or more tomograms.)

If no start reference is given, the user can manually select positions in the tomogram using aiClicker.

8. Reconstruction sub-volumes: aiRecSubvol

Whereas the goal in standard reconstruction is to obtain an evenly distributed power spectrum without any artificial amplifications and isotropic effects, the best way to reconstruct sub-volumes is to omit any processing prior to reconstruction, which means using weighted back-projection without weighting. Each individual particle in turn now has an over-weighted low-frequency range due to the overlapping tilts in the power-spectrum. But this will be compensated by the wedge-filter during sub-tomogram averaging – whereas signal that would be dampened too much due to a weighting filter is lost forever. It is not possible to use the SART algorithm as no complete forward projection can be computed only from sub-volumes. Hence only the WBP reconstruction is possible and we encourage the use of no filtering.

aiRecSubvol can operate in two different modes: particle or reference. In particle mode, each particle is reconstructed separately and stored in a file for each particle. The wedge information is stored in a metadata file for the first particle only as it is the same for all. In reference mode, for each reference in the motive list, the particles are summed up internally and an additional wedge file is created that contains the wedge normalization information. The two volumes are stored separated, as usually the references from multiple tomograms must be added prior to normalization. One thus adds the reference and wedge volumes for all tomograms and performs the wedge normalization with the summed-up volumes, both using aiToolbox.

9. Sub-tomogram averaging: aiAveraging

At this stage one has a motive list with the estimated positions of particles in the three-dimensional global coordinate system. In case of some apriori pre-alignment or after template matching also a rough orientation is known. aiAveraging is used to precisely determine the orientation and a local shift for each particle so that it is aligned to a common reference. The first reference can be a feature-less gauss blob or the random sum of all unaligned particles, in next iterations, the result of the previous iteration is used as reference.

The motive list stores the position, the shift and the orientation of each particle. The combination of shift and rotation is called a displacement and each particle can have a pre-reconstruction displacement and a post-reconstruction displacement. This means that a reconstructed sub-volume integrates the displacement given by the pre-reconstruction displacement. If an additional displacement is needed to perfectly fit to the reference, it is given in the post-reconstruction displacement values.

aiAveraging can handle multiple references at the same time, currently with the only restriction that the volume dimensions must be the same for all references. When performing a standard averaging with one common reference – as usually done in the first iterations together with a strong low-pass filter – a FSC after each iteration is calculated with a threshold of 0.5. When doing so-called gold-standard averaging, two separated reference volumes are used and particles that are associated with one reference never see the other reference. In this case, the resolution is estimated at the 0.143 threshold of the FSC.

Gold standard averaging thus needs two independent references and in order to encode multiple references, ARTIATOMI uses even-odd pairs: a gold-standard reference pair is given by two consecutive even-odd reference numbers like 0-1, 2-3, etc.

In case of multiple references in one tomogram, e.g. a protein of interest and some “supporting proteins” to allow for better shift and CTF-refinement, it is advised to only use even numbered reference numbers even if standard averaging is performed so that one can switch to gold-standard later on without the need to renumber any reference. (In this case the standard averaging must be performed with a low-pass filter or even better on binned sub-volumes so that high-resolution features are not shared across the gold-standard references.)

aiAveraging supports phase and cross-correlation to determine the particle orientations, the missing wedge is compensated at all times, either by a user-provided wedge-volume or a wedge shape is computed internally that also includes dose-weighting (including also missing frames from movie-stack alignation and any other weighting applied during reconstruction). To focus the averaging on a specific part of the reference, a mask is used and the mask can have any shape and should have soft edges. The area of interest has value 1 and the suppressed area has value 0, all mask values must be in the range [0..1].

10. 3D-classification: aiPCAClassification

As the name suggests, aiPCAClassification performs classification of sub-tomograms using PCA. The comparison is done in 3D on a masked area of the entire particle and compensates for the missing wedge so that the orientation of a sub-volume does not influence the classification results.

11. Refinement of the 2D shift in projection space: aiShiftRefine

aiShiftRefine estimates the individual shift for each particle/sub-volume so that individual particles are shift consistent. This means that the center of each particle in each projection is projected into the center of the sub-volume. Due to tilt-series alignment inaccuracies, this is not always the case.

aiShiftRefine offers two possible approaches: the correlation of the projected reference with the original micrograph or with the projection of the individual particle. As the signal of one single particle is not sufficient to overcome noise, we need to average the result of neighboring sub-volumes. In both cases, we can either average the CC-MAP or project multiple references (and if sub-volumes are used, sub-volumes) to a larger image that is then used for correlation.

Projecting only the sub-volumes effectively removes the entire background of the tomogram and allows to perform shift refinement even on thicker samples as only the particle of interest is correlated. One could also understand sub-volume-based shift refinement as sub-tomogram averaging where the shift is not determined globally for all tilts together but individually for each tilt.

12. CTF parameter refinement: aiCTFRefine

After shift refinement, one knows precisely the position of every sub-volume on each micrograph. This can be used to refine also the CTF parameters for each individual tilt in a tilt series. Remember that defocus is always defined as the defocus at Z=0 in the global coordinate system. Positioning the result of sub-tomogram averaging for each reference at their estimated point in 3D space of the corresponding sub-volume and projecting them in the direction of a tilt in the tilt series, one obtains a precise representation that ideally is identical to the original projection except all the surrounding elements in the tomogram outside of the sub-volumes.

If we now do the projection including the CTF-based phase-flipping, we obtain rings in the power spectrum with opposite phases, depending on the defocus used for projection. In order to precisely determine the correct CTF parameters, i.e. defocus and astigmatism for each tilt in the tilt series, we only have to project with varying parameters and scan for the best combination. As the CTF introduces a defocus dependent phase flip, phase correlation results in a strong peak once we get to the correct parameters and we simply store the best parameter for each tilt.

13. aiDoseBoosterTraining and aiDoseBooster

aiDoseBooster comes in two separated tools: aiDoseBooster and aiDoseBoosterTraining.

When recording tilt-series (or individual single particle-like frames) in dose-fractionation movie stacks for each tilt, a noise model can be trained based on the individual movie frames using the Noise2Noise-principle [https://proceedings.mlr.press/v80/lehtinen18a.html]. This noise model is encoded using a neural network in the shape of UNET and its parameters are trained using aiDoseBoosterTraining. As input, aiDoseBoosterTraining takes a list of movie-stacks (and potentially a gain reference) and the network is trained in multiple iterations. To avoid that the resulting network is biased for a specific defocus, a wide variety of defoci should be used for the training dataset. As a rule of thumb: As much images should be used for training as possible, the only limiting factor is likely the time that it takes to train the network. We experienced best results, when training the network specifically for each dataset to denoise.

Once the network is trained based on the entire dataset (or a subset), each individual tilt-series can then be denoised using aiDoseBooster. The result is then a second tilt-series using the same metadata as the original tilt-series but with lower noise levels. We also experienced that thon-rings might appear clearer in the denoised images, but this can also be artifacts due to training bias. It is hence better not to perform CTF detection / defocus determination in the denoised images. Finally, we only used so far the denoised images purely for visualization purposes. For processing we usually rely on the original tilt-series.

14. The swiss army knife of Artiatomi: aiToolbox

aiToolbox is a small command line based tool that provides a simple interface to most common tasks that occur in processing in cryo-electron tomography. aiToolbox allows to quickly set values in a motivelist, crop a volume, create a mask, compute the FSC of an average, import data and much more.