Structural data

The world wide Protein Data Bank (wwPDB) manages and archives information about the 3D structure of proteins. Its information can be accessed through several websites such as the Research Collaboratory for Structural Bioinformatics (RCSB).

When working with structure files it is important to keep a few things in mind such that mistakes can be avoided. One aspect to consider is that protein structure files come in two mainly used formats, namely the PDB file format (.pdb) and the mmCIF file format (.cif). The PDB file format is the legacy file format of the Protein Data Bank, which has been replaced with the newer mmCIF file format that stores data on biological macromolecules. All structures available in the Protein Data Bank are available as a mmCIF file, while not all are available as a PDB file. Mainly large protein complexes only contain mmCIF files.

The way both files store data and relate to each other can be quite complex. There are several resources explaining their relation and their content. The most important thing to know about these files for you as a user of protti is the fact that there are different naming conventions for elements of the structure.

PDB files solely use the naming provided by the author of the structure in order to match the identification used in the publication that describes the structure. Columns containing author information start with the prefix “auth_”. This includes the naming of chains (asym_id), atom IDs (atom_id) and residue IDs (seq_id).

mmCIF files use both the naming provided by the author and a naming following the standardised conventions for mmCIF files. Columns containing standardised information start with the prefix “label_”.

Additional information: The differences between author provided and standardised labeling may vary a lot depending on the structure and can be quite confusing. An example taken from the documentation of mmCIF files describes the relationship between author provided and standardised residue IDs as follows: The label_seq_id (standardised residue ID) column is required to be a sequential list of positive integers, while the auth_seq_id (author residue ID) is not necessarily a number and the values do not have to be positive. The author may assign values to auth_seq_id in any desired way. For instance, the values may be used to relate this structure to a numbering scheme in a homologous structure, including sequence gaps or insertion codes. Alternatively, a scheme may be used for a truncated polymer that maintains the numbering scheme of the full length polymer. In all cases, the scheme used here must match the scheme used in the publication that describes the structure.

For the use of the functions in protti it is important to know that these two different conventions exist. All functions usually take the author provided information as inputs. This has the reason that even though the PDB file format is a legacy format, many people still mainly use this format. With the find_peptide_in_structure() function we provide a convenient way for you to find both the standardised and the author provided residue numberings for a peptide, protein region or amino acid, returning their UniProt start and end position. In addition, we return the full stretch of residue annotations for the author provided information in the column auth_seq_id since they are not necessarily consecutive as previously mentioned.

Loading data

For this vignette we use a subset of two proteins from separate LiP-MS experiments that contain peptides significantly changing in abundance due to a structural change. We combined these two proteins into one data set to showcase that multiple proteins and structures can be analysed at once. The data was obtained from the paper “Dynamic 3D proteomes reveal protein functional alterations at high resolution in situ” (Cappelletti 2021). Both experiments were conducted on purified E. coli proteins spiked into E. coli lysates. The first protein is phosphoglycerate kinase 1 (pgk) which was treated with 25 mM 3-phosphoglyceric acid (3PG). The second protein is phosphoenolpyruvate-protein phosphotransferase (ptsI) which was treated with 25 mM fructose 1,6-bisphosphatase (FBP).

Based on structural changes observed in these proteins the authors conclude that FBP binds to the active site of ptsI and acts as a competitive inhibitor, which is confirmed by an activity assay. The structural change in pgk caused by 3PG is likely representative of substrate occupancy and correlates with metabolic flux changes.

The example data is included in protti and you can easily use it by calling the data() function. This will load the data into your R session.

utils::data("ptsi_pgk")

If you want to read your own data into R, we suggest using the read_protti() function. This function is a wrapper around the fast fread() function from the data.table package and the clean_names() function from the janitor package. This will allow you to not only load your data into R very fast, but also to clean up the column names into lower snake case. This will make it easier to remember them and to use them in your data analysis. Usually you would perform the following analysis after the analyses described in either the vignette for single dose treatments or the vignette for dose response data.

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

Preparing data

First we need to annotate the data with information required for structural data analysis. To find peptides in a protein structure we need to know their positions in the protein first. For this we can use the protti function find_peptide(). This function searches for the peptide within the protein sequence, which we need to fetch from UniProt. For this we can use the protti function fetch_uniprot(). In addition to the protein sequence we will also need information about PDB identifiers associated with the protein, which can also be accessed via UniProt.

# 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", "xref_pdb")
)

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

Fetching atomic structural information

In addition to the general structural information that can be retrieved with fetch_pdb() you can also access atom level information for each structure using the fetch_pdb_structure() function. This information is mainly useful for applications that rely on the coordinates of each atom in the structure. Later in this vignette we will demonstrate how this can be used in the create_structure_contact_map() function. Furthermore, fetch_pdb_structure() retrieves B-factors for each atom.

