Thermodynamic stability is a fundamental property of proteins that influences protein structure, function, expression, and solubility. Not surprisingly, prediction of the effect of point mutations on protein stability is a topic of great interest in many areas of biological research, including the field of protein post-translational modifications. Consequently, a wide range of strategies for estimating protein energetics has been developed. In the package ptm we have implemented two popular computational approaches for prediction of the effect of amino acid changes on protein stability.
On the one hand, foldx.mut() implements FoldX (buildmodel and positionscan methods), a computational approach that uses a force field method and although it has been proved to be satisfactorily accurate, it is also a time-consuming method. On the other hand, I-Mutant is a method based on machine-learning and it represents an alternative much faster to FoldX. The ptm function that implements this last approach is imutant().
Although using very different strategies, both functions assess
the thermodynamic stability effect of substituting a single amino acid residue in the protein of interest. To do that, both function perform an estimation of the change in Gibbs free energy, 
Where 
In the particular case of the function ddG.ptm(), a PDB model containing the modified residue is built and saved in the current directory. This, model is used to compute:
Where 
Currently, the target residues for the modifications implemented into ddG.ptm() are Ser, Thr and Tyr, which are phosphoacceptors and therefore, they are modified to pSer, pThr and pTyr. Another modification addressed by ddG.ptm() is the oxidation of Met to methionine sulfoxide (MetO). However, in the latter case the PDB model containing the modified residue will introduce, in the place of the target Met, either Gln or Thr instead of MetO (since MetO is not an option offered by FoldX). Oxidation of methionine adds polarity to an otherwise apolar side chain. More concretely, the side chain hydrophobicity index decreases from 0.738 (Met) to 0.238 (MetO) after the sulfoxidation reaction. Not surprisingly, polar amino acids such as Gln (hydrophobicity index of 0.251) and Thr (0.450) can sometimes mimic the sulfoxidized state of a protein. Furthermore, we also have evidence that nature may also employ this “trick”, but in reverse: evolving methionine sites from glutamine and threonine residues.
ddG.ptm dependency
The function ddG.ptm() requires computations that are carry out by an external program: FoldX. You can easily download the software (which requires free registration) here. After you have installed FoldX and added the binary location to PATH environment variable, you’re done! Nevertheless, if you are interested in the FoldX suite, you may wish to consult its manual.
ddG.ptm and Methionine Sulfoxidation
As a protein model we are going to use Alpha-1-antitrypsin (P01009), which presents 9 methionine residues in its mature form, the oxidation of two of which (M351 and M358) has been described to cause loss of anti-neutrophil elastase activity. We will assess the effect on stability of the modification of the 9 methionine residues:
seq <- get.seq('P01009', db = 'metosite')
met <- gregexpr('M', seq)[[1]]
meto <- lapply(met, function(x) ddG.ptm('3cwm', pos = x, 'A', ptm = "MetO-Q"))
meto <- as.numeric(unlist(meto))
barplot(meto, horiz = TRUE,
ylab = "Modified Met",
xlab = expression(paste(Delta, Delta, "G (kcal/mol)", sep = "")),
col = c(rep("red", 5), "blue", "blue", "red", "red"),
names.arg = met, cex.names = 0.6)
It is interesting to note that, while most of the methionine residues from alpha-1-antitrypsin destabilize the protein when oxidized, the effect is the opposite for the two methionine resides (M351 and M358) whose oxidation has been described to cause loss of anti-neutrophil elastase activity. We can observe that the modification of M351 and M358 has a neutral to slightly stabilizing effect.
ddG.ptm and Phosphorylation
Now, to illustrate the use of ddG.ptm() in the context of phosphorylation, we are going to use the human glyceraldehyde-3-phosphate dehydrogenase (GAPDH, PDB ID: 1U8F). To start, we will find out all the phosphorylatable sites in this protein:
sites <- ptm.scan("P04406")
psites <- sites[which(sites$p == TRUE), c(2,3,13)]
psites
## n aa reg ## 25 25 S NA ## 42 42 Y TRUE ## 45 45 Y NA ## 49 49 Y NA ## 51 51 S NA ## 52 52 T NA ## 59 59 T NA ## 75 75 T NA ## 83 83 S NA ## 94 94 Y TRUE ## 98 98 S TRUE ## 99 99 T TRUE ## 103 103 T NA ## 104 104 T NA ## 122 122 S TRUE ## 125 125 S NA ## 140 140 Y NA ## 143 143 S NA ## 148 148 S NA ## 151 151 S NA ## 153 153 T NA ## 154 154 T NA ## 176 176 T NA ## 177 177 T NA ## 182 182 T NA ## 184 184 T NA ## 187 187 T NA ## 192 192 S NA ## 210 210 S NA ## 211 211 T NA ## 229 229 T NA ## 237 237 T TRUE ## 241 241 S NA ## 246 246 T TRUE ## 255 255 Y NA ## 266 266 S NA ## 276 276 Y NA ## 284 284 S NA ## 312 312 S NA ## 314 314 Y NA ## 320 320 Y NA ## 321 321 S NA ## 333 333 S NA
As it can be observed in the table above, many of the phosphorylation sites described for GAPDH have been dectected in high-throughput experiments and therefore there is no evidence of their involvement in regulatory processes. To emphisize this fact, we are going to add a variable regf taking two possible values: Y (for regulatory sites) and N (for non-regulatory sites)
psites$regf <- NA
psites$regf[which(psites$reg == TRUE)] <- "Y"
psites$regf[is.na(psites$reg)] <- "N"
Now, with the help of ddG.ptm() we shall test the null hypothesis that the phosphorylation of regulatory and non-regulatory sites have the same effect on protein stability. So, we have to compute ΔΔG for all the PTMs:
ddG <- data.frame(pos = psites$n, res = psites$aa, ptm = NA, DDG = NA, reg = psites$reg, regf = psites$regf)
for (i in 1:nrow(psites)){
t <- psites$aa[i]
if (t == "S"){
ddG$ptm[i] <- 'pSer'
} else if (t == "Y"){
ddG$ptm[i] <- 'pTyr'
} else if (t == "T"){
ddG$ptm[i] <- 'pThr'
}
# ddG$DDG[i] <- ddG.ptm('./1u8f_Repair.pdb', 'O', pos = ddG$pos[i], ptm = ddG$ptm[i])
ddG$DDG[i] <- ddG.ptm('1U8F', 'O', pos = ddG$pos[i], ptm = ddG$ptm[i])
}
When the distribution of the ΔΔG values obtained is plotted:
boxplot(as.numeric(ddG$DDG) ~ ddG$regf, ylab = expression(paste(Delta, Delta, "G (kcal/mol)", sep ="")),
xlab = "Regulatory site?")
two conspicuous conclusions can be made: (i) most of the PTMs lead to destabilization, and (ii) it seems that regulatory phosphorylations are more destabilizing than their non-regulatory counterpart. In order to gain support for this last conclusion, let’s carry out a t-student test:
t.test(as.numeric(ddG$DDG) ~ ddG$regf, mu = 0, alt = "less", conf = 0.95, var.eq = F, paired = F)
## ## Welch Two Sample t-test ## ## data: as.numeric(ddG$DDG) by ddG$regf ## t = -2.38, df = 15.131, p-value = 0.01544 ## alternative hypothesis: true difference in means is less than 0 ## 95 percent confidence interval: ## -Inf -0.9430228 ## sample estimates: ## mean in group N mean in group Y ## 5.616452 9.190741
It seems that we can reject the null hypothesis that regulatory and non-regulatory have the same effect on the protein stability. Since, ddG.ptm() has built and saved the PDB models for each modification we have carried out, let’s examine some of them. To this end, let’s focus on the regulatory phosphorylation:
ddG[which(ddG$reg), ]
## pos res ptm DDG reg regf ## 2 42 Y pTyr 4.23529 TRUE Y ## 10 94 Y pTyr 11.3747 TRUE Y ## 11 98 S pSer 11.5082 TRUE Y ## 12 99 T pThr 8.21505 TRUE Y ## 15 122 S pSer 11.0302 TRUE Y ## 32 237 T pThr 5.94235 TRUE Y ## 34 246 T pThr 12.0294 TRUE Y
Note how phosphorylation of Thr246 in a single subunit is strongly destabilizing (ΔΔG = 12.03 kcal/mol). It has been reported that phosphorylation of Thr246 by PKCδ inhibits GAPDH-driven mitophagy. If we whish to explore the effect of this modification in a second subunit, we can repit the computation implemented by ddG.ptm(), but now passing as an argument the PDB file generated previously (the one incorporating pThr246 in one subunit):
ppT246OP <- ddG.ptm(pdb = './1u8f_pThr246.pdb', ch = 'P', pos = 246, ptm = 'pThr')
## [1] "THE FILE 1u8f_pThr246_pThr246.pdb HAS BEEN SAVED IN THE CURRENT DIRECTORY"
ppT246OP
## [1] "25.5512" ## attr(,"wild-type") ## [1] "THR" ## attr(,"position") ## [1] 246 ## attr(,"PTM") ## [1] "pThr" ## attr(,"units") ## [1] "kcal/mol"
Thus, the change in thermodynamic stability after phosphorylation of two subunits would be of ΔΔG = 12.03 + 25.55 = 37.58 kcal/mol). The figure below helps to rationalize this result.
A short note before finishing. Sometimes phosphoproteins can be crystallized and resolved (see Table 1 from the review Structural Basis for Control by Phosphorylation). Thus, if you have a PDB file containing phosphoresidues, you may be able to assess the effect of dephosphorylation on the protein stability:
rev <- ddG.ptm(pdb = '1JST', ch = 'A', pos = 160, ptm = 'pThr', dir = 'b')
## [1] "THE FILE 1JST_pThr160.pdb HAS BEEN SAVED IN THE CURRENT DIRECTORY"
rev
## [1] "2.92483" ## attr(,"wild-type") ## [1] "TPO" ## attr(,"position") ## [1] 160 ## attr(,"PTM") ## [1] "pThr" ## attr(,"units") ## [1] "kcal/mol"
ptm 