Mercurial > repos > marie-tremblay-metatoul > nmr_alignment
diff nmr_alignement/NmrAlignment_wrapper.R @ 1:58eecef626da draft
Uploaded
author | marie-tremblay-metatoul |
---|---|
date | Thu, 02 Mar 2017 10:23:37 -0500 |
parents | |
children | 908e1345d7ca |
line wrap: on
line diff
--- /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 +## marie.tremblay-franco@toulouse.inra.fr + +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())