view nmr_alignement/NmrAlignment_wrapper.R @ 4:265ee538e120 draft default tip

Uploaded

#!/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"]])){
	fileType="zip"
	zipfile= argLs[["zipfile"]]
	directory=unzip(zipfile, list=F)
	directory=paste(getwd(),strsplit(directory[1],"/")[[1]][2],sep="/")
} else if (!is.null(argLs[["tsvfile"]])){
	fileType="tsv"
	directory <- read.table(argLs[["tsvfile"]],check.names=FALSE,header=TRUE,sep="\t")
}


		## 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"]]

  ## Graphical zone(s)
graphicalZones <- argLs[["zone_graphical_choices.choice"]]
graphicalZonesBorders <- NULL
if (!is.null(argLs$zone_exclusion_left))
{
  for(i in which(names(argLs)=="zone_graphical_left"))
  {
    graphicalZonesBorders <- c(graphicalZonesBorders,list(c(argLs[[i]],argLs[[i+1]])))
  }
}


	# Outputs
logOut <- argLs[["logOut"]]
alignedSpectra <- argLs[["alignedSpectra"]]
graphOut <- argLs[["graphOut"]]


## Checking R packages
##--------------------
sink(logOut)
cat("\tPACKAGE INFO\n")
#pkgs=c("xcms","batch")
pkgs=c("batch","speaq")
for(pkg in pkgs) {
    suppressPackageStartupMessages( stopifnot( library(pkg, quietly=TRUE, logical.return=TRUE, character.only=TRUE)))
    cat(pkg,"\t",as.character(packageVersion(pkg)),"\n",sep="")
}

## Checking arguments
##-------------------
error.stock <- "\n"
if(length(error.stock) > 1)
  stop(error.stock)


## Computation
##------------
directory.alignement <- nmr.alignment(fileType=fileType,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")


## Graphical output: overlay of raw and estimated spectra
pdf(graphOut,onefile=TRUE)

graphical.zone.length <- length(graphicalZonesBorders)
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="")

if (graphical.zone.length != 0)
{
  #   par(mfrow=c(2,1))
  for (g in 1:graphical.zone.length)
  {
    print(paste(g, graphicalZonesBorders[[g]][1], graphicalZonesBorders[[g]][2]))
    drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == graphicalZonesBorders[[g]][1])[1]):(which(round(as.numeric(colnames(raw.spectra)),2) == graphicalZonesBorders[[g]][2])[1])],xlab="", main="")
    drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == graphicalZonesBorders[[g]][1])[1]):(which(round(as.numeric(colnames(aligned.spectra)),2) == graphicalZonesBorders[[g]][2])[1])],xlab="", main="")
  }
}
dev.off()


## Ending
##---------------------
cat("\nEnd of 'NMR alignment' Galaxy module call: ", as.character(Sys.time()), sep = "")
sink()
options(stringsAsFactors = strAsFacL)
rm(list = ls())