relion_particle_subtract resets rlnGroupNumber in (optional) --data star file but uses rlnGroupScaleCorrection from optimiser star.

[huwjenkins]

There is a bug in relion_particle_subtract when the optional particle star file supplied with the --data option contains particles from fewer micrographs than the run_model.star referenced from the run_optimiser.star star file supplied with --i.

rlnGroupNumber is reset based on the contents of the --data star file and these rlnGroupNumbers are used to select the scale factor from the run_model.star file referenced in the run_optimiser.star file:

relion/src/particle_subtractor.cpp

Lines 742 to 750 in b978db9

	if (opt.do_scale_correction)
	{
	int group_id = opt.mydata.getGroupId(part_id);
	RFLOAT myscale = opt.mymodel.scale_correction[group_id];
	FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Fsubtract)
	{
	DIRECT_MULTIDIM_ELEM(Fsubtract, n) *= myscale;
	}
	}

To demonstrate I made particle star files from Refine3D/job029/run_data.star from the tutorial data (relion50_tutorial_precalculated_results.tar.gz). These files are here:
relion_particle_subtract_bug.tar.gz

selected_particles.star contains 1 particle from every micrograph except that corresponding to rlnGroupNumber 2 which has 3 particles. selected_particles_no_rlnGroupNumber_2.star is that file with all particles from rlnGroupNumber 2 removed.

To make it clearer to visualise I also modified run_model.star to set rlnGroupScaleCorrection for rlnGroupNumber 2 to an extreme value:

loop_
_rlnGroupNumber #1
_rlnGroupName #2
_rlnGroupNrParticles #3
_rlnGroupScaleCorrection #4
           1 Movies/20170629_00021_frameImage.mrc           74     0.915112
           2 Movies/20170629_00022_frameImage.mrc           83     999999.9
           3 Movies/20170629_00023_frameImage.mrc           83     0.959629
...

Running relion_particle_subtract with selected_particles.star:

`which relion_particle_subtract_mpi` --i Refine3D/job029/run_optimiser.star --mask MaskCreate/job020/mask.mrc --data Refine3D/job029/selected_particles.star --o Subtract/job036/ --float16

produces subtracted particles where the 2nd, 3rd and 4th particle have crazy stats as expected (output shortened to make it easier to read):

