Introduction

#> [1] "!!! IMPORTANT !!!"
#> [1] "This Vignette has not been built completely on CRAN due to size limitations."
#> [1] "Please check the correct version here: "
#> [1] "https://jpquast.github.io/protti/articles/protein_structure_workflow.html"

This vignette will take you through all important tools protti provides for structural data analysis. These include functions that allow you to fetch information about protein structures both on an experimental as well as structural level. Additionally, structure contact maps can be created from distances between atoms. Lastly, B-factors* of structures can be manipulated to contain custom values. This gives you the ability to colour regions on structures as well as to apply continuous values associated with each amino acid as a colour gradient using biomolecular visualisation software such as PyMOL or ChimeraX or the R package r3dmol.

*B-factors or temperature factors are used in structural biology to describe the atomic displacement in the structure. In low resolution structures they might correspond to the flexibility of regions, while in high resolution structures they should not be used to infer information regarding flexibility.

In structural proteomics, protein structural changes can be determined using methods such as limited proteolysis coupled to mass spectrometry (LiP-MS), cross-linking mass spectrometry (XL-MS) and hydrogen/deuterium exchange mass spectrometry (HDX-MS). In case of LiP-MS, structural changes can be identified with peptide or precursor* centric proteomics in the form of peptide or precursor fold changes. Relating this information to positions within the 3D space of a protein can be very helpful to infer potential functional roles of those structural changes. If structural changes occur close to the active site of an enzyme, kinetics might be affected, if they occur at a protein-protein binding interface, the corresponding protein complex might have dissociated or assembled.

*Peptide precursors are the different charge states of a peptide (+/- modifications) and the actual molecular unit that is detected in the mass spectrometer.

Getting started

Before starting the structural analysis of your data you need to load protti and additional packages used for the analysis. As described on the main page of protti, the tidyverse package collection works well together with the rest of the package, by design. You can load packages with the library() function.

# Load packages
library(protti)
library(dplyr)
library(magrittr)
library(stringr)
library(tidyr)
library(ggplot2)

We will demonstrate the use of protti’s structure analysis functions based on a data set of peptides that were significantly changing in a LiP-MS experiment and thus point towards a structural change in the protein.

The functions used in this analysis can also be used in a different context than a structural proteomics experiment. Instead of peptides you can investigate the positioning of individual amino acids, or whole domains within a protein.

utils::data("ptsi_pgk")

# Load data
your_data <- read_protti("your_differential_abundance_data.csv")

# Input UniProt IDs
uniprot_ids <- unique(ptsi_pgk$pg_protein_accessions)

# Fetch UniProt information
uniprot_information <- fetch_uniprot(uniprot_ids = uniprot_ids,
                                     columns = c("sequence", "database(PDB)"))

# Add UniProt information and find peptide positions
ptsi_pgk_annotated <- ptsi_pgk %>% 
  left_join(uniprot_information, by = c("pg_protein_accessions" = "id")) %>% 
  find_peptide(protein_sequence = sequence, peptide_sequence = pep_stripped_sequence)

Fetching structural information

You can access structural information with protti in order to see if there are protein structures available for your protein, how they were acquired and if they fit your quality requirements.

In the previous step we already retrieved information about the availability of protein structures for your protein of interest using the fetch_uniprot() function. The column database_pdb contains PDB identifiers if there is a structure available. Since the identifiers are pasted together into one string we need to separate them into individual identifiers to provide them to the fetch_pdb() function. This can be achieved through a combination of the stringr function str_split() and the tidyr function unnest().

# Extract PDB IDs from UniProt information
ptsi_pgk_pdb_ids <- ptsi_pgk_annotated %>% 
  distinct(pg_protein_accessions, database_pdb) %>% 
  mutate(pdb_id = str_split(database_pdb, pattern = ";")) %>% 
  unnest(pdb_id) %>% 
  filter(pdb_id != "")

# Fetch pdb information
ptsi_pgk_pdb_information <- fetch_pdb(pdb_ids = unique(ptsi_pgk_pdb_ids$pdb_id))

