Computational biology ; protein structure prediction ; Model quality assessment programs ; Boltzmann distribution ; Decoys `R' Us ; Annsen's thermodynamic hypothesis ; Finite difference PoissonBoltzmann (FDPB) ; APBS ; statistical potentials ; protein sidechain ; decoy sets ; template based modeling ; ab initio protein structure prediction ;
We have modified our discussion to reflect the point we had missed previously regarding the pKa values in accordance with Dr Kammerlin's suggestion.
See the authors' detailed response to the review by Patricia C Weber
The challenge of deriving the native structure of a protein from its sequence has intrigued researchers for decades1. Methods that predict putative structures from sequences are based either on features from databases of known structures (template-based methods)2–4 or use first principles of atomic interactions (ab initio or de novo methods)5–7. Typically, these methods yield multiple, closely-ranked possibilities. Model quality assessment programs (MQAP) that validate accuracy of these predicted structures are used to select the best candidate from the set of predicted structures.
MQAPs can be classified as energy, consensus or knowledge based. Two major sources of errors in energy based methods used for refining or discriminating protein structures are inaccuracies in the force field due to the inherent approximations in equations that model multi-atomic configurations, and inadequate sampling of the conformational space8–12. Consensus based methods are based on the principle that structural features that are frequently observed in a population of structures are more likely to be present in the native structure13–16. These clustering methods outperform other MQAP methods14 and are “very useful for structural meta-predictors17”. However, they are prone to be computationally intensive due structure-to-structure comparison of all models16, and are of limited use when the number of possible structures is small18. Knowledge based methods proceed by deriving an empirical potential (also known as statistical potential) from the frequency of residue contacts in the known structures of native proteins19,20. For a system in thermodynamic equilibrium, statistical physics hypothesizes that the accessible states are populated with a frequency which depends on the free energy of the state and is given by the Boltzmann distribution. The Boltzmann hypothesis states that if the database of known native protein structures is assumed to be a statistical system in thermodynamic equilibrium, specific structural features would be populated based on the free energy of the protein conformational state. Applying a converse logic, Sippl reasoned that the frequencies of occurrence of structural features such as interatomic distances in the database of known protein structures could be used to assign a free energy (potential of mean force) for a given protein conformation21,22. Furthermore, this statistical potential can be used to discriminate the native structure23–27. The proper characterization of the reference state is a critical aspect in applying statistical potentials23. In spite of their popularity, the application of such empirical energy functions to predict and assess protein structures are vigorously debated28,29. Many MQAP programs perform better when multiple statistical metrics are combined30–33. The paramount importance of obtaining high quality protein structures from sequences using in silico methods can be estimated by the effort invested by researchers every two years34 to evaluate both structure prediction tools35 and MQAPs17,34,36.
Here, we propose a novel statistical potential to assess the quality of protein structures based on the electrostatic potential difference (EPD) of Cβ atoms in consecutive residues - EPD profile of sidechain atoms used in assessment of protein structures (ESCAPIST). Previously, we have established that the EPD is conserved in cognate pairs of active site residues in proteins with the same function37–40. The ability of finite difference methods to quickly obtain consistent electrostatic properties from peptide structures provides an invaluable tool for investigating other innate properties of protein structures41. We plot the EPD profiles for different atom types (Cα atoms, Cβ atoms and the C-N bond) in consecutive residues from a set of non-homologous protein structures obtained from the PISCES database (http://dunbrack.fccc.edu/PISCES.php)42. We proceed to show that the EPD between Cβ atoms in consecutive residues can be used to generate a scoring function that assesses the quality of protein structures. This EPD scoring function is then applied to standard decoy sets from the Decoys ‘R’ Us database (http://dd.compbio.washington.edu) to establish the validity of our method43.
To extract feature values we chose a set of 1000 proteins from the PISCES database with percentage identity cutoff of 20%, resolution cutoff of 1.6 Å and a R-factor cutoff of 0.25 (SI Table 1).
While these pairs have for a low standard deviation (SD) like all other pairs, the absolute value of their mean is different (higher) than any pair that does not include a proline. This also highlights the unique nature of proline in protein structures.
Adaptive Poisson-Boltzmann Solve (APBS) writes out the electrostatic potential in dimensionless units of kT/e where k is Boltzmann’s constant, T is the temperature in K and e is the charge of an electron. The units of EPD are same as that of the electrostatic potential. The EPD of the C-N peptide bond has a Gaussian distribution with mean = 420 EPD units and SD = 55 EPD units (Figure 1). In the probability distribution for four pairs of amino acids the mean of all pairs of amino acids are the same (Figure 1a). Figure 1b shows the scatter plot for the mean and standard deviation (SD). Thus, the amino acids are indistinguishable using the profile of the EPD of the C-N peptide bond across all protein structures since they have identical mean values and a large variance (SD=~50).
AA: Alanine/Alanine, AC: Alanine/Cysteine, HS: Histidine/Serine and DF: Aspartic-acid/Phenylalanine. (a) Probability distribution for four pairs of amino acids. (b) Scatter plot for all pairs of amino acids. It can be seen that the mean and SD for all pairs of amino acids are the same. Further, the variance is large (SD=~50), indicating that this feature is not tightly constrained in peptide structures.
The probability distribution for four pairs of amino acids for the EPD between the Cα atoms of consecutive residues (Figure 2a) have means that are slightly more varied than those for the C-N bond (Figure 1a). In the scatter plot for the mean and SD of all pairs (Figure 2b) the outliers are pairs that include proline, which have a higher mean, although the magnitude of SD is the same (Table 1).
A: Alanine/Alanine, AC: Alanine/Cysteine, HS: Histidine/Serine, DF: Aspartic-acid/Phenylalanine. (a) Probability distribution for four pairs of amino acids. (b) Scatter plot for all pairs of amino acids. It is seen that pairs of amino acids which include proline have a higher mean, although the magnitude of SD is the same.
In contrast to the results described above, the EPD between the Cβ atoms in consecutive residues in the peptide structure can be used to discriminate different amino acid pairs in the protein structure. The mean EPD of all amino acid pairs are much more varied (Figure 3a). These pairs do not include glycine, which lacks a sidechain. In the scatter plot for the mean and SD, the outliers are pairs that include cysteine (Figure 3b), which have a much higher SD (=~90) as compared to other pairs (SD=~35) (Table 2), and thus cannot be used for discriminatory purposes.
AA: Alanine/Alanine, AD: Alanine/Aspartic-acid, AE: Alanine/Glutamic-acid, DF: Aspartic-acid/Phenylalanine, DY - Aspartic-acid/Tyrosine, HT: Histidine/Threonine, HS: Histidine/Serine. (a) Probability distribution for seven pairs of amino acids. (b) Scatter plot for all pairs of amino acids. The pairs which include cysteine have a high standard deviation. It is seen that the mean is much more varied than the electrostatic potential difference (EPD) for Cα and the C-N peptide bond.
These pairs have a random values for the mean and a high standard deviation (SD), with the exception of the pair ‘CC’ (not the disulfide bond) which has a low mean value and SD. Consequently, these values can not discriminate between pairs of amino acids.
These values are used as a discriminator when choosing the native structure from a set of possible candidates (Table 3). To establish the non-triviality of these values, we also show that the variance of the EPD between these pairs increases with increasing sequence distance. Thus, the EPD between the pairs ‘DF’ and ‘HS’ has lesser correlation as the sequence distance between them increases (sample size for each sequence distance is > 30) (Figure 4). The SD for distance 1 (i.e. consecutive residues) is 29.8 EPD units and 31.8 EPD units for ‘DF’ and ‘HS’, respectively - and rises to around 60 EPD units with increasing sequence distance.
These pairs are used for discriminating predicted structures in order to obtain the native structure. The complete set is available at https://github.com/sanchak/mqap.
Each sequence distance has at least 30 sample points. DF: Aspartic-acid/Phenylalanine, HS: Histidine/Serine. As expected, there is lesser correlation in the EPD values between the shown amino acid pairs ‘DF’ and ‘HS’ as the sequence distance between the residues increases. The SD for distance 1 (i.e. consecutive residues) is 29.8 EPD units and 31.8 EPD units for ‘DF’ and ‘HS’, respectively - and rises to around 60 EPD units with increasing sequence distance.
We obtained the score (PDScore) of any given protein structure by comparing the electrostatics of the Cβ atoms based on Table 3. To benchmark model quality assessment programs, we used decoy sets from the Decoys ‘R’ Us database43. We detail our results from some of these decoy sets. Each set has several structures that are supposed to be ranked worse than the native structure.
The misfold decoy set has ~20 protein structures, each of which has a correct and an incorrect structure specified (three structures have two incorrect structures: we randomly chose the first)44. The PDScore of these proteins were computed (Table 4). Barring three structures (PDBids: 1CBH, 1FDX and 2SSI), the PDScore of the incorrect structure is higher than that of the correct structures.
This decoy set has ~20 protein structures - each of which has a correct and an incorrect structure specified. The PDBs are sorted based on the number of residues in the structure (NRes). Three of the structures (1CBH, 1FDX and 2SSI) have a lower PDScore for the incorrect structure.
The hg_structal set has about ~30 proteins. Each protein has 30 structures (including the native structure). Table 5A shows specificity obtained for structures in this decoy set. The average specificity obtained for this decoy set is 0.91 (Table 5A). The decoy set 4state_reduced has ~600 structures for each of the seven proteins. We obtain an average specificity of 0.94 for this decoy set (Table 5B). Similarly, for the ig_structal decoy set we obtain a specificity of 0.93 (Supplementary Table 1).
The PDBs are sorted based on specificity. (A) The hg_structal decoy set has ~30 protein structures - each of which has 30 structures. The average specificity obtained for the set is 0.91. (B) The 4state_reduced decoy set has 7 protein structures - each of which has ~600 structures. The average specificity obtained for the set is 0.94. (C) The fisa set has 4 protein structures - each of which has 500 structures. The electrostatic discriminator has low specificities in this case. We have previously demnostrated that this decoy set can be discriminated by a distance based criterion. It consists of physically nonviable structures, thus rendering an electrostatic analysis meaningless. NRes = number of residues, NStructures = number of structures in the decoy set.
The functional characterization of a protein from its sequence using in silico methods based on the ‘sequence to structure to function’ paradigm is contingent upon the availability of its 3D-structure. The rapidly developing field of next generation sequencing has exacerbated the bottleneck of obtaining structural data using crystallization techniques45. This ever-widening gap has been filled by methods that predict structures from sequences46, based either on features from databases of known structures2–4 or from first principles of atomic interactions5,6.
The various sources of error in protein structure prediction have been previously discussed in detail47. An incorrect model of a protein structure will result in an inaccurate analysis of its properties48. For example, continuum models49 that compute potential differences and pKa values from charge interactions in proteins50 are sensitive to the spatial arrangement of the atoms in the structure. It must be pointed out that other detailed methods using quantum mechanical/molecular mechanical (QM/MM) techniques and and doing extensive conformational sampling have been able to determine the side chain pKa values with high accuracy51. Accurate structural information is indispensable for in silico methods that extract the electrostatic profile of atoms in the peptide structure41,52, and for other methods widely used in pharmacology for drug discovery53. Model quality assessment programs (MQAP) that validate the accuracy of predicted structures are thus a critical aspect in the process of modeling a protein structure from its sequence. MQAPs can be classified as energy8–12, consensus13–16 or knowledge based (statistical potential)21–27. The state of the art methods for predicting structures35 and MQAPs17,34,36 are evaluated by researchers every two years.
Previously, we hypothesized and demonstrated that the electrostatic potential difference (EPD) in cognate pairs in the active site are conserved in proteins with the same functionality37,40,54, even when evolution has converged to the same catalytic from completely different sequences55. This similarity is observed in structures solved independently over many years and establishes the reliability of the APBS and PDB2PQR implementations41,56. We focused on unraveling other electrostatic features that are innate to protein structures. Here, we first demonstrate that the EPD between the C-N peptide bond and the Cα atoms of consecutive residues are independent of the amino acid type. This is expected, since the distance between these atoms are almost invariant across all structures. The EPD of the C-N bond has a high variance, implying that the backbone accommodates relatively large variations while seeking energetically minimized structures.
The true source of the chemical and structural diversity in protein structures is the side chain atoms. Every amino acid, except glycine, has a Cβ atom that hosts a unique moiety of atoms. Although the reactive groups are different for amino acids, we show that this difference is encapsulated in the backbone Cβ atoms. We first show that different pairs of amino acids have significantly different mean EPD values in side chain Cβ atoms (Figure 3), unlike the EPD of the C-N peptide bond (Figure 1) or the EPD between consecutive Cα atoms (Figure 2). Further, the variance is much less than in the EPD of the C-N bond. These facts suggested that the EPD between Cβ atoms of consecutive residues in the peptide structure might act as a discriminator. Our hypothesis is based on the insightful Boltzmann law that lower energy conformations are disproportionately sampled, on the thermodynamic hypothesis57 that the native structure has minimal energy, and the hypothesis that statistical derived features in known protein structures have a Gaussian distribution21. We apply our discriminator to standard decoy sets from the Decoys ‘R’ Us database to establish the validity of the method43.
Our work also highlights the unique properties of proline in the protein structure58. This is evident from the different magnitude of EPD in consecutive Cα atoms involving proline (Table 1). Another noteworthy aspect is the large variation in EPD in consecutive Cβ atoms involving cysteine (Table 2), demonstrating the unique role cysteine plays in providing flexibility to protein structures, a critical element in the evolution of complex organisms59. The discrimination of Cβ atoms also provides a uniform basis for methods that require a single-atom representation of a residue. Such methods depend on a correct parameterization of the reactive atoms37, a task further complicated by amino acids such as histidine which has two reactive atoms. For example, the EPD between the negatively charged E and D with respect to the aromatic phenylalanine is -108 and -93 EPD units, in spite of the difference in their reactive atom. Similarly, the EPD between alanine and the three aromatic amino acids (F, W and Y) are -67, -66 and -63 EPD units respectively.
We achieved an average specificity of 0.91, 0.94 and 0.93 on hg_structal, 4state_reduced and ig_structal decoy sets, respectively, taken from the Decoys ‘R’ Us database. We have previously demonstrated that the fisa decoy set can be discriminated by a distance based discriminator60. ESCAPIST does not discriminate the native structure in this decoy set (Table 5C). The physical implication of ESCAPIST results on the fisa decoy set, which has significant RMSD for backbone Cα atoms, needs further elaboration. The input to a finite difference Poisson-Boltzmann (FDPB) analysis is a charge distribution that might be unfeasible due to energy functions other than electrostatics. For example, van der Waals force or the elastic bond length force components might prevent two atoms from being in close proximity. However, if such a physically impossible configuration were presented to a FDPB-based analysis tool, such as APBS41, it would still generate an electrostatic potential landscape. Inferences based on this potential landscape would be incorrect due to its physical non-viability. Thus, before invoking the EPD constraints specified here, it is imperative that other spatial constraints that are rarely violated in structures are checked. Possibly for this same reason, MQAPs that combine many features in their scoring functions are superior. Moreover, it should be kept in mind that decoy sets, like most benchmarking sets, are prone to biases61 and possible errors31. In fact, the fisa decoy set has been shown to violate the van der Waals term61. To summarize, we present a novel discriminating feature in protein structures based on the electrostatic properties of the side chain atoms. We validated this discrimination in several decoy sets taken from the Decoys ‘R’ Us database, and achieved high specificities in most decoy sets.
Our proposed method has two phases. In the learning phase, we process multiple structures to extract the feature values - mean values of electrostatic potential difference (EPD) for each amino acid pair. These feature values are applied on query proteins to obtain a score (PDscore) that indicates the deviation of the feature values in the given structure from the ‘ideal’ values. Thus, a lower PDscore indicates a better candidate.
Algorithm 1 shows the procedure LearnFeatures() that extracts the feature values from a set of proteins ΦLearningPhaseproteins (Equation 1). We ignore the first x=IgnoreNTerm and last y=IgnoreCTerm pairs of residues in the protein structure to exclude the terminals. For every consecutive pair of residues in the structure, we calculate the EPD (see below for method) between two specified atoms (atomP and atomQ). Both atomP and atomQ are set to Cβ to obtain EPD values for Cβ atoms, while we set atomP to ‘C’ and atomQ to ‘N’ in order to obtain the C-N peptide bond EPD values. The mean (Mean learnt value - MLV) and standard deviation (SD) are computed for each amino acid type pair (AAType) in protein (Equation 2), and the mean computed for the globals set of proteins (MLV(AATypex, AATypey)) for each pair of amino acid types (Equation 3). Pairs whose EPD have a SD greater than a threshold value (sdThresh, 50 by default) are ignored. Means for all significant pairs (ϕpairMean) are noted to a file, which is the input to the quality assessment phase. The EPD between a pair of amino acid is order-independent - for example, the EPD statistics for the pair ‘AC’ (alanine-cysteine) includes the EPD of both ‘AC’ and ‘CA’ (with the sign reversed).
Algorithm 2 shows the function AssessEPDQuality() that generates the PDscore for a given protein from the template file generated by the learning phase. The set of proteins consists of the native structure P1 and N-1 decoys structures (Equation 4). Once again, barring x=IgnoreNTerm and y=IgnoreCTerm number of residues from the N and C terminals, the pairwise EPD for consecutive residues are computed. The absolute value of the difference of these values from their corresponding means, if they exist, in the template file is added to generate the absolute score (Equation 5). This score is normalized with the number of residues that have been compared to obtain the final PDscore. In summary, the PDscore encapsulates the average distance of the EPD for the given atom pairs (it may be Cβ, Cα or the C-N bond) of consecutive residues from their mean values. We hypothesize that in the native or a near native structure, the PDscore will be minimized for the EPD of Cβ atoms of consecutive residues, i.e. given a set of proteins Φproteins consists of the native structure P1 and N-1 decoys structures, P1 will have the minimum PDscore (Equation 6).
The top level procedure ESCAPIST() is shown in Algorithm 3. It invokes the function LearnFeatures() once, and applies the learnt values to assess the quality of structures based on the feature values obtained.
Adaptive Poisson-Boltzmann Solver41 (APBS) and the PDB2PQR package56 package was used to calculate the potential difference between the reactive atoms of the corresponding proteins. The APBS parameters are set as follows - solute dielectric: 2, solvent dielectric: 78, solvent probe radius: 1.4 Å, Temperature: 298 K and 0 ionic strength. APBS writes out the electrostatic potential in dimensionless units of kT/e where k is Boltzmann’s constant, T is the temperature in K and e is the charge of an electron.
Input: ϕproteins = {P1 ··· PM} : M Proteins in the learning set
Input: IgnoreNTerm: Ignore this number of residues in the N Terminal
Input: IgnoreCTerm: Ignore this number of residues in the C Terminal
Input: atomP: Atom type in first residue
Input: atomQ: Atom type in second residue
Input: sdThresh: Threshfold for standard deviation of the EPD
Output: ϕpairMean = {meanPDC1 ··· meanPDCK}: Mean values of EPD between specified atoms
of successive residues, there being K such significant pairs
begin
/*K pairs of amino acid type (sorted: AC and CA are equivalent)*/
/*Each set is initialize to be the null set*/
ϕpair = {ϕ1PDC ··· ϕKPDC} :
foreach Pi in ϕproteins do
N = NumberOfResidues(Pi);
for p ← IgnoreNTerm to (N − IgnoreCTerm) do
q = p + 1 ;
/* Amino acid pairs are order independent */
ResiduePairTypeString = ResidueTypeString(p) + ResidueTypeString(q);
ResiduePairTypeStringSorted = Sort(ResiduePairTypeString;
/* Reverse sign of potential difference accordingly */
multfactor = 1 ;
if (ResiduePairTypeStringSorted != ResiduePairTypeString) then
multfactor = -1 ;
end
PD = PotentialDifference(p, q, atomP, atomQ) * multfactor ;
/* Let the amino acid pair be the kth in the set ϕpair */
InsertInSet(PD, ϕkPDC);
end
end
/* Compute Mean and SD of each set - ignore pairs with SD greater than sdThresh*/
ϕpairMean = ∅;
foreach ϕipair in ϕpair do
(Meani, SDi) = MeanAndSD(ϕipair);
if (SDi > sdThresh) then
Add(Meani, ϕpairMean);
end
end
return (ϕpairMean);
end
Input: P1 : Protein under consideration
Input: IgnoreNTerm: Ignore this number of residues in the N Terminal
Input: IgnoreCTerm: Ignore this number of residues in the C Terminal
Input: atomP: Atom type in first residue
Input: atomQ: Atom type in second residue
Input: ϕpairMean = {meanPDC1 ··· meanPDCM}: Mean values of EPD between specified atoms of
successive residues
Output: PDscore: Score indicating the normalized distance of the observed values from the (mean)
learnt values from native structures
begin
PDscore = 0 ; NumberCompared = 0 ; N = NumberOfResidues(P1);
for p ← IgnoreNTerm to (N − IgnoreCT erm) do
q = p + 1 ;
/* Amino acid pairs are order independent */
ResiduePairTypeString = ResidueTypeString(p) + ResidueTypeString(q);
ResiduePairTypeStringSorted = Sort(ResiduePairTypeString;
/* Reverse sign of potential difference accordingly */
multfactor = 1;
if (ResiduePairTypeStringSorted != ResiduePairTypeString) then
multfactor = -1 ;
end
/* Let the amino acid pair be the kth in the set ϕpair */
PD = PotentialDifference(p, q, atomP, atomQ) * multfactor ;
if (∃meanPDCk) then
NumberCompared = NumberCompared + 1 ;
diff = absolute(PD − meanPDCk);
PDscore = PDscore + diff;
end
end
/* Normalize */
PDscore = PDscore/NumberCompared;
return (PDscore);
end
Input: ϕproteins: Learning set
Input: P1: Protein to be scored
Input: IgnoreNTerm: Ignore this number of residues in the N Terminal
Input: IgnoreCTerm: Ignore this number of residues in the C Terminal
Input: atomP: Atom type in first residue
Input: atomQ: Atom type in second residue
Input: sdThresh: Threshfold for standard deviation of the EPD
Output: PDscore: Score indicating the normalized distance of the observed values from the (mean)
learnt values from native structures
begin
/* This is invoked once*/ ϕpairMean =
LearnFeatures(ϕproteins, IgnoreNTerm, IgnoreCTerm, atomP, atomQ, sdT hresh);;
PDscore = AssessEPDQuality(P1, IgnoreNTerm, IgnoreCTerm, atomP, atomQ, ϕpairMean);
return (PDscore);
end
| Views | Downloads | |
|---|---|---|
| F1000Research | - | - |
|
PubMed Central
Data from PMC are received and updated monthly.
|
- | - |
Provide sufficient details of any financial or non-financial competing interests to enable users to assess whether your comments might lead a reasonable person to question your impartiality. Consider the following examples, but note that this is not an exhaustive list:
Sign up for content alerts and receive a weekly or monthly email with all newly published articles
Already registered? Sign in
The email address should be the one you originally registered with F1000.
You registered with F1000 via Google, so we cannot reset your password.
To sign in, please click here.
If you still need help with your Google account password, please click here.
You registered with F1000 via Facebook, so we cannot reset your password.
To sign in, please click here.
If you still need help with your Facebook account password, please click here.
If your email address is registered with us, we will email you instructions to reset your password.
If you think you should have received this email but it has not arrived, please check your spam filters and/or contact for further assistance.
Comments on this article Comments (0)