$ relion_image_handler --i Subtract/job036/particles_subtracted.star --stats
NOTE: the input (--i) is a STAR file but the output (--o) does not have .mrcs extension. The output is treated as a suffix, not a path.
000001@Subtract/job036/Particles/sub... ; avg= 0.00432972 stddev= 0.989898 minval= -4.48047 maxval= 4; angpix = 1.24484
000002@Subtract/job036/Particles/sub... ; avg= -nan stddev= nan minval= -inf maxval= inf; angpix = 1.24484
000003@Subtract/job036/Particles/sub... ; avg= inf stddev= nan minval= -65504 maxval= inf; angpix = 1.24484
000004@Subtract/job036/Particles/sub... ; avg= -nan stddev= nan minval= -inf maxval= inf; angpix = 1.24484
000005@Subtract/job036/Particles/sub... ; avg= 0.00468104 stddev= 0.996155 minval= -4.37891 maxval= 4.23828; angpix = 1.24484
000006@Subtract/job036/Particles/sub... ; avg= 0.00455536 stddev= 0.994194 minval= -4.66016 maxval= 4.00391; angpix = 1.24484
000007@Subtract/job036/Particles/sub... ; avg= 0.00454263 stddev= 0.98476 minval= -4.62891 maxval= 4; angpix = 1.24484
000001@Subtract/job036/Particles/sub... ; avg= 0.00485805 stddev= 0.99733 minval= -4.875 maxval= 4.03906; angpix = 1.24484
000002@Subtract/job036/Particles/sub... ; avg= 0.00469817 stddev= 1.00279 minval= -4.07422 maxval= 4; angpix = 1.24484
000003@Subtract/job036/Particles/sub... ; avg= 0.00479792 stddev= 0.998805 minval= -4.67969 maxval= 4.17578; angpix = 1.24484
000004@Subtract/job036/Particles/sub... ; avg= 0.00473028 stddev= 1.0032 minval= -4.50391 maxval= 4.07031; angpix = 1.24484
000005@Subtract/job036/Particles/sub... ; avg= 0.00470503 stddev= 0.999442 minval= -4.65234 maxval= 4.04297; angpix = 1.24484
000006@Subtract/job036/Particles/sub... ; avg= 0.00469012 stddev= 0.992939 minval= -4.58594 maxval= 3.64648; angpix = 1.24484
000007@Subtract/job036/Particles/sub... ; avg= 0.00442474 stddev= 0.999112 minval= -4.78125 maxval= 4; angpix = 1.24484
000001@Subtract/job036/Particles/sub... ; avg= 0.0044675 stddev= 1.01018 minval= -4.29688 maxval= 4.13281; angpix = 1.24484
000002@Subtract/job036/Particles/sub... ; avg= 0.00487376 stddev= 0.997607 minval= -4.51953 maxval= 4.00391; angpix = 1.24484
000003@Subtract/job036/Particles/sub... ; avg= 0.00495144 stddev= 0.991863 minval= -4.43359 maxval= 4.05078; angpix = 1.24484
000004@Subtract/job036/Particles/sub... ; avg= 0.00520154 stddev= 0.994975 minval= -4.20703 maxval= 4.46875; angpix = 1.24484
000005@Subtract/job036/Particles/sub... ; avg= 0.00487905 stddev= 1.00278 minval= -4.28125 maxval= 4; angpix = 1.24484
000006@Subtract/job036/Particles/sub... ; avg= 0.00477604 stddev= 1.00249 minval= -4.98438 maxval= 4.08203; angpix = 1.24484
000001@Subtract/job036/Particles/sub... ; avg= 0.00506635 stddev= 1.00124 minval= -4.36328 maxval= 4.02734; angpix = 1.24484
000002@Subtract/job036/Particles/sub... ; avg= 0.00496975 stddev= 0.998067 minval= -4.88281 maxval= 5.04297; angpix = 1.24484
000003@Subtract/job036/Particles/sub... ; avg= 0.00465052 stddev= 0.99753 minval= -4.46094 maxval= 4; angpix = 1.24484
000004@Subtract/job036/Particles/sub... ; avg= 0.00500923 stddev= 1.00828 minval= -4.61328 maxval= 4.05469; angpix = 1.24484
000005@Subtract/job036/Particles/sub... ; avg= 0.00492513 stddev= 1.0017 minval= -4.38281 maxval= 4.65234; angpix = 1.24484
000006@Subtract/job036/Particles/sub... ; avg= 0.00492848 stddev= 1.00095 minval= -4.32031 maxval= 4; angpix = 1.24484

However, running the same command with selected_particles_no_rlnGroupNumber_2.star:

`which relion_particle_subtract_mpi` --i Refine3D/job029/run_optimiser.star --mask MaskCreate/job020/mask.mrc --data Refine3D/job029/selected_particles_no_rlnGroupNumber_2.star --o Subtract/job037/ --float16

produces subtracted particles where only the 2nd particle has crazy stats and all subsequent particles have different stats from the previous example (again output shortened to make it easier to read):

