MolProbity: More and better reference data for improved all‐atom structure validation

Improving the worldwide database

The primary aim of structure validation, as we see it, is not just to identify mistakes, but rather to help fix them. We calibrate our success at that goal by annually tracking the “clashscore” (all‐atom steric clash overlaps ≥0.4Å per thousand atoms) for new worldwide PDB depositions in the resolution range 1.8‐2.2Å. All‐atom clashes are an especially sensitive and powerful indicator of local fitting problems, and still are provided only by MolProbity. Before the advent of MolProbity in 2002, clashscores were constant, but since then they have steadily improved, now by about a factor of 3, as shown in Figure 1(A). The change seems to be leveling off, as it must, since it cannot go below zero. In fact, this measure is sometimes interpreted too stringently: the goal is few clashes, as in the best reference data, not zero clashes.16, 17

The wwPDB has adopted four MolProbity criteria (clashscore, Ramachandran, rotamer, and RNA backbone) along with other validation types, at deposition, as PDF reports for reviewers, and as visual “sliders” showing relative percentile scores on each structure summary page at the RCSB, PDBe, and PDBj websites [Fig. 1(B)]. They are now used for X‐ray crystallography18, 19 and for nuclear magnetic resonance,20 and may also be adopted for cryo electron microscopy. Complete MolProbity validation and its ongoing updates are integrated into the Phenix structural biology software suite, in the graphical user interface as well as for automated procedures and command‐line use.21, 22, 23

CCTBX and python

The full range of MolProbity validation has now been incorporated into the Phenix crystallography suite21 to provide frequent user feedback, and specific MolProbity criteria are also used directly within the automated Phenix workflows. To accomplish this, MolProbity's mid‐level Java utilities have been re‐implemented in Python and use the same open‐source CCTBX (Computational Crystallography Toolbox) codebase24 that underlies Phenix. For Phenix users, MolProbity validation is accessible through the Phenix GUI, with real‐time links to outlier locations for fixup in Coot, or through the command‐line as phenix.molprobity. The individual components of MolProbity, including Reduce, Probe, Clashscore, Ramalyze, Rotalyze, Omegalyze, CaBLAM, Cbetadev, and Suitename are also separately available through the Phenix command‐line. The cctbx_project portion of Phenix is open source and freely available at https://github.com/cctbx/cctbx_project.

The specific validation functions on the MolProbity website call the same Python utilities and other CCTBX functions, so there is now only a single codebase for this part of the validation. Building MolProbity on the CCTBX project assures that validation on the MolProbity webserver and within Phenix stay synchronized. It also takes advantage of the robust support that the CCTBX project offers for access to the diffraction data and for evolving file types like mmCIF and other new code requirements.

Manchester mirror

