From 97e6d99d037abcb864bcd0a24546bc63e720b42b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 12 Jul 2026 21:34:52 +0200 Subject: [PATCH] docs: document French-Wilson amplitudes, R-free flags and model-based maps MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Update CPU_DATA_ANALYSIS.md for the new merge/analysis capabilities: - §10.7 R-free test-set flags: Friedel-merged Laue-ASU key (Bijvoet-safe), deterministic splitmix64 selection, resolution stratification; - §10.8 French-Wilson amplitudes: posterior-mean |F| from I and sigma with the Wilson prior (centric/acentric, epsilon multiplicity, strong-reflection short-circuit), written as MTZ F/SIGF and mmCIF F_meas_au; - §14 model-based validation (rugnux --model): Fcalc, flat bulk solvent + overall and per-shell scaling, R-work/R-free vs a model, and 2Fo-Fc / Fo-Fc maps. Also refresh the scope list and references, and correct the old "intensities only" note (the output now carries amplitudes too). Co-Authored-By: Claude Opus 4.8 (1M context) --- docs/CPU_DATA_ANALYSIS.md | 73 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 71 insertions(+), 2 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 3e753c03..cb4b9f67 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -14,7 +14,9 @@ This document describes the crystallographic algorithms implemented in Jungfrauj 8. Bragg integration by either 2D box summation or profile fitting (Kabsch, reference-free), 9. scaling and merging, 10. merge-level error modelling and outlier rejection, -11. auxiliary statistics (Wilson plot, ⟨I/σ(I)⟩, CC1/2, CCref). +11. auxiliary statistics (Wilson plot, ⟨I/σ(I)⟩, CC1/2, CCref), +12. amplitude estimation (French–Wilson) and R-free test-set flagging, +13. optional model-based validation: R-free against a supplied model and 2Fo−Fc / Fo−Fc electron-density maps. ## References @@ -25,6 +27,9 @@ The methods are inspired and reuising solutions implemented in: - T. A. White et al., CrystFEL method papers (spot finding, three‑ring integration, serial/still diffraction processing concepts). - J. Kieffer & J. P. Wright, "PyFAI: a Python library for high performance azimuthal integration on GPU", *Powder Diffraction* **28** (2013), S339-S350 (detector geometry definition, azimuthal integration) - H. Powell, "The Rossmann Fourier autoindexing algorithm in MOSFLM", *Acta Cryst.* **D55** (1999), 1690-1695 (FFT indexing) +- S. French & K. Wilson, "On the treatment of negative intensity observations", *Acta Cryst.* **A34** (1978), 517-525 (Bayesian amplitude estimation from intensities). +- A. T. Brünger, "Free R value: a novel statistical quantity for assessing the accuracy of crystal structures", *Nature* **355** (1992), 472-475 (R-free cross-validation). +- M. Wojdyr, "GEMMI: A library for structural biology", *J. Open Source Softw.* **7** (2022), 4200 (model / structure-factor / map machinery used in §14). (list is not exhaustive) ## 1. Geometry, reciprocal-space mapping, and basic quantities @@ -561,6 +566,34 @@ After scale-fulls, two **optional correction surfaces** can be fitted on the com Both surfaces are **cross-validated**: fitted on even-numbered frames and kept only if they improve the held-out odd-frame symmetry-equivalent agreement by a clear margin (and vice versa). A surface fitted to noise where its systematic is absent therefore does not generalize and is discarded — an opt-in correction never adds scatter. +### 10.7 R-free test-set flags + +A fraction of the unique reflections (`rfree_fraction`, default 0.05) is flagged as a **free (test) set**, written to the output (MTZ `FreeR_flag`, mmCIF `_refln.status_free`, a text-HKL column) for model validation (§14) and for downstream refinement. The flag is a pure function of the reflection's **Friedel-merged (Laue) ASU index**, which gives three properties: + +- all symmetry- and Friedel-equivalent reflections share one flag — in particular a Bijvoet pair $I(+)/I(-)$, kept as two separate merged rows in anomalous mode, is **never split** across the work and free sets (which would bias R-free); +- the free/work decision is a deterministic hash of that key, so the same reflection always lands in the same set — reproducible run-to-run and independent of the order in which observations were merged; +- the set is **stratified by resolution**: within each of 20 shells (equal-volume in $1/d^2$), exactly $\mathrm{round}(\text{fraction}\cdot n)$ of the distinct reflections — those with the smallest hash — are flagged free, so ~5 % of the data is free in every shell. + +### 10.8 French–Wilson amplitudes + +The last step of the merge estimates a Bayesian structure-factor amplitude $|F|$ for each unique reflection from its intensity $I$ and error $\sigma$, so the output carries amplitudes alongside intensities (a naïve $\sqrt{\max(I,0)}$ turns every weak or negative measurement into a biased — or zero — amplitude). With the Wilson prior for the true intensity $J\ge 0$ at that resolution, + +$ +P_\mathrm{acentric}(J) \propto e^{-J/\Sigma},\qquad +P_\mathrm{centric}(J) \propto J^{-1/2}\,e^{-J/2\Sigma}, +$ + +and a Gaussian likelihood $\mathcal{N}(I;J,\sigma^2)$, the posterior mean amplitude and its uncertainty are + +$ +\langle |F|\rangle = \frac{\int_0^\infty \sqrt{J}\,\mathcal{N}(I;J,\sigma^2)\,P(J)\,\mathrm{d}J}{\int_0^\infty \mathcal{N}(I;J,\sigma^2)\,P(J)\,\mathrm{d}J},\qquad +\sigma_F = \sqrt{\langle J\rangle - \langle|F|\rangle^2}. +$ + +The prior mean is $\Sigma = \varepsilon\,\langle I/\varepsilon\rangle_\mathrm{shell}$, where $\varepsilon$ is the reflection's epsilon (symmetry-enhancement) multiplicity and $\langle I/\varepsilon\rangle$ is the Wilson mean in its resolution shell (so reflections on symmetry elements, and each shell, are treated correctly). Strong reflections ($I>4\sigma$) short-circuit to $|F|=\sqrt{I}$, where the French–Wilson bias is negligible; a reflection with an unusable $I/\sigma$ falls back to $\sqrt{\max(I,0)}$. The integral is evaluated numerically with a log-shift for stability. + +Amplitudes are written as MTZ `F`/`SIGF`, mmCIF `_refln.F_meas_au`/`F_meas_sigma_au`, and appended to the text HKL, alongside the intensity columns. The **same** $|F|$ feed the model-validation step (§14), so the reflection file and the maps use one consistent set of amplitudes. + --- ## 11. Mosaicity and “profile radius” monitoring @@ -613,4 +646,40 @@ A linear regression of $\log\langle I\rangle$ vs $1/d^2$ provides an estimate of - **Twinning check.** A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ (taken per resolution shell with noise-only shells skipped and Wilson outliers rejected, so a single strong reflection in a collapsed-mean shell cannot skew it) are written to the merged mmCIF as a twinning diagnostic. Twinning is only flagged in Laue classes where a merohedral twin law can exist; the holohedral high-symmetry classes ($4/mmm$, $6/mmm$, $m\bar{3}m$, and $\bar{3}m$ on a rhombohedral lattice) are exempt, so a low $\langle|L|\rangle$ there is reported as a statistical artefact rather than twinning. - **Outlier rejection.** Merging applies an optional per-observation median-based $N\sigma$ cut (default 6σ for `rot3d`) and an optional per-crystal $\Delta\mathrm{CC}_{1/2}$ image rejection (`--reject-delta-cchalf`, CrystFEL-style, off by default). The same $N\sigma$ cut is fed back into the error model: after an initial $a,b$ fit the parameters are re-fit once on the reflections that survive rejection (dropping any whose squared deviation exceeds $N\sigma^2\,[a\,\sigma^2 + (b\,\langle I\rangle)^2]$), so the calibrated errors describe the reflections that actually enter the merge rather than the pre-rejection pool. - **Automatic resolution cutoff.** By default the reported/written high-resolution limit is trimmed where $\mathrm{CC}_{1/2}$ falls off (logistic, target 0.30); `--scaling-high-resolution` overrides it and `--resolution-cutoff off` disables it. -- **Intensities only.** The merged output carries intensities (mmCIF `intensity_meas`, MTZ `IMEAN`/`SIGIMEAN`); it does not convert to amplitudes $|F|$ (no French–Wilson / truncate step) — do that downstream. +- **Amplitudes and intensities.** The merged output carries both intensities (mmCIF `intensity_meas`, MTZ `IMEAN`/`SIGIMEAN`) and French–Wilson amplitudes (mmCIF `F_meas_au`, MTZ `F`/`SIGF`; §10.8), so a downstream program can refine against either. + +--- + +## 14. Model-based validation: R-free against a model and electron-density maps + +Offline (`rugnux --model model.pdb`) the merged data can be scored against a supplied atomic model and **initial** electron-density maps computed — enough to confirm that a model fits the data and to inspect the density, not a substitute for refinement. **The structure itself is not refined**; the model is only re-fractionalized into the data unit cell (a rigid cell adjustment, so a deposited model with a slightly different cell still lines up), and the observed amplitudes are the French–Wilson $|F|$ from §10.8, so the R-free and the maps use exactly the same amplitudes as the written reflection file. The model, structure-factor, bulk-solvent and FFT machinery is provided by GEMMI. + +### 14.1 Model structure factors + +The model electron density is sampled on a grid (IT92 X-ray form factors, with a Refmac-compatible Gaussian blur chosen for the grid spacing) and Fourier-transformed to structure factors $F_\mathrm{calc}(hkl)$ up to the data resolution. + +### 14.2 Bulk solvent and scaling + +A flat bulk-solvent mask around the model is transformed to $F_\mathrm{mask}$, and the model is scaled to the observed amplitudes in two stages: + +1. an overall least-squares fit of a scale $k$, an anisotropic $B$, and the flat-solvent parameters $k_\mathrm{sol}, B_\mathrm{sol}$, giving + + $ + F_\mathrm{model} = k\,e^{-\mathbf{h}^\top \mathbf{B}\,\mathbf{h}/4}\left(F_\mathrm{calc} + k_\mathrm{sol}\,e^{-B_\mathrm{sol}\,s^2}\,F_\mathrm{mask}\right),\quad s^2 = 1/4d^2; + $ + +2. a smooth **per-resolution-shell scale** $K(1/d^2)$ (least squares $\sum F_o F_c/\sum F_c^2$ per shell, fit on the work reflections only) applied on top, which absorbs the residual radial $F_o/F_\mathrm{model}$ mismatch a single overall $B$ cannot. + +### 14.3 R-work and R-free + +Crystallographic R-factors are reported over the work and free sets (the §10.7 flags): + +$ +R = \frac{\sum \big|\,|F_o| - |F_\mathrm{model}|\,\big|}{\sum |F_o|}, +$ + +with R-free the same sum restricted to the free set — an unbiased measure of how well the model explains data it was not scaled against. + +### 14.4 Electron-density maps + +Two maps are formed with the model phases $\varphi_\mathrm{model}$: a $2F_o-F_c$ map, coefficients $(2|F_o|-|F_\mathrm{model}|)\,e^{i\varphi_\mathrm{model}}$, and an $F_o-F_c$ difference map, $(|F_o|-|F_\mathrm{model}|)\,e^{i\varphi_\mathrm{model}}$, each inverse-Fourier-transformed to a real-space CCP4 map (`_2fofc.ccp4`, `_fofc.ccp4`). A map-coefficient MTZ (`_maps.mtz`: `FP`, `FC`, `PHIC`, `FWT`/`PHWT`, `DELFWT`/`PHDELWT`, `FREE`) is written alongside so the maps can be reopened or rebuilt in Coot / PyMOL. These are unweighted difference coefficients (no $\sigma_A$ / figure-of-merit weighting), which is why they are described as *initial* maps.