## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

ggplot2::theme_set(
  ggplot2::theme_minimal()
)

library(mspms)
library(dplyr)

## ----eval = FALSE-------------------------------------------------------------
# mspms::generate_report(
#   prepared_data = mspms::peaks_prepared_data,
#   outdir = "../Desktop/mspms_report"
# )

## -----------------------------------------------------------------------------
colData <- readr::read_csv(system.file("extdata/colData.csv",
  package = "mspms"
))

head(colData)

## -----------------------------------------------------------------------------

library(dplyr)
library(mspms)

# File path of peaks lfq file
lfq_filepath <- system.file("extdata/peaks_protein-peptides-lfq.csv",
  package = "mspms"
)

colData_filepath <- system.file("extdata/colData.csv", package = "mspms")

# Prepare the data
peaks_prepared_data <- mspms::prepare_peaks(lfq_filepath,
  colData_filepath,
  quality_threshold = 0.3,
  n_residues = 4
)

peaks_prepared_data


## -----------------------------------------------------------------------------
combined_peptide_filepath <- system.file("extdata/fragpipe_combined_peptide.tsv",
  package = "mspms"
)

colData_filepath <- system.file("extdata/colData.csv", package = "mspms")


fragpipe_prepared_data <- prepare_fragpipe(
  combined_peptide_filepath,
  colData_filepath
)

fragpipe_prepared_data

## -----------------------------------------------------------------------------
peptide_groups_filepath <- system.file(
  "extdata/proteome_discoverer_PeptideGroups.txt",
  package = "mspms"
)

colData_filepath <- system.file("extdata/proteome_discover_colData.csv",
  package = "mspms"
)

prepared_pd_data <- prepare_pd(peptide_groups_filepath, colData_filepath)

prepared_pd_data

## -----------------------------------------------------------------------------
set.seed(2)
processed_qf <- process_qf(peaks_prepared_data)
processed_qf

## -----------------------------------------------------------------------------
log2fc_t_test_data <- mspms::log2fc_t_test(processed_qf = processed_qf)

## -----------------------------------------------------------------------------
plot_qc_check(processed_qf,
  full_length_threshold = 10,
  cleavage_product_threshold = 5
)

## -----------------------------------------------------------------------------
plot_nd_peptides(processed_qf)

## -----------------------------------------------------------------------------
mspms_tidy_data <- mspms_tidy(processed_qf)

## ----eval = FALSE-------------------------------------------------------------
# plot_heatmap(mspms_tidy_data)

## -----------------------------------------------------------------------------
plot_pca(mspms_tidy_data, value_colname = "peptides_norm")

## -----------------------------------------------------------------------------
plot_volcano(log2fc_t_test_data)

## -----------------------------------------------------------------------------
sig_cleavage_data <- log2fc_t_test_data %>%
  dplyr::filter(p.adj <= 0.05, log2fc > 3)

p1 <- mspms::plot_cleavages_per_pos(sig_cleavage_data)

p1

## -----------------------------------------------------------------------------
catA_sig_cleavages <- log2fc_t_test_data %>%
  dplyr::filter(p.adj <= 0.05, log2fc > 3) %>%
  dplyr::filter(condition == "CatA") %>%
  dplyr::pull(cleavage_seq) %>%
  unique()

## -----------------------------------------------------------------------------
all_possible_8mers_from_228_library <- calculate_all_cleavages(
  mspms::peptide_library$library_real_sequence,
  n_AA_after_cleavage = 4
)

## -----------------------------------------------------------------------------
plot_icelogo(catA_sig_cleavages,
  background_universe = all_possible_8mers_from_228_library
)

## -----------------------------------------------------------------------------
sig_cleavage_data <- log2fc_t_test_data %>%
  dplyr::filter(p.adj <= 0.05, log2fc > 3)

plot_all_icelogos(sig_cleavage_data)

## -----------------------------------------------------------------------------
max_log2fc_pep <- log2fc_t_test_data %>%
  dplyr::filter(p.adj <= 0.05, log2fc > 3) %>%
  dplyr::filter(log2fc == max(log2fc)) %>%
  pull(peptide)

## -----------------------------------------------------------------------------
p1 <- mspms_tidy_data %>%
  dplyr::filter(peptide == max_log2fc_pep) %>%
  plot_time_course()

p1

## -----------------------------------------------------------------------------
sessionInfo()