SWATH2stats_vignette.R
In SWATH2stats: Transform and Filter SWATH Data for Statistical Packages

### 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-57 (eval = FALSE)
###################################################
## if (!require("BiocManager"))
##     install.packages("BiocManager")
## BiocManager::install("SWATH2stats")


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


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

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


###################################################
### code chunk number 5: SWATH2stats_vignette.Rnw:81-92 (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:97-101 (eval = FALSE)
###################################################
## # consult the manual page.
## help(import_data)
## # rename the columns
## data <- import_data(data)


###################################################
### code chunk number 7: SWATH2stats_vignette.Rnw:108-113
###################################################
# 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:119-126 (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:131-140
###################################################
# 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:149-150
###################################################
count_analytes(data.annotated)


###################################################
### code chunk number 11: SWATH2stats_vignette.Rnw:156-163
###################################################
# 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:169-178
###################################################
# 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:184-188
###################################################
# 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:198-201 (eval = FALSE)
###################################################
## protein_matrix <- write_matrix_proteins(data,
##                       filename = "SWATH2stats_overview_matrix_proteinlevel",
##                       rm.decoy = FALSE)


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


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


###################################################
### code chunk number 17: SWATH2stats_vignette.Rnw:239-247
###################################################
# 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:252-255
###################################################
# create plots from fdr_table
plot(overall_fdr_table, output = "Rconsole",
               filename = "FDR_report_overall")


###################################################
### code chunk number 19: SWATH2stats_vignette.Rnw:260-268
###################################################
# 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:273-276
###################################################
# create plots from fdr_table
plot(byrun_fdr_cube, output = "Rconsole",
              filename = "FDR_report_overall")


###################################################
### code chunk number 21: SWATH2stats_vignette.Rnw:284-287
###################################################
# 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:292-295
###################################################
# 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:302-305
###################################################
# 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:317-318
###################################################
data.filtered.mscore <- filter_mscore(data.annotated, 0.01)


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


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


###################################################
### code chunk number 27: SWATH2stats_vignette.Rnw:338-341
###################################################
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:347-349
###################################################
data <- filter_proteotypic_peptides(data.filtered.mscore)
data.all <- filter_all_peptides(data.filtered.mscore)


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


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


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


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


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


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


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


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


###################################################
### code chunk number 37: SWATH2stats_vignette.Rnw:420-432 (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 38: SWATH2stats_vignette.Rnw:437-439
###################################################
mapDIA.input <- convert4mapDIA(data.transition, RT =TRUE)
head(mapDIA.input)


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


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


###################################################
### code chunk number 41: SWATH2stats_vignette.Rnw:456-462 (eval = FALSE)
###################################################
## 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")


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