$ relion_image_handler --i Subtract/job037/particles_subtracted.star --stats
NOTE: the input (--i) is a STAR file but the output (--o) does not have .mrcs extension. The output is treated as a suffix, not a path.
000001@Subtract/job037/Particles/sub... ; avg= 0.00432972 stddev= 0.989898 minval= -4.48047 maxval= 4; angpix = 1.24484
000002@Subtract/job037/Particles/sub... ; avg= -nan stddev= nan minval= -inf maxval= inf; angpix = 1.24484
000003@Subtract/job037/Particles/sub... ; avg= 0.0046192 stddev= 0.994307 minval= -4.66016 maxval= 4; angpix = 1.24484
000004@Subtract/job037/Particles/sub... ; avg= 0.00459128 stddev= 0.984935 minval= -4.625 maxval= 4; angpix = 1.24484
000005@Subtract/job037/Particles/sub... ; avg= 0.00478292 stddev= 0.997296 minval= -4.875 maxval= 4.03906; angpix = 1.24484
000006@Subtract/job037/Particles/sub... ; avg= 0.00466403 stddev= 1.00278 minval= -4.07422 maxval= 4; angpix = 1.24484
000001@Subtract/job037/Particles/sub... ; avg= 0.00477285 stddev= 0.998801 minval= -4.67969 maxval= 4.17578; angpix = 1.24484
000002@Subtract/job037/Particles/sub... ; avg= 0.004719 stddev= 1.00321 minval= -4.50391 maxval= 4.07031; angpix = 1.24484
000003@Subtract/job037/Particles/sub... ; avg= 0.00477376 stddev= 0.999394 minval= -4.65234 maxval= 4.04297; angpix = 1.24484
000004@Subtract/job037/Particles/sub... ; avg= 0.00474791 stddev= 0.993044 minval= -4.58594 maxval= 4; angpix = 1.24484
000005@Subtract/job037/Particles/sub... ; avg= 0.00478911 stddev= 0.999216 minval= -4.78125 maxval= 4; angpix = 1.24484
000006@Subtract/job037/Particles/sub... ; avg= 0.00428254 stddev= 1.01014 minval= -4.29297 maxval= 4.13281; angpix = 1.24484
000001@Subtract/job037/Particles/sub... ; avg= 0.00463858 stddev= 0.997741 minval= -4.51953 maxval= 4.00391; angpix = 1.24484
000002@Subtract/job037/Particles/sub... ; avg= 0.00490564 stddev= 0.992162 minval= -4.42969 maxval= 4.05078; angpix = 1.24484
000003@Subtract/job037/Particles/sub... ; avg= 0.00501422 stddev= 0.9947 minval= -4.21094 maxval= 4.46875; angpix = 1.24484
000004@Subtract/job037/Particles/sub... ; avg= 0.00506702 stddev= 1.00284 minval= -4.28125 maxval= 4; angpix = 1.24484
000005@Subtract/job037/Particles/sub... ; avg= 0.00489512 stddev= 1.00221 minval= -4.98438 maxval= 4.08203; angpix = 1.24484
000006@Subtract/job037/Particles/sub... ; avg= 0.0048224 stddev= 1.00146 minval= -4.36328 maxval= 4.02734; angpix = 1.24484
000001@Subtract/job037/Particles/sub... ; avg= 0.00493526 stddev= 0.998076 minval= -4.88281 maxval= 5.04297; angpix = 1.24484
000002@Subtract/job037/Particles/sub... ; avg= 0.0048767 stddev= 0.99748 minval= -4.45703 maxval= 4; angpix = 1.24484
000003@Subtract/job037/Particles/sub... ; avg= 0.00488182 stddev= 1.00826 minval= -4.61328 maxval= 4.05469; angpix = 1.24484
000004@Subtract/job037/Particles/sub... ; avg= 0.00498734 stddev= 1.00136 minval= -4.38281 maxval= 4.65234; angpix = 1.24484
000005@Subtract/job037/Particles/sub... ; avg= 0.00492791 stddev= 1.00095 minval= -4.32031 maxval= 4; angpix = 1.24484

To further confirm running the command:

grep frameImage  Subtract/job0NN/particles_subtracted.star |awk '{print $17, $7}'

confirms rlnGroupNumber corresponds to different micrographs in the 2 files:

Subtract/job036/particles_subtracted.star:

1 MotionCorr/job002/Movies/20170629_00021_frameImage.mrc
2 MotionCorr/job002/Movies/20170629_00022_frameImage.mrc
2 MotionCorr/job002/Movies/20170629_00022_frameImage.mrc
2 MotionCorr/job002/Movies/20170629_00022_frameImage.mrc
3 MotionCorr/job002/Movies/20170629_00023_frameImage.mrc
4 MotionCorr/job002/Movies/20170629_00024_frameImage.mrc
5 MotionCorr/job002/Movies/20170629_00025_frameImage.mrc
6 MotionCorr/job002/Movies/20170629_00026_frameImage.mrc
7 MotionCorr/job002/Movies/20170629_00027_frameImage.mrc
8 MotionCorr/job002/Movies/20170629_00028_frameImage.mrc
9 MotionCorr/job002/Movies/20170629_00029_frameImage.mrc
10 MotionCorr/job002/Movies/20170629_00030_frameImage.mrc
11 MotionCorr/job002/Movies/20170629_00031_frameImage.mrc
12 MotionCorr/job002/Movies/20170629_00035_frameImage.mrc
13 MotionCorr/job002/Movies/20170629_00036_frameImage.mrc
14 MotionCorr/job002/Movies/20170629_00037_frameImage.mrc
15 MotionCorr/job002/Movies/20170629_00039_frameImage.mrc
16 MotionCorr/job002/Movies/20170629_00040_frameImage.mrc
17 MotionCorr/job002/Movies/20170629_00042_frameImage.mrc
18 MotionCorr/job002/Movies/20170629_00043_frameImage.mrc
19 MotionCorr/job002/Movies/20170629_00044_frameImage.mrc
20 MotionCorr/job002/Movies/20170629_00045_frameImage.mrc
21 MotionCorr/job002/Movies/20170629_00046_frameImage.mrc
22 MotionCorr/job002/Movies/20170629_00047_frameImage.mrc
23 MotionCorr/job002/Movies/20170629_00048_frameImage.mrc
24 MotionCorr/job002/Movies/20170629_00049_frameImage.mrc