The result of fetch_pdb() contains a multitude of information. These include for example the experimental method that was used to determine the structure (experimental_method and structure_method). Depending on the method, there are columns specific to them. For example the column ph_crystallisation contains information about the pH of the solution used for crystallisation if X-ray crystallography was used as a method, while the column ph_nmr contains equivalent information if the method was NMR. A structure resolution (resolution_combined) is for example only included if available, which is not the case for NMR structures.

Generally, the molecular information about the structure is divided into “polymer” and “nonpolymer”. While nonpolymer information concerns binding partners such as metal-ions, metabolites, drugs or other small molecules, the polymer information concerns the protein or polynucleotide.

You could for example only consider structures for your analysis that were determined by X-ray crystallography and that have a resolution of below 3 Å. In addition, you could select only the one structure that contains the longest sequence stretch for each of our proteins (the structure with the maximum length). This can be achieved in the following way:

filtered_structures <- ptsi_pgk_pdb_information %>% 
  filter(experimental_method == "X-ray",
         resolution_combined <= 3) %>% 
  group_by(reference_database_accession) %>% 
  filter(length == max(length)) %>% 
  ungroup()

Generally, you can use the information retrieved by fetch_pdb() in any way that fits your specific research question. However, you should note that by using very strict filtering criteria you might remove some of the proteins that actually have structural information. This could be because their structures do not have a high enough resolution or they are not solved using the method of your choice. Filtering by length should also be done with care since for some proteins there might not be any full length structure but two structures of different domains that would make up the whole structural information if used in combination.

# Input PDB IDs
pdb_ids <- unique(filtered_structures$pdb_ids) # "1ZMR", "2HWG"

# Fetch atom level structural information
ptsi_pgk_structure_information <- fetch_pdb_structure(pdb_ids = pdb_ids,
                                                      return_data_frame = TRUE)

# Fetch atom level structural prediction information from AlphaFold
ptsi_pgk_prediction_information <- fetch_alphafold_prediction(uniprot_ids = uniprot_ids,
                                                              return_data_frame = TRUE)

# Example for fetching all predictions for Methanocaldococcus jannaschii
# mj_predictions <- fetch_alphafold_prediction(organism_name = "Methanocaldococcus jannaschii")

Generation of structural contact maps

A very useful structural analysis function provided by protti is create_structure_contact_map(), which uses atom level structural information to calculate distances between atoms or amino acid residues. The generated contact map can be used in order to visually assess which regions are in close proximity to each other. Furthermore, you can use this function to identify all atoms that are in a certain proximity to a ligand or other small molecules. Additionally, protein-protein interfaces can be easily identified in structures that contain protein complexes.

As mentioned in the introduction it is important to find start and end positions of e.g. peptides in a specific structure file based on the UniProt positions. How UniProt and structure positions relate to each other can be extracted from data fetched through fetch_pdb(). We provide the convenient find_peptide_in_structure() function to actually perform the conversion of UniProt positions to structure specific positions. It either internally uses fetch_pdb() to obtain the necessary information or it takes the previously fetched output from fetch_pdb() as an input to its pdb_data argument. In the latter case it is then also possible to first filter the data to only perform the mapping on relevant structures, which will also speed up the process.

Note: If AlphaFold predictions are used for the contact map this step can be skipped and UniProt positions can be used directly since the AlphaFold prediction numbering is based on the UniProt positions.

In the following we are going to identify structure specific peptide start and end positions for our annotated example data (ptsi_pgk_annotated). By providing the previously created filtered_structures data frame as the input to the pdb_data argument, we will only determine the positions of our peptides with regards to the structures contained in the filtered data.

ptsi_pgk_peptide_structure_positions <- find_peptide_in_structure(
   peptide_data = ptsi_pgk_annotated,
   peptide = pep_stripped_sequence,
   start = start,
   end = end,
   uniprot_id = pg_protein_accessions,
   pdb_data = filtered_structures,
   retain_columns = c(eg_precursor_id, diff, adj_pval))

Note: The peptide argument in this function does not necessarily require the peptide sequence but can also take any other unique peptide identifier as input. This is possible because peptides are not matched to the structure sequence by their sequence but are identified based on their start and end positions in the protein.

