# ddG.ptm()

### Description

Builds a PDB model of the modified protein and computes the corresponding change in stability.

### Arguments

`pdb` the 4-letter identifier of a PDB structure or the path to a PDB file.

`ch `a letter identifying the chain of interest.

`pos` the position, in the primary structure, of the residue to be modified.

`ptm` the post-translational modification to be considered. It should be one among: ‘pSer’, ‘pThr’, ‘pTyr’, ‘MetO-Q’, ‘MetO-T’.

`dir` indicates the direction of the PTM reaction: either forward (‘f’), or backward (‘b’).

`pH` a numeric value between 0 and 14.

### Value

The function computes and returns the ΔΔG (kcal/mol) for the requested modification, defined as ΔΔG = ΔGmodified – ΔGunmodified, where ΔG is the Gibbs free energy for the folding of the protein from its unfolded state. Thus, a positive value means a destabilizing effect, and vice versa. A PDB model containing the modified target is saved in the current directory.

### References

Schymkowitz et al (2005) Nucl. Ac. Res. 33:W382-W388.

### Details

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, ΔΔG (kcal/mol), for the folding process. That is,

$\Delta \Delta G = \Delta G_{mut} - \Delta G_{wt}$

Where ΔG is the change in free energy corresponding to the folding of the protein from its unfolded state, either in the wild type (wt) or mutated (mut) version of the protein. Thus, a ΔΔG > 0 implies a decrease in stability (destabilizing mutation) while a ΔΔG < 0 is interpreted as an increase in stability (stabilizing mutation).

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:

$\Delta \Delta G = \Delta G_{modified} - \Delta G_{unmodified}$

Where $\Delta G_{modified}$ is the folding free energy of the protein containing the modified target residue. When interpreting results arising from the use of foldx.mut() or ddG.ptm() (both based on the FoldX suit) a practical rule of thumb may be as follows:

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.

For instance, Alpha-1-antitrypsin (P01009) 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. Next, 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"))
```
```## [1] "THE FILE 3cwm_MetO-Q63.pdb HAS BEEN SAVED IN THE CURRENT DIRECTORY"
```
```## [1] "THE FILE 3cwm_MetO-Q220.pdb HAS BEEN SAVED IN THE CURRENT DIRECTORY"
```
```## [1] "THE FILE 3cwm_MetO-Q221.pdb HAS BEEN SAVED IN THE CURRENT DIRECTORY"
```
```## [1] "THE FILE 3cwm_MetO-Q226.pdb HAS BEEN SAVED IN THE CURRENT DIRECTORY"
```
```## [1] "THE FILE 3cwm_MetO-Q242.pdb HAS BEEN SAVED IN THE CURRENT DIRECTORY"
```
```## [1] "THE FILE 3cwm_MetO-Q351.pdb HAS BEEN SAVED IN THE CURRENT DIRECTORY"
```
```## [1] "THE FILE 3cwm_MetO-Q358.pdb HAS BEEN SAVED IN THE CURRENT DIRECTORY"
```
```## [1] "THE FILE 3cwm_MetO-Q374.pdb HAS BEEN SAVED IN THE CURRENT DIRECTORY"
```
```## [1] "THE FILE 3cwm_MetO-Q385.pdb HAS BEEN SAVED IN THE CURRENT DIRECTORY"
```
```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)

```

We can observe that the modification of most of the residues from alpha-1-antitrypsin has a destabilizing effect, with the exception of M351 and M358, whose modification has a neutral to slightly stabilizing effect.

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])
}
```

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(phos\$DG) ~ phos\$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 protein stability. Since ddG.ptm() has built and saved PDB models for each modification, 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.2  TRUE    Y
## 10  94   Y  pTyr  11.4 TRUE    Y
## 11  98   S  pSer  11.5 TRUE    Y
## 12  99   T  pThr  8.2  TRUE    Y
## 15 122   S  pSer  11.0 TRUE    Y
## 32 237   T  pThr  5.9  TRUE    Y
## 34 246   T  pThr  12.0 TRUE    Y
```

Note how phosphorylation of Thr246 in a single subunit is strongly destabilizing (ΔΔG = 12.0 kcal/mol). It has been reported that phosphorylation of Thr246 by PKCδ inhibits GAPDH-driven mitophagy. If we wish to explore the effect of this modification in a second subunit, we can repeat 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_Repair.pdb', 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"
```