Mercurial > repos > marie-tremblay-metatoul > nmr_bucketing
diff NmrBucketing_wrapper.R @ 9:62c62e31fc80 draft
planemo upload for repository https://github.com/workflow4metabolomics/nmr_bucketing commit fc2be0f9fa66f830d592d74d14b47a63e761647b
author | lecorguille |
---|---|
date | Fri, 21 Apr 2017 08:53:40 -0400 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/NmrBucketing_wrapper.R Fri Apr 21 08:53:40 2017 -0400 @@ -0,0 +1,295 @@ +#!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file + +## 070115_NmrBucketing2galaxy_v1.R +## Marie Tremblay-Franco +## MetaboHUB: The French Infrastructure for Metabolomics and Fluxomics +## www.metabohub.fr/en +## marie.tremblay-franco@toulouse.inra.fr + +runExampleL <- FALSE + +if(runExampleL) { +##------------------------------ +## Example of arguments +##------------------------------ +argLs <- list(StudyDir = "Tlse_BPASourisCerveau", + upper = "10.0", + lower = "0.50", + bucket.width = "0.01", + exclusion = "TRUE", + exclusion.zone = list(c(6.5,4.5)), + graph="Overlay") + +argLs <- c(argLs, + list(dataMatrixOut = paste(directory,"_NmrBucketing_dataMatrix.tsv",sep=""), + sampleMetadataOut = paste(directory,"_NmrBucketing_sampleMetadata.tsv",sep=""), + variableMetadataOut = paste(directory,"_NmrBucketing_variableMetadata.tsv",sep=""), + graphOut = paste(directory,"_NmrBucketing_graph.pdf",sep=""), + logOut = paste(directory,"_NmrBucketing_log.txt",sep=""))) +} + +##------------------------------ +## Options +##------------------------------ +strAsFacL <- options()$stringsAsFactors +options(stringsAsFactors = FALSE) + + +##------------------------------ +## Libraries laoding +##------------------------------ +# For parseCommandArgs function +library(batch) +# For cumtrapz function +library(pracma) + +# 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("NmrBucketing_script.R") +source_local("DrawSpec.R") + +##------------------------------ +## Errors ????????????????????? +##------------------------------ + + +##------------------------------ +## Constants +##------------------------------ +topEnvC <- environment() +flagC <- "\n" + + +##------------------------------ +## Script +##------------------------------ +if(!runExampleL) + argLs <- parseCommandArgs(evaluate=FALSE) + + +## Parameters Loading +##------------------- + # Inputs +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") +} + +leftBorder <- argLs[["left_border"]] +rightBorder <- argLs[["right_border"]] +bucketSize <- argLs[["bucket_width"]] +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]]))) + } +} + +graphique <- argLs[["graphType"]] + + # Outputs +nomGraphe <- argLs[["graphOut"]] +dataMatrixOut <- argLs[["dataMatrixOut"]] +logFile <- argLs[["logOut"]] +if (fileType=="zip") +{ + sampleMetadataOut <- argLs[["sampleOut"]] + variableMetadataOut <- argLs[["variableOut"]] +} + +## Checking arguments +##------------------- +error.stock <- "\n" + +if(length(error.stock) > 1) + stop(error.stock) + + +## Computation +##------------ +outputs <- NmrBucketing(fileType=fileType, fileName=directory, leftBorder=leftBorder, rightBorder=rightBorder, bucketSize=bucketSize, + exclusionZones=exclusionZones, exclusionZonesBorders=exclusionZonesBorders, graph=graphique, nomFichier=nomGraphe, + savLog.txtC=logFile) +data_bucket <- outputs[[1]] +data_sample <- outputs[[2]] +data_variable <- outputs[[3]] +ppm <- outputs[[4]] +ppm <- round(ppm,2) + +## Graphical outputs +##------------------ +if (graphique != "None") +{ + excludedZone <- NULL + for (c in 1:length(exclusionZonesBorders)) + { + excludedZone <- c(excludedZone,exclusionZonesBorders[[c]]) + excludedZone <- sort(excludedZone) + } + nbZones <- length(excludedZone)/2 + n <- length(excludedZone) + + # Graphic Device opening + pdf(nomGraphe,onefile=TRUE) + + if (graphique == "Overlay") + { + # Global spectral window + spectra <- data.frame(t(data_bucket)) + drawSpec(spectra,xlab="", ylab="Intensity", main="") + + ## Zoomed spectral window depending on exclusion zone(s) + if (nbZones != 0) + { + BInf <- excludedZone[n] + if (round(BInf,1) == BInf) + { + BInf <- BInf+0.01 + } + spectra <- data.frame(t(data_bucket[1:(which(ppm == BInf)[[1]]),])) + drawSpec(spectra,xlab="", ylab="Intensity", main="") + n <- n - 1 + + while (n >= nbZones & nbZones > 1) + { + BInf <- excludedZone[n-1] + if (round(BInf,1) > BInf) + { + BInf <- BInf+0.01 + } + spectra <- data.frame(t(data_bucket[(which(ppm == excludedZone[n])[[1]]):(which(ppm == BInf)[[1]]),])) + drawSpec(spectra,xlab="", ylab="Intensity", main="") + n <- n - 2 + } + + BInf <- excludedZone[1] + if (round(BInf,1) <= BInf) + { + BInf <- BInf+0.01 + } + spectra <- data.frame(t(data_bucket[(which(ppm == BInf)[[1]]):nrow(data_bucket),])) + drawSpec(spectra,xlab="", ylab="Intensity", main="") + } + } + else + { + for (i in 1:ncol(data_bucket)) + { + par(mfrow=c((nbZones+2),1)) + n <- length(excludedZone) + spectra <- t(data_bucket[,i]) + names(spectra) <- rownames(data_bucket) + plot(1:length(spectra), spectra, type='l', xlab="", ylab="Intensity", main=colnames(data_bucket)[i], xaxt = "n") + xPos <- 1 + nAxisPos <- 4 + startP <- length(nAxisPos) + endP <- nrow(data_bucket) + GraphRange <- c(startP:endP) + tempVal = trunc(length(GraphRange)/nAxisPos) + xPos = c(0:nAxisPos) * tempVal + axis(1, at = xPos, labels = rownames(data_bucket)[xPos + startP]) + + ## Zoomed spectral window depending on exclusion zone(s) + if (nbZones != 0) + { + BInf <- excludedZone[n] + if (round(BInf,1) == BInf) + { + BInf <- BInf+0.01 + } + spectra <- t(data_bucket[1:(which(ppm == BInf)[[1]]),i]) + names(spectra) <- rownames(data_bucket)[1:(which(ppm == BInf)[[1]])] + plot(1:length(spectra), spectra, type='l',xlab="", ylab="Intensity", main="", xaxt = "n") + xPos <- 1 + nAxisPos <- 4 + startP <- length(nAxisPos) + endP <- length(spectra) + GraphRange <- c(startP:endP) + tempVal = trunc(length(GraphRange)/nAxisPos) + xPos = c(0:nAxisPos) * tempVal + axis(1, at = xPos, labels = rownames(data_bucket)[xPos + startP]) + n <- n - 1 + + while (n >= nbZones & nbZones > 1) + { + BInf <- excludedZone[n-1] + if (round(BInf,1) > BInf) + { + BInf <- BInf+0.01 + } + spectra <- t(data_bucket[(which(ppm == excludedZone[n])[[1]]):(which(ppm == BInf)[[1]]),i]) + names(spectra) <- rownames(data_bucket)[(which(ppm == excludedZone[n])[[1]]):(which(ppm == BInf)[[1]])] + plot(1:length(spectra), spectra, type='l',xlab="", ylab="Intensity", main="", xaxt = "n") + xPos <- 1 + nAxisPos <- 4 + startP <- length(nAxisPos) + endP <- length(spectra) + GraphRange <- c(startP:endP) + tempVal = trunc(length(GraphRange)/nAxisPos) + xPos = c(0:nAxisPos) * tempVal + axis(1, at = xPos, labels = rownames(data_bucket)[xPos + startP]) + n <- n - 2 + } + + BInf <- excludedZone[1] + if (round(BInf,1) <= BInf) + { + BInf <- BInf+0.01 + } + spectra <- t(data_bucket[(which(ppm == BInf)[[1]]):nrow(data_bucket),i]) + names(spectra) <- rownames(data_bucket)[(which(ppm == BInf)[[1]]):nrow(data_bucket)] + plot(1:length(spectra), spectra, type='l',xlab="", ylab="Intensity", main="", xaxt = "n") + xPos <- 1 + nAxisPos <- 4 + startP <- length(nAxisPos) + endP <- length(spectra) + GraphRange <- c(startP:endP) + tempVal = trunc(length(GraphRange)/nAxisPos) + xPos = c(0:nAxisPos) * tempVal + axis(1, at = xPos, labels = rownames(data_bucket)[xPos + startP]) + } + } + } + dev.off() +} +## Saving +##------- + # Data +data_bucket <- cbind(rownames(data_bucket),data_bucket) +colnames(data_bucket) <- c("Bucket",colnames(data_bucket)[-1]) +write.table(data_bucket,file=argLs$dataMatrixOut,quote=FALSE,row.names=FALSE,sep="\t") + # Sample +data_sample <- cbind(rownames(data_sample),data_sample) +colnames(data_sample) <- c("Sample",colnames(data_sample)[-1]) +write.table(data_sample,file=argLs$sampleOut,quote=FALSE,row.names=FALSE,sep="\t") + # Variable +data_variable <- cbind(rownames(data_variable),data_variable) +colnames(data_variable) <- c("Bucket",colnames(data_variable)[-1]) +write.table(data_variable,file=argLs$variableOut,quote=FALSE,row.names=FALSE,sep="\t") + + +## Ending +##--------------------- + +cat("\nEnd of 'NMR bucketing' Galaxy module call: ", as.character(Sys.time()), sep = "") + +## sink(NULL) + +options(stringsAsFactors = strAsFacL) + +rm(list = ls())