Before we calculate the structure contact map we are going to filter our input data to only contain significant peptides. Those are the peptides that meet our cutoff criteria for the fold change (log2(2)) and adjusted p-values (0.01). The contact map will give us an insight into regions that are in close proximity to our significant peptides. In this way we will be able to assess whether significant peptides are in close proximity to each other.

As mentioned previously, there are multiple conventions for positions in structure files. For the create_structure_contact_map() function it is important to provide the author defined positions and chain names which are marked with the auth prefix. The distance_cutoff argument is set to 10 Angstrom and will only keep distances below that value, which will reduce the size of our output. Since contact maps can easily reach a size of several hundred megabytes (per map!) it is advisable to only keep the information strictly required. Another argument that is by default set to TRUE is return_min_residue_distance. This will return a contact map that contains minimal residue distances instead of atom distances, based on which the map is calculated initially.

By default contact maps are always created for the provided selection with regards to the whole structure. If you provide an additional data frame to the data2 argument it is possible to compare your selection provided to the data argument to this other subsection. This will reduce the size of contact maps if only a subset of distances should be calculated.

Note: You can also create contact maps for AlphaFold predictions. In that case you can provide the UniProt ID to the id argument.

Note: If you want to create a contact map of an unpublished or simply downloaded structure that you have the .pdb or .cif file of, you can provide the path to that structure to the structure_file argument. In this case the id argument should contain a column that only has the content “my_structure”.

Note: The input to the create_structure_contact_map() function can be a data frame containing specific regions (e.g. significant peptides) of a structure that you want to create the contact map for (as explained above and demonstrated below). It is however also possible to create contact maps for certain chains of the structure. In that case you only need to provide the structure ID and chain ID. Furthermore, it is possible to generate a complete contact map of each residue in the structure. In that case you will only need to provide the structure ID. Each of these options will increase the size of your map and also the time it takes to create it. Therefore, you should use these options with caution.

# Filter data for significant peptides.
significant_peptides <- ptsi_pgk_peptide_structure_positions %>% 
  filter(abs(diff) > 2, adj_pval <= 0.01)

# Create a structure contact maps 
contact_map <- create_structure_contact_map(
  data = significant_peptides,
  id = pdb_ids,
  chain = auth_asym_id,
  auth_seq_id = auth_seq_id,
  distance_cutoff = 10,
  pdb_model_number_selection = c(0, 1),
  return_min_residue_distance = TRUE
)

# This is a helper function for the plot.
# It allows the display of integers on the axis.
integer_breaks <- function(n = 5, ...) {
  fxn <- function(x) {
    breaks <- floor(pretty(x, n, ...))
    names(breaks) <- attr(breaks, "labels")
    breaks
  }
  return(fxn)
}

# Plot structure contact maps
# 1ZMR
contact_map[["1ZMR"]] %>% # Extract data frame from list
  mutate(chain_combinations = paste0("chain_", label_asym_id_var1, "_vs_chain_", label_asym_id_var2)) %>% 
  ggplot(aes(x = label_seq_id_var1, y = label_seq_id_var2, fill = min_distance_residue)) +
  geom_tile() +
  scale_y_continuous(breaks = integer_breaks()) +
  scale_x_continuous(breaks = integer_breaks()) +
  facet_wrap(~chain_combinations, scale = "free") +
  labs(title = "Structure contact map 1ZMR") +
  theme_bw()
  
# 2HWG
contact_map[["2HWG"]] %>% # Extract data frame from list
  mutate(chain_combinations = paste0("chain_", label_asym_id_var1, "_vs_chain_", label_asym_id_var2)) %>% 
  ggplot(aes(x = label_seq_id_var1, y = label_seq_id_var2, fill = min_distance_residue)) +
  geom_tile() +
  scale_y_continuous(breaks = integer_breaks()) +
  scale_x_continuous(breaks = integer_breaks()) +
  labs(title = "Structure contact map 2HWG") +
  facet_wrap(~chain_combinations, scale = "free") +
  theme_bw()

The contact map for each structure provided as input will be saved as a separate element of a list, which is returned as the output of the create_structure_contact_map() function. These elements can be easily accessed using two square brackets ([[]]) right behind the name of the list. The brackets should either contain the name of the element or the index. You can also access each element in the list in a loop using the map() function from the purrr package.

