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