inst/doc/data_analysis_single_dose_treatment_workflow.R

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

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

## ----setup, eval = test_protti, message = FALSE, warning = FALSE--------------
# library(protti)
# library(magrittr)
# library(dplyr)

## ----eval=FALSE---------------------------------------------------------------
# # To read in your own data you can use read_protti()
# your_data <- read_protti(filename = "mydata/data.csv")

## ----load_data, eval = test_protti--------------------------------------------
# data("rapamycin_10uM")

## ----median_normalisation_filtering, eval = test_protti, message = FALSE, warning = FALSE----
# data_normalised <- rapamycin_10uM %>%
#   filter(eg_is_decoy == FALSE) %>%
#   mutate(intensity_log2 = log2(fg_quantity)) %>%
#   normalise(
#     sample = r_file_name,
#     intensity_log2 = intensity_log2,
#     method = "median"
#   )
# 
# data_filtered <- data_normalised %>%
#   filter_cv(
#     grouping = eg_precursor_id,
#     condition = r_condition,
#     log2_intensity = intensity_log2,
#     cv_limit = 0.25,
#     min_conditions = 1
#   )

## ----data_preparation_prot_pep, eval = test_protti, message = FALSE, warning = FALSE----
# data_filtered_proteotypic <- data_filtered %>%
#   filter(pep_is_proteotypic == TRUE)

## ----fetching_database_information, eval = test_protti, message = FALSE, warning = FALSE----
# uniprot_ids <- unique(data_filtered_proteotypic$pg_protein_accessions)
# 
# uniprot <-
#   fetch_uniprot(
#     uniprot_ids = uniprot_ids,
#     columns = c(
#       "protein_name",
#       "gene_names",
#       "go_f",
#       "xref_string",
#       "cc_interaction",
#       "ft_act_site",
#       "ft_binding",
#       "xref_pdb",
#       "length",
#       "sequence"
#     )
#   ) %>%
#   rename(pg_protein_accessions = accession)
# 
# data_filtered_uniprot <- data_filtered_proteotypic %>%
#   left_join(
#     y = uniprot,
#     by = "pg_protein_accessions"
#   ) %>%
#   find_peptide(
#     protein_sequence = sequence,
#     peptide_sequence = pep_stripped_sequence
#   ) %>%
#   assign_peptide_type(
#     aa_before = aa_before,
#     last_aa = last_aa,
#     aa_after = aa_after,
#     protein_id = pg_protein_accessions
#   ) %>%
#   calculate_sequence_coverage(
#     protein_sequence = sequence,
#     peptides = pep_stripped_sequence
#   )

## ----coverage_plot, eval = test_protti, fig.align= "center", fig.width = 6, fig.height = 5----
# qc_sequence_coverage(
#   data = data_filtered_uniprot,
#   protein_identifier = pg_protein_accessions,
#   coverage = coverage
# )

## ----t_test, eval = test_protti, message = FALSE, warning = FALSE-------------
# diff_abundance_data <- data_filtered_uniprot %>%
#   assign_missingness(
#     sample = r_file_name,
#     condition = r_condition,
#     grouping = eg_precursor_id,
#     intensity = normalised_intensity_log2,
#     ref_condition = "control",
#     completeness_MAR = 0.7,
#     completeness_MNAR = 0.25,
#     retain_columns = c(
#       pg_protein_accessions,
#       go_f,
#       xref_string,
#       start,
#       end,
#       length,
#       coverage
#     )
#   ) %>%
#   calculate_diff_abundance(
#     sample = r_file_name,
#     condition = r_condition,
#     grouping = eg_precursor_id,
#     intensity_log2 = normalised_intensity_log2,
#     missingness = missingness,
#     comparison = comparison,
#     method = "moderated_t-test",
#     retain_columns = c(
#       pg_protein_accessions,
#       go_f,
#       xref_string,
#       start,
#       end,
#       length,
#       coverage
#     )
#   )

