inst/doc/data_analysis_dose_response_workflow.R

In protti: Bottom-Up Proteomics and LiP-MS Quality Control and Data Analysis Tools

## ----include = FALSE----------------------------------------------------------
build_vignette_on_cran <- identical(Sys.getenv("BUILD_VIGNETTE"), "true")
test_protti <- identical(Sys.getenv("TEST_PROTTI"), "true") & build_vignette_on_cran
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----CRAN_comment, message=FALSE, warning=FALSE, echo=FALSE-------------------
if (build_vignette_on_cran == FALSE) {
  print("!!! IMPORTANT !!!")
  print("This Vignette has not been built completely on CRAN due to size limitations.")
  print("Please check the correct version here: ")
  print("https://jpquast.github.io/protti/articles/https://jpquast.github.io/protti/articles/data_analysis_dose_response_workflow.html")
}

## ----load_packages, eval = test_protti, warning = FALSE, message = FALSE------
# # Load packages
# library(protti)
# library(dplyr)
# library(magrittr)
# library(ggplot2)

## ----use_read_protti, eval=FALSE----------------------------------------------
# # Load data
# rapamycin_dose_response <- read_protti("your_data.csv")

## ----load_data, eval = test_protti, warning = FALSE---------------------------
# utils::data("rapamycin_dose_response")

## ----filter_transform_normalise, eval = test_protti---------------------------
# # 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)

## ----intensity_distribution, eval = test_protti, fig.align = "center", fig.width = 7, fig.height = 5----
# qc_intensity_distribution(
#   data = data_normalised,
#   grouping = eg_precursor_id,
#   intensity = normalised_intensity_log2,
#   plot_style = "histogram"
# )

## ----intensity_filtering, eval = test_protti----------------------------------
# data_normalised <- data_normalised %>%
#   filter(normalised_intensity_log2 > 5)

## ----sample_correlation, eval = test_protti, fig.align = "center", fig.width = 7, fig.height = 5----
# qc_sample_correlation(
#   data = data_normalised,
#   sample = r_file_name,
#   grouping = eg_precursor_id,
#   intensity_log2 = normalised_intensity_log2,
#   condition = r_condition
# )

## ----sample_pca, eval = test_protti, fig.align = "center", fig.width = 7, fig.height = 5, warning = FALSE, message = FALSE----
# qc_pca(
#   data = data_normalised,
#   sample = r_file_name,
#   grouping = eg_precursor_id,
#   intensity = normalised_intensity_log2,
#   condition = r_condition
# )

## ----model_fit, eval = test_protti--------------------------------------------
# fit <- data_normalised %>%
#   fit_drc_4p(
#     sample = r_file_name,
#     grouping = eg_precursor_id,
#     response = normalised_intensity_log2,
#     dose = r_condition,
#     filter = "post",
#     n_replicate_completeness = 2, # specifiy this based on your data and cutoff criteria
#     n_condition_completeness = 4, # specifiy this based on your data and cutoff criteria
#     retain_columns = c(pg_protein_accessions)
#   )
# 
# # make sure to retain columns that you need later but that are not part of the function

## ----parallel_model_fit, eval = FALSE-----------------------------------------
# # 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,
#     n_replicate_completeness = 2, # specifiy this based on your data and cutoff criteria
#     n_condition_completeness = 4, # specifiy this based on your data and cutoff criteria
#   )
# 
# # remove workers again after you are done
# future::plan(future::sequential)

## ----result_analysis, eval = test_protti, echo=FALSE, results='asis'----------
# fit %>%
#   filter(rank <= 20) %>%
#   select(
#     rank,
#     score,
#     eg_precursor_id,
#     pg_protein_accessions,
#     anova_adj_pval,
#     correlation,
#     ec_50
#   ) %>%
#   mutate(
#     anova_adj_pval = format(anova_adj_pval, digits = 3),
#     correlation = format(correlation, digits = 3),
#     ec_50 = format(ec_50, digits = 2),
#     score = format(score, digits = 3)
#   ) %>%
#   knitr::kable(caption = "All hits")

## ----model_plot, eval = test_protti, fig.align = "center", fig.width = 7, fig.height = 5, message = FALSE, warning = FALSE----
# # Model plotting
# drc_4p_plot(fit,
#   grouping = eg_precursor_id,
#   dose = r_condition,
#   response = normalised_intensity_log2,
#   targets = "_VFDVELLKLE_.2",
#   unit = "pM",
#   x_axis_limits = c(10, 1e+08),
#   export = FALSE
# )

## ----uniprot_annotation, eval = test_protti-----------------------------------
# # 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")

## ----post_analysis, eval = FALSE----------------------------------------------
# ### 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
#   )