### R code from vignette source 'SWATH2stats_vignette.Rnw'

###################################################
### code chunk number 1: SWATH2stats_vignette.Rnw:48-49
###################################################
options(width=70)


###################################################
### code chunk number 2: SWATH2stats_vignette.Rnw:54-56 (eval = FALSE)
###################################################
## source("http://www.bioconductor.org/biocLite.R")
## biocLite("SWATH2stats")


###################################################
### code chunk number 3: SWATH2stats_vignette.Rnw:60-61
###################################################
library(SWATH2stats)


###################################################
### code chunk number 4: SWATH2stats_vignette.Rnw:70-75
###################################################
data('OpenSWATH_data', package='SWATH2stats')
data <- OpenSWATH_data

data('Study_design', package='SWATH2stats')
head(Study_design)


###################################################
### code chunk number 5: SWATH2stats_vignette.Rnw:80-88 (eval = FALSE)
###################################################
## # set working directory
## setwd('~/Documents/MyWorkingDirectory/')
## 
## # Input data file (openSWATH output)
## file.name <- 'OpenSWATH_output_file.txt'
## 
## # File name for annotation file
## annotation.file <- 'Study_design_file.txt'


###################################################
### code chunk number 6: SWATH2stats_vignette.Rnw:93-95 (eval = FALSE)
###################################################
## # load data
## data <- data.frame(fread(file.name, sep='\t', header=TRUE))


###################################################
### code chunk number 7: SWATH2stats_vignette.Rnw:100-104 (eval = FALSE)
###################################################
## # consult the manual page.
## help(import_data)
## # rename the columns
## data <- import_data(data)


###################################################
### code chunk number 8: SWATH2stats_vignette.Rnw:111-116
###################################################
# reduce number of columns
data <- reduce_OpenSWATH_output(data)

# remove the iRT peptides (or other proteins)
data <- data[grep('iRT', data$ProteinName, invert=TRUE),]


###################################################
### code chunk number 9: SWATH2stats_vignette.Rnw:122-129 (eval = FALSE)
###################################################
## # list number and different Files present
## nlevels(factor(data$align_origfilename))
## levels(factor(data$align_origfilename))
## 
## # load the study design table from the indicated file
## Study_design <- read.delim2(annotation.file,
##                             dec='.', sep='\t', header=TRUE)


###################################################
### code chunk number 10: SWATH2stats_vignette.Rnw:134-143
###################################################
# annotate data
data.annotated <- sample_annotation(data, Study_design)

head(unique(data$ProteinName))
# OPTIONAL: for human, shorten Protein Name to remove non-unique information
#(sp|Q9GZL7|WDR12_HUMAN --> Q9GZL7)
data$ProteinName <- gsub('sp\\|([[:alnum:]]+)\\|[[:alnum:]]*_HUMAN',
                         '\\1', data$ProteinName)
head(unique(data$ProteinName))


###################################################
### code chunk number 11: SWATH2stats_vignette.Rnw:152-153
###################################################
count_analytes(data.annotated)


###################################################
### code chunk number 12: SWATH2stats_vignette.Rnw:159-166
###################################################
# Plot correlation of intensity
correlation <- plot_correlation_between_samples(data.annotated, column.values = "Intensity")
head(correlation)

# Plot correlation of retention times
correlation <- plot_correlation_between_samples(data.annotated, column.values = "RT")



###################################################
### code chunk number 13: SWATH2stats_vignette.Rnw:172-181
###################################################
# plot variation of transitions for each condition across replicates
variation <- plot_variation(data.annotated)
head(variation[[2]])

# plot variation of summed signal for Peptides for each condition across replicates
variation <- plot_variation(data.annotated,
               Comparison = FullPeptideName + Condition ~ BioReplicate,
               fun.aggregate = sum)



###################################################
### code chunk number 14: SWATH2stats_vignette.Rnw:187-191
###################################################
# plot variation of transitions for each condition within replicates
# compared to total variation
variation <- plot_variation_vs_total(data.annotated)
head(variation[[2]])


###################################################
### code chunk number 15: SWATH2stats_vignette.Rnw:201-204 (eval = FALSE)
###################################################
## protein_matrix <- write_matrix_proteins(data,
##                       filename = "SWATH2stats_overview_matrix_proteinlevel",
##                       rm.decoy = FALSE)


###################################################
### code chunk number 16: SWATH2stats_vignette.Rnw:211-214 (eval = FALSE)
###################################################
## peptide_matrix <- write_matrix_peptides(data,
##                       filename = "SWATH2stats_overview_matrix_peptidelevel",
##                       rm.decoy = FALSE)


###################################################
### code chunk number 17: SWATH2stats_vignette.Rnw:236-237
###################################################
assess_decoy_rate(data)


###################################################
### code chunk number 18: SWATH2stats_vignette.Rnw:242-250
###################################################
# count decoys and targets on assay, peptide and protein level
# and report FDR at a range of m_score cutoffs
assess_fdr_overall(data, FFT = 0.7, output = "pdf_csv", plot = TRUE,
                   filename='assess_fdr_overall_testrun')

# The results can be reported back to R for further calculations
overall_fdr_table <- assess_fdr_overall(data, FFT = 0.7,
                                        output = "Rconsole")


###################################################
### code chunk number 19: SWATH2stats_vignette.Rnw:255-258
###################################################
# create plots from fdr_table
plot(overall_fdr_table, output = "Rconsole",
               filename = "FDR_report_overall")


###################################################
### code chunk number 20: SWATH2stats_vignette.Rnw:263-271
###################################################
# count decoys and targets on assay, peptide and protein level per run
# and report FDR at a range of m_score cutoffs
assess_fdr_byrun(data, FFT = 0.7, output = "pdf_csv", plot = TRUE,
                 filename='assess_fdr_byrun_testrun')

