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(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.1        <-c('A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L','M','N',
                'P','Q','R','S','T','V','W','Y')
mw <-c(89,121,133,146,165,75,155,131,146,131,149,132,115,147,174,105,119,117,204,181)
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)

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

Display the finished dataframe with pander::pander()

# show table with pander
pander(habk.df2)
aa.1 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(x ="mw",        
          y ="pI" , add = "reg.line",cor.coef = TRUE,ellipse = TRUE,conf.int=TRUE,
          label ="aa.1" , 
          ylab ="pI" , #TODO
          xlab ="MW" , #TODO
          title ="Prediction of pI from MW",
          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.064
pyrrolysine.MW<-255.313

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
lobf <- lm(pI ~ mw, data = habk.df2)
lobf_slope_int <- coef(lobf)
lobf_slope_int
## (Intercept)          mw 
##  4.50516096  0.01131509

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

#get coefficients / parameters for equation
slope<-lobf_slope_int[2]
intercept<-lobf_slope_int[1]

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

# Estimate pI for Se and Py
Se.pI<-slope*selenocysteine.MW+intercept
Py.pI<-slope*pyrrolysine.MW+intercept

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 an d 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
U_O_summary<-data.frame(amino.acid=c("Selenocystein","Pyrrolysine."),
                       three.letter=c("Sec","Pyl"),
                       one.letter=c("U","O"),
                       molecular.weight=c(selenocysteine.MW,pyrrolysine.MW),
                       pI.estimate=c(Se.pI,Py.pI)
                       
                       )

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

# make table

pander(U_O_summary
       )
amino.acid three.letter one.letter molecular.weight pI.estimate
Selenocystein Sec U 168.1 6.407
Pyrrolysine. Pyl O 255.3 7.394

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?

# correlation coefficient
slope
##         mw 
## 0.01131509

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)
# made dataframe
habk.df <- data.frame( weight=mw,
                      Vol = vol, Bulk = bulk, Polarity = pol,
                      pI = isoelec, Hyd.1 = H2Opho.34,
                      Hyd.2 = H2Opho.35, Surface_area = saaH2O,
                      Fract_area = faal.fold
                      )

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

# make correlation matrix
corr.mat=cor(habk.df)
corr.mat
##                  weight         Vol        Bulk   Polarity          pI
## weight        1.0000000  0.93373477  0.55416842  0.2904128  0.19743827
## Vol           0.9337348  1.00000000  0.72741358  0.2358472  0.36251434
## Bulk          0.5541684  0.72741358  1.00000000 -0.2022089  0.08225406
## Polarity      0.2904128  0.23584720 -0.20220894  1.0000000  0.26567728
## pI            0.1974383  0.36251434  0.08225406  0.2656773  1.00000000
## Hyd.1        -0.2714765 -0.08022302  0.44390869 -0.6687590 -0.20412730
## Hyd.2        -0.2452480 -0.15858483  0.31852492 -0.8535320 -0.26083174
## Surface_area  0.9676090  0.98694808  0.63576635  0.2883267  0.34675188
## Fract_area    0.1090068  0.18249961  0.48621057 -0.5292383 -0.17986683
##                    Hyd.1      Hyd.2 Surface_area Fract_area
## weight       -0.27147648 -0.2452480    0.9676090  0.1090068
## Vol          -0.08022302 -0.1585848    0.9869481  0.1824996
## Bulk          0.44390869  0.3185249    0.6357664  0.4862106
## Polarity     -0.66875899 -0.8535320    0.2883267 -0.5292383
## pI           -0.20412730 -0.2608317    0.3467519 -0.1798668
## Hyd.1         1.00000000  0.8495286   -0.1813011  0.8409746
## Hyd.2         0.84952863  1.0000000   -0.2258995  0.7875690
## Surface_area -0.18130107 -0.2258995    1.0000000  0.1234664
## Fract_area    0.84097463  0.7875690    0.1234664  1.0000000

Round off the correlation matrix to 1 decimal [lace]

# round off
corr.mat=round(corr.mat,1)

Display the matrix with pander::pander()

# display with pander()
pander(corr.mat)
Table continues below
  weight Vol Bulk Polarity pI Hyd.1 Hyd.2
weight 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
weight 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).

the Vol has the strongest correlation with pI

# find maximum
#max(corr.mat$pI)