Contact maps can be plotted using the ggplot2 package if they are not too large. Contact maps are based on positions of atoms or residues in the structure. While atom IDs are unique, residues are only unique within a certain chain. That means residue position 1 is likely present in each of the chains. To correctly display contact map plots for residues they should be faceted by the unique combinations of chains. If e.g. the x-axis displays residues from chain A and the y-axis residues of chain B this should be displayed in one specific facet.

Note: For the plot we recommend using the data base defined chain and residue identifiers in contrast to the author defined ones that were used as input to the create_structure_contact_map() function. The reason is that in the data base definitions each molecular entity has their own chain. If a certain protein is for example associated with a metal ion, these might be stored in the same chain in the author defined identifiers but separated into different chains for the data base identifiers. The problem is that in the author defined identifiers the residue numbers might not be consecutive and there could be a large gap in the numbering between the protein and the metal residue.

The two contact maps that were generated based on significant peptides reveal some interesting information about the positioning of the peptides.

The “1ZMR” structure contains only one protein chain (A). In the contact map of significant peptides that are all located in chain A in relation to all other residues in chain A we can see that even though some of the peptides are quite far apart in the sequence, some are within 10 Angstrom distance in 3D space. Chain B and chain D are a calcium atom, respectively. In case of chain B the atom is actually in close proximity to one of the significant peptides. Chain E stores information about water molecules in the structure. The contact map is not very informative in this case, since the y-axis residue does not indicate any specific relationship of the water molecules with each other. Here it would be interesting to look at the distances of the water molecules to the peptide. If there are no water molecules in close proximity, this might indicate that the peptide is not surface exposed.

The “2HWG” structure contains two protein chains (A and B). Within each protein the significant peptides seem to be in close proximitiy to each other. In this case we can also observe that the peptide located around residue 555 is in close proximity to the other protein chain, which indicates that it might be located within a protein-protein interaction interface. Chain C and E contain a magnesium ion that seems to be in the proximity of two peptides. Chain D and F contain an oxalate ion which is known to inhibit ptsI. The oxalate ion is in very close proximity to to the same peptide that is close to the magnesium ion. This is likely the case because the oxalate ion and magnesium ion are located in the same site. Chain H and G contain water molecules.

Mapping of amino acids, peptides or regions on 3D protein structures and predictions

Contact maps provide quantitative information about the distances of residues to each other. We can use them to perform further high throughput analysis. In contrast to the rather complex contact map plots it is easier to look at the actual protein structures and the respective peptide positions to get a better overview of the proximity of peptides to each other and in relation to interesting sites. While this analysis is more subjective, it can make certain relationships in the data more obvious.

For the mapping of peptides, protein regions or amino acids on PDB structures or AlphaFold predictions, we provide the map_peptide_on_structure() function. Mapping is accomplished based on the replacement of B-factor information in the structure file with peptide specific values. When the structure is coloured by B-factor in PyMOL or ChimeraX, peptides will get highlighted. This function can on the one hand just display categorical information such as if the peptide is significant, present or undetected. On the other hand, continuous information can be displayed such as a score associated with each individual amino acid.

The mapping is performed based on positions and not peptide sequence matching. If you want to perform the mapping onto a PDB structure, you should run the find_peptide_in_structure() beforehand to obtain structure specific positions. If you want to map onto an AlphaFold prediction it is not necessary to carry out this step. However, the output of the find_peptide_in_structure() function will contain all necessary information to perform also AlphaFold predictions in case there is no available structure file.

For this function it is important to provide the author defined positions and chain names which are marked with the auth prefix. The function performs structural mapping based on this naming convention.

Note: If you perform mapping onto an AlphaFold prediction the content of the pdb_id provided column should be NA.

Note: If you want to map peptides to an unpublished or downloaded structure that you have the .pdb or .cif file of, you can provide the path to that structure to the structure_file argument. In this case you do not need to provide the pdb_id argument. It is very important that your input data frame peptide_data only contains information for the protein in the provided structure file and not for other unrelated proteins, otherwise mapping is not performed correctly.