Since 2016 there has been an identical MolProbity mirror website at the University of Manchester, UK (http://molprobity.manchester.ac.uk), with a link from the main site at Duke. It is hosted by Simon Lovell, who worked on early development of the validation criteria.9, 10, 11 This provides redundancy for rare downtimes and a closer site for users in the UK or Europe.

GitHub

GitHub now serves MolProbity both for version control and for distribution. The MolProbity source code and its key dependencies are all open source and freely available from the Richardson Lab GitHub repositories at https://github.com/rlabduke. MolProbity, Reduce, Probe, Suitename, and KiNG (as javadev) are each available as their own repositories. The reference_data repository contains the Top8000 dataset versions and the Ramachandran, rotamer, and CaBLAM distributions derived from it. This availability is used for Phenix nightly builds, and it allows users to install their own local MolProbity server, if they have limited internet access, confidentiality concerns, or a need to script many MolProbity runs using the command‐line tools; CCTBX is archived on GitHub, also open source. The GitHub pull request interface has also allowed community members outside the lab to submit suggested code improvements, which we have implemented.

Top8000 and RNA11 datasets

Many MolProbity validations rely on statistical expectations for macromolecular structure. These expectations are drawn from data of high‐quality residues in high‐resolution protein and nucleic acid structures. As the number of depositions in the PDB has grown, we have increased both the size and the accuracy of our reference datasets to better capture the depth and diversity of real molecular structure.

The latest iterations of our reference datasets are the Top8000 for proteins and RNA11 for RNA (see Methods for details of their construction). The “standard” Top8000 used for Ramachandran and CaBLAM is filtered at the 70% homology level and contains 7957 protein chains, up by an order of magnitude from the previous Top500.11 RNA11 contains 311 RNA structures (including proteins if present, to allow study of the interactions), up from 171 in the previous RNA05 dataset used in defining the community‐consensus RNA suite backbone conformers,25 with classification by content and function.

Chain‐level filtering is an important first step in preparing a high‐quality dataset. However, even high‐resolution structures almost always contain regions of local disorder. To guard against inclusion of the resulting local modeling errors, residue‐level filtering is a necessity. In preparation of Ramachandran and CaBLAM contours, we excluded residues having any mainchain atom with a B‐factor > 30. This simple filter worked quite well at excluding poor conformations and producing quite clean, reproducible contours. The enforced deposition of structure‐factor data along with each new entry in the worldwide PDB lets us use local electron‐density criteria in addition to B‐factors, for a more accurate and complete residue‐level filter (see Methods for specific details). All future validations and revisions will take advantage of this improved residue‐level filtering, applied across the atom types appropriate for backbone or for sidechain criteria.

Methyl orientations

Rotation of all methyl groups was originally written into Reduce but was deprecated almost immediately, partly because of the computational expense but mainly because steric clashes of methyl H atoms were nearly always caused by misplaced parent C atoms, so that curing them by methyl rotation was not the correct approach. Most methyls rotate no more than about 15° off stagger, but terminal Met methyls (with a longer S‐C bond) can rotate as much as 40–50°, so for many years we allowed rotation of Met methyls. However, in a survey of sub‐1Å resolution crystal structures we discovered that even for those structures such freedom produced the wrong answer more often than the right one, so we now keep Met methyls staggered also.

Methyl groups attached to planar aromatic rings (on groups such as hemes, thymines, or modified bases) are an interestingly different case, since the joining of a planar sp2 to a tetrahedral sp3 atom produces a flatter energy profile for rotation. Instead of a 3‐fold stagger, there is a preference for one of the two conformations with one H perpendicular to the ring plane and the other two 30° from it, confirmed by a survey of very high‐resolution H difference peaks. Figure 2 shows an example methyl on a heme at 0.88Å, illustrating clear positive difference peaks (blue) for all three H atoms in one of those preferred orientations, the incorrect positioning of our previous default (white bonds, off by 30°), and the resulting small, incorrect hydrogen clash (hotpink spikes). The new system solves such problems by trying out the two preferred orientations and choosing the better of the two (green lines) by all‐atom contact criteria, using a procedure analogous to the determination of Asn/Gln/His flips but not requiring an extended interaction network.

All‐atom contacts and OneDotEach

The standard MolProbity all‐atom contact representation,8 both for scoring and for visualization, consists of paired patches of dots (or overlap spikes) on the surface of each atom of the contacting pair. All hydrogen atoms are included. The intuition here is that atoms can be treated as interacting at the effective surfaces of their electron clouds (approximated by van der Waals radii if nuclear positions are being used). Thus contact is characterized by surface‐to‐surface interaction, as distinct from overlap volume or from pairwise center‐to‐center distance. That surface interaction is attractive until overlap, then repulsive unless defined to be a hydrogen bond. Since surface‐to‐surface approach defines interaction, a third intervening atom occludes pairwise interaction of the original two. Scalings are tuned by overall score match with observed occurrence in very well built experimental crystallographic models.

The many‐dots Probe protocol starts by placing at each atom position a sphere of surface dots approximately uniformly distributed (16 dots per Ų by default). The scoring algorithm uses nested loops measuring unoccluded line‐of‐sight distances ≤0.5Å (by default) between each dot on a source atom and all dots of neighboring, not‐covalently‐bonded target atoms. The result is better than pairwise‐from‐centers algorithms at giving an intuitive understanding of packing inside and between molecules, and is more powerful for diagnosing and guiding corrections of fitting errors in the model.

However, this calculation both produces verbose output and is inherently slow, suitable for one‐time evaluations but not usable in a fast computational inner loop. As a first step in correcting that problem, Probe now includes an optional calculation called OneDotEach. It solves the verbosity and helps the speed problem without losing the non‐pairwise aspect, but it does not yet incorporate the smooth derivatives necessary for dynamic calculations such as refinement.

OneDotEach starts from the simple center‐to‐center vector between each non‐covalently‐bonded pair of atoms tested (those closer than van der Waals contact plus 0.5Å, by default). It then tests whether that vector passes within any other nearby atom's van der Waals radius of its center: this is our definition of a third atom occluding a contact. If so, the potential contact pair is rejected. Of course, in compensation there will be an all‐atom contact between the occluding atom and one of the original pair. The logic was simplified by the realization that, except in a direly bad‐geometry model, the only type of third atom that can occlude a pair within 0.5Å of each other is a hydrogen covalently bonded to one of the pair atoms. The output of OneDotEach is a single dot for each atom in an accepted pair, at the position where the interatomic vector intersects its surface. The main application of OneDotEach so far has been in calculating the distributions of nearest‐neighbor distances used in defining new van der Waals radii for the H parameter update (see below).

Better‐idealized output coordinates from NQH flips

In order to perform all‐atom contact analysis, MolProbity uses the C++ program Reduce to add and optimize hydrogen atoms, including analysis of each complete H‐bond network, rotation of OH, SH, NH₃ (but not methyls), and consideration of first‐layer waters.9 As part of that process, Reduce optionally performs automated “flips” on Asn, Gln, and His residues to correct a common error where one of these sidechains (near‐symmetric in electron density) is modeled with its terminal χ angle 180° from the correct position. Such an error usually produces a pattern of clashes and missing H‐bonds that can be recognized automatically. When Reduce recognizes such a pattern, it will recommend a “flip” of that sidechain. The automated assignments are very reliable, but 3D “flipkin” kinemages show views with animation to compare the pattern of contacts between alternatives, allowing the user to see the evidence behind each potential flip and make an informed personal decision if they wish.

The NQH flip was previously performed by simply renaming the terminal N and O atoms of Asn and Gln or the C and N ring atoms of His, without changing any of the modeled coordinates. This method has the important advantage of preserving the fit‐to‐density of those atoms perfectly, and of showing the two competing sets of all‐atom contacts most correctly. Therefore the scores, flip decisions, and flipkins are still produced this way. However, the bond geometry around the His ring and between the C‐O and C‐N bonds is not quite symmetrical. Performing the flip by renaming atoms generated aberrant bond lengths and angles around the renamed atoms, by up to 6σ (0.1Å and 5° for Gln). This geometry would be easily corrected by another round of refinement after performing flips, but enough users have deposited structures without refining again that the PDB now contains a statistically significant population of aberrant bonds resulting from Reduce. To protect the integrity of the database, we developed a new method for producing the output coordinates from NQH flips.

The goal of the new method is to closely maintain the optimized fit of the relevant sidechain atoms in their electron density (and thus also their contacts), but without distorting covalent geometry. The intuitive, simple 180° rotation is not acceptable, because asymmetry of the sidechain head (the His ring or the terminal C/N/O of Asn and Gln) changes atom positions significantly. We chose a three‐step docking procedure. First, the sidechain head is rotated 180°, as intuition dictates. Second, the head group is hinged back into the plane of the head group in the original model, compensating for cases in which the head group is not in plane with its stem [Fig. 3(A,B)]. Third and finally, the whole sidechain is three‐point docked as a rigid body [Fig. 3(B,C)]. The Cα position is held constant and the two terminal H‐bonding atoms in the sidechain are docked as closely as possible onto their switched‐identity original positions.

This new flip method still affects sidechain geometry, but much less so. It changes the bond angles between the mainchain and the sidechain at the Cα, but usually by less than one σ. No bond lengths are affected. Users should bear in mind that the output coordinates will now not precisely match what is seen in the flipkins. The new output coordinates correct the database pollution, but should serve as a reminder that all local structure corrections can generate geometry problems where they rejoin the rest of the model, and refinement after correction is always necessary to settle these details.

Re‐definition of H‐atom parameter sets

Explicit hydrogen atoms have increasingly become an important part of both experimental and computational methods for structural biology. In contrast to the accurately determined geometrical parameters for heavier atoms, the various bond‐length and van der Waals parameter sets in current use for hydrogen were derived many decades ago from limited data. Their specific values differ by as much as 20% between libraries, including within our now‐integrated system of MolProbity and Phenix. That is a big discrepancy for an effective bond length, and it can fairly often mean the difference between an acceptable steric contact and a serious steric clash of the van der Waals spheres (>0.4Å overlap counts as serious in MolProbity). Around that threshold of 0.4Å, we aim to optimize the balance between diagnosing serious conformational errors and raising false alarms. Packing analysis and validation both depend on the total system of hydrogen x‐H distances, vdW radii, and the 0.4Å threshold defined for clashes. Although there may not be a single right answer at all resolution levels, several factors convinced us that our previous system was slightly too strict – primarily because even for the best structures, clashscore bottomed out at about 5 clashes per 1000 atoms rather than at zero.

For the above reasons, we set out to provide the scientific basis for two updated sets of H parameters, one specific for the electron‐cloud‐center positions suitable for x‐ray crystallography and one for the nuclear positions used in most other methods, and as accurate as currently feasible across bonded‐atom types and geometries within each set. Many information sources were utilized, including a search of the older literature, spherical‐patch fitting to quantum‐calculated (QM) electron‐density contours, small‐molecule neutron and x‐ray coordinates from the Cambridge Structural Database (CSD26), H difference peaks in <1Å protein crystal structures from the Protein Data Bank and, most decisive of all, a combination of H atom coordinates and electron‐density difference peaks in small‐molecule x‐ray structures from the Crystallography Open Database (COD27). New H van der Waals radii were tuned for use with the new H positions, by pairwise nearest‐neighbor atom‐atom distance distributions in the Top8000 quality‐filtered dataset of protein chains. The extensive methodological details are described in the Supplementary Information. Along the way, several interesting categories of H atom contacts were clarified. Very short carboxyl O‐O pair H‐bonds form the short side of a cleanly bimodal distribution and are unequivocally real, even outside the transition states where they are usually studied;28, 29 however, they require a narrowly specific relative geometry.30 Shortened CH…O interactions are fairly common in β sheet, but for good reference data they show distances within the tail of the overall distribution and no preference for biologically functional sites.

Historically, Reduce used nuclear hydrogen positions because it was originally developed for packing calculations rather than for crystallography, and because the best‐documented set of van der Waals radii was tuned to those positions.31, 32 However, for the major current uses that is the wrong choice, because the electron cloud is what diffracts X‐rays and its outer region is where the atoms actually interact – repelling when separate electron clouds overlap.33 Phenix, appropriately for X‐ray refinement, has placed hydrogens at the shorter distances of the electron‐cloud centroids. The Phenix x‐H distance values were adopted from ShelX34 by way of the CCP4 monomer library,35 but the original procedures were not explicitly documented and there was no value for S‐H. Therefore, the work described here has combined evidence from a variety of sources in order to define optimal contemporary sets of electron‐cloud and of nuclear x‐H bond lengths and corresponding van der Waals parameters. We also re‐examined the nuclear positions for H and D, finding no significant difference, and we carefully proofread large libraries to correct the inevitable few typographical errors. A number of software modifications were made, both for accomplishing the underlying research and for implementing its results.

Heavier atoms such as carbon have electron‐cloud center and nuclear position essentially coincident, but for hydrogens it has long been known that they differ quite significantly.36, 37, 38 The single electron of a hydrogen atom must provide its share of electron occupancy in the covalent bond to its parent heavier atom, which systematically contracts the electron cloud inward from the nuclear position. In addition, the H electron cloud can be shifted sidewise by steric bumps or H‐bonding with surrounding atoms.39 As an example of both these effects, the Hɛ1 hydrogen on Trp 37 of PDB 1yk4 at 0.69Å resolution in Figure 4 was placed by Reduce at the nuclear position and planar to the ring. However, the clear Fo‐Fc difference peak (blue) for Hɛ1 is shifted both inward toward the ring N and also upward toward a line to the H‐bond acceptor, a carboxyl oxygen. This work aims to correctly reflect the systematic shifts along the covalent bond direction, but does not attempt to model second‐order shifts due to local environment. The following sections summarize the database and literature results and the overall decisions, organized by parent‐atom element type.

Upload: interpretation, conversions, and cleaning

The first step in analyzing a structure on the MolProbity webservice is either to upload a coordinate file from your computer or to fetch one from the PDB or other database. About 85% of MolProbity use is by upload (presumably structural biologists analyzing their own models during the process of structure solution) and 15% by fetch (presumably biological, biomedical, bioinformatic, educational, and other end‐users).

Input coordinate files are parsed to provide user feedback about their interpreted contents, such as number of chains, alternate conformations, presence and type of hydrogens, etc. The user should check that the displayed interpretation seems correct, since it is easy to confuse O vs 0 or l vs 1 in a PDB code and fetch the wrong file. Some format problems now generate error messages, such as old PDB format (pre‐v3.0, now rare) which is converted, mispaired MODEL and ENDMDL records, or duplicate atoms (usually a result of missing or inconsistent chain or segid fields), which will fail in CCTBX. Files submitted in mmCIF format are now automatically converted to PDB format on upload or fetch, with the hybrid36 system used for files too large for standard PDB format. Hybrid36 uses 2‐character chain‐ids and numbers atoms normally through atom 99999, then uses a combination of letters and numerals starting with A0000; the equivalent operation for residue number starts after 9999. The hybrid36 format is supported by Phenix and by all MolProbity validations. Due to heightened security issues, we also now check for content that seems to be executable code or not to be an interpretable coordinate file.

Hydrogen addition and NQH flips

On the main page there are edit options, such as dropping extra chain copies, and file‐choice options, such as choosing to keep the H atoms of your own input file rather than have Reduce optimize them. However, the first validation step is nearly always H addition, absolutely necessary to take advantage of the all‐atom contact analysis that is MolProbity's most unique and powerful feature. Roughly half the atoms in a given macromolecule are hydrogens, and the vast majority of interatomic contacts in a macromolecule are between H atoms. However, due to their single, weakly scattering electron, hydrogens are rarely visible except in difference density at ultra‐high resolution. As a result, hydrogens have seldom been modeled historically, and MolProbity must add and optimize H atoms for most structures before full validation can be performed. It uses Reduce for that task. The default is now to place those hydrogens at the electron‐cloud center positions suitable for x‐ray crystallography (see above), but the user can choose to use nuclear positions. In either case, they are optimized across complete H‐bond networks, including optimization of Asn/Gln/His (NQH) “flips” unless turned off by the user. The NQH flips are very reliable and the easiest type of automatic structure correction, so we recommend their use each time your structure has changed substantially. However, they are not perfect, and if you are already near completion and have decided to reject some flips, then by all means turn the process off. Note that you must download the resulting file to take advantage of the improved flips, and that their coordinates will now be better idealized than in previous versions of MolProbity (see above).

Better treatment of ensembles

We have improved the functionality of MolProbity for multi‐model ensemble PDB files. NMR structures are often deposited as ensembles of models, with potentially over a hundred models within a single PDB file. Additionally, ensemble crystallography is becoming more common, as a way of expressing either the modeling uncertainty or the likely conformational heterogeneity within a crystal. Previously, if a user uploaded an ensemble file to MolProbity only a limited set of analysis options were available, primarily an option to generate a multi‐model multi‐criterion kinemage, and a multi‐model Ramachandran plot PDF file. In order to obtain the results of the entire MolProbity validation suite for an ensemble, users were forced to manually select each model and run the analysis on each model one at a time.

In MolProbity 4.4, we have updated the ensemble analysis so that users can run our entire suite of analysis options on the complete ensemble file. Summary charts for each model are presented under separate tabs within the results page, allowing users to easily click between different models and compare their validation results. Additionally, the analysis options have been updated to include all of the improvements documented in this article. In order to maintain a reasonable runtime, we limit the ensemble analysis to the first 80 models, which covers the majority of ensembles currently deposited in the PDB and still provides a good sense of quality for a larger ensemble.

Finally, note that biological‐unit files from the PDB are expressed as multi‐model files, but they have a very different logic, since the “models” are actually separate instances of the same chain within a molecule, rather than alternative models for the same chain in the same place. MolProbity deals with these biological‐unit multi‐model files by attempting to remap the chains in the different models into a single model, giving each chain a new unique chain ID. This code has been updated to use 2‐character chain IDs, which allows analysis of larger biological‐unit files.

Primary validation run and summary

After adding hydrogens, the next step is to run “Analyze all‐atom contacts and geometry”. In that setup, MolProbity turns on appropriate options dynamically depending on the contents, size, and resolution of the structure being analyzed, but you can of course change them. Help explanations are linked from each option to help you decide.

The report of results starts with a summary table coded in stoplight colors, plus full details in graphics and in chart form. The following sections describe each individual validation type, in the order they appear on the summary.

All‐atom contact analysis

The first line in the summary table reports the all‐atom clashscore (number of bad atom‐atom overlaps ≥ 0.4Å per thousand atoms). Clashes are the single most powerful diagnostic for many kinds of local fitting problems. For instance, the backward His in Figure 7(A) has the wrong protonation and N placement as well as clashes, and can be corrected automatically as in panel B. The clashing “water” in Figure 7(C) is really a positive ion, as shown in panel D. If you are the structural biologist, start fixups on the largest clashes first when you return to model rebuilding (you can sort any column in the html chart by severity). Clashes are directional as well as local, so that, in context of the electron density, it is usually possible to figure out the underlying problem. For a severe clash, some group is usually turned around backward or misidentified; just pushing the two atoms apart is almost never a good answer. If you are an end‐user, zoom in on the part that interests you most. If that region is clear of outliers, then the details there are probably quite reliable even if the overall scores are poor.

All‐atom contact analysis also evaluates presence or absence of H‐bonds and favorable van der Waals contacts. Those add further useful information, intuitively understandable in the graphics, such as the improbably missing H‐bonds between an Arg sidechain and the adjacent RNA phosphate shown in Figure 7(E) and corrected in Figure 7(F). Although not considered in standard refinement or model building, the absence of H‐bonding in such an arrangement is actually a strong signal of probable misfitting.

Sidechain rotamers

The new Top8000 dataset also allowed us to revisit and improve our empirical distributions for rotamer validation.44 Nearly a million (983,574) non‐Gly, non‐Ala residues passed the new electron‐density and other quality filters to be included in the final reference dataset from which the individual multi‐dimensional distributions are made. This increased size and accuracy allowed us to define 3 validation regions for rotamers as has long been done for Ramachandran‐plot criteria, with 98% of observed reference rotamers in the favored region, 0.3% in the outlier region, and the remainder falling in the allowed region between those cutoffs. Overall outlier frequency in general data has remained about the same, although some rotamer centers have shifted slightly or even divided, and additional rotamer peaks have reached above the outlier level for the long, multi‐dihedral sidechains. See the rotamer paper44 for methods, a complete list of rotamer names, χ values and reference‐data frequencies, and an explanation of why we can call this MolProbity's “ultimate” rotamer library.

Ramachandran backbone criteria

In Ramachandran analysis, for each residue the backbone ϕ and ψ dihedrals are calculated and compared against outer contours for the expected distribution of those angles. Those outer contours match very closely for most amino acid types (general case), but differ strongly for the other 5 types, since they have idiosyncratic configurations near the Cα (see18 supplement). The Top8000 dataset allowed us to add unique distributions for the cis‐proline and the branched Cβ (Ile/Val) cases and to update our existing Gly, Pro, prePro, and general cases (Fig. 8). The cutoffs for Favored vs Allowed vs Outlier are unchanged, with 98% of observed reference datapoints in the favored region and 0.5% (1 in 2000) in the outlier region. The only changes from older scores commonly seen in new Ramachandran validation results are for Ile/Val residues that were acceptable by the old but not the new specific criterion, or general‐case residues newly acceptable which lie along the more continuous region that runs down the positive side of ϕ.

MolProbity score

At the frequent request of users, we developed the MolProbity score as a combined single indicator of model quality.13 It uses a weighted function of clashes, Ramachandran favored, and rotamer outliers, scaled and normalized so that its value approximates the resolution at which that score would be average. For MolProbity score and for clashscore, the summary table includes a percentile relative to structures near the same resolution. MolProbity scores and percentiles give a quick, rough guide for end‐users to compare entries at different resolutions for their molecule of interest, and structural biologists should aim to improve their model if scores are below average.

Covalent geometry

MolProbity performs covalent bond geometry analysis for both mainchain and sidechain atoms in protein and nucleic acids, now using the Phenix geometry libraries. Bond lengths or angles more than 4σ from the expected value are considered outliers, and are flagged in both chart and graphical forms. Another validation unique to MolProbity is the Cβ deviation,11 a combination of covalent angles and chirality around the Cα that flags geometry problems there much more effectively than simply analyzing individual variables. Large Cβ deviations usually mean that either the sidechain or the backbone has been misfit at that residue.

By default, the expected values for bond geometry are drawn from a single‐value library derived from Engh and Huber.45 However, single values, even with standard deviations, do not fully capture covalent bond geometry in the complex environment of macromolecules, since bond geometry is dependent on local backbone conformation. For this, a Conformation‐Dependent Library (CDL) has been developed46, 47 and implemented in Phenix48 for protein refinement. The CDL relates the expected covalent bond geometry to local backbone Ramachandran conformation. Because the expected bond geometry values in the CDL differ from those in the single‐value library (especially for the N‐Cα‐C τ angle), MolProbity validation now uses the CDL values for structures refined with the CDL, as detected from the REMARK 3 information of a submitted file. Similarly, for RNA, geometry targets are dependent on ribose pucker.

Cis or twisted non‐trans peptides

The peptide bond that joins adjacent amino‐acid residues in a protein has partial double‐bond character and therefore assumes a trans, or more rarely a cis, configuration. The cis configuration is significantly more common preceding a proline and results in a unique Ramachandran distribution for cis‐proline. To maintain this special relationship, we associate peptide bonds with their following residue. About 5% of prolines are cis, while only about 0.03% of all non‐proline residues are genuinely cis.

Recently, we were alerted to a surprising and improbable increase in the number of cis non‐proline peptide bonds being modeled,49 as shown in the plot (updated) of Figure 9(A). These are due to model‐building without consideration of prior probabilities, but also in part due to the lack of validation that flagged cis‐nonPro peptides, in MolProbity or other systems. We have therefore implemented a new validation and visual markup for non‐trans peptides. Matching the PDB definition, we define a cis peptide as one with an ω angle between −30° and +30°, and a trans peptide as one with an ω angle > +150° or < −150°. We add an additional definition of “twisted peptides” for ω angles that are more than 30° from planar. Justifiable twisted peptides are even rarer than non‐proline cis,50 and twisted peptides should virtually always be considered modeling errors.

MolProbity reports on non‐trans peptides by providing counts of cis prolines, cis non‐prolines, and twisted peptides. Counts of cis non‐prolines or twisted peptides that constitute a suspiciously high percentage of the structure are flagged with yellow or red in the summary statistics chart. In the multi‐criterion chart that reports on each residue individually, each non‐trans residue is marked with its category (cis Pro, cis nonPro, twisted Pro, twisted nonPro) and the measured value of its omega peptide dihedral. In the multi‐criterion kinemage, each non‐trans peptide is marked with a surface that fills in the trapezoidal shape between the backbone trace and the Cα trace of a cis peptide [Fig. 9(B)]. These trapezoids are offset slightly from the model for ease of reading and are color‐coded by severity. Cis prolines are marked in gentle sea green, as they are relatively common and expected. Cis non‐prolines are marked in a more aggressive lime green, as they are likely to be errors. The vanishingly rare twisted peptides are marked in a warning yellow. Additionally, the interior angle between the two triangles of the trapezoidal shape indicates the severity of the twist.

Genuine cis‐nonPro peptides, like most very rare conformations, are nearly always found at functional sites, because evolution does not conserve unfavorable conformations unless they are biologically useful. Validation, and subsequent correction, of the incorrect ones is valuable because it improves the signal‐to‐noise for identifying the important, genuine cases. Following the discussions, papers, newsletters [e.g. Williams and Richardson (2015)23], and screaming markups in MolProbity and in Coot,51 it seems from the last few points in Figure 9(A) that this epidemic of overuse is now on the way to being cured.

CaBLAM for lower resolutions

We have introduced another new validation targeted at low‐resolution crystal structures, the Cα Based Low‐resolution Annotation Method, or CaBLAM.52 Low‐resolution structures pose two particular challenges to validation. First, they tend to contain many modeling errors, making it difficult to choose where to start. Second, sensitive validations, which work well in high‐resolution structures, may be confounded by the compound errors common in low‐resolution models. Ramachandran validation is a particular case in point ‐ Ramachandran validation is sensitive to small dihedral changes and is a powerful tool at high resolution. However, very large distortions of the Ramachandran dihedrals, especially multiple distortions in series, can mislead Ramachandran analysis.

CaBLAM is a validation of protein backbone that is more robust in the 2.5–4.0Å regime than Ramachandran analysis. CaBLAM takes advantage of the phenomenon that the overall Cα trace is relatively well represented in low‐resolution electron density and is relatively well modeled by humans and programs even when other details of the model are not. CaBLAM therefore uses Cα geometry to determine the intended structure of a region. It then checks the details of the model against this assessment of intended structure using contours derived from the Top8000, and declares outliers where there is a mismatch, as for most of the black datapoints in Figure 10(A,B).

The measures CaBLAM uses are, for each residue, two Cα pseudo dihedrals we call μ_in and μ_out and a dihedral that relates carbonyl oxygen orientation across the residue we call ν. Taken together, μ_in and μ_out indicate the intended structure. α helix, 3₁₀ helix, and β strands can be identified from μ_in and μ_out. The ν dihedral is sensitive to incorrect orientations of the peptide plane, which constitute a variety of common modeling errors at low resolution, especially misplacement of the carbonyl oxygen due to ambiguous or truncated sidechain density.53 Together these measures both identify modeling errors and suggest probable secondary structure elements, as shown on a local region in Figure 10(C,D). In contrast, Ramachandran and DSSP analyses are unhelpful and even misleading in cases such as this.

CaBLAM is most useful for validating secondary structure elements in low‐resolution models. In a high‐quality model, it provides little information beyond Ramachandran validation. Because loop structure is highly varied, CaBLAM's training set cannot recognize all valid loop conformations, and it tends to generate false positives in loops. Nevertheless, CaBLAM provides uniquely useful validation in a resolution regime where other methods struggle. For these reasons, MolProbity runs CaBLAM by default only for structures at 2.5Å or worse from either X‐ray or cryoEM.

RNA pucker and backbone conformers

In addition to the all‐atom contact and covalent geometry criteria applicable to nucleic acids as well as proteins, MolProbity provides conformational criteria tailored for RNA. C3'‐endo vs C2′‐endo ribose pucker cannot be seen directly at typical RNA resolutions, but it can be deduced from features which are observable and reliably modeled: the spherical phosphate density and the direction of the glycosidic bond that connects the base blob to the ribose blob.25, 54 The “Pperp” criterion (see Methods) is easily approximated visually while fitting, as how far the 3' P is out from the plane of the previous base,55 and it has enabled pucker‐specific torsion and geometry targets during Phenix refinement. A pucker outlier by this criterion is nearly always wrong.

In a community consensus study,25 42 clearly valid and several “wannabe” RNA backbone conformers were defined and named. This process was aided by the “suite” parsing of RNA backbone from sugar to sugar rather than the phosphate‐to‐phosphate nucleotide, which relates adjacent bases and within which the dihedrals show higher correlation.56 These backbone conformers, with a few updates, (see Methods) are a validation criterion in MolProbity. They are reported in the chart by name and “suiteness”: 1.0 at the 7‐dihedral cluster mean and 0 at the edge, as calculated by the Suitename program.25 They cover most but not all valid conformations, so each outlier should be examined but may be valid, especially if it is in an extended arrangement.

Multi‐criterion chart and graphics

The multi‐criterion kinemage produced by MolProbity or phenix.kinemage is the most powerful feedback we provide for a detailed exploration of a structure and its challenges. The visual markup employed in the kinemage encodes not just the location of outliers, but also their severity. Larger steric overlaps result in larger clash spikes, the sphere that marks a Cβ deviation has a radius equal to the deviation distance, and the “fans” that mark bond angle outliers extend between the modeled angle and the ideal, as in Figure 11(A). Additionally, the visual density of outlier markup is key to identifying problem regions in a model. We strongly recommend studying the multi‐criterion kinemage as linked either in the Phenix GUI or on‐line in MolProbity. If Java is impossible in your browser, download the kinemage and view it on your computer, using the KiNG program57 (available standalone on GitHub or packaged within Phenix both in the GUI and as a command‐line tool) and Java. We are investigating the addition of MolProbity markup in other software, for alternative future on‐line viewing (see Methods/Software).

The multi‐criterion chart (an html page) gives detailed, sortable information on every residue, or only on each residue with an outlier if that option was specified. It now uses an enhanced coloring scheme to reproduce some of that visual intuition for outlier severity and problem regions in chart form [Fig. 11(B)]. Previously, any table cell representing an outlier would be colored hot pink, and all other cells would be uncolored. Now there are three colors ‐ light pink, hot pink, and bright red. These colors were selected to be distinguishable in gray‐scale as well.

Light pink is used for less‐favored conformations and minor outliers. During our CASP8 assessment,58 we found that it was often useful to overlook small clashes (with overlaps of 0.4–0.5Å) in favor of larger problems in particularly challenging predicted models. The same proved true for experimental models solved at low resolution. Small clashes are now colored light pink in the chart [Fig. 1(B)], where previously they had been hot pink. Residues in the Ramachandran “Allowed” region were previously only identified by text in the multi‐criterion chart. Now these less favored, but non‐outlier conformations are marked with light pink, as are the new “Allowed” rotamer conformations. Similarly, the new CaBLAM validation uses light pink for residues in its “Disfavored” region. Residues colored light pink in the chart should be considered worthy of attention, but are not necessarily outliers.

Bright red indicates particularly severe outliers. Clashes with an overlap of >0.9Å, bond length and angle outliers of >10σ, Cβ deviations >0.7Å, twisted peptide bonds >45° from planar, and Cα geometry outliers identified by CaBLAM are marked in the chart with bright red [Fig. 11(B)]. Ramachandran and rotamer validations do not currently produce a “severe outlier” designation suitable for marking in this way. The general outliers marked with hot pink are sometimes justified as valid by strong density, hydrogen bonding, structural homology, or other factors. The cutoffs for the severe outliers are set such that any outlier marked with bright red is almost impossible to justify.

Website service

MolProbity is a research, software, and service project that deals with large volumes of complex data. The service component performs comprehensive validation on individual macromolecular structures, where it pays special attention to the local anomalies which are usually errors but sometimes valid and biologically important. After producing effective, user‐friendly software that creates user demand, the over‐riding requirement for service is near‐24/7 uptime of the website. That both requires long‐term management and equipment upgrades to maintain capacity, run speed, and security, and also requires constant short‐term attention: automatic monitoring where feasible, attention to user alerts, and frequent checks of the server state for possible hardware failures, persistent attack trials, or hung jobs. Our team, especially the system manager, almost always notice and fix problems within a few hours, nearly all of which is set up so it can be done remotely.

The website also produces error messages, with an email bug report option. We respond promptly to problems and queries, fix actual bugs as soon as feasible, and prioritize requests for new features if they are of fairly broad use and we have staff that can provide them. File format problems are the most common cause of bug reports, although those have decreased since we now diagnose the most frequent ones at the input stage (either fetch or upload).

Software

Our programs where speed or complexity is an issue are in C or C++. As described in the Infrastructure section, our connecting scripts and utilities formerly in Java have been rewritten in Python, as are new ones. The website is primarily controlled by PHP, especially for its dynamic responses to input, user choice, and intermediate results. Our online 3D interactive molecular display is done by the Java KiNG program.12, 57 Its structure modeling and corrections, docking construction, image and movie creation, and of course interactive display can still be done offline by us and others, but Java is nearly dead for the important online display in MolProbity. We have established collaborative initial systems to read and display 3D MolProbity validation markup on the molecule in ChimeraX60 with Tom Goddard and Tristan Croll, and in the RCSB PDB's new Javascript NGL viewer61 with Alex Rose. We will continue to enhance those capabilities and spread them to as many different viewers as feasible.

As described in Infrastructure, we use Git and GitHub as both our software version control for development and as our open‐source distribution system; it also provides a second backup. As Phenix developers, our components participate in their unit and regression test system. Most importantly, we ourselves both use in our own research and also set up deliberate tests of the detailed results from our software and look critically for unanticipated problems.

Reference datasets

The research component of MolProbity selects and quality filters high‐resolution, non‐redundant structures to obtain our basic chain‐level reference datasets. For each distinct use, those chains are quality filtered at the residue level, eliminating most errors by an optimal compromise between total number of residues and reliability of each one of them. For this reference data, our policy accepts a moderate rate of false negatives (i.e., correct residues not included) in order to get a low rate of false positives, provided that the rejection criteria are not logically connected to the criterion we will be testing and evaluating. In order to ensure satisfying the low false‐positive rate, we examine several dozen individual examples near the proposed thresholds to check for distortions, unconvincing electron density, or other circumstances which could allow acceptance of a problematic local model.

Conformational and validation parameters for each residue of the reference data are stored in, and queried from, a database. Earlier work used MySQL,62 and the recent rotamer work44 set up a Mongo database http://github.com/mongodb/mongo. Results from such structural‐bioinformatics queries are converted into multi‐dimensional data distributions smoothed and contoured by iterative kernel‐density methods.11 The validation measures eventually resulting from this process are characterized in stringency by the percent of outliers in the residue‐filtered, high‐resolution data; currently this varies from 0.03% for cis non‐proline peptides to 0.3% for rotamers to 1% for CaBLAM outliers.

The two primary reference datasets used in the work described here are the Top8000 for proteins and RNA11 for RNA. To construct the Top8000, we considered all protein chains solved by x‐ray crystallography and released by the PDB as of March 25, 2011. To be eligible for inclusion in the dataset, a chain had to have been solved at a resolution better than 2.0Å and to be of good structural quality. General structural quality was enforced by requiring chains to have a MolProbity score (a combination of clash, rotamer, and Ramachandran measures13) better than 2.0, and overall geometric quality was enforced by requiring chains to have ≤5% of residues with Cβ deviations, ≤5% of residues with bond angle outliers, and ≤5% of residues with bond length outliers. Chains that passed these filters were grouped according to PDB homology clusters, separately at 90%, 70%, and 50% sequence‐identity levels, and separately with and without requiring deposited diffraction data. In each homology cluster, the highest quality chain was selected for inclusion in the Top8000, where quality was determined by the average of a chain's resolution and MolProbity score. The “standard” Top8000 used for Ramachandran and CaBLAM is filtered at the 70% homology level and contains 7957 protein chains. In preparation of Ramachandran and CaBLAM contours, we excluded residues having any mainchain atom with a B‐factor > 30. This simple filter provides a proxy for local model fit to the electron density, and it indeed produced quite clean, reproducible contours.

For RNA, high‐resolution data is still relatively sparse, so criteria must be more forgiving, but more hand selection is feasible. For the RNA11 reference dataset, all X‐ray structures as of 11/11/11 containing at least a 3‐nucleotide RNA chain at ≤3.0Å resolution were hand‐selected for homology, in order to allow for inclusion of more than one complex or condition of the same sequence showing a significantly different RNA conformation (such as 5S rRNA alone versus in the ribosome). That process produced a chain‐level set of 311 structures, up from the RNA05 dataset of 171.25

The enforced deposition of structure‐factor data along with each new entry in the worldwide Protein Data Bank (wwPDB15) provides an opportunity for a more accurate and complete residue‐level filter. To this end, we revised the Top8000 in 2015 to require each chain to have deposited structure factors in addition to the other criteria. We then developed a new residue‐level filter combining local real‐space correlation coefficient, 2mFo‐DFc σ value at each atom, and B‐factor.44 Optimization produced threshold criteria using all 3 terms: real‐space correlation coefficient ≥0.7, 2mFo‐DFc σ ≥ 1.1, and atomic B‐factor <40, for all relevant atoms in the residue. Future validations and revisions will take advantage of this improved residue‐level filtering.

Electron‐cloud hydrogen positions

This study involved multiple methods and the work of over half the authors over several years. The method description is therefore very lengthy and is found in the Supplement.

CaBLAM

CaBLAM was developed to be a robust validation even in models where many atoms are placed incorrectly.53 As a result it seeks to use a minimal set of atoms ‐ Cαs and COs ‐ to define its four parameters. The main parameters are Cα pseudo dihedrals μ_in (defined as Cα_i‐2‐Cα_i‐1‐Cα_i‐Cα_i+1 for residue i) and μ_out (defined as Cα_i‐1‐Cα_i‐Cα_i+1‐Cα_i+2 for residue i), as shown in Figure 12. The third dihedral ν is defined using Cα_i‐1, Cα_i, Cα_i+1, O_i‐1, and O_i for residue i. Additional pseudoatom points are constructed: X_i‐1 on the Cα_i‐1‐Cα_i line at the point closest to O_i‐1, and X_i on the Cα_i‐Cα_i+1 line at the point closest to O_i. The ν dihedral is then defined as O_i‐1‐X_i‐1‐X_i‐O_i (pink in Fig. 12). Finally, the Cα virtual angle for residue i is defined as Cα_i‐1‐Cα_i‐Cα_i+1.

These four parameters were calculated for each protein residue in the Top8000. Residues for which any of the atoms used in the calculation were missing or had B‐factor >30 were excluded. Thus only residues from fully modeled and confidently modeled regions were included in the final training set.

Three‐dimensional contours were generated using two different combinations of these parameters: contours in the μ_in‐μ_out‐ν parameter space define expected protein backbone behavior, and are used to identify the majority of CaBLAM outliers; contours in the μ_in‐μ_out‐Cα virtual angle parameter space define expected Cα trace behavior and are used to identify severe Cα trace modeling errors. As in Ramachandran evaluation, proline residues populate a restricted portion of the CaBLAM parameter spaces and glycine residues populate a less restricted portion than the general case. Therefore, separate contours are defined for the proline, glycine, and general cases of residue type.

A further set of two‐dimensional contours were generated to define the observed behavior of secondary structure in the μ_in‐μ_out parameter space. Residues in the training set were identified as α helix, 3₁₀ helix, or β strand based on hydrogen bonding pattern. A residue i was identified as α helix if residues i‐1, i, and i + 1 all participated in n to n + 4 hydrogen bonding or if residues i‐1, i, and i + 1 all participated in n‐4 to n hydrogen bonding. These paired, one‐sided definitions of helix allowed the correct identification of helix residues near the ends of helices. A residue i was identified as 3₁₀ if residues i‐2, i‐1, and i all participated in n to n + 3 hydrogen bonding and residues i, i + 1, and i + 2 all participated in n‐3 to n hydrogen bonding. A residue i was identified as β strand if it was on a middle strand of a β sheet, and a full cycle of β‐pattern hydrogen bonding (parallel, antiparallel, or a combination of the two) both preceded and followed that residue. This definition of β structure is somewhat restrictive, but was necessary to generate clean contours. This provided a mapping between Cα geometry and hydrogen bonding patterns.

Setting cutoff values for the contours was done heuristically by manual inspection of structures with known modeling errors. For the μ_in‐μ_out‐ν CaBLAM space, the outlier cutoff was set at the bottom 1%, and a second “Disfavored” cutoff roughly analogous to the Ramachandran “Allowed” cutoff was set at the bottom 5%. The top 95% are considered favored conformations. For the Cα‐only μ_in‐μ_out‐Cα virtual angle space, a single outlier cutoff was set at the bottom 0.5%. This space did not receive an “Allowed” analogue cutoff, as it is intended primarily as a sanity check on Cα trace interpretation. The lower bounds for α and 3₁₀ helix contours were set at 0.1%, and for β strand at 0.01%. Additionally, for a residue to be identified as secondary structure, that residue's preceding and following residues are also required to score well on that secondary structure's contours. All these definitions strive to achieve a balance between a generous interpretation of structure appropriate to addressing low‐resolution structures and a requirement for structure continuity that prevents false identifications.

When CaBLAM validation is run on a structure, the four CaBLAM parameters ‐ μ_in, μ_out, ν, and Cα virtual angle ‐ are calculated for each protein residue. Much as Ramachandran validation is not possible for the first or last residue in a chain, CaBLAM assessment is not possible for the first two or last two residues because not all atoms necessary to calculate μ_in and μ_out are present. For each residue with a complete parameterization, its parameters are compared to the μ_in‐μ_out‐ν contours to identify CaBLAM outliers. Each residue is also compared to the μ_in‐μ_out‐Cα virtual angle contours to identify Cα geometry outliers. Cα geometry outliers are usually excluded from the following secondary structure identification step on grounds that secondary structure identification is dependent on a reliable Cα trace. Residues are then compared to the μ_in‐μ_out contours for secondary structure. Individual residues identified as matching the appropriate μ_in‐μ_out geometry are assembled into secondary structure elements ‐ helices and strands ‐ based on adjacency to other residues of like geometry. Residues that are part of these assembled secondary structure elements are reported as probable secondary structure, with their % contour as a score; residues with isolated secondary‐structure‐like geometry are not reported.

CaBLAM validation results are available in several forms. On the Phenix commandline, phenix.cablam produces colon‐delimited text output by default. On the MolProbity website, this same information is available in the structure summary and residue‐level tables. CaBLAM also provides visual markup for outliers in kinemage format. CaBLAM outliers and disfavored residues are marked with lines that follow the ν dihedral of that residues, colored hotpink for outliers and purple for disfavored [see Fig. 9(A)]. Cα geometry outliers are also marked, using red lines that follow the Cα trace along the Cα virtual angle of the residue.

CaBLAM is a unique tool for validating protein backbone, and especially secondary structure, at low resolution or in otherwise difficult models. Because CaBLAM's outlier cutoff is unusually punishing for a validation ‐ excluding a full 1% of observed structure ‐ a significant number of residues are expected to present as outliers even in good structures. As a result, CaBLAM performs best on regular regions of structure and tends to offer false positives in loop regions due to their high irregularity. Ramachandran evaluation is therefore still recommended as the backbone validation of choice for high‐resolution structures. Nevertheless, CaBLAM extends the reach of protein backbone validation into a resolution range at which Ramachandran evaluation is not reliable.

RNA pucker and backbone

Ribose pucker outliers (C3'‐endo when they should be C2′‐endo or vice versa) are checked by the “Pperp” versus δ dihedral criterion from Richardson et al.25 A perpendicular is dropped from the 3′ P atom to the vector direction of the glycosidic bond; its length is >2.9Å for C3′ pucker and ≤2.9Å for C2′. Backbone suite (sugar‐to‐sugar) conformers are defined, both originally and here, by 7‐dimensional suite‐torsion datapoint clustering, shown and analyzed in the Mage63 or KiNG57 display programs. A named conformer requires ≥5 examples with B‐factor <60, or ≥3 if one is at <2.0Å resolution.25 Their validation for new structures is calculated by the Suitename program.

In present work with RNA11, we made use of the ERRASER program,64, 65, 66 which uses MolProbity, Phenix, and Rosetta to perform and evaluate exhaustive backbone‐conformer search, for problematic dinucleotides. It does not use the information of our defined suite conformers, but almost always reproduces them in its output. The most frequent apparently valid RNA backbone outliers we see are for bulged or junction suites with widely separated bases, which have high variability and therefore few examples per conformer. To test whether such a case with good electron density might be a valid new suite conformer, we ran ERRASER on it and related structures to see if a consistent consensus conformer was found. For instance, Figure 13 shows a spread‐out junction at 2.9Å in the 3gx5 riboswitch, previously an unrecognized and thus outlier (!!) conformation, but with completely unambiguous electron‐density support. ERRASER left it as is, and corrected several less well‐fit examples in related structures to match. This conformer has now been named 3h and added to the suite conformer list.

Supplementary material

MolProbity4_Supplement.pdf describes in detail the methods used in the revision of hydrogen distance parameters and lists all of the structures used.