diff nmr_preprocessing/NmrPreprocessing_wrapper.R @ 2:7304ec2c9ab7 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/nmr_preprocessing/NmrPreprocessing_wrapper.R	Mon Jul 30 10:33:03 2018 -0400
@@ -0,0 +1,411 @@
+#!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file
+
+## 170116_NmrPreprocessing.R
+## Manon Martin and Marie Tremblay-Franco
+
+##======================================================
+##======================================================
+# Preamble
+##======================================================
+##======================================================
+
+runExampleL <- FALSE
+
+
+##------------------------------
+## Options
+##------------------------------
+strAsFacL <- options()$stringsAsFactors
+options(stringsAsFactors = FALSE)
+
+##------------------------------
+## Libraries laoding
+##------------------------------
+library(batch)
+library(ptw)
+library(Matrix)
+library(ggplot2)
+library(gridExtra)
+library(reshape2)
+
+
+# R script call
+source_local <- function(fname)
+{
+	argv <- commandArgs(trailingOnly = FALSE)
+	base_dir <- dirname(substring(argv[grep("--file=", argv)], 8))
+	source(paste(base_dir, fname, sep="/"))
+}
+#Import the different functions
+source_local("NmrPreprocessing_script.R")
+source_local("DrawFunctions.R")
+
+##------------------------------
+## Script
+##------------------------------
+runExampleL <- FALSE
+
+
+if(!runExampleL)
+  argLs <- parseCommandArgs(evaluate=FALSE)
+
+sink(argLs$logOut)
+
+
+##------------------------------
+## Errors ?????????????????????
+##------------------------------
+
+
+##------------------------------
+## Constants
+##------------------------------
+topEnvC <- environment()
+flagC <- "\n"
+
+
+
+
+# log file
+# print(argLs[["logOut"]])
+
+## Starting
+cat("\nStart of 'Preprocessing' Galaxy module call: ", as.character(Sys.time()), "\n", sep = "")
+
+
+##======================================================
+##======================================================
+## Parameters Loading
+##======================================================
+##======================================================
+
+# graphical inputs
+FirstOPCGraph <- argLs[["FirstOPCGraph"]]
+SSGraph <- argLs[["SSGraph"]]
+ApodGraph <- argLs[["ApodGraph"]]
+FTGraph <- argLs[["FTGraph"]]
+SRGraph <- argLs[["SRGraph"]]
+ZeroOPCGraph <- argLs[["ZeroOPCGraph"]]
+BCGraph <- argLs[["BCGraph"]]
+FinalGraph <- argLs[["FinalGraph"]]
+
+
+# 1rst order phase correction ------------------------
+  # Inputs
+	## Data matrix
+Fid_data0 <- read.table(argLs[["dataMatrixFid"]],header=TRUE, check.names=FALSE, sep='\t')
+# Fid_data0 <- Fid_data0[,-1]
+Fid_data0 <- as.matrix(Fid_data0)
+
+	## Samplemetadata
+samplemetadataFid <- read.table(argLs[["sampleMetadataFid"]],check.names=FALSE,header=TRUE,sep="\t")
+samplemetadataFid <- as.matrix(samplemetadataFid)
+
+
+# water and solvent(s) correction ------------------------
+  # Inputs
+lambda <- argLs[["lambda"]]
+
+
+
+# apodization -----------------------------------------
+  # Inputs
+phase=0
+rectRatio=1/2
+gaussLB=1
+expLB=1
+apodization <- argLs[["apodizationMethod"]]
+
+if (apodization=='exp'){
+  expLB <- argLs[["expLB"]]
+  } else if (apodization=='cos2'){
+  phase <- argLs[["phase"]]
+  } else if (apodization=='hanning'){
+  phase <- argLs[["phase"]]
+  } else if (apodization=='hamming'){
+  phase <- argLs[["phase"]]
+  } else if (apodization=='blockexp'){
+  rectRatio <- argLs[["rectRatio"]]
+  expLB <- argLs[["expLB"]]
+  } else if (apodization=='blockcos2'){
+  rectRatio <- argLs[["rectRatio"]]
+  } else if (apodization=='gauss'){
+  rectRatio <- argLs[["rectRatio"]]
+  gaussLB <- argLs[["gaussLB"]]
+  }		
+
+
+# Fourier transform ----------------------------------
+  # Inputs
+
+
+# Zero Order Phase Correction -------------------------------
+  # Inputs
+
+angle = NULL
+excludeZOPC = NULL
+
+
+zeroOrderPhaseMethod <- argLs[["zeroOrderPhaseMethod"]]
+										   
+if (zeroOrderPhaseMethod=='manual'){
+  angle <- argLs[["angle"]]
+}
+
+excludeZoneZeroPhase <- argLs[["excludeZoneZeroPhase.choice"]]
+if (excludeZoneZeroPhase == 'YES') {
+  excludeZoneZeroPhaseList <- list()
+  for(i in which(names(argLs)=="excludeZoneZeroPhase_left")) {
+    excludeZoneZeroPhaseLeft <- argLs[[i]]
+    excludeZoneZeroPhaseRight <- argLs[[i+1]]
+    excludeZoneZeroPhaseList <- c(excludeZoneZeroPhaseList,list(c(excludeZoneZeroPhaseLeft,excludeZoneZeroPhaseRight)))
+  }
+  excludeZOPC <- excludeZoneZeroPhaseList
+}
+
+
+# Internal referencering ----------------------------------
+# Inputs
+shiftTreshold = 2 # c
+ppm = TRUE
+shiftReferencingRangeList = NULL  # fromto.RC
+pctNearValue = 0.02 # pc 
+rowindex_graph = NULL
+ppm_ref = 0 # ppm.ref
+
+# 
+# shiftReferencing <- argLs[["shiftReferencing"]]
+# print(shiftReferencing)
+# 
+# if (shiftReferencing=="YES")
+# {
+#   
+# shiftReferencingMethod <- argLs[["shiftReferencingMethod"]]
+# 
+# if (shiftReferencingMethod == "thres")	{
+# 	shiftTreshold <- argLs[["shiftTreshold"]]
+# }
+
+shiftReferencingRange <- argLs[["shiftReferencingRange"]]
+
+if (shiftReferencingRange == "near0"){
+  pctNearValue <- argLs[["pctNearValue"]]
+}
+
+if (shiftReferencingRange == "window"){
+  shiftReferencingRangeList <- list()
+  for(i in which(names(argLs)=="shiftReferencingRangeLeft")) 
+  {
+    shiftReferencingRangeLeft <- argLs[[i]]
+    shiftReferencingRangeRight <- argLs[[i+1]]
+    shiftReferencingRangeList <- c(shiftReferencingRangeList,list(c(shiftReferencingRangeLeft,shiftReferencingRangeRight)))
+  }
+}
+
+shiftHandling <- argLs[["shiftHandling"]]
+
+ppmvalue <- argLs[["ppmvalue"]]
+
+
+
+# }
+
+
+# Baseline Correction -------------------------------
+  # Inputs
+lambdaBc <- argLs[["lambdaBc"]] 
+pBc <- argLs[["pBc"]] 
+epsilon <- argLs[["epsilon"]] 
+
+excludeBC = NULL
+
+excludeZoneBC <- argLs[["excludeZoneBC.choice"]]
+if (excludeZoneBC == 'YES') {
+  excludeZoneBCList <- list()
+  for(i in which(names(argLs)=="excludeZoneBC_left")) {
+    excludeZoneBCLeft <- argLs[[i]]
+    excludeZoneBCRight <- argLs[[i+1]]
+    excludeZoneBCList <- c(excludeZoneBCList,list(c(excludeZoneBCLeft,excludeZoneBCRight)))
+  }
+  excludeBC <- excludeZoneBCList
+}
+
+# transformation of negative values -------------------------------
+  # Inputs
+NegativetoZero <- argLs[["NegativetoZero"]]
+
+
+  # Outputs
+nomGraphe <- argLs[["graphOut"]]
+# dataMatrixOut <- argLs[["dataMatrixOut"]]
+log <- argLs[["logOut"]]
+
+
+
+## Checking arguments
+##-------------------
+error.stock <- "\n"
+
+if(length(error.stock) > 1)
+  stop(error.stock)
+  
+
+##======================================================
+##======================================================  
+## Computation
+##======================================================
+##======================================================
+
+pdf(nomGraphe, onefile = TRUE, width = 13, height = 13)
+
+# FirstOrderPhaseCorrection ---------------------------------
+Fid_data <- GroupDelayCorrection(Fid_data0, Fid_info = samplemetadataFid, group_delay = NULL)
+
+if (FirstOPCGraph == "YES") {
+  title = "FIDs after Group Delay Correction"
+  DrawSignal(Fid_data, subtype = "stacked",
+             ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, 
+             xlab = "Frequency", num.stacked = 4, 
+             main = title, createWindow=FALSE)
+}
+
+# SolventSuppression ---------------------------------
+Fid_data <- SolventSuppression(Fid_data, lambda.ss = lambda, ptw.ss = TRUE, plotSolvent = F, returnSolvent = F)
+	
+if (SSGraph == "YES") {
+  title = "FIDs after Solvent Suppression "
+  DrawSignal(Fid_data, subtype = "stacked",
+             ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, 
+             xlab = "Frequency", num.stacked = 4, 
+             main = title, createWindow=FALSE)
+}
+
+
+# Apodization ---------------------------------	
+Fid_data <- Apodization(Fid_data, Fid_info = samplemetadataFid, DT = NULL, 
+                         type.apod = apodization, phase = phase, rectRatio = rectRatio, gaussLB = gaussLB, expLB = expLB, plotWindow = F, returnFactor = F)
+
+if (ApodGraph == "YES") {
+  title = "FIDs after Apodization"
+  DrawSignal(Fid_data, subtype = "stacked",
+             ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, 
+             xlab = "Frequency", num.stacked = 4, 
+             main = title, createWindow=FALSE)
+}
+
+
+# FourierTransform ---------------------------------
+Spectrum_data <- FourierTransform(Fid_data, Fid_info = samplemetadataFid, reverse.axis = TRUE)
+
+
+if (FTGraph == "YES") {
+  title = "Fourier transformed spectra"
+  DrawSignal(Spectrum_data, subtype = "stacked",
+             ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, 
+             xlab = "Frequency", num.stacked = 4, 
+             main = title, createWindow=FALSE)
+}
+
+
+
+# ZeroOrderPhaseCorrection ---------------------------------
+Spectrum_data  <- ZeroOrderPhaseCorrection(Spectrum_data, type.zopc = zeroOrderPhaseMethod,
+                                           plot_rms = NULL, returnAngle = FALSE,
+                                           createWindow = TRUE,angle = angle,
+                                           plot_spectra = FALSE,
+                                           ppm.zopc = TRUE, exclude.zopc = excludeZOPC)
+
+
+# InternalReferencing ---------------------------------
+# if (shiftReferencing=="YES") {
+Spectrum_data <- InternalReferencing(Spectrum_data, samplemetadataFid, method = "max", range = shiftReferencingRange,
+                                     ppm.value = ppmvalue, shiftHandling = shiftHandling, ppm.ir = TRUE,
+                                     fromto.RC = shiftReferencingRangeList, pc = pctNearValue)
+
+if (SRGraph == "YES") {
+  title = "Spectra after Shift Referencing"
+  DrawSignal(Spectrum_data, subtype = "stacked",
+             ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, 
+             xlab = "Frequency", num.stacked = 4, 
+             main = title, createWindow=FALSE)
+}
+
+# }
+
+if (ZeroOPCGraph == "YES") {
+title = "Spectra after Zero Order Phase Correction"
+DrawSignal(Spectrum_data, subtype = "stacked",
+           ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, 
+           xlab = "Frequency", num.stacked = 4, 
+           main = title, createWindow=FALSE)
+}
+
+
+# BaselineCorrection ---------------------------------									 
+Spectrum_data <- BaselineCorrection(Spectrum_data, ptw.bc = TRUE, lambda.bc = lambdaBc, 
+                                    p.bc = pBc, eps = epsilon, ppm.bc = TRUE, 
+                                    exclude.bc = excludeBC,
+                                    returnBaseline = F) 
+
+
+
+if (BCGraph == "YES") {
+title = "Spectra after Baseline Correction"
+DrawSignal(Spectrum_data, subtype = "stacked",
+           ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, 
+           xlab = "Frequency", num.stacked = 4, 
+           main = title, createWindow=FALSE)
+}
+
+
+# NegativeValuesZeroing ---------------------------------
+if (NegativetoZero=="YES") {
+  Spectrum_data <- NegativeValuesZeroing(Spectrum_data)
+}
+
+if (FinalGraph == "YES") {
+  title = "Final preprocessed spectra"
+  DrawSignal(Spectrum_data, subtype = "stacked",
+             ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, 
+             xlab = "Frequency", num.stacked = 4, 
+             main = title, createWindow=FALSE)
+}
+
+invisible(dev.off())
+
+
+data_variable <- matrix(NA, nrow = 1, ncol = dim(Spectrum_data)[2], dimnames = list("ID", NULL)) 
+colnames(data_variable) <- colnames(Spectrum_data)
+data_variable[1,] <- colnames(data_variable)
+
+
+##======================================================
+##======================================================
+## Saving
+##======================================================
+##======================================================
+
+# Data Matrix
+write.table(t(Re(Spectrum_data)),file=argLs$dataMatrix, quote=FALSE, row.names=TRUE, sep="\t", col.names=TRUE)
+
+# Variable metadata
+write.table(data_variable,file=argLs$variableMetadata, quote=FALSE, row.names=TRUE, sep="\t", col.names=TRUE)
+
+# log file
+# write.table(t(data.frame(argLs)), file = argLs$logOut, col.names = FALSE, quote=FALSE)
+
+# input arguments
+cat("\n INPUT and OUTPUT ARGUMENTS :\n")
+
+argLs
+
+
+## Ending
+
+cat("\nEnd of 'Preprocessing' Galaxy module call: ", as.character(Sys.time()), sep = "")
+
+sink()
+
+options(stringsAsFactors = strAsFacL)
+
+rm(list = ls())