## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=8, 
  fig.height=8
)

## ---- eval = FALSE------------------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly = TRUE))
#      install.packages("BiocManager")
#  
#  BiocManager::install("MSstatsPTM")

## ---- message=FALSE, warning=FALSE--------------------------------------------
library(MSstatsPTM)
library(MSstats)

## ----raw_data_ptm, eval = FALSE-----------------------------------------------
#  
#  annotation = data.frame('Condition' = c('Control', 'Control', 'Control',
#                                           'Treatment', 'Treatment', 'Treatment'),
#                           'BioReplicate' = c(1,2,3,4,5,6),
#                           'Run' = c('prot_run_1', 'prot_run_2', 'prot_run_3',
#                                     'phos_run_1', 'phos_run_2', 'phos_run_3'),
#                           'Type' = c("Protein", "Protein", "Protein", "PTM",
#                                      "PTM", "PTM"))
#  
#  # Run MSstatsPTM converter with modified and unmodified datasets.
#  raw.input = ProgenesistoMSstatsPTMFormat(raw_ptm_df, annotation,
#                                            raw_protein_df, fasta_path)

## ----raw_data, eval = FALSE---------------------------------------------------
#  # Add site into ProteinName column
#  raw_ptm_df$ProteinName = paste(raw_ptm_df$ProteinName,
#                                  raw_ptm_df$Site, sep = "_")
#  
#  # Run MSstats Converters
#  PTM.data = MSstats::ProgenesistoMSstatsFormat(raw_ptm_df, annotation)
#  PROTEIN.data = MSstats::ProgenesistoMSstatsFormat(raw_protein_df, annotation)
#  
#  # Combine into one list
#  raw.input = list(PTM = PTM.data,
#                    PROTEIN = PROTEIN.data)

## -----------------------------------------------------------------------------
head(raw.input$PTM)
head(raw.input$PROTEIN)

## ----summarize, message=FALSE, warning=FALSE----------------------------------

MSstatsPTM.summary = dataSummarizationPTM(raw.input, verbose = FALSE, 
                                          use_log_file = FALSE, append = FALSE)

head(MSstatsPTM.summary$PTM$ProteinLevelData)
head(MSstatsPTM.summary$PROTEIN$ProteinLevelData)

## ----qcplot, message=FALSE, warning=FALSE-------------------------------------
dataProcessPlotsPTM(MSstatsPTM.summary,
                    type = 'QCPLOT',
                    which.PTM = "allonly",
                    address = FALSE)

## ----profileplot, message=FALSE, warning=FALSE--------------------------------

dataProcessPlotsPTM(MSstatsPTM.summary,
                    type = 'ProfilePlot',
                    which.Protein = "Q9Y6C9",
                    address = FALSE)

## ----model, message=FALSE, warning=FALSE--------------------------------------

# Specify contrast matrix
comparison = matrix(c(-1,0,1,0),nrow=1)
row.names(comparison) = "CCCP-Ctrl"
colnames(comparison) = c("CCCP", "Combo", "Ctrl", "USP30_OE")

MSstatsPTM.model = groupComparisonPTM(MSstatsPTM.summary, 
                                      data.type = "LabelFree",
                                      contrast.matrix = comparison,
                                      use_log_file = FALSE, append = FALSE,
                                      verbose = FALSE)
head(MSstatsPTM.model$PTM.Model)
head(MSstatsPTM.model$PROTEIN.Model)
head(MSstatsPTM.model$ADJUSTED.Model)

## ----volcano, message=FALSE, warning=FALSE------------------------------------
groupComparisonPlotsPTM(data = MSstatsPTM.model,
                        type = "VolcanoPlot",
                        FCcutoff= 2,
                        logBase.pvalue = 2,
                        address=FALSE)

## ----heatmap, message=FALSE, warning=FALSE------------------------------------
groupComparisonPlotsPTM(data = MSstatsPTM.model,
                        type = "Heatmap",
                        which.PTM = 1:30,
                        address=FALSE)

## ----sample_size, message=FALSE, warning=FALSE--------------------------------

# Specify contrast matrix
sample_size = designSampleSizePTM(MSstatsPTM.model, c(2.0, 2.75), FDR = 0.05, 
                                  numSample = TRUE, power = 0.8)

head(sample_size)

## ----sample_size_plot, message=FALSE, warning=FALSE---------------------------

MSstats::designSampleSizePlots(sample_size)