## ----pval_distribution, eval = test_protti, message = FALSE, warning = FALSE, fig.align= "center", fig.width = 6, fig.height = 5----
# pval_distribution_plot(
#   data = diff_abundance_data,
#   grouping = eg_precursor_id,
#   pval = pval
# )

## ----volcano_plot, eval = test_protti, fig.align= "center", fig.width = 6, fig.height = 5, message = FALSE, warning = FALSE----
# volcano_plot(
#   data = diff_abundance_data,
#   grouping = eg_precursor_id,
#   log2FC = diff,
#   significance = pval,
#   method = "target",
#   target_column = pg_protein_accessions,
#   target = "P62942",
#   x_axis_label = "log2(fold change) Rapamycin treated vs. untreated",
#   significance_cutoff = c(0.05, "adj_pval")
# )
# 
# # The significance_cutoff argument can also just be used for a
# # regular cutoff line by just providing the cutoff value, e.g.
# # signficiance_cutoff = 0.05

## ----barcode_plot, eval = test_protti, fig.align = "center", fig.width = 6, message = FALSE, warning = FALSE----
# FKBP12 <- diff_abundance_data %>%
#   filter(pg_protein_accessions == "P62942")
# 
# barcode_plot(
#   data = FKBP12,
#   start_position = start,
#   end_position = end,
#   protein_length = length,
#   coverage = coverage,
#   colouring = diff,
#   cutoffs = c(diff = 1, adj_pval = 0.05),
#   protein_id = pg_protein_accessions
# )

## ----woods_plot, eval = test_protti, fig.align = "center", fig.width = 6, message = FALSE, warning = FALSE----
# FKBP12 <- FKBP12 %>%
#   mutate(significant = ifelse(adj_pval < 0.01, TRUE, FALSE))
# 
# woods_plot(
#   data = FKBP12,
#   fold_change = diff,
#   start_position = start,
#   end_position = end,
#   protein_length = length,
#   coverage = coverage,
#   colouring = adj_pval,
#   protein_id = pg_protein_accessions,
#   facet = FALSE,
#   fold_change_cutoff = 1,
#   highlight = significant
# )

## ----protile_plot, eval = test_protti, fig.align = "center", fig.width = 20, fig.height = 6, message = FALSE, warning = FALSE----
# FKBP12_intensity <- data_filtered_uniprot %>%
#   filter(pg_protein_accessions == "P62942")
# 
# peptide_profile_plot(
#   data = FKBP12_intensity,
#   sample = r_file_name,
#   peptide = eg_precursor_id,
#   intensity_log2 = normalised_intensity_log2,
#   grouping = pg_protein_accessions,
#   targets = "P62942",
#   protein_abundance_plot = FALSE
# )

## ----additional_functions, eval=FALSE-----------------------------------------
# diff_abundance_significant <- diff_abundance_data %>%
#   # mark significant peptides
#   mutate(is_significant = ifelse((adj_pval < 0.01 & abs(diff) > 1), TRUE, FALSE)) %>%
#   # mark true positive hits
#   mutate(binds_treatment = pg_protein_accessions == "P62942")
# 
# ### GO enrichment using "molecular function" annotation from UniProt
# 
# calculate_go_enrichment(
#   data = diff_abundance_significant,
#   protein_id = pg_protein_accessions,
#   is_significant = is_significant,
#   go_annotations_uniprot = go_f
# )
# 
# ### Network analysis
# 
# network_input <- diff_abundance_significant %>%
#   filter(is_significant == TRUE)
# 
# analyse_functional_network(
#   data = network_input,
#   protein_id = pg_protein_accessions,
#   string_id = xref_string,
#   binds_treatment = binds_treatment,
#   organism_id = 9606
# )
# 
# ### 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
# 
# diff_abundance_significant %>%
#   # 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 = is_significant,
#     pathway_id = pathway_id,
#     pathway_name = pathway_name
#   )
# 
# ### Treatment enrichment analysis
# 
# calculate_treatment_enrichment(diff_abundance_significant,
#   protein_id = pg_protein_accessions,
#   is_significant = is_significant,
#   binds_treatment = binds_treatment,
#   treatment_name = "Rapamycin"
# )