DIALS σ(I) estimates are considerably smaller than those from XDS

[dagewa]

I had a case where I could not solve a small molecule 3D ED structure with shelxt 2018/2, but it solved fine with data from XDS. The cause appears to be because in that version of sheltx:

The mean I/sigma(I) for systematic absences (where at least 10 of them were recorded) is now reported and space groups rejected for which it is greater than 5.0.

The "absences" were slightly present for this data set, and due to rather low error estimates from DIALS, the correct space group was rejected.

This ties in with anecdotal evidence that DIALS error estimates generally seem too small, or equivalently, reported I/σ(I) seems too high.

All this sent me on an investigation into error estimation in DIALS compared to XDS, using the following script:

#!/bin/bash
set -e

# Check script input
if [ "$#" -ne 1 ]; then
    echo "You must supply the location of a diffraction image from " \
"a single scan"
    exit 1
fi

IMAGE=$(realpath "$1")

mkdir -p xia2-dials && cd xia2-dials
xia2 image="$IMAGE" xia2.settings.resolution.d_min=2
cd ..

mkdir -p xia2-xds && cd xia2-xds
xia2 pipeline=3dii image="$IMAGE" xia2.settings.resolution.d_min=2
cd ..

curl -O https://gist.githubusercontent.com/dagewa/16c9caa5841443c7611b6d6b70d2c4f9/raw/compare_errors.py
chmod +x compare_errors.py
time ./compare_errors.py

1. This processes a single sweep with xia2 pipeline=dials and then with xia2 pipeline=3dii
2. The scaled reflections from both DIALS and XDS are then read in by the compare_errors.py script
3. The reflections are matched
4. The DIALS intensities and errors are scaled to the XDS ones using a robust linear fit between the DIALS and XDS intensities
5. The XDS vs DIALS intensities are plotted with the robust fit line (i.e. with a slope of 1 after the rescaling)
6. The XDS vs DIALS intensity sigmas are plotted. A robust fit line is calculated on the stronger reflections and added to the plot to give an idea of the overall scale between XDS and DIALS errors

What follows are results from some of the data sets available in DIALS data. In summary, the XDS error estimates are always higher than the DIALS ones, but this varies between 1.2 times for thaumatin_i04 and 4.2 times for x4wide. No claim is made as to which is more "correct", but the discrepancies and some of the features in the plots seem surprising.

I would appreciate if anyone could double-check my analysis and comment on the observations.

fumarase

compare_errors.sh $(dials.data get -q fumarase)/12287_1_E1_001.img

insulin

compare_errors.sh $(dials.data get -q insulin)/insulin_1_001.img

mpro_x0692

compare_errors.sh $(dials.data get -q mpro_x0692)/Mpro-x0692_1_0001.cbf

spring8_ccp4_2018

compare_errors.sh $(dials.data get -q spring8_ccp4_2018)/ccp4school2018_bl41xu/05/data03/data03_master.h5

small_molecule_example

compare_errors.sh $(dials.data get -q small_molecule_example)/x3_1_0001.cbf.gz

NB xia2 pipeline=3dii failed at pointless, but we had got far enough anyway

thaumatin_i04

compare_errors.sh $(dials.data get -q thaumatin_i04)/th_8_2_0001.cbf

x4wide

compare_errors.sh $(dials.data get -q x4wide)/X4_wide_M1S4_2_0001.cbf

[dagewa]

Simulated data

As a counter example to the trend above, when I run this analysis on one of James Holton's simulated data sets, then the DIALS error estimates are actually larger than XDS. In this case the scale factors are essentially flat.

curl -O https://bl831.als.lbl.gov/~jamesh/benchmarks/get_test_data.com
chmod +x get_test_data.com
./get_test_data.com
../compare_errors.sh data/core_001.cbf

(in this case there was not enough memory to run the robust linear regression, so that was replaced with np.polyfit(x, y, deg=1))

[huwjenkins]

This is very interesting! I looked at the thaumatin_i04 data with a slightly different approach.

1. Integrated data as usual with DIALS.
2. Scaled with dials.scale using:

dials.scale symmetrized.expt symmetrized.refl d_min=1.35 intensity_choice=profile

high res cutoff much more generous than you used as I/σ(I) was ~20 at 2 Å for this dataset.

3. Exported only scaled reflections with:

dials.filter_reflections scaled.refl flag_expression="scaled&integrated_prf"
dials.export scaled.expt filtered.refl format=xds_ascii intensity=profile

4. Scaled DIALS.HKL with XSCALE.
5. Ran your script (I also had to use np.polyfit(x, y, deg=1))
   
   

Some observations:

I have no strange behaviour for reflections where I ~ 0 from DIALS. Does SCALED_NATIVE_SWEEP1.HKL contain σ(I) < 0? the CORRECT step of XDS uses this to mark misfits that should be removed - your script doesn't appear to check for these. Ignore - these are reflections for which I from DIALS = 0.

