Key vocab

Key functions / packages

Predict pI for an Selenocysteine and Pyrrolysine

Amino acids can be described by a number of chemical properties. This information is fairly easy to locate for the 20 standard proteinogenic amino acids coded by the standard codon table. For other amino acids it can be very difficult to find this information.

A key question in the study of the origin of the universal genetic code is: of the hundreds of amino acids that occur on earth, why does the codon table code for the 20 that it does? That is, why isn’t life based on a different set of 20 amino acids?

Many studies compare and contrast the chemical properties of the 20 proteinogenic amino acids with other amino acids, such as the amino acids that occur meteorites. Unfortunately, these studies rarely publish their data. Moreover, it seems like there is not always experimentally derived data on the chemical properties of non-proteinogenic amino acids. Authors’ therefore use chemistry modeling software to predict the chemical attributes non-standard amino acids.

In this portfolio assignment you will build a simple regression model (line of best fit) using the lm() function, then use the molecular weight of non-standard amino acids to predict other chemical values.

You’ll then assess whether this model is likely to be a very good one for predicting pI.

The assignment consists only of instructions - no code! To complete this assignment, you will need to gather the necessary data and code from recent assignments to make a self-sufficient script to carry out the following analysis.

Step 0: Load the necessary packages

# packages 
library(ggplot2)
library(ggpubr)
library(scatterplot3d)
library(pander)

Step 1: Make a dataframe

First, we need a dataframe to hold data on the 20 proteinogenic amino acids. Make a dataframe with the following columns:

  1. amino acid codes
  2. molecular weight
  3. pI (the vector is not called this though in my code code - make sure you set up your code properly)

Call the dataframe habk.df2.

# vectors
aa.codes <- c('A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L','M','N',
                'P','Q','R','S','T','V','W','Y')
aa.mw <- c(89, 121, 133, 146, 165, 75, 155, 131, 146, 131,
               149, 132, 115, 147, 174, 105, 119, 117, 204, 181)
pI <- c(6,5.07,2.77,3.22,5.48,5.97,7.59,6.02,9.74,5.98, 5.74, 5.41, 6.3, 
            5.65, 10.76, 5.68, 6.16, 5.96, 5.89, 5.66)


# make dataframe
habk.df2 <- data.frame(aa.codes = aa.codes,
                       aa.mw = aa.mw,
                       pI = pI)

Display the finished dataframe with pander::pander()

# show table with pander
pander(habk.df2)
aa.codes aa.mw pI
A 89 6
C 121 5.07
D 133 2.77
E 146 3.22
F 165 5.48
G 75 5.97
H 155 7.59
I 131 6.02
K 146 9.74
L 131 5.98
M 149 5.74
N 132 5.41
P 115 6.3
Q 147 5.65
R 174 10.76
S 105 5.68
T 119 6.16
V 117 5.96
W 204 5.89
Y 181 5.66

Step 2: Plot pI versus molecular weight

Use ggpubr to make a plot of the pH at the isoelectric point (pI) versus molecular mass.

  1. Plot pI on one axis and molecular weight on the other.
  2. Include regression line (line of best fit)
  3. Include other information, including a confidence intervals (CI), confidence ellipse, and the correlation coefficient
  4. Have each point labeled with the 1-letter amino acid code.
  5. Include relevant x-axis and y-axis labels and a title.
# make plot
ggscatter( y = "pI", x = "aa.mw",
          xlab = "Molecular Mass (Daltons)",
          ylab = "pI (pH at Isoelectric Point)",
          add = "reg.line", ellipse = TRUE,
          cor.coef = TRUE, label = "aa.codes", conf.int = TRUE, 
          data = habk.df2)
## `geom_smooth()` using formula 'y ~ x'

Step 3: Determine MW for Selenocysteine and Pyrrolysine

Look up the molecular mass of the non-standard amino acids Selenocysteine and Pyrrolysine from https://en.wikipedia.org/wiki/Amino_acid

Place this information in objects called selenocysteine.MW and pyrrolysine.MW

# molecular weights
selenocysteine.MW <- 168.1
pyrrolisyine.MW <- 255.3

Selenocysteine and Pyrrolysine are amino acids that can be added to proteins by re-coding stop codons. (On an exam you should recognize these as amino acids that occur in proteins via this mechanism and be able to answer simple multiple choice questions about them.)

For more information see:

Scitable: https://www.nature.com/scitable/topicpage/an-evolutionary-perspective-on-amino-acids-14568445/

Wikipedia: https://en.wikipedia.org/wiki/Selenocysteine https://en.wikipedia.org/wiki/Pyrrolysine

Step 4: Estimate pI for Selenocysteine (U) and Pyrrolysine (O)

