Skewed angle distributions

[daniel-s-d-larsson]

Several times I've noticed that the distribution of angles during 3D classification and 3D refinement jobs looked suspicious. It is especially noticeable when applying symmetries, with lots of particle ending up with angles on the edge of the asymmetric unit and under sampling on the "other side". But it is also noticeable for asymmetric (C1) reconstructions as "bolding" patches. It seems to affect both the CPU and the GPU versions of the code, but the result is more noticeable in my GPU runs.

My sample is approximately round, so there should not be any reason for orientational biasing from the grid. Attached is an example where I used icosahedral symmetry (I4) to perform a 3D refinement of an isometric virus using relion_refine_mpi both with and without --gpu. The GPU job ran on a single node with 12*2 hyper threaded cores and 4 Titan Black GPUs. The CPU job ran on 8 nodes. The relion version used was 2.0.4 (commit 76f90ea in the public repo, "repaired critical bug auto-refine non-GPU" by Sjors Scheres).

Submit commands:

mpirun -np $((1*6+1)) --oversubscribe $(which relion_refine_mpi) --o Refine3D/job071/run --auto_refine --split_random_halves --i Class3D/job070/run_ct27_it050_data.star --ref Refine3D/job037/run_ct21_class001_4.11A_130px.mrc --ini_high 15 --dont_combine_weights_via_disc --pool 20 --ctf --ctf_corrected_ref --particle_diameter 380 --flatten_solvent --oversampling 1 --healpix_order 2 --auto_local_healpix_order 4 --offset_range 2 --offset_step 1 --sym I4 --low_resol_join_halves 40 --norm --scale --j 4 --gpu ""

mpirun -np $((8*6+1)) --oversubscribe $(which relion_refine_mpi) --o Refine3D/job071b/run --auto_refine --split_random_halves --i Class3D/job070/run_ct27_it050_data.star --ref Refine3D/job037/run_ct21_class001_4.11A_130px.mrc --ini_high 15 --dont_combine_weights_via_disc --pool 20 --ctf --ctf_corrected_ref --particle_diameter 380 --flatten_solvent --oversampling 1 --healpix_order 2 --auto_local_healpix_order 4 --offset_range 2 --offset_step 1 --sym I4 --low_resol_join_halves 40 --norm --scale --j 2

Angles from the GPU run:

Angles from the CPU run:

[daniel-s-d-larsson]

Here are the angle distributions plotted for each iteration as animated GIFs:

GPU run:

CPU run:

[bforsbe]

Thanks for the detailed report!

A couple of follow-up Qs;

1. Do both runs reach the same resolution?
2. Do both runs reach the same sampling density?
3. What happens if you run the GPU run with "--healpix_order 4 --auto_local_healpix_order 5" ?
4. What happens if you run any of the runs with another icosahedral symmetry, like I2?

Thanks!

[daniel-s-d-larsson]

1. Do both runs reach the same resolution?
2. Do both runs reach the same sampling density?

These are the output from the last step of each job:

GPU job:

