Bayesian Polishing questions.

hgxy15 commented 3 years ago

Dear developers of RELION,

I'm going through the bayesian polishing related stuff in the src directory and came across a question.

In the file motioncorr_runner.cpp, the global motion seems to be calculated in terms of original movie pixel size, the local motion model however, seems to be calculated in terms of binned micrograph pixel size as the local_xshifts, local_yshifts are calculated from Ipatches which is already binned, unlike the global motion which is multiplied by the prebinning factor, there seems to be no prebinning factor correction for local motion. The function I'm looking at is the following one:

    bool MotioncorrRunner::executeOwnMotionCorrection(Micrograph &mic);

However, in the following getShiftAt function from micrograph_model.cpp:

    int Micrograph::getShiftAt(RFLOAT frame, RFLOAT x, RFLOAT y, RFLOAT &shiftx, RFLOAT &shifty, bool use_local, bool normalise) const;

When loading inital tracks for bayesian polishing, both local and global motion seems to be calculated in terms of binned micrograph pixel size. In my understanding, this would lead to decreased precision of initial track estimation when early binning is used. Is there anything I'm missing here or is this something we need to consider when using motion model star files generated by relion's own motioncor2 implementation? (e.g. Should I manually multiply all coefficients by the binning factor I used ?)

Many thanks in advance!

Best, Gaoxing

biochem-fan commented 3 years ago

I don't understand your question. The shifts are in floating point numbers, not integers. So working in the original pixel size or the binned pixel size does not really matter, as long as they are consistent.

hgxy15 commented 3 years ago

from my understanding, the bayesian polishing workflow works as follows: First, the motioncor from relion's own implementation calculates global and local residual shifts and store these information in the motion star file of each micrographs. After 2D, 3D classifications and 3D refinement, the motion_refine program extracts the initial track information from motion star files using global shift stored in the global shift section and local motion model coefficients. For this to achieve the best outcome, all units used to describe the motion through out the whole process should be consistent. My question is that after going through the codes, it seems that the global shifts stored in the motion star files are in units of movie pixel size. When extracting initial tracks, the motion_refine program also treats all extracted information as in units of unbinned moviel pixels. However, the global motion for each frame is corrected by the prescaling factor in motioncorr_runner and stored in the final motion star file, the local residual motions are not, which means that the local residual motions are stored as binned pixels (reflected in the scale of the coefficients) but interpreted as unbinned movie pixels in the motion_refine program, if this is the case, maybe this would lead to some inaccuracies in the initial tracks?

hgxy15 commented 3 years ago

This is merely something that's been troubling me for some time, its most likely that I have missed something in the code. Just wanted to confirm this with you for my own peace of mind :)

biochem-fan commented 3 years ago

Are you aware that data_local_shift table from MotionCorr is never read by Polish but data_local_motion_model is used?

hgxy15 commented 3 years ago

Yes, which is why I specifically checked the coefficient estimation part of the motioncorr_runner, I think the vecX and vecY fed into the final svd are in units of binned pixel, not movie pixel, so the numbers restored by these coefficients should also be in units of binned micrograph pixels, but I think they are now being interpreted as in units of original movie pixels.

biochem-fan commented 3 years ago

At the moment I don't have time to investigate this carefully so below is my random guess.

Patch shifts are normalized to the image size, e.g. const RFLOAT x = patch_xs[i] / nx - 0.5; in MotionCorr before fitting. Doesn't this solve the issue?

hgxy15 commented 3 years ago

The final SVD takes three arguments. matA * coeffX = vecX; I thinks the code you mentioned takes the unit out of matA and makes matA essentially a matrix with no unit. coeffX and vecX then should have the same unit, the unit in vecX is still binned micrograph pixel, so coeffX should also be in units of binned micrograph pixel.

biochem-fan commented 3 years ago

Aha, OK, this indeed sounds suspicious. I will investigate this later.

Meanwhile, please don't worry too much. Bayesian Polishing is not very sensitive to initial values. When you use UCSF MotionCors, Polish sees only global trajectories but it works just fine.

hgxy15 commented 3 years ago

Thanks a lot! looking forward to your reply :)

hgxy15 commented 3 years ago

Also on a related topic, I noticed that when I throw away the first 2 frames in motioncorr and specify first_frame to be 3 and last frame to be default in motion_refine, the output tracks.star file only contains motion information for 28 frames (32 frames collected in total), when I specify first frame as 1, the last frame using default (which should be the last frame collected), the output tracks.star file contains 32 frames, if I specify first frame to be 3 and last frame to be 32, the output tracks.star file contains 30 frames as I want it to be (since 2 frames are thrown away at the time of motion correction).

In short, if you have thrown away frames at the time of motion correction, in order to obtain correct number of output frames in bayesian polishing, one should specify both the --first_frame and --last_frame option and leave none to default. I'm leaving this comment here just in case someone else runs into the same problem, although I think it would be more reassuring if the next release would allow us to leave the last frame or first frame to default when throwing away frames at the start or end of a movie stack.

biochem-fan commented 3 years ago

Brief notes on my investigation:

Indeed SVD is called in binned pixels when early_binning is ON (i.e. the default).
But if we change that, we have to fix realSpaceInterpolation as well.
Polish expects Micrograph::getShiftAt to return in the original unbinned pixels (micrograph_handler.cpp).
Micrograph::getShiftAt says (shiftx, shifty) are always UNBINNED pixels in the original movie.

Plan:

Fix SVD and realSpaceInterpolation.

P.S. I won't fix --first_frame and --last_frame issues.

3dem / relion

Bayesian Polishing questions. #771