### 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-91 (eval = FALSE)
###################################################
## # set working directory
## setwd('~/Documents/MyWorkingDirectory/')
## 
## # Define input data file (e.g. OpenSWATH_output_file.txt)
## file.name <- 'OpenSWATH_output_file.txt'
## 
## # File name for annotation file
## annotation.file <- 'Study_design_file.txt'
## 
## # load data
## data <- data.frame(fread(file.name, sep='\t', header=TRUE))


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


###################################################
### code chunk number 7: SWATH2stats_vignette.Rnw:107-112
###################################################
# 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 8: SWATH2stats_vignette.Rnw:118-125 (eval = FALSE)
###################################################
## # list number and different Files present
## nlevels(factor(data$filename))
## levels(factor(data$filename))
## 
## # load the study design table from the indicated file
## Study_design <- read.delim2(annotation.file,
##                             dec='.', sep='\t', header=TRUE)


###################################################
### code chunk number 9: SWATH2stats_vignette.Rnw:130-139
###################################################
# 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 10: SWATH2stats_vignette.Rnw:148-149
###################################################
count_analytes(data.annotated)


###################################################
### code chunk number 11: SWATH2stats_vignette.Rnw:155-162
###################################################
# 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 12: SWATH2stats_vignette.Rnw:168-177
###################################################
# 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 13: SWATH2stats_vignette.Rnw:183-187
###################################################
# 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 14: SWATH2stats_vignette.Rnw:197-200 (eval = FALSE)
###################################################
## protein_matrix <- write_matrix_proteins(data,
##                       filename = "SWATH2stats_overview_matrix_proteinlevel",
##                       rm.decoy = FALSE)


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


###################################################
### code chunk number 16: SWATH2stats_vignette.Rnw:232-233
###################################################
assess_decoy_rate(data)


###################################################
### code chunk number 17: SWATH2stats_vignette.Rnw:238-246
###################################################
# 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 18: SWATH2stats_vignette.Rnw:251-254
###################################################
# create plots from fdr_table
plot(overall_fdr_table, output = "Rconsole",
               filename = "FDR_report_overall")


###################################################
### code chunk number 19: SWATH2stats_vignette.Rnw:259-267
###################################################
# 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 20: SWATH2stats_vignette.Rnw:272-275
###################################################
# create plots from fdr_table
plot(byrun_fdr_cube, output = "Rconsole",
              filename = "FDR_report_overall")


###################################################
### code chunk number 21: SWATH2stats_vignette.Rnw:283-286
###################################################
# 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 22: SWATH2stats_vignette.Rnw:291-294
###################################################
# 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 23: SWATH2stats_vignette.Rnw:301-304
###################################################
# 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 24: SWATH2stats_vignette.Rnw:316-317
###################################################
data.filtered.mscore <- filter_mscore(data.annotated, 0.01)


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


###################################################
### code chunk number 26: SWATH2stats_vignette.Rnw:329-330 (eval = FALSE)
###################################################
## data.filtered.mscore <- filter_mscore_condition(data.annotated, 0.01, 3)


###################################################
### code chunk number 27: SWATH2stats_vignette.Rnw:337-340
###################################################
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 28: SWATH2stats_vignette.Rnw:346-348
###################################################
data <- filter_proteotypic_peptides(data.filtered.mscore)
data.all <- filter_all_peptides(data.filtered.mscore)


###################################################
### code chunk number 29: SWATH2stats_vignette.Rnw:354-355
###################################################
data.filtered.max <- filter_on_max_peptides(data.filtered.mscore, 5)


###################################################
### code chunk number 30: SWATH2stats_vignette.Rnw:360-361
###################################################
data.filtered.max.min <- filter_on_min_peptides(data.filtered.max, 2)


###################################################
### code chunk number 31: SWATH2stats_vignette.Rnw:373-375
###################################################
data.transition <- disaggregate(data)
head(data.transition)


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


###################################################
### code chunk number 33: SWATH2stats_vignette.Rnw:386-387
###################################################
data.python <- convert4pythonscript(data)


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


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


###################################################
### code chunk number 36: SWATH2stats_vignette.Rnw:408-412 (eval = FALSE)
###################################################
## MSstats.input <- convert4MSstats(data.transition)
## 
## library(MSstats)
## quantData <- dataProcess(MSstats.input)


###################################################
### code chunk number 37: SWATH2stats_vignette.Rnw:419-431
###################################################
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 38: SWATH2stats_vignette.Rnw:436-438
###################################################
mapDIA.input <- convert4mapDIA(data.transition, RT =TRUE)
head(mapDIA.input)


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


###################################################
### code chunk number 40: SWATH2stats_vignette.Rnw:448-450
###################################################
PECA.input <- convert4PECA(data)
head(PECA.input)


###################################################
### code chunk number 41: SWATH2stats_vignette.Rnw:455-462
###################################################
library(PECA)
group1 <- c("Hela_Control_1", "Hela_Control_2", "Hela_Control_3")
group2 <- c("Hela_Treatment_1", "Hela_Treatment_2", "Hela_Treatment_3")

# PECA_df
results <- PECA_df(PECA.input, group1, group2, id="ProteinName", test = "rots")
head(results)


###################################################
### code chunk number 42: SWATH2stats_vignette.Rnw:468-475 (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)