Similar to fetch_pdb() PDB identifiers are provided as input to fetch_pdb_structure(). In the example below we provide the two PDB IDs that were left after the filtering from the previous step was applied to the data. If not too many structures (~<20) are fetched at once it is possible to combine all information in a data frame by setting the return_data_frame argument to TRUE. This is by default FALSE, which will return a list instead with each element representing one structure.

# 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
)

Fetching atomic structure prediction information from AlphaFold

If a protein does not have any available structure you can alternatively fetch protein structure predictions from AlphaFold using the fetch_alphafold_prediction() function. The function works similar to fetch_pdb_structure() but it takes UniProt identifiers instead of PDB identifiers as input. Furthermore, there is the option to fetch all predictions for an organism as a whole. The B-factor information in AlphaFold predictions represents the quality score of the atom position prediction.

It is of course also possible to retrieve structure prediction information for proteins that have a known protein structure. In the example below we will retrieve prediction information for our two proteins ptsi (P08839) and pgk (P0A799). The output column score_quality annotates the score from the prediction_score column based on the official AlphaFold score quality annotations and cutoffs.

# 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")

AlphaFold domain predictions using the predicted aligned error (PAE)

AlphaFold predictions come with two kinds of error metrics, the per-residue confidence score (pLDDT) that is saved in the b-factor information of structure files and the predicted aligned error (PAE) that can be downloaded separately from the EBI database. According to the description on the EBI website, the PAE matrix at position (x, y) indicates AlphaFold’s expected position error at residue x, when the predicted and true structures are aligned on residue y. This is useful for assessing inter-domain accuracy. One interesting use case for the PAE is the prediction of domains in AlphaFold predictions. Tristan Croll created a python script that uses a graph-based community clustering algorithm on PAEs in order to infer pseudo-rigid protein domains. These domain predictions can be for example used in order to check the validity of distance calculations between residues of a predicted structure. If residues are on two different domains distances should be used with caution.

With protti we provide the function fetch_alphafold_aligned_error() that allows you to retrieve the predicted aligned error for AlphaFold predictions directly into your R session. In addition, we adapted the python function from Tristan and created the protti function predict_alphafold_domain() that predicts protein domains for AlphaFold predictions. In the example below we will retrieve PAEs for our two proteins ptsi (P08839) and pgk (P0A799) and subsequently predict protein domains. The error_cutoff argument of the fetch_alphafold_aligned_error() function defines the maximum error that should be saved. All higher values will be discarded in order to reduce the size of the returned list. For the predict_alphafold_domain() function we try out two different graph_resolution values which will determine the size of the clusters and thus how strict domains are created.

# Fetch aligned errors
aligned_error <- fetch_alphafold_aligned_error(
  uniprot_ids = uniprot_ids,
  error_cutoff = 4
)

# Predict protein domains with graph_resolution of 1
af_domains_res_1 <- predict_alphafold_domain(
  pae_list = aligned_error,
  return_data_frame = TRUE,
  graph_resolution = 1
) # Default

# Predict protein domains with graph_resolution of 0.5
af_domains_res_05 <- predict_alphafold_domain(
  pae_list = aligned_error,
  return_data_frame = TRUE,
  graph_resolution = 0.5
)

3D structure mapping in R using r3dmol

It is possible to use the r3dmol package to visualize protein structures directly in R. The package is still under development and provides support for 3Dmol.js. You can find extensive documentation and examples on their GitHub page. Below we provide a specific example of how to use the package in order to visualise structures created by map_peptides_on_structure() in R colouring by B-factors. Later we provide an example of how to use a colour gradient. In fact the same function that is used later can also be used in this case.

If you run the code snippets on your own make sure you replace paste0(tempdir(), "/1ZMR_P0A799.pdb") with the location and file you want to visualise.

r3dmol facilitates direct interaction with the structure if this vignette is exported as an html document. You can for example zoom, rotate or move the structure.

# 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()

PyMOL

The following command colours structures by B-factor. You can specify a gradient with predefined or custom colours separating them by underscores. The gradient can also contain more than 3 different colour steps. The minimum and maximum values of the gradient should also be specified. These are based on the previously mentioned scaling. In the below gradient all non-covered regions will be gray. Covered by non-significant regions will be “palegreen” and signficiant regions will be “deepsalmon”.

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

Custom colours can be defined with the following command and subsequently used in the command above.

Note: Custom names should not contain underscores since this confuses PyMOL when the gradient function is used.

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

ChimeraX

In Chimera you have similar options like in PyMOL for colouring structures by B-factors. You can use the command below to create a colour gradient over the structure.

color bfactor palette darkgrey:darkseagreen:salmon

If you only want to colour proteins by B-factors while leaving other entities such as nucleic acids or small molecules in their default colour you can use the following command instead.

color byattribute bfactor protein palette darkgrey:darkseagreen:salmon

Result

The figures below are the results of performing the colouring steps on the two structures of ptsI and pgk with significantly changing peptides highlighted in red. Areas in gray were not covered by any detected peptides.

