Accurately estimating the thermodynamic cost of a mutation is a building block of protein engineering and design. This task is especially tricky for membrane proteins because the calculations must account for the lipid bilayer. In this tutorial, we will walk through the protocol for estimating the ∆∆G for lipid facing positions using RosettaMP and the franklin2019 energy function.
As an example, we will examine mutations in the integral membrane enzyme PagP. PagP is a beta-barrel protein that transfers a palmitoyl group fron the sn-1 position of a glycerophospholipid to the endotoxin of lipopolysacharide (LPS). The enzyme provides bacterial resistance to host immune defenses such as antimicrobial pepetides (Guo et al. 1998; Kawasaki et al. 2004). Recently, Marx & Fleming measured the energetic cost of making mutations at the V111 position on PagP (Marx & Fleming, 2017). Here, we will perform the same set of mutations with Rosetta and compare with the experimental values.
The first step is to initialize PyRosetta and load the protein of interest. In this tutorial, we will use PagP (PDB 3GP6). The starting structure is from the Orientations of Proteins in Membranes database (https://opm.phar.umich.edu/) which provides spatial arrangements of membrane proteins in the lipid bilayer.
!pip install pyrosettacolabsetup
import pyrosettacolabsetup; pyrosettacolabsetup.install_pyrosetta()
import pyrosetta; pyrosetta.init()
from pyrosetta import *
init( extra_options="-mp:lipids:has_pore false")
Make sure you are in the right directory for accessing the .pdb
files:
cd google_drive/MyDrive/student-notebooks/
#cd google_drive/MyDrive/student-notebooks/
pose = pose_from_pdb( "inputs/3gp6_A.pdb")
Then, initialize the protein in the membrane using AddMembraneMover.
Here, the protein is already oriented in the bilayer so we can estimate the transmembrane spans from the structure and orientation. Thus, we use the from_structure
option to initialize the spanning topology.
from pyrosetta.rosetta.protocols.membrane import *
add_memb = AddMembraneMover("from_structure")
add_memb.apply(pose)
Next, we will compute the ∆∆G for several point mutations in PagP. In the Marx & Fleming experiment, position V111 was first mutated to alanine. Therefore, we will create this variant first. We will use the mutate_residue
function from the predict_ddG
PyRosetta module included in this package. In this tutorial, we will use a repack radius of 8.0 Å.
An important note - Pyrosetta residue numbering may differ from the PDB numbering because PyRosetta requires continuous numbering for calculations. Here, the PyRosetta residue number for V111 is 104.
from additional_scripts import predict_ddG
# Create a franklin2019 energy function
sfxn = create_score_function("franklin2019")
# Repack and score the native conformation
reference_pose = predict_ddG.mutate_residue(pose, 104, "A", 8.0, sfxn)
To demonstrate the ΔΔG calculation, we will now compute the energetic cost of mutating alanine to tryptophan.
# Score the alanine reference pose
score_A111 = sfxn.score(reference_pose)
# Repack and score the L111 conformation
pose_W111 = predict_ddG.mutate_residue(pose, 104, "W", 8.0, sfxn)
score_W111 = sfxn.score(pose_W111)
# Compute the ddG of mutation as mutant_score - native_score (final-initial
ddG = score_W111 - score_A111
print(ddG)
The ΔΔG for mutating alanine to tryptophan at position 111 is -1.84 Rosetta Energy Units (REU). A Rosetta Energy Unit is an arbitrary unit for the Rosetta energy function. Next, we would like to compute the ΔΔG for mutating alanine to all 19 canonical amino acids. To do so, we will generalize the code above into a function for easy calculations of multiple single point mutations.
def compute_ddG(pose, native_res, site_no, mutant_res, sfxn):
"""A function for computing the ddG of single point mutations
Example:
$ compute_ddG(pose, "V", 49, "A", sfxn)
Arguments:
- pose = Object containing the coordinates for the biomolecular system
- native_res = Native amino acid
- site_no = Host site amino acid position
- mutant_res = Mutant amino acid
= sfxn = Score function object
"""
repacked_native = predict_ddG.mutate_residue(pose, site_no, native_res, 8.0, sfxn)
native_score = sfxn.score(repacked_native)
repacked_mutant = predict_ddG.mutate_residue(pose, site_no, mutant_res, 8.0, sfxn)
mutant_score = sfxn.score(repacked_mutant)
ddG = mutant_score - native_score
return ddG
Now, we will write a loop that computes the ΔΔG for all canonical amino acids and store the results in a python dictionary.
# List of canonical amino acid one-letter codes
amino_acids = [ 'A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y' ]
# Initialize an empty dictionary to store the data
ddG_data = {}
# Loop through all amino acids
for aa in amino_acids:
ddG = compute_ddG(reference_pose, 'A', 104, aa, sfxn)
ddG_data[ aa ] = ddG
The next step is to compare the ΔΔG predictions to the experimentally measured values. The experimental data for Marx & Fleming are located in inputs
in a file called PagP_Marx_Fleming_set.dat
. We will parse the file and then import these values into a |dictionary.
# Read contents of file into a list (file format is space delimited)
with open( 'inputs/PagP_Marx_Fleming_set.dat', 'rt' ) as f:
data = f.readlines()
data = [ x.strip() for x in data ]
data = [ x.split(' ') for x in data ]
# Convert the list into a dictionary
exp_ddG_data = {}
for i in range(1, len(data)):
exp_ddG_data[ data[i][2] ] = float(data[i][3])
We now convert the dictionary format to numpy arrays that are compatible with analysis. Here, we will compute the correlation coefficient and make a scatterplot of the experimentally measured vs. predicted values.
import numpy as np
mutations = np.asarray( ddG_data.keys() )
ddG_values = np.asarray( list(ddG_data.values()) )
exp_values = np.asarray( list(exp_ddG_data.values()) )
# Compute the correlation coefficient
corr = np.corrcoef( exp_values, ddG_values )
print(corr[0,1])
We initially find that the correlation coefficient is low (0.376). We will want to find any outliers in the dataset that are lowering this value.
def find_outliers(x):
upper_quartile = np.percentile(x, 75)
lower_quartile = np.percentile(x, 25)
IQR = (upper_quartile - lower_quartile)
quartile_set = (lower_quartile - IQR, upper_quartile + IQR)
for y in x.tolist():
if (y < quartile_set[0]) or (y > quartile_set[1]):
print(y)
find_outliers(ddG_values)
Here, the ΔΔG value for proline is an outlier. We will investigate this more later. For now, we will remove it from the set and then recompute the correlation coefficient.
# Proline is the 13th, amino acid of 20
exp_data_no_P = []
pred_data_no_P = []
for i in range(0, 20):
if ( i != 12 ):
exp_data_no_P.append( list(exp_ddG_data.values())[i] )
pred_data_no_P.append( list(ddG_data.values())[i] )
exp_data_no_P = np.asarray( exp_data_no_P )
pred_data_no_P = np.asarray( pred_data_no_P )
corr = np.corrcoef( exp_data_no_P, pred_data_no_P )
print(corr[0,1])
The new value of R = 0.692 is much more encouraging! Next, we will visualized the predicted vs. experimentally measured values with a scatterplot.
import matplotlib.pyplot as plt
plt.axhline(y=0, color='k', linestyle='-')
plt.axvline(x=0, color='k', linestyle="-")
plt.ylim([-6,6])
plt.xlim([-6,6])
# compute the best fit line
from numpy.polynomial.polynomial import polyfit
b, m = polyfit(exp_data_no_P, pred_data_no_P, 1)
x = np.linspace(-6, 6, num=50)
plt.plot(x, b + m * x, color='r', linestyle='-')
# plot the data
plt.scatter(exp_data_no_P, pred_data_no_P)
Finally, we would like to use the models to learn why some mutations stabilized PagP, whereas other side chains did not. Of course, we need a metric for identifying the most confident predictions, especially since the correlation coefficient is not perfect. To do so, we will compute the residuals from the line of best fit and set an empirical cutoff of 1.5 REU.
import seaborn as sns
sns.set(style="whitegrid")
resid = sns.residplot(x=exp_data_no_P, y=pred_data_no_P, color="b")
resid.set_ylabel("Residual")
resid.set_xlabel("Exp (kcal/mol)")
print(exp_data_no_P)
Here, we find five ∆∆G values that are predicted outside of the cutoff: Glycine, Leucine, Valine, Tryptophan, and Threonine. Next, we will use this information to hypothesize a mechanism for a reasonable prediction (lysine) and rationalize incorrect predictions for proline and leucine. The first step is to quantify which energy components make the largest contribution to the overall ∆∆G of mutation. To do so, we will write a function that can extrapolate this information from the energy function.
# Store models of mutated PagP proteins
mutant_tyr = predict_ddG.mutate_residue(pose, 104, "Y", 8.0, sfxn)
mutant_lys = predict_ddG.mutate_residue(pose, 104, "K", 8.0, sfxn)
mutant_leu = predict_ddG.mutate_residue(pose, 104, "L", 8.0, sfxn)
def get_energy_components( native_pose, mutated_pose, sfxn):
# Extract & parse scores
tmp_native = native_pose.energies().total_energies().weighted_string_of( sfxn.weights() )
tmp_mutant = mutated_pose.energies().total_energies().weighted_string_of( sfxn.weights() )
array_native = list(filter( None, tmp_native.split(' ') ))
array_mutant = list(filter( None, tmp_mutant.split(' ') ))
# Pull out only the scores from these arrays
native_scores = []
for i in range( len(array_native) ):
if ( i % 2 != 0 ):
native_scores.append( float( array_native[i] ) )
mutant_scores = []
for i in range( len(array_mutant) ):
if ( i % 2 != 0 ):
mutant_scores.append( float( array_mutant[i] ) )
# Calculate ddG of individual components
ddGs = []
for i in range( len( mutant_scores ) ):
ddG_component = mutant_scores[i] - native_scores[i]
ddGs.append( round( ddG_component, 3 ) )
# Get labels
labels = []
for i in range( len(array_native) ):
if ( i % 2 == 0 ):
labels.append( array_native[i].translate(':').strip(":") )
return labels, ddGs
# Compute the ddG breakdown
labels, tyr_ddGs = get_energy_components( reference_pose, mutant_tyr, sfxn )
labels, lys_ddGs = get_energy_components( reference_pose, mutant_lys, sfxn )
labels, leu_ddGs = get_energy_components( reference_pose, mutant_leu, sfxn )
print(labels)
Next, we will make a bar graph to visualize the contributions of each energy component for the ddG of these three single point mutations.
x = np.arange(len(labels))
width = 1 # the width of the bars
# Plot for tyrosine
fig, ax = plt.subplots()
rects1 = ax.bar(x, tyr_ddGs, width, label='Y')
ax.set_ylabel('∆∆G (REU)')
ax.set_title('Contributions to the ∆∆G of Mutation for A104Y')
ax.set_xticks(x)
ax.set_xticklabels(labels, rotation='vertical')
ax.legend()
fig.tight_layout()
# Plot for lysine
fig, ax = plt.subplots()
rects2 = ax.bar(x, lys_ddGs, width, label='K')
ax.set_ylabel('∆∆G (REU)')
ax.set_title('Contributions to the ∆∆G of Mutation for A104K')
ax.set_xticks(x)
ax.set_xticklabels(labels, rotation='vertical')
ax.legend()
fig.tight_layout()
# Plot for lysine
fig, ax = plt.subplots()
rects3 = ax.bar(x, leu_ddGs, width, label='L' )
ax.set_ylabel('∆∆G (REU)')
ax.set_title('Contributions to the ∆∆G of Mutation for A104L')
ax.set_xticks(x)
ax.set_xticklabels(labels, rotation='vertical')
ax.legend()
fig.tight_layout()
Finally, we will export the model files so that we can visualize them in PyMOL.
reference_pose.dump_pdb( "PagP_ala_ref.pdb" )
mutant_tyr.dump_pdb( "PagP_A104Y.pdb" )
mutant_lys.dump_pdb( "PagP_A104K.pdb" )
mutant_leu.dump_pdb( "PagP_A104L.pdb" )