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/https://jpquast.github.io/protti/articles/data_analysis_dose_response_workflow.html"

This vignette will take you through all important steps for the analysis of dose-response experiments with protti. For experiments with treatment titrations (i.e. samples treated with different amounts of the same compound, protein, metabolite, element, ion, etc.) you can fit sigmoidal log-logistic regression curves to the data and determine significantly changing precursors, peptides or proteins based on Pearson correlation and ANOVA q-value.

Before analysing your data please ensure that it is of sufficient quality by using protti’s quality control functions (see Quality control (QC) workflow) If you would like to compare a control condition with one treatment condition (or multiple unrelated treatments) you can check out the data analysis vignette for single dose treatment data. For help with data input please click here.

A typical dose-response experiment contains multiple samples that were treated with different amounts of e.g. a drug. Replicates of samples treated with the same dose make up a condition. Commonly, the first concentration is 0 (i.e. the control, in which treatment is with the solvent of the drug). Dose-response treatments require a minimal number of treatments to fit curves with sufficient quality. For analysis with protti at least 5 different conditions should be present. Another consideration for your experiment is the range of concentrations. They should not be too close together since effects are usually best identified over a larger concentration range. But you should make sure not to space them out too far. It is generally advised to space them evenly on a logarithmic scale of the base 10 or \(e\) (Euler’s number). You can also include steps in between e.g. 100, 500, 1000 or 100, 200, 1000 for a log10 scale. It is advisable to use a rather broad concentration range for an experiment in which you do not know what to expect.

protti fits four-parameter log-logistic dose-response models to your data. It utilizes the drm() and LL.4() functions from the drc package for this. You can also select a non logarithmic model using the L.4() function in case your data does not follow a log-logistic but a logistic regression.

For limited proteolysis-coupled to mass spectrometry (LiP-MS) data, dose-response curves have been used previously for the identification of drug binding sites in complex proteomes (Piazza 2020). Since LiP-MS data is analysed on the peptide or precursor* level, using additional information about peptide behaviour from multiple conditions reduces false discovery rate.

_*A peptide precursor is the actual molecular unit that was detected on the mass spectrometer. This is a peptide with one specific charge state and its modification(s)._

Getting started

Before starting the 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 after you installed them with the library() function.

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

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

utils::data("rapamycin_dose_response")

# Filter, log2 transform and normalise data
data_normalised <- rapamycin_dose_response %>%
  filter(eg_is_decoy == FALSE) %>%
  mutate(intensity_log2 = log2(fg_quantity)) %>%
  normalise(
    sample = r_file_name,
    intensity_log2 = intensity_log2,
    method = "median"
  ) %>%
  filter(pep_is_proteotypic == TRUE)

qc_intensity_distribution(
  data = data_normalised,
  grouping = eg_precursor_id,
  intensity = normalised_intensity_log2,
  plot_style = "histogram"
)

data_normalised <- data_normalised %>%
  filter(normalised_intensity_log2 > 5)

qc_sample_correlation(
  data = data_normalised,
  sample = r_file_name,
  grouping = eg_precursor_id,
  intensity_log2 = normalised_intensity_log2,
  condition = r_condition
)

qc_pca(
  data = data_normalised,
  sample = r_file_name,
  grouping = eg_precursor_id,
  intensity = normalised_intensity_log2,
  condition = r_condition
)

Fitting dose-response curves

The main part of a dose-response data analysis is the fit_drc_4p() function. It fits four-parameter log-logistic dose-response curves for every precursor (or peptide/protein) with the following equation:

\[y = d + \frac{a - d}{1 + (\frac{x}{c})^{-b}}\] \(x\) is the dependent variable
\(y\) is the independen variable
\(a\) is the minimum value that can be obtained
\(d\) is the maximum value that can be obtained
\(c\) is the point of inflection
\(b\) is the Hill’s coefficient (e.g. the negative slope at the inflection point)

The output of fit_drc_4p() provides extensive information on the goodness of fit. Furthermore, the function filters and ranks fits based on several criteria which include completeness of data and a significance cutoff based on an adjusted p-value obtained from ANOVA. For details about the filtering steps you can read the function documentation by calling ?fit_drc_4p. If you only care about your potential hits (and exclude precursors/peptides/proteins that do not pass the filtering e.g. due to too few observations), you can choose filter = "pre", which filters the data before model fitting. This speeds up the process because less models need to be fit. If you want to perform enrichment analysis later on you should keep the default filter = "post". This will fit all possible models and then only annotate models based on whether they passed or failed the filtering step (passed_filter).

fit <- data_normalised %>%
  fit_drc_4p(
    sample = r_file_name,
    grouping = eg_precursor_id,
    response = normalised_intensity_log2,
    dose = r_condition,
    filter = "post",
    retain_columns = c(pg_protein_accessions)
  )

# make sure to retain columns that you need later but that are not part of the function

