diff filterFragSpectra.R @ 0:f52287a06c02 draft

"planemo upload for repository https://github.com/computational-metabolomics/mspurity-galaxy commit cb903cd93f9378cfb5eeb68512a54178dcea7bbc-dirty"

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/filterFragSpectra.R	Wed Nov 27 13:44:58 2019 -0500
@@ -0,0 +1,100 @@
+library(optparse)
+library(msPurity)
+library(xcms)
+print(sessionInfo())
+
+
+option_list <- list(
+  make_option("--out_rdata", type="character"),
+  make_option("--out_peaklist_prec", type="character"),
+  make_option("--out_peaklist_frag", type="character"),
+  make_option("--pa", type="character"),
+
+  make_option("--ilim", default=0.0),
+  make_option("--plim", default=0.0),
+
+  make_option("--ra", default=0.0),  
+  make_option("--snr", default=0.0),
+
+  make_option("--rmp", action="store_true"),
+  make_option("--snmeth", default="median", type="character")
+)
+
+opt <- parse_args(OptionParser(option_list=option_list))
+print(opt)
+
+
+loadRData <- function(rdata_path, name){
+    #loads an RData file, and returns the named xset object if it is there
+    load(rdata_path)
+    return(get(ls()[ls() %in% name]))
+}
+
+# Requires
+pa <- loadRData(opt$pa, 'pa')
+
+if(is.null(opt$rmp)){
+  opt$rmp = FALSE
+}else{
+  opt$rmp = TRUE
+}
+
+pa <- filterFragSpectra(pa, 
+                        ilim=opt$ilim, 
+                        plim=opt$plim,
+                        ra=opt$ra,
+                        snr=opt$snr,
+                        rmp=opt$rmp,
+                        snmeth=opt$snmeth)
+
+print(pa)
+save(pa, file=opt$out_rdata)
+
+# get the msms data for grpid from the purityA object
+msmsgrp <- function(grpid, pa){
+  msms <- pa@grped_ms2[grpid]
+  
+  grpinfo <- pa@grped_df[pa@grped_df$grpid==grpid,]
+  
+  grpinfo$subsetid <- 1:nrow(grpinfo)
+  result <- plyr::ddply(grpinfo, ~subsetid, setid, msms=msms)
+  return(result)
+}
+
+# Set the relevant details 
+setid <- function(grpinfo_i, msms){
+  msms_i <- msms[[1]][[grpinfo_i$subsetid]]
+  n <- nrow(msms_i)
+  msms_i <- data.frame(msms_i)
+  colnames(msms_i)[1:2] <- c('mz', 'i')
+  m <- cbind('grpid'=rep(grpinfo_i$grpid,n), 'pid'=rep(grpinfo_i$pid,n), 'fileid'=rep(grpinfo_i$fileid,n), msms_i)
+  return(m)
+}
+
+
+
+if (length(pa)>0){
+
+  if (length(pa@grped_ms2)==0){
+    message('No spectra available')
+  } else{
+
+    # get group ids
+    grpids <- unique(as.character(pa@grped_df$grpid))
+
+    # loop through all the group ids
+    df_fragments = plyr::adply(grpids, 1, msmsgrp, pa=pa)
+    df_fragments = merge(df_fragments, pa@puritydf[,c("pid", "acquisitionNum", "precursorScanNum")], by="pid")
+    df_fragments = df_fragments[order(df_fragments$grpid, df_fragments$pid, df_fragments$mz),]
+    #select and reorder columns
+    df_fragments = df_fragments[,c("grpid", "pid", "precursorScanNum", "acquisitionNum", "fileid", "mz", "i", "snr", "ra", "purity_pass_flag", "intensity_pass_flag", "ra_pass_flag", "snr_pass_flag", "pass_flag")]
+
+    pa@grped_df$filename = sapply(pa@grped_df$fileid, function(x) names(pa@fileList)[as.integer(x)])
+
+    print(head(pa@grped_df))    
+    write.table(pa@grped_df, opt$out_peaklist_prec, row.names=FALSE, sep='\t')
+    print(head(df_fragments))
+    write.table(df_fragments, opt$out_peaklist_frag, row.names=FALSE, sep='\t')
+  }
+}
+