The plot of σ(I) has a slope of 1.182 which exactly matches the estimated I/σ(I) asymptotic limits reported by dials.scale and XSCALE

Most interesting are the error model parameters...

Error model details:
  Type: basic
  Parameters: a = 0.67441, b = 0.02787

I think you need to take square root of both a and b from XDS to match dials.scale as the formulae are different:

v(I) = a[v₀(I) + bI²] vs σ'² = a²(σ² + (bI)²)

so:

       XDS      dials.scale
a    0.80932    0.67441
b    0.02742    0.02787
ISa  45.06      53.196

[huwjenkins]

So of course I did:

dials.scale symmetrized.expt symmetrized.refl d_min=1.35 intensity_choice=profile a=0.80932 b=0.02742

This gives:

Error model details:
  Type: basic
  Parameters: a = 0.80932, b = 0.02742
  Error model formula: σ'² = a²(σ² + (bI)²)
  estimated I/sigma asymptotic limit: 45.062

The estimated I/σ(I) asymptotic limit is now identical to that reported by XSCALE and the graphs are as you'd expect from this

So I would concentrate the investigation on the differences in error model refinement between XDS and dials.scale.

[graeme-winter]

Some thoughts:

• collecting "the same" data set a number of times is something I have done a few times in the past, which could allow some semi-independent way of determining uncertainties on reflections (at least, as a reference)
• we can do a little "shuffling" on that data to generate statistically valid if not strictly true real random samples from a population: take every count on pixel[k,j,i] and sum them across n data sets before re-distributing them
• looking at the σ(I) estimates from very weak / absent reflections as commented above can also be explored by taking a clean data set from e.g. a centred sample (ohai insulin with I213) and processing with a primitive lattice, thus generating a lot which should be 0+/-σ(I)

I will dig out my scripts for the second of these: when I discussed this with Randy at a GRC many years ago he did suggest that it wasn't really horrible - I think I called it "semi-synthetic" or something at the time

[biochem-fan]

