### R code from vignette source 'SWATH2stats_vignette.Rnw' ################################################### ### code chunk number 1: SWATH2stats_vignette.Rnw:49-50 ################################################### options(width=70) ################################################### ### code chunk number 2: SWATH2stats_vignette.Rnw:55-57 (eval = FALSE) ################################################### ## source("http://www.bioconductor.org/biocLite.R") ## biocLite("SWATH2stats") ################################################### ### code chunk number 3: SWATH2stats_vignette.Rnw:61-63 ################################################### rm(list=ls()) 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) ################################################### ## n# 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:101-106 ################################################### # 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:112-119 (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(file.path(getwd(), annotation.file), ## dec='.', sep='\t', header=TRUE) ################################################### ### code chunk number 9: SWATH2stats_vignette.Rnw:124-133 ################################################### # 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:152-155 ################################################### data('OpenSWATH_data_FDR', package='SWATH2stats') data.FDR<-sample_annotation(OpenSWATH_data_FDR, Study_design) assess_decoy_rate(data.FDR) ################################################### ### code chunk number 11: SWATH2stats_vignette.Rnw:160-168 ################################################### # count decoys and targets on assay, peptide and protein level # and report FDR at a range of m_score cutoffs assess_fdr_overall(data.FDR, 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.FDR, FFT = 0.7, output = "Rconsole") ################################################### ### code chunk number 12: SWATH2stats_vignette.Rnw:173-176 ################################################### # create plots from fdr_table plot(overall_fdr_table, output = "Rconsole", filename = "FDR_report_overall") ################################################### ### code chunk number 13: SWATH2stats_vignette.Rnw:181-189 ################################################### # 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.FDR, 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.FDR, FFT = 0.7, output = "Rconsole") ################################################### ### code chunk number 14: SWATH2stats_vignette.Rnw:194-197 ################################################### # create plots from fdr_table plot(byrun_fdr_cube, output = "Rconsole", filename = "FDR_report_overall") ################################################### ### code chunk number 15: SWATH2stats_vignette.Rnw:205-208 ################################################### # select and return a useful m_score cutoff in order # to achieve the desired FDR quality for the entire table mscore4assayfdr(data.FDR, FFT = 0.7, fdr_target=0.02) ################################################### ### code chunk number 16: SWATH2stats_vignette.Rnw:213-216 ################################################### # select and return a useful m_score cutoff # in order to achieve the desired FDR quality for the entire table mscore4pepfdr(data.FDR, FFT = 0.7, fdr_target=0.02) ################################################### ### code chunk number 17: SWATH2stats_vignette.Rnw:223-226 ################################################### # select and return a useful m_score cutoff in order # to achieve the desired FDR quality for the entire table mscore4protfdr(data.FDR, FFT = 0.7, fdr_target=0.02) ################################################### ### code chunk number 18: SWATH2stats_vignette.Rnw:238-239 ################################################### data.filtered.mscore <- filter_mscore(data.annotated, 0.01) ################################################### ### code chunk number 19: SWATH2stats_vignette.Rnw:244-246 ################################################### data.filtered.mscore <- filter_mscore_requant(data.annotated, 0.01, 0.8, rm.decoy=FALSE) ################################################### ### code chunk number 20: SWATH2stats_vignette.Rnw:251-252 (eval = FALSE) ################################################### ## data.filtered.mscore <- filter_mscore_condition(data.annotated, 0.01, 3) ################################################### ### code chunk number 21: SWATH2stats_vignette.Rnw:259-262 ################################################### data.filtered.fdr <- filter_mscore_fdr(data.FDR, FFT=0.7, overall_protein_fdr_target = 0.03, upper_overall_peptide_fdr_limit = 0.05) ################################################### ### code chunk number 22: SWATH2stats_vignette.Rnw:268-270 ################################################### data <- filter_proteotypic_peptides(data.filtered.mscore) data.all <- filter_all_peptides(data.filtered.mscore) ################################################### ### code chunk number 23: SWATH2stats_vignette.Rnw:276-277 ################################################### data.filtered.max <- filter_on_max_peptides(data.filtered.mscore, 5) ################################################### ### code chunk number 24: SWATH2stats_vignette.Rnw:282-283 ################################################### data.filtered.max.min <- filter_on_min_peptides(data.filtered.max, 2) ################################################### ### code chunk number 25: SWATH2stats_vignette.Rnw:297-299 (eval = FALSE) ################################################### ## write_matrix_proteins(data, filename = "SWATH2stats_overview_matrix_proteinlevel", ## rm.decoy = FALSE) ################################################### ### code chunk number 26: SWATH2stats_vignette.Rnw:305-308 ################################################### data.peptide <- data data.peptide$aggr_Fragment_Annotation <- NULL data.peptide$aggr_Peak_Area <- NULL ################################################### ### code chunk number 27: SWATH2stats_vignette.Rnw:311-313 (eval = FALSE) ################################################### ## write.csv(data.peptide, file='peptide_level_output.csv', ## row.names=FALSE, quote=FALSE) ################################################### ### code chunk number 28: SWATH2stats_vignette.Rnw:318-321 (eval = FALSE) ################################################### ## write_matrix_peptides(data, ## filename = "SWATH2stats_overview_matrix_peptidelevel", ## rm.decoy = FALSE) ################################################### ### code chunk number 29: SWATH2stats_vignette.Rnw:327-328 ################################################### data.transition <- disaggregate(data) ################################################### ### code chunk number 30: SWATH2stats_vignette.Rnw:331-333 (eval = FALSE) ################################################### ## write.csv(data.transition, file='transition_level_output.csv', ## row.names=FALSE, quote=FALSE) ################################################### ### code chunk number 31: SWATH2stats_vignette.Rnw:339-342 (eval = FALSE) ################################################### ## data <- convert4pythonscript(data) ## write.table(data, file="input.tsv", sep="\t", row.names=FALSE, quote=FALSE) ## head(data) ################################################### ### code chunk number 32: SWATH2stats_vignette.Rnw:350-352 (eval = FALSE) ################################################### ## data.transition <- data.frame(fread('output.csv', ## sep=',', header=TRUE)) ################################################### ### code chunk number 33: SWATH2stats_vignette.Rnw:359-363 ################################################### MSstats.input <- convert4MSstats(data.transition) library(MSstats) quantData <- dataProcess(MSstats.input) ################################################### ### code chunk number 34: SWATH2stats_vignette.Rnw:370-382 ################################################### 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 35: SWATH2stats_vignette.Rnw:387-389 ################################################### mapDIA.input <- convert4mapDIA(data.transition, RT =TRUE) head(mapDIA.input) ################################################### ### code chunk number 36: SWATH2stats_vignette.Rnw:392-394 (eval = FALSE) ################################################### ## write.table(mapDIA.input, file='mapDIA.txt', quote=FALSE, ## row.names=FALSE, sep='\t') ################################################### ### code chunk number 37: SWATH2stats_vignette.Rnw:401-408 (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)