For ptsI which was treated with fructose 1,6-bisphosphatase (FBP) the significantly changing peptides are all located in the same region of the protein. They are in close proximity to the active site that contains the magnesium atom coloured in blue and the oxalate ion coloured in purple. The contact maps showed a significant peptide in close proximity to the protein-protein interaction surface of the homodimer. This peptide is indicated by black arrows. It is not buried in the binding site but rather located at the edge of it. Its residues are, however, directly participating in the formation of the interaction surface through hydrogen bonds and likely a pi-pi stacking between a proline and histidine residue. The fact that there are no significant peptides buried in the interaction surface likely indicates that the complex does not dissociate into its subunits upon FBP treatment. It is more likely that FBP treatment induces a general structural change observable in the vicinity of the active site due to its likely binding to this site. This observation is also in line with the fact that FBP is an inhibitor of ptsI (Cappelletti 2021).

Mapping of significantly changing peptides onto the pgk structure reveals that most of the protein seems to be structurally affected by the treatment with 3-phosphoglyceric acid (3PG). Based on UniProt annotations and the similarity of the protein to its homologues in other organisms, the 3PG binding site is located on the left side indicated with a black arrow. None of the residues involved in the interaction with 3PG are found in a significantly changing peptide. Due to the compact size of the protein and the even distribution of significant peptides it is hard to draw any conclusions about potential structural effects of the treatment.

Mapping of amino acid scores

In order to narrow down and prioritize structurally affected protein regions, we have developed the amino acid scoring function calculate_aa_scores(). For each amino acid position in the protein the function calculates the average product of the -log10(adjusted p-value) and the absolute log2(fold change) per peptide, covering the amino acid position. In other words, the information of all the peptides in a specific region is integrated, which strengthens the confidence of a structural change being true since false positive peptides are averaged out. One downside of this approach is that it relies on very good sequence coverage since otherwise false positive peptides with high scores that solely cover a specific region are not distinguishable from real hits. In our example we are lucky that both proteins have a very high coverage.

To map the amino acid scores, the amino acid positions within the protein structure need to be correctly assigned first using find_peptide_in_structure(), with the residue column as input for peptide, the start and the end position. This step is required after the amino acid score calculation, since calculate_aa_score() needs a numeric input for start and end positions, and therefore, auth_seq_id (a character column) cannot be used as an input.

The map_peptides_on_structure() function is capable of mapping continuous values to a structure. In this case values are again scaled between 50 and 100. It is likely that values fall in between 50 and 100 which is displayed as a colour gradient step in between the colour in the middle of the gradient (50) and the colour at the end of the gradient (100).

Important: If you run this Vignette locally, make sure to save the structures generated next in a different folder than the previous ones. Otherwise the last structures will be overwritten.

Note: If multiple structures are mapped at the same time, by default each structure scales its map values between 50 and 100 based on its minimum and maximum score. It is also possible to use one scale for all structures that are mapped if the scale_per_structure argument of the map_peptides_on_structure() function is set to FALSE. In this case only one of the structures (the one containing the highest score) will have the 100 value. The highest score of the other structures will be lower than 100. This means that the scale of colours between both structures is comparable.

# 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()
)

While the previous method of displaying only significantly changing peptides that have a minimum fold change is useful for the identification of regions that change the most, it is not resistant to false positive peptides. The score on the other hand will only be high if an amino acid is consistently part of significant and differentially abundant peptides. The problem of the score is that it is hard and visually even impossible to interpret what a certain score means. This is due to the fact that the same score can be caused by two completely different scenarios. On the one hand an amino acid might be associated with peptides that have a low fold change but are highly significant. On the other hand it can be associated with peptides that have a high fold change but are of low significance. Both scenarios could yield the same amino acid score. Therefore, the score should not be misinterpreted as describing the extent of the change but should rather be seen as a probability that a certain amino acid is really structurally affected.

You can visualise the score as a colour gradient using the r3dmol package in R. For that you need to first create a gradient using the colorRampPalette() function, which creates a function capable of creating a colour gradient over the specified colours. Instead of providing only three colours to this function you can provide additional colours to make a gradient with multiple substeps. It is important to always use 101 colours to get the full range for the structure. The colour gradient needs to be concatenated using paste0() and provided to the colorfunc argument within the shown javascript function context.

# 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()

Applying the score to both structures refines the region that is most robustly affected by the treatment. For ptsI it seems that a region of amino acids directly in the active center are consistently affected. In addition we see a small region below, which was previously not seen due to its low fold changes that did not meet the cutoff. The score tells us that this region is indeed consistently affected and that there is likely a structural change occurring even though the fold change of peptides is rather low. The region in the binding interface previously indicated with a black arrow is now not highlighted anymore, indicating that this was likely a false positive hit and that other peptides in the same region do not exhibit a strong and significant response.

For ptg we see that the most robustly changing region is very close to the binding site of 3PG, which the protein was treated with. In this case, the score helped to refine the binding site of the compound since previously almost the whole protein was covered with significantly changing peptides.