EDSTATS

NAME

SYNOPSIS

DESCRIPTION

• The real-space R factor (RSR) is defined (Brändén & Jones, 1990; Jones et al., 1991) as:
  RSR = Σ |ρ_obs - ρ_calc| / Σ |ρ_obs + ρ_calc|

• The real-space correlation coefficient (RSCC) is defined as:
  RSCC = cov(ρ_obs,ρ_calc) / sqrt(var(ρ_obs) var(ρ_calc))

  where cov(.,.) and var(.) are the sample covariance and variance (i.e. calculated with respect to the sample means of ρ_obs and ρ_calc).

• EDSTATS computes two real-space correlation coefficients: the 'sample' correlation coefficient defined above, and the 'population' correlation coefficient, i.e. with respect to the population (overall) means, which will be zero if the F(000) terms were not included in the map calculation (OVERLAPMAP uses only the sample means).  The RSCC based on the population means seems to be better at detecting weak correlations.

• The real-space Z_obs metric (RSZO) is defined (Tickle, 2011) as:
  RSZO = mean(ρ_obs) / σ(Δρ)

  where σ(Δρ) is the standard uncertainty of the difference Fourier map.  Note that this is the standard uncertainty of the 'Fo-Fc' map, NOT the RMS value of the '2Fo-Fc' map, which bears no relationship whatsoever to the uncertainty!

• The real-space Z_diff metrics (RSZD- and RSZD+) are defined (Tickle, 2011) as follows for the sets of negative and positive values respectively of Δρ at the grid points that are covered by the group of main- or side-chain atoms under consideration:

  1. Order the values in each set in increasing numerical value (i.e. ignoring the sign).

  2. For each of N subsets of size 1, 2, ..., N-1, N of the numerically highest values of the original set of size N, compute the cumulative probability of chi-square (χ² = Σ (Δρ/σ(Δρ))²) for the subset.  So the subset of size 1 is simply the numerically highest value ('maximum order statistic') in the original set, the subset of size 2 consists of the 2 highest values of the set, the subset of size N-1 excludes the lowest value, and the subset of size N is just the set itself.

  3. In practice this χ² cumulative probability is very difficult to compute (even by stochastic numerical integration) for subsets other than those of size 1 and N (it involves integrals up to dimension N where N may be anything from 10 to 1000).  Note that the standard χ² cumulative probability assumes that the sample is selected randomly, whereas here we are selecting the highest values.  Therefore we approximate it as the product of two components: the standard cumulative probability of χ² for a randomly selected subset, and a correction, the Dunn-Šidák correction (Sokal & Rohlf, 1995; Gibbons & Chakraborti, 2003), in this case the cumulative probability of the order statistic, for the fact that we are selecting the highest values.

  4. Take the highest cumulative probability over all subsets, and convert this to the corresponding normal Z-score, making the Z-score negative for the set of negative values; this is the final RSZD- or RSZD+ score.  The program also computes a combined RSZD score which is simply the maximum of |RSZD-| and RSZD+.

• The real-space Z-scores RSZO, RSZD- and RSZD+ require estimates of the standard uncertainty σ(Δρ) and offset of the 'Fo-Fc' map (the offset arises from omission of the F(000) term, which may differ from zero since the model is not necessarily complete).  The recommended procedure is to use as an initial estimate the value of the σ(Δρ) in the map header, with zero as the offset, and then rescale σ(Δρ) and the offset separately for each chain and the bulk solvent.  Bulk solvent is assigned the chain ID '%' for this purpose and ordered waters are considered to belong to a chain with ID '0' whatever their actual chain IDs in the PDB file.

• The sample size correction above arises because the greater the sample size the more likely it is that high values will occur purely by chance.  This correction takes into account the fact that the number of grid points is not the same for all residues, because obviously different residue types contain different numbers of atoms, and also different limiting atom radii will enclose different numbers of grid points, because the radius varies with atom type and B_iso.  The correction is therefore necessary to make the metrics comparable between different residues and to be able to apply a common threshold to the metric for all residues.  Note that the RSR and RSCC metrics do not apply a sample size correction: it is assumed that all sample points contribute equally to the metrics independent of the sample size.

• The number of grid points referred to above is the number of statistically independent grid points covering the atoms; this is the actual number of grid points with an over-sampling correction factor.  According to the Nyquist-Shannon sampling theorem, the grid spacing required for statistical independence is 1/2 the high resolution cut-off (d_min), so e.g. if a grid spacing of d_min/4 is used then the effective number of grid points is the actual number / 2³.