CurrentResolution= 8.34844 Angstroms, which requires orientationSampling of at least 2.51748 degrees for a particle of diameter 380 Angstroms
Oversampling= 0 NrHiddenVariableSamplingPoints= 1305
OrientationalSampling= 0.46875 NrOrientations= 145
TranslationalSampling= 0.228 NrTranslations= 9
=============================
Oversampling= 1 NrHiddenVariableSamplingPoints= 41760
OrientationalSampling= 0.234375 NrOrientations= 1160
TranslationalSampling= 0.114 NrTranslations= 36
=============================
Expectation iteration 18
13.35/13.35 min ............................................................~~(,_,">
Averaging half-reconstructions up to 40 Angstrom resolution to prevent diverging orientations ...
Note that only for higher resolutions the FSC-values are according to the gold-standard!
Calculating gold-standard FSC ...
Maximization ...
9.65/9.65 min ............................................................~~(,_,">
Auto-refine: Refinement has converged, stopping now... 
Auto-refine: + Final reconstruction from all particles is saved as: Refine3D/job071/run_class001.mrc
Auto-refine: + Final model parameters are stored in: Refine3D/job071/run_model.star
Auto-refine: + Final data parameters are stored in: Refine3D/job071/run_data.star
Auto-refine: + Final resolution (without masking) is: 8.34844
Auto-refine: + But you may want to run relion_postprocess to mask the unfil.mrc maps and calculate a higher resolution FSC

CPU job:

CurrentResolution= 8.34844 Angstroms, which requires orientationSampling of at least 2.51748 degrees for a particle of diameter 380 Angstroms
 Oversampling= 0 NrHiddenVariableSamplingPoints= 1305
 OrientationalSampling= 0.46875 NrOrientations= 145
 TranslationalSampling= 0.225 NrTranslations= 9
=============================
 Oversampling= 1 NrHiddenVariableSamplingPoints= 41760
 OrientationalSampling= 0.234375 NrOrientations= 1160
 TranslationalSampling= 0.1125 NrTranslations= 36
=============================
 Expectation iteration 19
10.12/10.12 min ............................................................~~(,_,">
 Averaging half-reconstructions up to 40 Angstrom resolution to prevent diverging orientations ...
 Note that only for higher resolutions the FSC-values are according to the gold-standard!
 Calculating gold-standard FSC ...
 Maximization ...
15.72/15.72 min ............................................................~~(,_,">
 Auto-refine: Refinement has converged, stopping now... 
 Auto-refine: + Final reconstruction from all particles is saved as: Refine3D/job071b/run_class001.mrc
 Auto-refine: + Final model parameters are stored in: Refine3D/job071b/run_model.star
 Auto-refine: + Final data parameters are stored in: Refine3D/job071b/run_data.star
 Auto-refine: + Final resolution (without masking) is: 8.34844
 Auto-refine: + But you may want to run relion_postprocess to mask the unfil.mrc maps and calculate a higher resolution FSC

[bforsbe]

Looks identical from a sampling point-of-view. There's something strange about the sampling in the CPU-run though, judging from the GIFs. The sampling points at iterations 2-6 (or thereabouts) are all extended along a line... curious.

You said those were angular distributions, but are they darker for higher intensities, or are they plotted points with coordinates given bu those in the _data.star? I'm guessing the latter?

[daniel-s-d-larsson]

You said those were angular distributions, but are they darker for higher intensities, or are they plotted points with coordinates given bu those in the _data.star? I'm guessing the latter?

Yes, the points are plotted on top of each other.

[daniel-s-d-larsson]

3. What happens if you run the GPU run with "--healpix_order 4 --auto_local_healpix_order 5" ?

With these settings the angles look better. The reconstruction converged to the same resolution as the other two runs. The areas without sampling are smaller and at different locations. But they are still there and with corresponding clumping up of particles on the other side of the asymmetric-unit boundary. I guess I could try to tweak these two parameters further, but I think also there could be problems in the code. Will the increased initial sampling affect convergence for difficult to align data-sets?

 CurrentResolution= 8.34844 Angstroms, which requires orientationSampling of at least 2.51748 degrees for a particle of diameter 380 Angstroms
 Oversampling= 0 NrHiddenVariableSamplingPoints= 1260
 OrientationalSampling= 0.46875 NrOrientations= 140
 TranslationalSampling= 0.225 NrTranslations= 9
=============================
 Oversampling= 1 NrHiddenVariableSamplingPoints= 40320
 OrientationalSampling= 0.234375 NrOrientations= 1120
 TranslationalSampling= 0.1125 NrTranslations= 36
=============================
 Expectation iteration 12
13.65/13.65 min ............................................................~~(,_,">
 Averaging half-reconstructions up to 40 Angstrom resolution to prevent diverging orientations ...
 Note that only for higher resolutions the FSC-values are according to the gold-standard!
 Calculating gold-standard FSC ...
 Maximization ...
2.72/2.72 min ............................................................~~(,_,">
 Auto-refine: Refinement has converged, stopping now... 
 Auto-refine: + Final reconstruction from all particles is saved as: Refine3D/job071c/run_class001.mrc
 Auto-refine: + Final model parameters are stored in: Refine3D/job071c/run_model.star
 Auto-refine: + Final data parameters are stored in: Refine3D/job071c/run_data.star
 Auto-refine: + Final resolution (without masking) is: 8.34844
 Auto-refine: + But you may want to run relion_postprocess to mask the unfil.mrc maps and calculate a higher resolution FSC

[bforsbe]

Your issue resulted in a very interesting discussion, and I believe we can offer some insight.

Relion uses the initial --healpix_order to determine the most basic tessellation of the Tilt/Rot-space, according to the healPix library (see this article, and especially Table 1 for details). Any healPix pixel-centers which lie outside the Tilt/Rot region corresponding to the specified symmetry (in this case icosahedral), gets thrown away. As you go higher in resolution and sampling, Relion over-samples the Tilt/Rot region by increasing the healpix order. Although from your plots it seems as though the prior used in Relion is only ever assigned to the most basic, original sampling.

This means that when you use a very low --healpix_order compared to the symmetry, the prior is trying to cover a very small Tilt/Rot region with very big healpix Pixels. In your original run, the large pieces of excluded space thus correspond to the original healPix pixels. Only when Relion starts doing local searches do some particles find it possible to migrate to those regions.

In your most recent run, your initial sampling is much higher as a result of --healpix_order 4, meaning your starting healPix pixels are smaller. Thus, only very small regions around the border (about 4 regions by the looks of it) get excluded, and they are much smaller as well.

Is a higher --healpix_order really a "solution"

As you use higher symmetries, you can afford (computationally) to go up in healpix sampling. However Relion should perhaps assure that your prior always completely covers the subspace determined by the symmetry at the sampling used for the present iteration. We will direct some attention to this, since it is in principle not very complicated to fix (famous last words).

Sidenote: The only reason this arose at all is because healPix is not compatible with arbitrary rotational symmetries; that was not something they kept in mind when developing it, because they had astronomy in mind, and the night sky does not often display e.g. icosahedral symmetry.

Thank you for a very illuminating issue report!

P.S.

I'm still not sure why you get higher population on the border, but it seems to be much reduced in the latest run. Will think more about that, so don't close the issue just yet!

[bforsbe]

I should perhaps also clarify and emphasize that this will only ever be an effect when you use symmetries, and more of an effect he higher the symmetry. Higher-symmetry usually means "easier to align" (in some sense), so while a detrimental effect, I don't think it's going to mean night and day differences. Your resolution for instance, didn't suffer to much, but I think you could judge best from your densities whether there was any significant improvement.

[daniel-s-d-larsson]

In terms of compute time, increased initial sampling almost cut the number of iterations in half. Wall-time the original GPU run took 3 hrs 40 min and the increased sampling GPU run took 2 hrs 40 min.

[bforsbe]

Convergence doing its thing! Nice result. Any noticeable differences in the densities?

[daniel-s-d-larsson]

Any noticeable differences in the densities?

No. The two maps looks very similar. A difference map has maximum values of 0.024.

[daniel-s-d-larsson]

To illustrate the bug, I've rotated the distribution plot around the two-fold and three-fold axes. The aggregation of particles along the rim of the asymmetric unit perfectly line up with the white patches. So there seems to be something that prevents particles to optimize across the asymmetric boundary, and hence particles that by chance end up on the wrong side will get stuck.

[bforsbe]

That confirms or at least strongly indicates that the increased density at certain points along the symmetry unit edges is due to the same excluded-area effect. I'll try to get this patched during next week. Thanks!

[daniel-s-d-larsson]

Here is another interesting/pathogenic angle plot. The particle angles bunch up in discrete patches on a lattice. Increasing the initial angular sampling did not help in this case. Rather it caused some of the "spots" to completely disappear.

This is for the C1 reconstruction of the same data set as above. Despite the fact that the sample has strong inherent symmetry, I still don't think this is sample-derived.

The above projection does not show it very clearly, but the sampling at the "north" and "south" poles is also less (i.e. "ozone holes"). Perhaps that is related to the bug with the boundaries of symmeterized reconstructions.

[daniel-s-d-larsson]

Here is a follow-up on the issue with the auto-refine C1 reconstruction in my last comment with an additional GIF showing a detailed view of the angles:

Parameters used in this job were: --oversampling 1 --healpix_order 2 --auto_local_healpix_order 4 (first four iterations ran with 4 GPUs, then restarted without GPU support)

Iteration 1-4 the angular sampling was 15° (7.5° with oversampling). Iteration 5-9 the angular sampling was 7.5° (3.75° with oversampling). Iteration 10-12 the angular sampling was 3.75° (1.875° with oversampling). Iteration 13-14 the angular sampling was 1.875° (0.9375° with oversampling). Iteration 15-16 the angular sampling was 0.9375° (0.46875° with oversampling). Iteration 17-20 the angular sampling was 0.46875° (0.234375 with oversampling)

Observations (some already mentioned above):

1. Higher healpix levels do not seem to sample the sphere evenly. In particular from iteration 10 and onwards, "allowed" islands are separated by non-sampled regions.
2. The lattice seems to drift between consecutive iterations at lower healpix levels, even when the same sampling rate remains (iteration 1-9). Is that a property of healpix or a bug?
3. During the first phase (iteration 1-9) the pixels are "smeared" out into streaks. This might be related to the drift of the lattice. As mentioned in another post above, this seems more pronounced in CPU runs, although the first four iterations in this run were with GPU support. Another thought is that it is related to misalignment between the two independent reconstructions, since I plot particles from both the models on top of each other.
4. Smearing was not observed in another run with the same parameters with GPU support throughout the entire run. See this animated GIF:

I am now trying with --auto_local_healpix_order set to 5 (job not converged yet). This gives a different spacing of the "islands" and the "forbidden zones":

EDIT: Updated the second animated GIF to show the same field of view as the other two GIFs.

[bforsbe]

These are results of a C1-run on something with high symmetry, right?

P.S. It much easier to refer to your observations if you make it a numbered list rather than bullet:ed.

[daniel-s-d-larsson]

Yes, this is an attempt of an asymmetric reconstruction of a virus particle where the capsid has icosahedral symmetry. (Changed the list into a numbered one.)

[daniel-s-d-larsson]

I can also report back on a test I did for the issue with the icosahedral averaged reconstruction above (job071c).

To test whether I could heal the distribution, I wrote a script that moved all particles assigned to pixels with a large population to the "other side" of the asymmetric unit boundary. (Achieved by flipping the sign of the rotAngle parameter in the _data.star file). I also tried to move it a bit further in on the other side (see figure below). Neither attempt worked. Instead I got back approximately the same distribution as I started with.

Figure: Heat map of angle distributions (in jet color scale using matplotlib's hexbin). (left) Distribution of iteration 11 from job071c. (middle) Moved bins with many particles. (right) Final, converged distribution in iteration 13.

Very puzzling. Just to verify, I looked at the particles in large classes to see if they looked strange (e.g. very low contrast, bad ice etc), but they looked perfectly fine to me and not different from other particles.

[bforsbe]

The "drift" you mention (observation 2) is intended measure which RELION does actively to avoid issues arising from the finite sampling grid and lingering bias. It's a perturbation which you can turn off by using --perturb 0. You should only turn this off for testing, any results you get with --perturb 0 might have bias or be self-reinforced ("over-fitting").

Also, I'd be very cautious to not jump to conclusions about results when you stretch the scope of the algorithm. When you run C1 on something highly symmetrical you are posing a very strange question for relion to answer. Relion uses orientation-distributions as statistics for prior belief in later iterations. It does this by mapping the assigned orientation of each image to the sampling grid. Now imagine you align a symmetrical object to a symmetrical reference, and imagine the fit to be a heatmap on a sphere. That heatmap is going to be symmetrical, meaning that you'll find (in your case) 60 points which are "the best". But since we sample finitely, there will be some grid-points which are closer to one of those "best" points than others. In just the same way, you are mapping that "best assignment" to a specific point on the sampling grid for the next iterations prior belief.

What I am trying to describe is that you are implicitly asking relion to take the icosahedral asymmetric unit and split it up over the entire sphere. Your are expecting that this should smear it all over, but due to finite sampling and the fact that there is a 60-fold redundancy in the fitting landscape, relion is rather mapping that unit to isolated regions. This is not a problem, unless you overlap this distribution in all icosahedral symmetry mapping and still find regions which are not sampled. My gut feeling is that you won't.

[daniel-s-d-larsson]

This is indeed a tricky problem. The symmetry will definitely interfere with the orientation determination. My initial model is actually asymmetric (produced by doing a symmetry expansion and 3D classification, as you suggested to me before). When I tried to use a symmetric initial structure, the reconstruction converged with all angles ending up on a ring or at two opposing "poles".

[bforsbe]

I'm not sure there is a "problem", at least not a major one. It might be a problem for what you are specifically trying to do, if you are trying to detect small deviations away from something symmetrical, but that's rather something relion has not been designed to do. Your initial concerns regarded whether or not the angle-distributions were legitimately ok, and apart from when using fairly coarse sampling at high symmetry, they are. At least, I have not seen anything to indicate anything else.