I remember Andrew Leslie (by Garib's suggestion?) performed the third analysis on MOSFLM by doubling the unit cell lengths.

[graeme-winter]

I remember Andrew Leslie (by Garib's suggestion?) performed the third analysis on MOSFLM by doubling the unit cell lengths.

Yes, you can do that with DIALS as well with dials.reindex space_group=P1 change_of_basis_op=2a,2b.2c straight out of dials.index - the space_group= is important as cctbx will "remember" otherwise

[huwjenkins]

The reflections from dials with I=0 appear to be a xia2 thing. I cannot reproduce with 'normal' dials workflow. They have I=0 σ(I)=0:

>>> from dials.array_family import flex
>>> refl = flex.reflection_table.from_file("xia2-dials/DataFiles/AUTOMATIC_DEFAULT_scaled.refl")
>>> sel = (refl['intensity.scale.value'] == 0 and refl['intensity.scale.variance'] == 0)
>>> sel.count(True)
15914

processed using xia2 image=../../th_8_2_0001.cbf xia2.settings.resolution.d_min=1.35

@dagewa suggest you filter these out in your script. Culprit seems to be output.include_bad_reference=True option in xia2's dials.integrate command. These reflections correctly do not have scaled flag set and so should be excluded in this analysis.

$ dials.filter_reflections AUTOMATIC_DEFAULT_scaled.refl

...
No filter specified. Performing analysis instead.
+--------------------------------+--------+
| flag                           |   nref |
|--------------------------------+--------|
| background_includes_bad_pixels |  35400 |
| bad_for_refinement             |      2 |
| bad_reference                  |  15914 | <-----
| centroid_outlier               |  14312 |
...

[dagewa]

Thanks @huwjenkins, I will update the script. @graeme-winter investigating this with a "semi-synthetic" dataset (or a few) sounds like an excellent idea.

As for the systematic absences, I still need to really dig into it for the data set mentioned at the top of the thread. I suspect the absence-breaking to be "true", as this is ED on an inorganic crystal. Olex2 produces a plot of the intensity distribution of absence-breaking reflections, which looks like this:

[dagewa]

For XDS the plot for the same data set looks like this. Really quite different!

[dagewa]

Hmm, absence-breaking does seem quite significant in this case. If I've done this right, I selected systematic absences only with this script:

#!/bin/env dials.python

from dxtbx.model.experiment_list import ExperimentList
from dials.array_family import flex

expt = ExperimentList.from_file("scaled.expt")
refl = flex.reflection_table.from_file("scaled.refl")

sel_abs = flex.bool(len(refl), False)

for i in range(len(expt)):
    i_sel = (refl["id"] == i).iselection()
    i_refl = refl.select(i_sel)
    i_expt = expt[i]
    abs_flag = i_refl.as_miller_array(i_expt).sys_absent_flags().data()
    sel_abs.set_selected(i_sel, abs_flag)

refl = refl.select(sel_abs)
refl.as_file("absences.refl")

then looked at the 364 remaining reflections in the image viewer. Some of them have very clear spots, like this:

The space group is $F d d 2$

[huwjenkins]

For the small molecule ED dataset what does the table

+--------------------------+----------+------------------------+----------------------+
| Intensity range (<Ih>)   |   n_refl |   Uncorrected variance |   Corrected variance |
|--------------------------+----------+------------------------+----------------------|

in the dials.scale logfile look like? For small molecule ED datasets I've often found unless you're scaling datasets from multiple crystals together you need to reduce min_Ih to get a reasonable number of reflections included in the refinement of the error model parameters.

[huwjenkins]

Also how are you going from XDS_ASCII -> SHELX HKLF format - it looks from the plot from Olex2 that the XDS intensities are much larger than the ones from DIALS?

[dagewa]

Thanks Huw, yes I reduced min_Ih to 10 in this case, which did increase the sigmas a little. This is merged from 3 crystals though. I also increased the gain of the detector (Timepix) to 1.35 as it looked like it might be a little over 1 due to charge-sharing between pixels. That also increased sigmas a bit. Here's the table:

+--------------------------+----------+------------------------+----------------------+
| Intensity range (<Ih>)   |   n_refl |   Uncorrected variance |   Corrected variance |
|--------------------------+----------+------------------------+----------------------|
| 2257.50 - 358.69         |      126 |                 18.9   |                0.904 |
| 358.69 - 226.58          |      126 |                 12.661 |                1.2   |
| 226.58 - 172.80          |      126 |                 11.386 |                1.299 |
| 172.80 - 137.33          |      126 |                  8.102 |                0.992 |
| 137.33 - 112.71          |      126 |                  7.377 |                0.946 |
| 112.71 - 87.41           |      161 |                  9.543 |                1.296 |
| 87.41 - 50.84            |      454 |                  6.918 |                0.998 |
| 50.84 - 29.57            |      655 |                  7.187 |                1.092 |
| 29.57 - 17.20            |      899 |                  6.335 |                0.985 |
| 17.20 - 9.99             |      988 |                  5.241 |                0.822 |
+--------------------------+----------+------------------------+----------------------+

[dagewa]

Scaling and conversion for XDS:

cat <<EOF > XSCALE.INP
OUTPUT_FILE=temp.ahkl
INPUT_FILE=../Data1/XDS_ASCII.HKL
INPUT_FILE=../Data3/XDS_ASCII.HKL
INPUT_FILE=../Data4/XDS_ASCII.HKL
SPACE_GROUP_NUMBER= 22
UNIT_CELL_CONSTANTS= 18.374 18.669 6.767 90 90 90
RESOLUTION_SHELLS=4 3 2 1 0.61
EOF
xscale

cat <<EOF > XDSCONV.INP
INPUT_FILE=temp.ahkl
OUTPUT_FILE=xds.hkl  SHELX
FRIEDEL'S_LAW=FALSE
EOF
xdsconv

NB the space group is $Fdd2$ as correctly found by shelxt. For the DIALS sys abs results above that was set prior to scaling, which I didn't do here for XDS

[dagewa]

There does seem to be something odd going on. Here's the compare_errors.py plots for this data

[dagewa]

If I don't rescale the DIALS intensities and errors to match the XDS ones (because anyway the linear fit looks bad) then I get this.


So indeed, XDS intensities are much higher than DIALS here 🤔

[dagewa]

Maybe just an effect of the overall scale factor being arbitrary and XDS adjusting to best use a fixed-width format?

[huwjenkins]

In my experience XSCALE always results in very large intensities - there'll be a line in XSCALE.LP with something like

FACTOR TO PLACE ALL DATA SETS TO AN APPROXIMATE ABSOLUTE SCALE 0.636846E+04
 (ASSUMING A PROTEIN WITH 50% SOLVENT)

but I'm intrigued that these are getting through XDSCONV. This must have changed because I'm sure it used to scale so intensities fitted into the 8.2F format but now maybe it doesn't and takes advantage that SHELX programs should be able to read 9999999. from this format. What is maximum and minimum in xds.hkl

[dagewa]

For this case,

 FACTOR TO PLACE ALL DATA SETS TO AN APPROXIMATE ABSOLUTE SCALE 0.469060E+02
 (ASSUMING A PROTEIN WITH 50% SOLVENT)

The five-number summary of the I column in xds.hkl looks like this:
-18426.00 7010.28 17819.00 47216.00 999999.00

I guess it is scaled so that the maximum intensity is 999999.00?

[huwjenkins]

Probably I should change format=shelx option in dials.export to do this too - I think based on the various fixes to iotbx/shelx/hklf.py over the years I now understand exactly what this format can contain 🙂

[huwjenkins]

possibly the file contains 999999. which is the maximum that column can contain and still leave a space after l but yes it looks like both XSCALE and XDSCONV are scaling up increasing the the intensity and dials.export is presumably scaling it so maximum I is 9999.99 (at least I think that's the default) so I don't think the absolute values on the plots in Olex2 are anything to worry about. I was just interested because I thought I remembered XDSCONV behaviour being different but probably I hallucinated this as I can't see anything in XDS release notes.

The shape of those I(XDS) vs I(dials) plots are quite surprising though. For those are you using the XSCALE and scaled.refl output as for the other cases?

[huwjenkins]

@dagewa I think your script should really filter the reflection table by refl.get_flags(refl.flags.scaled). I just tested this on some of our ED data and this behaviour gave me quite a shock till I realised what was happening!

>>> from dials.array_family import flex
>>> refl = flex.reflection_table.from_file('scaled.refl')
>>> max(refl["intensity.scale.value"])
246564.43974757672
>>> min(refl["intensity.scale.value"])
-402839.77312233415
>>> refl = refl.select(refl.get_flags(refl.flags.scaled))
>>> max(refl["intensity.scale.value"])
278.9691218882821
>>> min(refl["intensity.scale.value"])
-4.726715612925455

[dagewa]

Yes, that seems sensible. Ok, I updated the script. I think maybe the robust fit line should only be calculated on the stronger intensities too. Plus sometimes calculating that fails when there are too many points. After I improve the script a bit more I may re-calculate plots for the examples above.

[graeme-winter]

XSCALE places the intensities on an absolute scale assuming 50% solvent

Output from CORRECT is on image scale

[huwjenkins]

@dagewa - I looked a bit more at this with our biotin ED datasets on zenodo. I think there are 2 different issues here. The MX cases where you are comparing one dataset processed with XDS and DIALS is different to the ED case with 3 datasets.

When you process the ED data are you scaling all three datasets together in DIALS and refining 1 set of error model parameters for all 3 datasets (error_model.grouping.combined) wheras in XDS/XSCALE you are scaling all datasets to one reference (by default the first INPUT_FILE) and refining error model parameters for all 3 datasets separately?

Does that explain the plots? Could you try adding id = flex.int(int(i) for i in xds_refl[9]) and xds_refl["id"] = id to the script so you can colour the points by the dataset they come from?

[huwjenkins]

Here's my analysis. The data are biotin_xtal{1..5}.tar.xz from here.

Process with this script.

Colour series is orange, blue, green, pink, yellow for datasets 0-4 (xtal1 to xtal5)

min_Ih=25:

Error model details:
  Type: basic
  Parameters: a = 4.30799, b = 0.00000
  Error model formula: σ'² = a²(σ² + (bI)²)
  estimated I/sigma asymptotic limit: 52482507136502743040.000

min_Ih=10:

Error model details:
  Type: basic
  Parameters: a = 3.49346, b = 0.01498
  Error model formula: σ'² = a²(σ² + (bI)²)
  estimated I/sigma asymptotic limit: 19.103

min_Ih=5:

Error model details:
  Type: basic
  Parameters: a = 2.79170, b = 0.02591
  Error model formula: σ'² = a²(σ² + (bI)²)
  estimated I/sigma asymptotic limit: 13.823

min_Ih=1:

Error model details:
  Type: basic
  Parameters: a = 1.49682, b = 0.07238
  Error model formula: σ'² = a²(σ² + (bI)²)
  estimated I/sigma asymptotic limit: 9.231

min_Ih=0:

Error model details:
  Type: basic
  Parameters: a = 1.29630, b = 0.08832
  Error model formula: σ'² = a²(σ² + (bI)²)
  estimated I/sigma asymptotic limit: 8.734

XSCALE:

     a        b          ISa    ISa0   INPUT DATA SET
 4.706E+00  4.609E-02    2.15  316.23 DIALS_0.HKL
 1.715E+01  4.047E-03    3.80  316.23 DIALS_1.HKL
 1.965E+01  1.840E-03    5.26  316.23 DIALS_2.HKL
 7.135E+00  1.563E-02    2.99  316.23 DIALS_3.HKL
 7.425E+00  3.587E-03    6.13  316.23 DIALS_4.HKL

Isn't this the issue discussed before with respect to applying a pedestal to the raw data - XDS/XSCALE fits the error model parameters using all reflections i.e. no Ih cutoff like in dials.scale

[dagewa]

Thanks @huwjenkins. I got caught up with some deadlines so haven't been able to look at this for a while, but I will pick it up again eventually.