Subtract/job037/particles_subtracted.star:

1 MotionCorr/job002/Movies/20170629_00021_frameImage.mrc
2 MotionCorr/job002/Movies/20170629_00023_frameImage.mrc
3 MotionCorr/job002/Movies/20170629_00024_frameImage.mrc
4 MotionCorr/job002/Movies/20170629_00025_frameImage.mrc
5 MotionCorr/job002/Movies/20170629_00026_frameImage.mrc
6 MotionCorr/job002/Movies/20170629_00027_frameImage.mrc
7 MotionCorr/job002/Movies/20170629_00028_frameImage.mrc
8 MotionCorr/job002/Movies/20170629_00029_frameImage.mrc
9 MotionCorr/job002/Movies/20170629_00030_frameImage.mrc
10 MotionCorr/job002/Movies/20170629_00031_frameImage.mrc
11 MotionCorr/job002/Movies/20170629_00035_frameImage.mrc
12 MotionCorr/job002/Movies/20170629_00036_frameImage.mrc
13 MotionCorr/job002/Movies/20170629_00037_frameImage.mrc
14 MotionCorr/job002/Movies/20170629_00039_frameImage.mrc
15 MotionCorr/job002/Movies/20170629_00040_frameImage.mrc
16 MotionCorr/job002/Movies/20170629_00042_frameImage.mrc
17 MotionCorr/job002/Movies/20170629_00043_frameImage.mrc
18 MotionCorr/job002/Movies/20170629_00044_frameImage.mrc
19 MotionCorr/job002/Movies/20170629_00045_frameImage.mrc
20 MotionCorr/job002/Movies/20170629_00046_frameImage.mrc
21 MotionCorr/job002/Movies/20170629_00047_frameImage.mrc
22 MotionCorr/job002/Movies/20170629_00048_frameImage.mrc
23 MotionCorr/job002/Movies/20170629_00049_frameImage.mrc

This probably doesn't have much effect as the rlnGroupScaleCorrection values are often close to 1 but it is incorrect. The relevant code seems to be here:

relion/src/exp_model.cpp

Lines 954 to 958 in b978db9

	// The group number is only set upon reading: it is not read from the STAR file itself,
	// there the only thing that matters is the order of the micrograph_names
	// Write igroup+1, to start numbering at one instead of at zero
	MDimg.setValue(EMDL_MLMODEL_GROUP_NO, group_id + 1, part_id);

[huwjenkins]

As mentioned on CCPEM mailing list, if the particle star file contains more group_ids than the run_model.star file referenced in the run_optimiser.star file then opt.mymodel.scale_correction[group_id] eventually becomes undefined.

To demonstrate I made 100 copies of the first particle and incremented rlnMicrographName from MotionCorr/job002/Movies/20170629_00001_frameImage.mrc
to
MotionCorr/job002/Movies/20170629_00100_frameImage.mrc.
Some of the subtracted particles then have nonsense pixel values:

$ relion_image_handler --i Subtract/job038/particles_subtracted.star --stats
NOTE: the input (--i) is a STAR file but the output (--o) does not have .mrcs extension. The output is treated as a suffix, not a path.
000001@Subtract/job038/Particles/sub... ; avg= 0.00442051 stddev= 0.989596 minval= -4.48153 maxval= 3.99965; angpix = 1.24484
...
000071@Subtract/job038/Particles/sub... ; avg= 4.81034e-05 stddev= 0.990762 minval= -4.52113 maxval= 4.00155; angpix = 1.24484
000072@Subtract/job038/Particles/sub... ; avg= nan stddev= -nan minval= 1.79769e+308 maxval= -1.79769e+308; angpix = 1.24484
000073@Subtract/job038/Particles/sub... ; avg= nan stddev= -nan minval= 1.79769e+308 maxval= -1.79769e+308; angpix = 1.24484
000074@Subtract/job038/Particles/sub... ; avg= 4.81034e-05 stddev= 0.990762 minval= -4.52113 maxval= 4.00155; angpix = 1.24484

[biochem-fan]

Thank you very much for your report. Indeed this is a bug but I myself do not have time to fix this. I am happy to review a pull request if someone sends one.