Build a regression model using lm() to determine the y = m*x + b equation for the line in your scatterplot.

A simple regression model has a slope (m) and an intercept and can be set up as an equation

y = m*x + b

In this situation the model would be

pI = m*MW.da + b

“b” is called the “intercept” in R output.

When you build a regression model you can access output that looks like this: (Intercept) MW.da 2.30536091 1.01131509

Which is interpreted as

(Intercept) MW.da b m

If the numbers above were accurate (they aren’t), then our equation would be pI = 1.01131509*MW.da + 2.30536091

First, build the regression model (aka linear model)

# Build model
lm(aa.mw ~ pI)
## 
## Call:
## lm(formula = aa.mw ~ pI)
## 
## Coefficients:
## (Intercept)           pI  
##     115.898        3.445

Next get the coefficients for the model (the m of m*x, and the b).

#get coefficients / parameters for equation
lm.b <- 115.898
lm.m <- 3.445

Use the equation y = m*x + b to estimate the pI of selenocystein an pyrrolysine.

# Estimate pI for Se and Py
estimate.pI.Se <- lm.m * selenocysteine.MW + lm.b
estimate.pI.Py <- lm.m * pyrrolisyine.MW + lm.b

Step 5: Make a table of summary information for U and O

Make a table of your results, including the following columns

  1. amino.acid: full amino acid names for U and O
  2. three.letter: 3-letter code for the amino acids (given on Wikipedia)
  3. one.letter: 1-letter code for the amino acids (given on Wikipedia)
  4. molecular.weight: from wikipedia
  5. pI.estimate: from your regression equation
# make table
amino.acid <- c("Selenocysteine", "Pyrrolysine")
three.letter <- c("Sec", "Pyl")
one.letter <- c("U", "O")
molecular.weight <- c(selenocysteine.MW, pyrrolisyine.MW)
pI.estimate <- c(estimate.pI.Se, estimate.pI.Py)

habk.df2.addition <- data.frame(amino.acid = amino.acid,
                               three.letter = three.letter,
                               one.letter = one.letter,
                               molecular.weight = molecular.weight,
                               pI.estimate = pI.estimate)

Display the dataframe with the function pander() from the pander package.

# make table

pander(habk.df2.addition)
amino.acid three.letter one.letter molecular.weight pI.estimate
Selenocysteine Sec U 168.1 695
Pyrrolysine Pyl O 255.3 995.4

Step 6: Do you think your predictions of pI will be very accurate?

Using R, determine the correlation coefficient for molecular mass versus pI. Is this a very strong correlation coefficient? No this is not a very strong coefficient

# correlation coefficient
MW.pI.corr <- cor.test(habk.df2$pI, habk.df2$aa.mw)
MW.pI.corr
## 
##  Pearson's product-moment correlation
## 
## data:  habk.df2$pI and habk.df2$aa.mw
## t = 0.85448, df = 18, p-value = 0.4041
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2685455  0.5885382
## sample estimates:
##       cor 
## 0.1974383

Build a dataframe with all 8 variables used in Higgs and Attwood Chapter 2. Add molecular weight to this (which wasn’t included in their original table).

Determine the variable which has the highest correlation with pI.

First, make a datafrome with volume, bulk etc, along with molecular weight.

# data
vol <- c(67, 86, 91, 109, 135, 48, 118, 124, 135, 124, 124, 96, 90,
            114, 148, 73, 93, 105, 163, 141)

bulk  <-c(11.5, 13.46, 11.68, 13.57, 19.8, 3.4, 13.69, 21.4, 15.71, 21.4,
         16.25, 12.28,17.43, 14.45, 14.28, 9.47, 15.77, 21.57, 21.67,
         18.03)

pol <-c(0,1.48,49.7,49.9,0.35,0,51.6,0.13,49.5,0.13,
          1.43,3.38,1.58,3.53,52,1.67,1.66,0.13,2.1,1.61)

isoelec <-c(6,5.07,2.77,3.22,5.48,5.97,7.59,6.02,9.74,5.98, 5.74, 5.41,
            6.3, 5.65, 10.76, 5.68, 6.16, 5.96, 5.89, 5.66)
H2Opho.34 <-c(1.8,2.5,-3.5,-3.5,2.8,-0.4,-3.2,4.5, -3.9, 3.8, 1.9, -3.5,
              -1.6, -3.5, -4.5, -0.8, -0.7, 4.2, -0.9, -1.3)
H2Opho.35 <-c(1.6,2,-9.2,-8.2,3.7,1,-3,3.1,-8.8,2.8,3.4,-4.8,
              -0.2,-4.1,-12.3,0.6,1.2,2.6,1.9,-0.7)