# The results can be reported back to R for further calculations
byrun_fdr_cube <- assess_fdr_byrun(data, FFT = 0.7,
                                   output = "Rconsole")


###################################################
### code chunk number 21: SWATH2stats_vignette.Rnw:276-279
###################################################
# create plots from fdr_table
plot(byrun_fdr_cube, output = "Rconsole",
              filename = "FDR_report_overall")


###################################################
### code chunk number 22: SWATH2stats_vignette.Rnw:287-290
###################################################
# select and return a useful m_score cutoff in order
# to achieve the desired FDR quality for the entire table
mscore4assayfdr(data, FFT = 0.7, fdr_target=0.01)


###################################################
### code chunk number 23: SWATH2stats_vignette.Rnw:295-298
###################################################
# select and return a useful m_score cutoff
# in order to achieve the desired FDR quality for the entire table
mscore4pepfdr(data, FFT = 0.7, fdr_target=0.02)


###################################################
### code chunk number 24: SWATH2stats_vignette.Rnw:305-308
###################################################
# select and return a useful m_score cutoff in order
# to achieve the desired FDR quality for the entire table
mscore4protfdr(data, FFT = 0.7, fdr_target=0.02)


###################################################
### code chunk number 25: SWATH2stats_vignette.Rnw:320-321
###################################################
data.filtered.mscore <- filter_mscore(data.annotated, 0.01)


###################################################
### code chunk number 26: SWATH2stats_vignette.Rnw:326-328
###################################################
data.filtered.mscore <- filter_mscore_freqobs(data.annotated, 0.01, 0.8,
                                              rm.decoy=FALSE)


###################################################
### code chunk number 27: SWATH2stats_vignette.Rnw:333-334 (eval = FALSE)
###################################################
## data.filtered.mscore <- filter_mscore_condition(data.annotated, 0.01, 3)


###################################################
### code chunk number 28: SWATH2stats_vignette.Rnw:341-344
###################################################
data.filtered.fdr <- filter_mscore_fdr(data, FFT=0.7,
                                   overall_protein_fdr_target = 0.03,
                                   upper_overall_peptide_fdr_limit = 0.05)


###################################################
### code chunk number 29: SWATH2stats_vignette.Rnw:350-352
###################################################
data <- filter_proteotypic_peptides(data.filtered.mscore)
data.all <- filter_all_peptides(data.filtered.mscore)


###################################################
### code chunk number 30: SWATH2stats_vignette.Rnw:358-359
###################################################
data.filtered.max <- filter_on_max_peptides(data.filtered.mscore, 5)


###################################################
### code chunk number 31: SWATH2stats_vignette.Rnw:364-365
###################################################
data.filtered.max.min <- filter_on_min_peptides(data.filtered.max, 2)


###################################################
### code chunk number 32: SWATH2stats_vignette.Rnw:376-377
###################################################
data.transition <- disaggregate(data)


###################################################
### code chunk number 33: SWATH2stats_vignette.Rnw:380-382 (eval = FALSE)
###################################################
## write.csv(data.transition, file='transition_level_output.csv',
##           row.names=FALSE, quote=FALSE)


###################################################
### code chunk number 34: SWATH2stats_vignette.Rnw:388-390
###################################################
data <- convert4pythonscript(data)
head(data)


###################################################
### code chunk number 35: SWATH2stats_vignette.Rnw:393-394 (eval = FALSE)
###################################################
## write.table(data, file="input.tsv", sep="\t", row.names=FALSE, quote=FALSE)


###################################################
### code chunk number 36: SWATH2stats_vignette.Rnw:402-404 (eval = FALSE)
###################################################
## data.transition <- data.frame(fread('output.csv',
##                                     sep=',', header=TRUE))


###################################################
### code chunk number 37: SWATH2stats_vignette.Rnw:411-415 (eval = FALSE)
###################################################
## MSstats.input <- convert4MSstats(data.transition)
## 
## library(MSstats)
## quantData <- dataProcess(MSstats.input)


###################################################
### code chunk number 38: SWATH2stats_vignette.Rnw:422-434 (eval = FALSE)
###################################################
## aLFQ.input <- convert4aLFQ(data.transition)
## 
## library(aLFQ)
## prots <- ProteinInference(aLFQ.input, peptide_method = 'top',
##                           peptide_topx = 3,
##                           peptide_strictness = 'loose',
##                           peptide_summary = 'mean',
##                           transition_topx = 3,
##                           transition_strictness = 'loose',
##                           transition_summary = 'sum',
##                           fasta = NA, model = NA,
##                           combine_precursors = FALSE)


###################################################
### code chunk number 39: SWATH2stats_vignette.Rnw:439-441
###################################################
mapDIA.input <- convert4mapDIA(data.transition, RT =TRUE)
head(mapDIA.input)


###################################################
### code chunk number 40: SWATH2stats_vignette.Rnw:444-446 (eval = FALSE)
###################################################
## write.table(mapDIA.input, file='mapDIA.txt', quote=FALSE,
##           row.names=FALSE, sep='\t')


###################################################
### code chunk number 41: SWATH2stats_vignette.Rnw:453-460 (eval = FALSE)
###################################################
## data.annotated.full <- sample_annotation(OpenSWATH_data, Study_design)
## data.annotated.full <- filter_mscore(data.annotated.full,
##                                      mscore4assayfdr(data.annotated.full, 0.01))
## data.annotated.full$decoy <- 0 ### imsbInfer needs the decoy column
## 
## library(imsbInfer)
## specLib <- loadTransitonsMSExperiment(data.annotated.full)