If you have to fit many models this function may take about an hour or longer. If you have the computational resources you can speed this up by using parallel_fit_drc_4p(), a parallelised version of this function.

Note: Keep in mind that by spreading the model fitting over multiple cores, more memory than usual is needed. Furthermore, using the maximal number of available cores may not be the most efficient, because it takes time to copy the data to and from each core - this process will not be performed in parallel. Therefore, it is advised that the run time on each core is at least as long as the time to set up (data transfer, environment preparation etc.) all cores.

# setup of cores. Make sure you have the future package installed
future::plan(future::multisession, workers = 3)

# fit models in parallel
parallel_fit <- data_normalised %>%
  parallel_fit_drc_4p(
    sample = r_file_name,
    grouping = eg_precursor_id,
    response = normalised_intensity_log2,
    dose = r_condition,
    retain_columns = c(pg_protein_accessions),
    n_cores = 3
  )

# remove workers again after you are done
future::plan(future::sequential)

If we examine all precursors based on their rank (calculated from correlation and ANOVA q-value), we can see that as expected most of them are FKBP12 (P62942) peptides. One benefit of analysing LiP-MS data on the precursor level (rather than peptide level) is that multiple lines of evidence for a change in the specific peptide can be used (we have data for each different charge and modification state). Therefore, it is always good to check if there are other precursors of a good scoring peptide that do not show any regulation at all. This would mean that the reason for the observed response is not based on a biological effect.

# Model plotting
drc_4p_plot(fit,
  grouping = eg_precursor_id,
  dose = r_condition,
  response = normalised_intensity_log2,
  targets = "_VFDVELLKLE_.2",
  unit = "pM",
  export = FALSE
)

Further analysis

After models are fit, you can have a deeper look into the hits you obtained. protti provides several functions that help you to identify patterns in your data.

You can test your data set for gene ontology (GO) term enrichment (calculate_go_enrichment()). Alternatively, pathway annotations provided in the KEGG database can be used to test for pathway enrichments (calculate_kegg_enrichment()).

If you know which proteins bind or interact with your specific treatment, you can provide your own list of true positive hits and check if these are enriched in your significant hits by using protti’s calculate_treatment_enrichment() function. For our LiP-MS experiment using rapamycin we are probing direct interaction with proteins in contrast to functional effects. The only protein that rapamycin is known to bind to is FKBP12 and testing the significance of enrichment for a single or even a few proteins is not appropriate. However, testing for enrichment is especially useful if your treatment affects many proteins, since it can help you to reduce the complexity of your result.

The STRING database provides a good resource for the analysis of protein interaction networks. It is often very useful to check for interactions within your significant hits. For LiP-MS data this sometimes explains why proteins that do not directly interact with your treatment are still significantly affected. With analyse_functional_network(), protti provides a useful wrapper around some STRINGdb package functions.

# fetching of UniProt information
unis <- unique(fit$pg_protein_accessions)
uniprot <- fetch_uniprot(unis)

# annotation of fit data based on information from UniProt
fit_annotated <- fit %>%
  # columns containing proteins IDs are named differently
  left_join(uniprot, by = c("pg_protein_accessions" = "accession")) %>%
  # mark peptides that pass the filtering
  mutate(passed_filter = !is.na(rank)) %>%
  # create new column with prior knowledge about binding partners of treatment
  mutate(binds_treatment = pg_protein_accessions == "P62942")

### GO enrichment using "molecular function" annotation from UniProt

calculate_go_enrichment(fit_annotated,
  protein_id = pg_protein_accessions,
  is_significant = passed_filter,
  go_annotations_uniprot = go_f
) # column obtained from UniProt

### KEGG pathway enrichment

# First you need to load KEGG pathway annotations from the KEGG database
# for your specific organism of interest. In this case HeLa cells were
# used, therefore the organism of interest is homo sapiens (hsa)

kegg <- fetch_kegg(species = "hsa")

# Next we need to annotate our data with KEGG pathway IDs and perform enrichment analysis

fit %>%
  # columns containing proteins IDs are named differently
  left_join(kegg, by = c("pg_protein_accessions" = "uniprot_id")) %>%
  calculate_kegg_enrichment(
    protein_id = pg_protein_accessions,
    is_significant = passed_filter,
    pathway_id = pathway_id, # column name from kegg data frame
    pathway_name = pathway_name
  ) # column name from kegg data frame

### Treatment enrichment analysis

calculate_treatment_enrichment(fit_annotated,
  protein_id = pg_protein_accessions,
  is_significant = passed_filter,
  binds_treatment = binds_treatment,
  treatment_name = "Rapamycin"
)

### Network analysis

fit_annotated %>%
  filter(passed_filter == TRUE) %>% # only analyse hits that were significant
  analyse_functional_network(
    protein_id = pg_protein_accessions,
    string_id = xref_string, # column from UniProt containing STRING IDs
    organism_id = 9606,
    # tax ID can be found in function documentation or STRING database
    binds_treatment = binds_treatment,
    plot = TRUE
  )