• The advantage of the real-space Z-scores over the real-space R factor and correlation coefficient scores (including the 'population' CC metric) is that the former depend purely on model accuracy (RSZD) or model precision (RSZO), whereas RSR and RSCC depend on both (e.g. it's obvious from the plots that RSR and RSCC are at least partially correlated with the atomic B_isos); this means that it's impossible to say how much of the observed effect on the metric is due to lack of accuracy and how much to lack of precision.

  Note that model accuracy is related to the likelihood of the model (i.e. the consistency of the model with the data), and is what is improved by model building and refinement.  The difference Fourier density is obviously a measure of any discrepancy between the model and the data, so is a direct measure of model accuracy.

  Model precision is a property of the crystal and the data (assuming the refinement is done optimally), and is related to data quality and completeness, resolution, atom type (or atomic scattering factor), occupancy and atomic B_iso; hence model precision can only be improved by crystallizing in a different crystal form and/or collecting better (e.g. more precise and/or higher resolution) data.  The ρ_obs density, divided by its standard uncertainty (note: this is not the same as RMS(ρ_obs)), is a measure of model precision which incorporates all the above factors correlated with precision (e.g. the atomic B_iso is also a precision metric but it doesn't take account of the variation of precision with atom type and occupancy).

• The sums and min/max functions required to compute all residue or atom metrics are taken over all map grid points within a specified distance of each atom centre.  This distance limit is naturally a function of the atom type (via the atomic scattering factors computed from the 5-Gaussian approximation table in $CLIBD/atomsf.lib), the atomic B_iso values and the resolution limits, as shown in the following table of the distance limit r_max for an O atom.  Values used by SFALL are also shown for comparison: note that the latter depend only on B_iso and are independent of atom type and resolution:

  B/Å2:                   10    20    30    40    50    60    70    80    90
  dmin/Å
        rmax/Å (SFALL: all atoms)
   All                   2.35  2.67  2.95  3.21  3.45  3.67  3.88  4.08  4.27
  
        rmax/Å (EDSTATS: O atom)
   3.5                   1.72  1.78  1.83  1.89  1.95  2.02  2.08  2.15  2.22
   3.0                   1.51  1.58  1.65  1.72  1.80  1.88  1.97  2.06  2.14
   2.5                   1.31  1.39  1.49  1.59  1.70  1.80  1.91  2.02  2.12
   2.0                   1.12  1.24  1.38  1.52  1.66  1.79  1.91  2.02  2.13
   1.5                   0.96  1.16  1.35  1.52  1.66  1.79  1.91  2.02  2.13
  <=1.0                  0.91  1.16  1.35  1.52  1.66  1.79  1.91  2.02  2.13
  

  Note that the limiting high-resolution values of r_max are attained at ~ d_min = 1.5Å.

• The resolution-dependent distance limit is computed by first performing an analytical truncated Fourier transform of the atomic scattering factor f(s) to obtain the equation for the calculated electron density ρ(r) for data between specified resolution cut-offs, at distance r from the atom centre:

                                 smax
          ρ(r) = FT(f(s)) = (8/r) ∫ f(s) exp(-Bs2) sin(4πrs) s ds
                                 smin
  

  for specified limits s_min and s_max of sin(θ)/λ.

  Then the ratio of the radius integral of ρ(r) integrated out to the outer limit r_max relative to the radius integral integrated to infinite distance is:

                                 rmax          ∞
          Radius integral ratio = ∫ ρ(r) dr / ∫ ρ(r) dr
                                  0           0
  

  and this equation solved to obtain r_max for a radius integral ratio = 0.95 (i.e. 95% of the integral lies within distance r_max of the atom centre).  The integrals with respect to r can be obtained analytically; the integrals with respect to s in general have no analytical solution and must be computed numerically (using e.g. the QUADPACK library).  Note that ideally the volume integral of ρ(r):

                           rmax           rmax
          Volume integral = ∫ ρ(r) dV = 4π ∫ ρ(r) r2 dr
                            0              0
  

  should be used, but unfortunately this integral does not converge.

• For the RSZO metric EDSTATS uses the ρ_obs map with Fourier coefficients 2mF_o-DF_c for acentric reflections or mF_o for centrics (Main, 1979; Read, 1986); for the RSZD metric it uses the Δρ (difference Fourier) map with Fourier coefficients 2(mF_o-DF_c) for acentrics, or mF_o-DF_c for centrics.  For the RSR and RSCC metrics it uses the ρ_obs and ρ_calc maps. 

  However, for the latter, because we cannot rely on the correct Fourier coefficient for ρ_calc being present in the file of map coefficients, it is necessary to obtain it as the difference between the ρ_obs and Δρ coefficients.  Since we have:

          Δρ = ρobs - ρcalc
  

  or:

          ρcalc = ρobs - Δρ
  

  therefore for acentrics:

          ρcalc = F(2mFo-DFc) - F(2(mFo-DFc))
                = F(DFc)
  

  whereas for centrics:

          ρcalc = F(mFo) - F(mFo-DFc)
                = F(DFc)
  

  Hence the correct Fourier coefficient for ρ_calc is DF_c for all reflections.  Note that it is frequently stated that the coefficient for acentrics is mF_o-DF_c but if this were used it would give completely the wrong result for the ρ_calc coefficient (it would give mF_o !).

• EDSTATS also has options to output data for the histogram, 'P-P difference' and 'Q-Q difference' plots of the difference Fourier and observed Fourier maps.  Note that the 'P-P difference' and 'Q-Q difference' plots are functionally identical to the standard 'P-P plot ' (probability-probability) and 'Q-Q plot ' (quantile-quantile: 'quantile' is just another name for 'normalised deviate' or 'Z-score').  The distinction is purely one of presentation: whereas the standard 'P-P' or 'Q-Q' plot plots x vs. y, where x and y are respectively the normal expected and observed probabilitlies or quantiles, the 'P-P difference' or 'Q-Q difference' plot plots x vs. y-x.

The data file pdb-edstats.out, or a link to it, must be present in the current directory; alternatively set the environment variable PDB_EDSTATS to point to it.  The data were obtained by running edstats on ~ 600 supposedly 'good' structures (anonymous!) from PDB_REDO with Rfree < 0.175 and > 100 residues (protein only).  This is not ideal, since it would clearly be much better to bin the known structures by high resolution cut-off and compare your structure only with known structures at roughly the same resolution; however this will require a much larger database than I have the resources to set up in the short term.  Hopefully this feature will be developed and improved in a future release.

The columns in the data file pdb-edstats.out contain:

1. High resolution cut-off.

2. Resolution-weighted average B_iso (e.g. the effective average of B_isos 10 and 100 is not 55 but something much closer to 10, depending on the resolution cut-offs, since the atom with B_iso = 100 only contributes significantly to the scattering at low resolution).

3. Q-Q plot ZD- metric: this gives an overall indication of how much the distribution of all negative difference density in the asymmetric unit deviates from the expected normal distribution for purely random errors.  Significant negative density outliers giving a high numerical Q-Q plot ZD- metric probably indicates wrongly placed atoms, over-restrained B factors, problems with the bulk solvent parameters (e.g. due to low completeness at low resolution), or generally low data completeness.

4. Q-Q plot ZD+ metric: ditto for all positive difference density.  Low per-cent ranks (large values) for this metric are not as indicative of problems as are low ranks for the Q-Q plot ZD- metric above, because it can be difficult to interpret residual density outliers due to disorder, buffer ions, cryo-protectants and other additives, so uninterpreted (and uninterpretable) density is quite common in deposited structures.  Consequently a high value for the Q-Q plot ZD+ metric does not necessarily indicate a serious problem; it would be better to check the per-residue RSZD+ scores.

5. Percentage of residue RSZD- metrics numerically above the 3σ threshold (see also next section).

6. Percentage of residue RSZD+ metrics above the 3σ threshold.

The percent-rank.pl script prints out the per-cent ranks for metrics 3-6 above.

Examples of usage:

percent-rank.pl edstats.log
or
percent-rank.pl *.log

The RSR and RSCC (both 'sample' and 'population') metrics are tabulated for comparison but are correlated with both accuracy and precision, so they can be useful in some circumstances, but they don't always help with telling you whether adjustment of the model is required, or whether the problem is actually an intrinsic property of the structure, or lies with the data.  Note that the RSR and RSCC metrics vary with the program used, since they depend strongly on the radius cut-off, scaling algorithm and other variables which can vary a lot between programs.

REFERENCES

#!/bin/tcsh
# Fix up the map coefficients: FLABEL specifies the label for Fobs &
# σ(Fobs) (defaults are F/SIGF or FOSC/SIGFOSC).  Here, 'in.mtz'
# is the output reflection file from the refinement program in MTZ
# format.

rm -f fixed.mtz
mtzfix  FLABEL FP  HKLIN in.mtz  HKLOUT fixed.mtz  >mtzfix.log
if($?) exit $?

# Good idea to check the mtzfix output before proceeding!

less mtzfix.log

# If no fix-up was needed, use the original file.

if(! -e fixed.mtz)  ln -s  in.mtz fixed.mtz

# Compute the 2mFo-DFc map; you need to specify the correct labels for
# the F and phi columns: 'FWT' & 'PHWT' should work for Refmac.
# Note that EDSTATS needs only 1 asymmetric unit (but will also work
# with more).  Grid sampling must be at least 4.

echo 'labi F1=FWT PHI=PHWT\nxyzl asu\ngrid samp 4.5'  | fft  \
HKLIN fixed.mtz  MAPOUT fo.map
if($?) exit $?

# Compute the 2(mFo-DFc) map; again you need to specify the right
# labels.

echo 'labi F1=DELFWT PHI=PHDELWT\nxyzl asu\ngrid samp 4.5'  | fft  \
HKLIN fixed.mtz  MAPOUT df.map
if($?) exit $?

#!/bin/tcsh
# Q-Q difference plot & main- & side-chain residue statistics.

echo resl=50,resh=2.1  | edstats  XYZIN in.pdb  MAPIN1 fo.map  \
MAPIN2 df.map  QQDOUT q-q.out  OUT stats.out
if($?) exit $?

#!/bin/tcsh
# Main- & side-chain atom statistics, using chains A & I only & writing
# PDB file with per-atom Zdiff metrics.
echo mole=AI,resl=50,resh=2.1,main=atom,side=atom  | edstats  \
XYZIN in.pdb  MAPIN1 fo.map  MAPIN2 df.map  XYZOUT out.pdb  \
OUT stats.out
if($?) exit $?