CTFAPPLY

[sthennin]

We should replace CTFAPPLY by a proper CTF inversion routine, including Wiener filtering with correct weighting of frequency-dependent SNR, see RELION papers. Also, this should be done BEFORE and not AFTER unbending, to avoid truncation of the ripples of the PSFs. The higher we get in resolution (K2 images), the more important it becomes that CTF correction happens before unbending.

[sthennin]

This is implemented now in 2dx_ctfapply2k.for. There are now two options, to use CTF correction, or Wiener filtration. If choosing the latter, then the calculated CTF values are actually the inverse of the Wiener filter correction factor. These new "Wiener-CTF" values are sometimes quite large, ranging between zero and 100.0. So, we have to see if the other programs like CENTRIC and AVRAMPHS can correctly deal with this.

[sthennin]

If CTFAPPLY calculates a "Wiener-CTF", by which the AMP should be multiplied to get best weighting of signal for unbending and also on the same time for reconstruction, how would then the FOM have to be calculated? Should FOM now all be 100.0, since the AMP correction already took care of that? And how does lattice line fitting react to this?

[nikhilbiyani]

I was looking at the code.. Wiener filter is first calculated and then it's reciprocal is taken to cope up with the data flow. Buy if the reciprocal of the wiener filter is calculated at first place one could avoid few nans (/0) !

[sthennin]

Thinking about this, I don't think this is the way to go for electron crystallography. Lattice line fitting needs the amplitudes as precise as possible (without any down-weighting by SNR), plus some FOM, and then fits a least-square trace through those measurements.

However, single particle processing would need the Wiener filter approach.

Unbending requires the best SNR, which we get by multiplying the CTF-modulated image with the CTF, and the correlate with the CTF-free reference. This is what we do currently for non-tilted images (CTF correction option active), but not for highly tilted images (TTF correction option active).

[sthennin]

... extension to the above: Lattice line fitting needs the AMPlitudes as precise as possible, plus some PHAses and FOM values, and then fits a least-square trace through those measurements. For this, the AMP values should be CTF corrected. That is in the MRC software conventionally done by dividing by the CTF, which increases the AMP value. Conventionally, then there is a cutoff of 3 applied to the amplification by CTF correction, i.e., the smallest CTF value that is decided by is 0.333.

Here, we can do this better, by doing correct Wiener filtration. Instead of dividing by CTF, this should then be:

CTF / ( CTF**2 + 1 / SNR)

whereby SNR is AMP2/BCK2. This is then equal to

( CTF / BCK2 ) / ( CTF2 / BCK2 + 1 / AMP2 )

This is now implemented in line 400 of 2dx_ctfapplyk.for. BTW, the implementation checks for division by zero correctly, so it should not produce any NANs.

This should result in smaller error bars for the AMP values in the lattice line fitting.

[scherers]

Well, appreciate this method development discussion here.

It think we should record a test dataset of BR to benchmark all improvements of our software because:

• We already have large, well ordered 2D crystals of BR in the fridge and they diffracted perfectly for ED.
• Using a state of the art test dataset allows us to quantify improvements in an unbiased way and simplifies teaching.
• The known X-ray structure of BR allows trivial judgment of map quality at all resolution levels.

Thus I think it's worth to record and process 200-300 BR images to use them as sustainable benchmark for 2dx.

[sthennin]

The Wiener filter is Equation 2.45 in Penczek, P.A., Image Restoration in Cryo-Electron Microscopy, Methods in Enzymology, Volume 482, (2010) DOI: 10.1016/S0076-6879(10)82002-6 It still includes the envelope due to the B-factor weighting, which then should be frame-specific. It then becomes

(CTF * E) / ( (CTFE)*2 + 1 / SNR)

or, using SNR = AMP2 / BCK2,

( CTF * E / BCK2 ) / ( CTF2 * E2 / BCK2 + 1 / AMP**2 )

whereby E is E(v), as:

E(v) = exp ( -B / 4 * v**2 )

with v being the resolution in reciprocal Angstroems, and B being the B-factor (in A-2).

Sjors Scheres expands this E(v) in his eLife (2014) paper to

E(v) = exp ( -B / 4 * v**2 + C )

whereby C is a frequency-independent component of this envelope. B and C are frame-specific. Sjors normalizes over the sum of all E(v) for all frames, so that in the sum the different C values don't matter in their total amount. Only their relative ratios matter. My feeling is that C is a way to reduce or increase (C is mostly negative) the impact of the frequency-dependent B.

[sthennin]

Unbending before CTF correction means that we destroy the high-resolution details of the PSF. We need to do CTF correction before unbending.

Proposition: We have 40 frames of a tilted 2D crystal, called Fn. The drift-corrected average is called AVG.

We should process AVG the conventional way, with CTF correction and unbending. Thereby measure defocus and tilt precisely.

Apply CTF phase flipping to all frames Fn, in a tile-wise manner: Working on small tiles, cut 200x200 pixel tiles from Frame Fn, CTF-phase-flip this in Fourier space, crop in real space the central 100x100 pixels from each corrected tile, and write those back into the output image at the correct tile position. Using overlapping 200x200 tiles and writing back only 100x100 tiles will make sure we don't get edge effects between the titles.

These CTF-phase-flipped frames now have no delocalized point spread function, so that we can unbend such an image without cutting the PSF into pieces.

Now calculate a series of running averages over odd-numbered running sub-sets of frames, e.g. using 7 frames. This results in 40 running-average frames Fra.n Unbend each frame Fn, as developed by Sebastian, whereby Fra.n should be used to determine the unbending profile, which is then used to unbend Fn.

Evaluate AMP, PHS, BCK from each FFT of the unbent Fn. The PHS values are already phase-flipped.

Run CTFAPPLY, but ignore the sign of the CTF, so that the phase flipping isn't done a second time. This should write out AMP, the (not) corrected PHS, BCK, and CTF. This CTF value should be the result of the Wiener filter calculation above. However, this is easy for an image of a non-tilted crystal, as the defocus is constant. But for an image of a titled crystal, we need to calculate the summed CTF shape for each spot position: If on the tilt axis, the CTF oscillations are as for non-tilted samples, i.e., between 0 and 1. But for spots perpendicular to the tilt axis, the oscillations of the CTF may already have been reduced to 0.5, depending on the resolution of that spot.

[sthennin]

Now implemented.