rfree flag propogation

[cvazz]

Hi all,
I have been resorting to a hack for a little while now and think I should report this so that this can get fixed.

The following procedure "breaks" phenix, so that it does not recognize my file as xray data:

ds_mtz = rs.read_mtz(path)
meteor_map = rsmap.Map(ds_mtz)
meteor_map.write_mtz(loc)

So, what I have been using in the past, was to modify my meteor map and save it as extra columns to a ds_mtz2 and potentially delete unneeded columns.

So far I couldn't figure out why one works and the other doesn't, but if this is not just an issue on my machine, it definitely needs fixing!

[cvazz]

While I feel like this was not the only issue, I did notice that I also had problems when leaving out R-free flags. As of right now, it is not possible to include rfree flags in a meteor map, what do you think about allowing R-free flags to be part of rsmap?

[cvazz]

To clarify my thoughts and challenges:

1. As of right now, it is not possible to include rfree flags in a meteor map, but to me it feels like something maps should have (and keep)?

2. After combining a difference map with the dark map, some SIGF values were missing, which the difference map function was not able to compute. This creates headaches, for example when using the created data for refinement. This headache could be removed by having an approved arithmetic function, which could contain flags depending on what data should be returned (only nonnan rows for example).

As an example for how a meteor addition/subtraction function could look like, here my version. Not sure if i should relabel this thread a feature request now.

def adding_maps(map1: rsmap.Map, map2: rsmap.Map, *, factor1=1, factor2=1, suppress_warnings=False) -> rsmap.Map:
    finite_mask1 = np.isfinite(map1.amplitudes)
    finite_mask2 = np.isfinite(map2.amplitudes)
    map_1 = map1[finite_mask1]
    map_2 = map2[finite_mask2]
    if (~finite_mask1).any() or (~finite_mask2).any() and not suppress_warnings:
        logger.warning(
            f"Some Maps have non-overlapping finite amplitudes: map1 {np.sum(~finite_mask1)} NaNs, map2 {np.sum(~finite_mask2)} NaNs"
        )

    common_indices = map_1.index.intersection(map_2.index)

    structure_factors1 = map_1.loc[common_indices].to_structurefactor()
    structure_factors2 = map_2.loc[common_indices].to_structurefactor()
    added_structure_factors = (
        factor1 * structure_factors1 + factor2 * structure_factors2
    )

    sum_of_map  = rsmap.Map.from_structurefactor(
        added_structure_factors,
        index=common_indices,
        cell=map1.cell,
        spacegroup=map1.spacegroup,
    )
    logger.debug(f"High resolution limit after addition: {np.min(sum_of_map.compute_dHKL()):.2f}.")

    if  map1.has_uncertainties and map2.has_uncertainties:
        sigmaF = np.sqrt(
          (factor1 * map_1.uncertainties[common_indices]) ** 2
        + (factor2 * map_2.uncertainties[common_indices]) ** 2
        )
        sum_of_map.set_uncertainties(pd.Series(sigmaF, index=common_indices))
    return sum_of_map

[tjlane]

Hi @cvazz, two requests to help us:

1. Can you post the phenix error along with a complete set of code/input data needed to reproduce it?
2. Let's branch the Map operator discussion to a new issue, please!

[cvazz]

As I said above, the main issue is that it could be nice to carry around rfree-flags. I cannot save a meteor map back to file together with its R-free-flags. One needs r-free-flags in phenix. Without using the original r-free flags (or manipulating the original pdb file), one gets errors. For example, when creating this map

path = meteor_data_rsEGFP2
loc = "test.mtz"
ds_mtz = rs.read_mtz(path, )
meteor_map = rsmap.Map(ds_mtz,amplitude_column="F_on", uncertainty_column="SIGF_on", phase_column="PHIC_nochrom")
meteor_map.write_mtz(loc)

together with

phenix.refine test.mtz 8a6g-chromophore-removed.pdb xray_data.r_free_flags.generate=True

Together this leads to the following (resovable) error message.

As stated above, there are two fixes. One is to clean the pdb file. The other would be to save the rfree flags with the data (which currently means reverting to a rs.dataset.

A concrete To Do from here would be to to clean up the pdb file 8a6g-no-chrom by deleting the first 50 or lines of remarks. (Until the line that says "delete everything above if you want to submit to pdb"). This allows one to refine at least the present dataset.

Another, more long term, option would be to add an rfree column.

I also had issues when sharing files with Dmitrii. I have to verify if this was the only issue, however.

Generating a new array of R-free flags.

*******************************************************************************
*******************************************************************************

The MD5 checksum for the R-free flags array summarized above is:
  1c0bb84bc16dde71a563e1f99ca45f3d

The corresponding MD5 checksum in the PDB file summarized above is:
  9a82d61b30a74bac78dca939236526d2

These checksums should be identical but are in fact different. This is
because the R-free flags used at previous stages of refinement are
different from the R-free flags summarized above. As a consequence,
the values for R-free could be biased and misleading.

However, there is no problem if the R-free flags were just extended to
a higher resolution, or if some reflections with no data or that are
not part of the R-free set have been added or removed. In this case,
simply remove the

  REMARK r_free_flags.md5.hexdigest 9a82d61b30a74bac78dca939236526d2

record from the input PDB file to proceed with the refinement.

Otherwise it is best to recover the previously used R-free flags
and use them consistently throughout the refinement of the model.
Run this command again with the name of the file containing the
original flags as an additional input.

If the original R-free flags are unrecoverable, remove the REMARK
record as indicated above. In this case the values for R-free will
become meaningful only after many cycles of refinement.

*******************************************************************************
*******************************************************************************

Sorry: Please resolve the R-free flags mismatch.

[tjlane]

@cvazz great thanks for the more detailed report, makes sense now

I reclassified this as a feature request rather than a bug -- right now, meteor relies on the user to handle their own R-free flags.

Refinement with meteor maps is not supported just yet, since we need to fully support extrapolation: robust occupancy refinement, error propagation, etc. So let's organize this underneath a project aiming to enable refinement as standard practice!