saaH2O  <-c(113, 140,151,183,218,85,194,182,211,180, 204, 158, 143, 189, 
            241, 122, 146,160,259,229)
faal.fold <-c(0.74,0.91,0.62,0.62,0.88,0.72,0.78,0.88,0.52, 0.85, 0.85,
              0.63, 0.64,0.62,0.64,0.66,0.7, 0.86, 0.85, 0.76)

# Use aa.codes and aa.mw from above data

# make dataframe
habk.df3 <- data.frame(MW.da = aa.mw, Vol = vol, Bulk = bulk,
                       Polarity = pol, pI = isoelec,
                       Hyd.1 = H2Opho.34, Hyd.2 = H2Opho.35,
                       Surface_Area = saaH2O, Fract_area = faal.fold)

pander(habk.df3)
MW.da Vol Bulk Polarity pI Hyd.1 Hyd.2 Surface_Area Fract_area
89 67 11.5 0 6 1.8 1.6 113 0.74
121 86 13.46 1.48 5.07 2.5 2 140 0.91
133 91 11.68 49.7 2.77 -3.5 -9.2 151 0.62
146 109 13.57 49.9 3.22 -3.5 -8.2 183 0.62
165 135 19.8 0.35 5.48 2.8 3.7 218 0.88
75 48 3.4 0 5.97 -0.4 1 85 0.72
155 118 13.69 51.6 7.59 -3.2 -3 194 0.78
131 124 21.4 0.13 6.02 4.5 3.1 182 0.88
146 135 15.71 49.5 9.74 -3.9 -8.8 211 0.52
131 124 21.4 0.13 5.98 3.8 2.8 180 0.85
149 124 16.25 1.43 5.74 1.9 3.4 204 0.85
132 96 12.28 3.38 5.41 -3.5 -4.8 158 0.63
115 90 17.43 1.58 6.3 -1.6 -0.2 143 0.64
147 114 14.45 3.53 5.65 -3.5 -4.1 189 0.62
174 148 14.28 52 10.76 -4.5 -12.3 241 0.64
105 73 9.47 1.67 5.68 -0.8 0.6 122 0.66
119 93 15.77 1.66 6.16 -0.7 1.2 146 0.7
117 105 21.57 0.13 5.96 4.2 2.6 160 0.86
204 163 21.67 2.1 5.89 -0.9 1.9 259 0.85
181 141 18.03 1.61 5.66 -1.3 -0.7 229 0.76

Then make a correlation matrix from the data. Save it to an object called corr.mat

# make correlation matrix
corr.mat <- cor(habk.df3, habk.df3)

Round off the correlation matrix to 1 decimal [lace]

# round off
corr.mat.rounded <- round(corr.mat, digits = 1)

Display the matrix with pander::pander()

# display with pander()
pander(corr.mat.rounded)
Table continues below
  MW.da Vol Bulk Polarity pI Hyd.1 Hyd.2
MW.da 1 0.9 0.6 0.3 0.2 -0.3 -0.2
Vol 0.9 1 0.7 0.2 0.4 -0.1 -0.2
Bulk 0.6 0.7 1 -0.2 0.1 0.4 0.3
Polarity 0.3 0.2 -0.2 1 0.3 -0.7 -0.9
pI 0.2 0.4 0.1 0.3 1 -0.2 -0.3
Hyd.1 -0.3 -0.1 0.4 -0.7 -0.2 1 0.8
Hyd.2 -0.2 -0.2 0.3 -0.9 -0.3 0.8 1
Surface_Area 1 1 0.6 0.3 0.3 -0.2 -0.2
Fract_area 0.1 0.2 0.5 -0.5 -0.2 0.8 0.8
  Surface_Area Fract_area
MW.da 1 0.1
Vol 1 0.2
Bulk 0.6 0.5
Polarity 0.3 -0.5
pI 0.3 -0.2
Hyd.1 -0.2 0.8
Hyd.2 -0.2 0.8
Surface_Area 1 0.1
Fract_area 0.1 1

Which variable has the strongest (most positive OR most negative) correlation with pI? (if possible write this with code, otherwise just state what it is).

# find maximum
round(corr.mat.rounded[5, ], digit = 1)
##        MW.da          Vol         Bulk     Polarity           pI        Hyd.1 
##          0.2          0.4          0.1          0.3          1.0         -0.2 
##        Hyd.2 Surface_Area   Fract_area 
##         -0.3          0.3         -0.2
text <- "maximum correlation with pI is Volume with a correlation of 0.3"
print(text)
## [1] "maximum correlation with pI is Volume with a correlation of 0.3"