diff nmr_alignement/NmrAlignment_wrapper.R @ 1:58eecef626da draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/nmr_alignement/NmrAlignment_wrapper.R	Thu Mar 02 10:23:37 2017 -0500
@@ -0,0 +1,176 @@
+#!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file
+
+## 06102016_NmrAlignment_wrapper.R
+## Marie Tremblay-Franco
+## [e-mail]
+
+runExampleL <- FALSE
+
+##------------------------------
+## Options
+##------------------------------
+strAsFacL <- options()$stringsAsFactors
+options(stringsAsFactors = FALSE)
+
+
+##------------------------------
+## Libraries loading
+##------------------------------
+	# ParseCommandArgs function
+library(batch) 
+	# Alignment
+library(speaq)
+
+
+# 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="/"))
+}
+# Function import
+source_local("NmrAlignment_script.R")
+
+
+##------------------------------
+## Errors ?????????????????????
+##------------------------------
+
+
+##------------------------------
+## Constants
+##------------------------------
+topEnvC <- environment()
+flagC <- "\n"
+
+
+##------------------------------
+## Script
+##------------------------------
+if(!runExampleL)
+    argLs <- parseCommandArgs(evaluate=FALSE)
+
+	
+## Parameters Loading
+##-------------------
+	# Inputs
+		## Library of spectra to align
+if (!is.null(argLs[["zipfile"]])){
+	zipfile= argLs[["zipfile"]]
+	directory=unzip(zipfile, list=F)
+	directory=paste(getwd(),strsplit(directory[1],"/")[[1]][2],sep="/")
+} else if (!is.null(argLs[["library"]])){
+	directory=argLs[["library"]]
+    	if(!file.exists(directory)){
+		error_message=paste("Cannot access the directory :",directory,".Please verify if the directory exists or not.")
+		print(error_message)
+		stop(error_message)
+	}
+}
+
+
+		## Spectral width
+leftBorder <- argLs[["left_border"]]
+rightBorder <- argLs[["right_border"]]
+		
+		##Exclusion zone(s)
+exclusionZones <- argLs[["zone_exclusion_choices.choice"]]
+exclusionZonesBorders <- NULL
+if (!is.null(argLs$zone_exclusion_left))
+{
+   for(i in which(names(argLs)=="zone_exclusion_left"))
+   {
+     exclusionZonesBorders <- c(exclusionZonesBorders,list(c(argLs[[i]],argLs[[i+1]])))
+   }
+}
+		
+		## Reference spectrum
+reference <- argLs[["reference"]]
+
+		## Size of a small nDivRange
+nDivRange <- argLs[["nDivRange"]]
+		
+		## Intensity threshold for peak removal
+baselineThresh <- argLs[["baselineThresh"]]
+
+
+	# Outputs
+logOut <- argLs[["logOut"]]
+alignedSpectra <- argLs[["alignedSpectra"]]
+graphOut <- argLs[["graphOut"]]
+
+
+## Checking arguments
+##-------------------
+error.stock <- "\n"
+if(length(error.stock) > 1)
+  stop(error.stock)
+  
+  
+## Computation
+##------------
+directory.alignement <- nmr.alignment(directory=directory,leftBorder=leftBorder,rightBorder=rightBorder,exclusionZones=exclusionZones,
+                                  exclusionZonesBorders=exclusionZonesBorders, reference=reference, nDivRange=nDivRange, 
+                                  baselineThresh=baselineThresh, maxshift=50, verbose=FALSE)
+directory.raw <- directory.alignement[[1]]
+directory.aligned <- directory.alignement[[2]]
+
+## Saving
+##-------
+	# Aligned spectra
+t.directory.aligned <- t(directory.aligned)
+rownames(t.directory.aligned) <- colnames(directory.aligned)
+# colnames(t.directory.aligned) <- c("Bucket",colnames(t.directory.aligned))
+write.table(t.directory.aligned,file=alignedSpectra,row.names=TRUE,quote=FALSE,sep="\t") 
+
+
+excludedZone <- NULL
+for (c in 1:length(exclusionZonesBorders))
+{
+  excludedZone <- c(excludedZone,exclusionZonesBorders[[c]])
+  excludedZone <- sort(excludedZone)
+}
+
+## Graphical output: overlay of raw and estimated spectra
+pdf(graphOut,onefile=TRUE)
+par(mfrow=c(2,1))
+
+raw.spectra <- data.frame(directory.raw)
+colnames(raw.spectra) <- substr(colnames(raw.spectra),2,7)
+
+aligned.spectra <- data.frame(directory.aligned)
+colnames(aligned.spectra) <- substr(colnames(aligned.spectra),2,7)
+
+drawSpec(raw.spectra,xlab="", ylab="Raw spectra", main="")
+drawSpec(aligned.spectra,xlab="", ylab="Aligned spectra", main="")
+
+nbZones <- length(excludedZone)/2
+if (nbZones != 0)
+{
+  n <- length(excludedZone)
+  drawSpec(raw.spectra[,1:which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[n])[1]],xlab="", ylab="Raw spectra", main="")
+  drawSpec(aligned.spectra[,1:which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[n])[1]],xlab="", ylab="Aligned spectra", main="")
+
+  n <- n - 1
+  while (n >= nbZones & nbZones > 1)
+  {
+    drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[n])[1]):(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[n-1])[1])],xlab="", ylab="Raw spectra", main="")
+    drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[n])[1]):(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[n-1])[1])],xlab="", ylab="Aligned spectra", main="")
+    n <- n - 2
+  }
+  
+  drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[1])[1]):ncol(raw.spectra)],xlab="", ylab="Raw spectra", main="")
+  drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[1])[1]):ncol(aligned.spectra)],xlab="", ylab="Aligned spectra", main="")
+}
+drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == 2.4)[1]):(which(round(as.numeric(colnames(raw.spectra)),2) == 2.8)[1])],xlab="", ylab="Raw spectra", main="")
+drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == 2.4)[1]):(which(round(as.numeric(colnames(aligned.spectra)),2) == 2.8)[1])],xlab="", ylab="Aligned spectra", main="")
+
+dev.off()
+
+
+## Ending
+##---------------------
+cat("\nEnd of 'NMR alignment' Galaxy module call: ", as.character(Sys.time()), sep = "")
+options(stringsAsFactors = strAsFacL)
+rm(list = ls())