You provide information on colouring through the map_value argument. This should contain a numeric column. In the example below we insert the number 100 into this column if a peptide is significantly changing based on our previous cutoff criteria and the value 0 if it is not but was detected in the experiment. In this case you could simply insert any two numbers. However, the first number that colours your positive hits should be higher. The function will internally scale these values between 50 and 100. In this case 0 would become 50 and 100 stays 100. All regions of the structure that were not covered by any of you peptides receive the final value 0. This allows us to distinguish between covered (50-100) and non-covered (0) regions as well as covered non-significant (50) and significant regions (100). This is especially helpful if peptides are coloured according to a continuous score (see later in this vignette).

Lastly, we can specify an export location to indicate a place where the structure file should be saved. In this example the location is specified as a temporary file directory using the tempdir() function. You should change that location to any more accessible place. If you leave the argument empty, structures will be automatically saved in your working directory.

ptsi_pgk_peptide_structure_positions %>% 
  mutate(map_value = ifelse(eg_precursor_id %in% significant_peptides$eg_precursor_id, 
                            100, 
                            0)) %>% 
  map_peptides_on_structure(
   uniprot_id = pg_protein_accessions,
   pdb_id = pdb_ids,
   chain = auth_asym_id,
   auth_seq_id = auth_seq_id,
   map_value = map_value,
   file_format = ".pdb", 
   export_location = tempdir() # change to a location of your choice
 )

The function does not return anything in the R environment. Structures are directly saved in the desired format. The files contain the name of the structure and separated by “_” the name of all UniProt IDs present in the structure file.

To visualise coloured regions in your proteins you should open the file in either PyMOL or ChimeraX. There you have the option to colour structures by B-factor. Alternatively, we present an option within R using the r3dmol package to interactively display structures with the help of htmlwidgets.

# Install the r3dmol package if it is not installed
# install.packages("r3dmol")

# Load the r3dmol package
library(r3dmol) 

# Create structure
r3dmol() %>%
  m_add_model(data = paste0(tempdir(), "/1ZMR_P0A799.pdb"), format = "pdb") %>% 
  m_set_style(style = m_style_cartoon(
    colorfunc = "
        function(atom) {
          if (atom.b == 50) {return '#90EE90'};
          if (atom.b == 100) {return '#FF7276'};
          if (atom.b == 0) {return 'white'};
          return 'white';
        }"
  )) %>%
  m_zoom_to()

spectrum b, gray70_palegreen_deepsalmon, minimum = 0, maximum = 100

set_color cpink, [235, 52, 189]
set_color cblue, [52, 171, 235]

color bfactor palette darkgrey:darkseagreen:salmon

color byattribute bfactor protein palette darkgrey:darkseagreen:salmon

# Calculate the amino acid score
amino_acid_score <- calculate_aa_scores(
  data = ptsi_pgk_peptide_structure_positions,
  protein = pg_protein_accessions,
  diff = diff,
  adj_pval = adj_pval,
  start_position = start,
  end_position = end,
  retain_columns = c(pdb_ids, auth_asym_id)
)

# Find amino acid positions in the structure
ptsi_pgk_amino_acid_structure_positions <- find_peptide_in_structure(
   peptide_data = amino_acid_score,
   peptide = residue,
   start = residue,
   end = residue,
   uniprot_id = pg_protein_accessions,
   pdb_data = filtered_structures,
   retain_columns = c(amino_acid_score))

# Map the score on structure
map_peptides_on_structure(
  peptide_data = ptsi_pgk_amino_acid_structure_positions,
   uniprot_id = pg_protein_accessions,
   pdb_id = pdb_ids,
   chain = auth_asym_id,
   auth_seq_id = auth_seq_id,
   map_value = amino_acid_score,
   file_format = ".pdb",
   export_location = tempdir()
 )

# create a color gradient with 101 colors
color_gradient <- paste0('"', 
                         paste(colorRampPalette(c("white", "#90EE90", "#FF7276"))(101), 
                               collapse = '", "'), 
                         '"')

# create structure
r3dmol() %>%
  m_add_model(data = paste0(tempdir(), "/2HWG_P08839.pdb"), format = "pdb") %>%
  m_set_style(style = m_style_cartoon(
    colorfunc = paste0("
        function(atom) {
          const color = [", color_gradient,"]
          return color[Math.round(atom.b)]
        }")
  )) %>%
  m_zoom_to()