| Next changeset 1:cbea5e9fd0b4 (2017-05-23) |
|
Commit message:
Uploaded |
|
added:
nmr_preprocessing/.Rhistory nmr_preprocessing/Apodization_wrapper_Manon.R nmr_preprocessing/DrawFunctions.R nmr_preprocessing/Int_Funct.R nmr_preprocessing/NmrPreprocessing_script.R nmr_preprocessing/NmrPreprocessing_wrapper.R nmr_preprocessing/NmrPreprocessing_xml.xml nmr_preprocessing/ReadFids.R nmr_preprocessing/ReadFids_Manon.R nmr_preprocessing/ReadFids_wrapper.R nmr_preprocessing/ReadFids_xml.xml nmr_preprocessing/ptw.R nmr_preprocessing/static/images/NmrPreprocessing.png nmr_preprocessing/test-data/MTBLS1.zip nmr_preprocessing/test-data/MTBLS1_alignedSpectra.tabular |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/Apodization_wrapper_Manon.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nmr_preprocessing/Apodization_wrapper_Manon.R Fri May 05 07:13:02 2017 -0400 |
| [ |
| @@ -0,0 +1,187 @@ +#!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file + +## 12012017_Apodization_wrapper.R +## Manon Martin +## manon.martin@uclouvain.be + + +#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# Use of runExampleL ? +#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +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 +#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# stringsAsFactors = FALSE utilis?? o?? ?? +#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +options(stringsAsFactors = FALSE) + + + + +##------------------------------ +## Libraries loading +##------------------------------ +# For parseCommandArgs function +library(batch) + +# 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("Apodization_Manon.R") + +##------------------------------ +## Errors ????????????????????? +##------------------------------ + + +##------------------------------ +## Constants +##------------------------------ +#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# interviennent ou ? +#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +topEnvC <- environment() +flagC <- "\n" + + +##------------------------------ +## Script +##------------------------------ +if(!runExampleL) + argLs <- parseCommandArgs(evaluate=FALSE) # If FALSE valeurs par d??faut d??finies dans le software + + +## Parameters Loading +##------------------- +## Inputs +# data + +# --- Galaxy_preformatted data +# data <- read.table(argLs[["dataMatrix"]],check.names=FALSE,header=TRUE,sep="\t") +# rownames(data) <- data[,1] +# data <- data[,-1] +# --- + +# --- data from ReadFids +Fid_data <- read.table(argLs[["dataMatrixOut"]],check.names=FALSE,header=TRUE,sep="\t") +Fid_info <- read.table(argLs[["sampleOut"]],check.names=FALSE,header=TRUE,sep="\t") +# --- + +# other inputs (cf. XML wrapper) +DT = argLs[["DT"]] +type.apod = argLs[["ApodizationMethod"]] + +# set default values for optional arguments +phase=0 +rectRatio=1/2 +gaussLB=1 +expLB=1 + +# change the default values +if (type.apod %in% c("cos2", "hanning", "hamming")) { + phase = argLs[["phase"]] +} + +if (type.apod == "blockexp") { + rectRatio = argLs[["rectRatio"]] + expLB = argLs[["expLB"]] +} + +if (type.apod == "blockcos2") { + rectRatio = argLs[["rectRatio"]] +} + +if (type.apod == "gauss") { +gaussLB = argLs[["gaussLB"]] +} + +if (type.apod == "exp") { + expLB = argLs[["expLB"]] +} + +plotWindow = FALSE +returnFactor = FALSE + + + +#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# Utility of Outputs ?? +#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# Outputs +dataMatrixOut <- argLs[["dataMatrixOut"]] # Names from Saving + + + +## Checking arguments +##------------------- +error.stock <- "\n" + +#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# error.stock utilis?? ou ? +#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +if(length(error.stock) > 1) + stop(error.stock) + + +## Computation +##------------ +outputs <- Apodization(Fid_data = Fid_data, Fid_info = Fid_info, DT = DT, + type.apod = type.apod, phase = phase, rectRatio = rectRatio, + gaussLB = gaussLB, expLB = expLB, plotWindow = plotWindow, returnFactor = returnFactor) + +data_matrix <- outputs # Data matrix + + + +## Saving +##------- +# Data matrix +data_matrix <- cbind(rownames(data_matrix),data_matrix) +colnames(data_matrix) <- c("Sample",colnames(data_matrix)[-1]) +write.table(data_matrix,file=argLs$dataMatrixOut,quote=FALSE,row.names=FALSE,sep="\t") + + + +## Ending +##--------------------- + +cat("\nEnd of 'Apodization' Galaxy module call: ", as.character(Sys.time()), sep = "") + +## sink(NULL) + +options(stringsAsFactors = strAsFacL) + +rm(list = ls()) \ No newline at end of file |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/DrawFunctions.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nmr_preprocessing/DrawFunctions.R Fri May 05 07:13:02 2017 -0400 |
| [ |
| b'@@ -0,0 +1,257 @@\n+require(ggplot2)\n+require(gridExtra)\n+require(reshape2)\n+\n+\n+Draw <- function(Signal_data, type.draw = c("signal", "pca"), output = c("default", \n+ "window", "png", "pdf"), dirpath = ".", filename = "%003d", height = 480, \n+ width = 640, pdf.onefile = TRUE, ...) {\n+ \n+ # Data initialisation and checks ----------------------------------------------\n+ type.draw <- match.arg(type.draw)\n+ output <- match.arg(output)\n+ fullpath <- paste(file.path(dirpath, filename), output, sep = ".")\n+ createFile <- TRUE\n+ createWindow <- FALSE\n+ \n+ # Drawing --------------------------------------------------------------------\n+ # output\n+ switch(output, default = {\n+ createFile <- FALSE\n+ }, window = {\n+ createWindow <- TRUE\n+ createFile <- FALSE\n+ }, png = {\n+ grDevices::png(fullpath, width, height)\n+ }, pdf = {\n+ grDevices::pdf(fullpath, width = width/72, height = height/72, \n+ onefile = pdf.onefile)\n+ }, {\n+ stop("Unknown output type.")\n+ })\n+ \n+ # Drawing type (signal/spectrum or PCA)\n+ funs <- list(signal = DrawSignal, pca = DrawPCA)\n+ if (type.draw %in% names(funs)) {\n+ fun <- funs[[type.draw]]\n+ } else {\n+ stop(paste("Unknown type:", type.draw))\n+ }\n+ \n+ # Plot finalisation ----------------------------------------------\n+ if (is.vector(Signal_data)) {\n+ Signal_data <- vec2mat(Signal_data)\n+ }\n+ fun(Signal_data, createWindow = createWindow, ...)\n+ if (createFile) {\n+ grDevices::dev.off()\n+ }\n+}\n+\n+\n+\n+##### DrawSignal\n+\n+DrawSignal <- function(Signal_data, subtype = c("stacked", "together", \n+ "separate", "diffmean", "diffmedian", "diffwith"), \n+ ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, \n+ xlab = "rowname", RowNames = NULL, row = 1, num.stacked = 4, \n+ main = NULL, createWindow) {\n+ # nticks\n+ \n+ # Data initialisation and checks ----------------------------------------------\n+ \n+ subtype <- match.arg(subtype)\n+ \n+ \n+ n <- nrow(Signal_data)\n+ m <- ncol(Signal_data)\n+ \n+ \n+ scale <- colnames(Signal_data)\n+ \n+ num.plot <- sum(ReImModArg)\n+ \n+ Var <- rowname <- value <- NULL # only for R CMD check\n+ \n+ # Drawing array\n+ if (num.plot <= 0) {\n+ stop("Nothing selected in ReImModArg.")\n+ } else if (num.plot <= 2) {\n+ if (vertical) {\n+ nrow <- num.plot\n+ ncol <- 1\n+ } else {\n+ nrow <- 1\n+ ncol <- num.plot\n+ }\n+ } else {\n+ nrow <- 2\n+ ncol <- 2\n+ }\n+ \n+ # RowNames \n+ if (is.null(RowNames)) {\n+ RowNames <- rownames(Signal_data)\n+ if (is.null(RowNames)) {\n+ RowNames <- 1:n\n+ }\n+ } else {\n+ if (!is.vector(RowNames)) {\n+ stop("RowNames is not a vector")\n+ }\n+ if (length(RowNames) != n) {\n+ stop(paste("RowNames has length", length(RowNames), "and there are", n, "FIDs."))\n+ }\n+ }\n+ \n+ if (n == 1) {\n+ RowNames <- deparse(substitute(Signal_data))\n+ }\n+ \n+ elements <- list()\n+ if (ReImModArg[1]) {\n+ elements[["Re"]] <- Re(Signal_data)\n+ }\n+ if (ReImModArg[2]) {\n+ elements[["Im"]] <- Im(Signal_data)\n+ }\n+ if (ReImModArg[3]) {\n+ elements[["Mod"]] <- Mod(Signal_data)\n+ }\n+ if (ReImModArg[4]) {\n+ elements[["Arg"]] <- Arg(Signal_data)\n+ }\n+ \n+ \n+ \n+ \n+ # Drawing --------------------------------------------------------------------\n+ \n+ y = x = NULL # only for R CMD check\n+ \n+ \n+ # SEPARATE or STACKED ===============\n+ if (subtype == "separate" | subtype == "stacked") {\n+ \n+ i <- 1\n+ while (i <= n) {\n+ if (createWindow) {\n+ grDevices::dev.new(noRStudioGD = TRUE)\n+ }\n+ if (subtype == "separate") {\n+ # The other uses gridExtra to do that\n+ graphics::par(mfrow = c(nrow, ncol))\n+ }\n+ plots <- list()\n+ if (subtype == "separate") {\n+ last <- i\n+ '..b' } else {df <- data.frame(x = as.numeric(scale), y = elements[[name]][i, ])\n+ }\n+ \n+ plots[[name]] <- ggplot2::ggplot(data = df, ggplot2::aes(x = x, y = y)) + \n+ ggplot2::geom_line() + \n+ ggplot2::theme(legend.position = "none") + \n+ ggplot2::labs(x = xlab, y = name) +\n+ ggplot2::ggtitle(RowNames[i]) +\n+ ggplot2::theme_bw()\n+ \n+ if ((df[1, "x"] - df[(dim(df)[1]), "x"]) > 0) {\n+ plots[[name]] <- plots[[name]] + ggplot2::scale_x_reverse()\n+ }\n+ \n+ } else {\n+ \n+ if (n == 1 ) {\n+ melted <- data.frame(rowname = rep(name, m), \n+ Var = as.numeric(scale), value = elements[[name]][i,])\n+ } else if (last==i){ \n+ melted <- data.frame(rowname = rep(rownames(elements[[name]])[i], m), \n+ Var = as.numeric(scale), value = elements[[name]][i,])\n+ } else {melted <- reshape2::melt(elements[[name]][i:last, ], \n+ varnames = c("rowname", "Var"))\n+ }\n+ \n+ \n+ plots[[name]] <- ggplot2::ggplot(data = melted, ggplot2::aes(x = Var, y = value)) + \n+ ggplot2::geom_line() + \n+ ggplot2::facet_grid(rowname ~ ., scales = "free_y") + \n+ ggplot2::theme(legend.position = "none") + \n+ ggplot2::labs(x = xlab, y = name) +\n+ ggplot2::ggtitle(label = main) +\n+ ggplot2::theme_bw()\n+ \n+ if ((melted[1, "Var"] - melted[(dim(melted)[1]), "Var"]) > 0) {\n+ plots[[name]] <- plots[[name]] + ggplot2::scale_x_reverse()\n+ }\n+ }\n+ }\n+ \n+ if (subtype == "stacked") {\n+ do.call(gridExtra::grid.arrange, c(plots, list(nrow = nrow, ncol = ncol)))\n+ }\n+ i <- last + 1\n+ }\n+ } else if (subtype %in% c("together", "diffmean", "diffmedian", "diffwith")) {\n+ \n+ # TOGHETER or DIFFMEAN or DIFFMEDIAN or DIFFWITH ===============\n+ \n+ rainbow_colors <- grDevices::rainbow(n)\n+ \n+ if (createWindow) {\n+ grDevices::dev.new(noRStudioGD = TRUE)\n+ }\n+ graphics::par(mfrow = c(nrow, ncol))\n+ \n+ plots <- list()\n+ \n+ # Loop for Re, Im, Mod and Arg\n+ for (name in names(elements)) {\n+ # Get this part of the signal\n+ element <- elements[[name]]\n+ \n+ # Express the signal according to a reference if asked by `subtype\'\n+ if (subtype == "diffmean") {\n+ element <- sweep(element, MARGIN = 2, colMeans(element), `-`)\n+ } else if (subtype == "diffmedian") {\n+ element <- sweep(element, MARGIN = 2, matrixStats::colMedians(element), `-`)\n+ } else if (subtype == "diffwith") {\n+ element <- sweep(element, MARGIN = 2, element[row, ], `-`)\n+ if (row == 1 & n > 1) {\n+ # Since we use plot on the first row and lines on the following, the y\n+ # scale is calculated at the first row so if the first row is all 0, it\n+ # causes problems\n+ tmp <- element[1, ]\n+ element[1, ] <- element[2, ]\n+ element[2, ] <- tmp\n+ }\n+ }\n+ \n+ \n+ melted <- reshape2::melt(elements[[name]], varnames = c("rowname", "Var"))\n+ \n+ \n+ plots[[name]] <- ggplot2::ggplot(melted, ggplot2::aes(x = Var, \n+ y = value, group = rowname, colour = rowname)) + ggplot2::geom_line() + \n+ ggplot2::labs(x = xlab, y = name) + ggplot2::scale_colour_discrete(name = NULL) + \n+ ggplot2::ggtitle(main)\n+ \n+ if ((melted[1, "Var"] - melted[(dim(melted)[1]), "Var"]) > \n+ 0) {\n+ plots[[name]] <- plots[[name]] + ggplot2::scale_x_reverse()\n+ }\n+ \n+ do.call(gridExtra::grid.arrange, c(plots, list(nrow = nrow, \n+ ncol = ncol)))\n+ }\n+ }\n+ \n+ \n+}\n' |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/Int_Funct.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nmr_preprocessing/Int_Funct.R Fri May 05 07:13:02 2017 -0400 |
| [ |
| @@ -0,0 +1,236 @@ +beginTreatment <- function(name, Signal_data = NULL, Signal_info = NULL, + force.real = FALSE) { + + cat("Begin", name, "\n") + + + # Formatting the Signal_data and Signal_info ----------------------- + + vec <- is.vector(Signal_data) + if (vec) { + Signal_data <- vec2mat(Signal_data) + } + if (is.vector(Signal_info)) { + Signal_info <- vec2mat(Signal_info) + } + if (!is.null(Signal_data)) { + if (!is.matrix(Signal_data)) { + stop("Signal_data is not a matrix.") + } + if (!is.complex(Signal_data) && !is.numeric(Signal_data)) { + stop("Signal_data contains non-numerical values.") + } + } + if (!is.null(Signal_info) && !is.matrix(Signal_info)) { + stop("Signal_info is not a matrix.") + } + + + Original_data <- Signal_data + + # Extract the real part of the spectrum --------------------------- + + if (force.real) { + if (is.complex(Signal_data)) { + Signal_data <- Re(Signal_data) + } else { + # The signal is numeric Im(Signal_data) is zero anyway so let's avoid + # using complex(real=...,imaginary=0) which would give a complex signal + # in endTreatment() + force.real <- FALSE + } + } + + + # Return the formatted data and metadata entries -------------------- + + return(list(start = proc.time(), vec = vec, force.real = force.real, + Original_data = Original_data, Signal_data = Signal_data, Signal_info = Signal_info)) +} + +binarySearch <- function(a, target, lower = TRUE) { + # search the index i in a such that a[i] == target + # if it doesn't exists and lower, it searches the closer a[i] such that a[i] < target + # if !lower, it seraches the closer a[i] such that a[i] > target + # a should be monotone but can be increasing or decreasing + + # if a is increasing INVARIANT: a[amin] < target < a[amax] + N <- length(a) + if ((a[N] - target) * (a[N] - a[1]) <= 0) { + return(N) + } + if ((a[1] - target) * (a[N] - a[1]) >= 0) { + return(1) + } + amin <- 1 + amax <- N + while (amin + 1 < amax) { + amid <- floor((amin + amax)/2) + if ((a[amid] - target) * (a[amax] - a[amid]) < 0) { + amin <- amid + } else if ((a[amid] - target) * (a[amax] - a[amid]) > 0) { + amax <- amid + } else { + # a[amid] == a[amax] or a[amid] == target In both cases, a[amid] == + # target + return(amid) + } + } + if (xor(lower, a[amin] > a[amax])) { + # (lower && a[amin] < a[amax]) || (!lower && a[min] > a[max]) + # If increasing and we want the lower, we take amin + # If decreasing and we want the bigger, we take amin too + return(amin) + } else { + return(amax) + } +} + +checkArg <- function(arg, checks, can.be.null=FALSE) { + check.list <- list(bool=c(is.logical, "a boolean"), + int =c(function(x){x%%1==0}, "an integer"), + num =c(is.numeric, "a numeric"), + str =c(is.character, "a string"), + pos =c(function(x){x>0}, "positive"), + pos0=c(function(x){x>=0}, "positive or zero"), + l1 =c(function(x){length(x)==1}, "of length 1") + ) + if (is.null(arg)) { + if (!can.be.null) { + stop(deparse(substitute(arg)), " is null.") + } + } else { + if (is.matrix(arg)) { + stop(deparse(substitute(arg)), " is not scalar.") + } + for (c in checks) { + if (!check.list[[c]][[1]](arg)) { + stop(deparse(substitute(arg)), " is not ", check.list[[c]][[2]], ".") + } + } + } +} + +endTreatment <- function(name, begin_info, Signal_data) { + + # begin_info: object outputted from beginTreatment + + + # Formatting the entries and printing process time ----------------------- + end_time <- proc.time() # record it as soon as possible + start_time <- begin_info[["start"]] + delta_time <- end_time - start_time + delta <- delta_time[] + cat("End", name, "\n") + cat("It lasted", round(delta["user.self"], 3), "s user time,", round(delta["sys.self"],3), + "s system time and", round(delta["elapsed"], 3), "s elapsed time.\n") + + + if (begin_info[["force.real"]]) { + # The imaginary part is left untouched + i <- complex(real = 0, imaginary = 1) + Signal_data <- Signal_data + i * Im(begin_info[["Original_data"]]) + } + + if (begin_info[["vec"]]) { + Signal_data <- Signal_data[1, ] + } + + # Return the formatted data and metadata entries -------------------- + return(Signal_data) +} + +Interpol <- function(t, y) { + # y: sample + # t : warping function + + m <- length(y) + # t <= m-1 + # because if t > m-1, y[ti+1] will be NA when we compute g + valid <- 1 <= t & t <= m-1 # FIXME it was '<' in Bubble v2 + s <- (1:m)[valid] + ti <- floor(t[s]) + tr <- t[s] - ti + g <- y[ti + 1] - y[ti] + f <- y[ti] + tr * g + list(f=f, s=s, g=g) +} + +vec2mat <- function(vec) { + + return(matrix(vec, nrow = 1, dimnames = list(c(1), names(vec)))) + +} + +getArg <- function(arg, info, argname, can.be.absent=FALSE) { + if (is.null(arg)) { + start <- paste("impossible to get argument", argname, "it was not given directly and"); + if (!is.matrix(info)) { + stop(paste(start, "the info matrix was not given")) + } + if (!(argname %in% colnames(info))) { + if (can.be.absent) { + return(NULL) + } else { + stop(paste(start, "is not in the info matrix")) + } + } + if (nrow(info) < 1) { + stop(paste(start, "the info matrix has no row")) + } + arg <- info[1,argname] + if (is.na(arg)) { + stop(paste(start, "it is NA in the info matrix")) + } + } + return(arg) +} + +binarySearch <- function(a, target, lower = TRUE) { + # search the index i in a such that a[i] == target + # if it doesn't exists and lower, it searches the closer a[i] such that a[i] < target + # if !lower, it seraches the closer a[i] such that a[i] > target + # a should be monotone but can be increasing or decreasing + + # if a is increasing INVARIANT: a[amin] < target < a[amax] + N <- length(a) + if ((a[N] - target) * (a[N] - a[1]) <= 0) { + return(N) + } + if ((a[1] - target) * (a[N] - a[1]) >= 0) { + return(1) + } + amin <- 1 + amax <- N + while (amin + 1 < amax) { + amid <- floor((amin + amax)/2) + if ((a[amid] - target) * (a[amax] - a[amid]) < 0) { + amin <- amid + } else if ((a[amid] - target) * (a[amax] - a[amid]) > 0) { + amax <- amid + } else { + # a[amid] == a[amax] or a[amid] == target In both cases, a[amid] == + # target + return(amid) + } + } + if (xor(lower, a[amin] > a[amax])) { + # (lower && a[amin] < a[amax]) || (!lower && a[min] > a[max]) + # If increasing and we want the lower, we take amin + # If decreasing and we want the bigger, we take amin too + return(amin) + } else { + return(amax) + } +} + +# returns the discrete borders of the interval for a numeric vector a + +indexInterval <- function (a, from, to, inclusive=TRUE) { + # If inclusive and from <= to, we need to take the lower + # If not inclusive and from > to, we need to take the lower too + lowerFrom <- xor(inclusive, from > to) + fromIndex <- binarySearch(a, from, lowerFrom) + toIndex <- binarySearch(a, to, !lowerFrom) + return(fromIndex:toIndex) +} |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/NmrPreprocessing_script.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nmr_preprocessing/NmrPreprocessing_script.R Fri May 05 07:13:02 2017 -0400 |
| [ |
| b'@@ -0,0 +1,1194 @@\n+## ==========================\n+# Internal functions\n+## ==========================\n+\n+# beginTreatment \n+beginTreatment <- function(name, Signal_data = NULL, Signal_info = NULL, \n+ force.real = FALSE) {\n+ \n+ cat("Begin", name, "\\n")\n+ \n+ \n+ # Formatting the Signal_data and Signal_info -----------------------\n+ \n+ vec <- is.vector(Signal_data)\n+ if (vec) {\n+ Signal_data <- vec2mat(Signal_data)\n+ }\n+ if (is.vector(Signal_info)) {\n+ Signal_info <- vec2mat(Signal_info)\n+ }\n+ if (!is.null(Signal_data)) {\n+ if (!is.matrix(Signal_data)) {\n+ stop("Signal_data is not a matrix.")\n+ }\n+ if (!is.complex(Signal_data) && !is.numeric(Signal_data)) {\n+ stop("Signal_data contains non-numerical values.")\n+ }\n+ }\n+ if (!is.null(Signal_info) && !is.matrix(Signal_info)) {\n+ stop("Signal_info is not a matrix.")\n+ }\n+ \n+ \n+ Original_data <- Signal_data\n+ \n+ # Extract the real part of the spectrum ---------------------------\n+ \n+ if (force.real) {\n+ if (is.complex(Signal_data)) {\n+ Signal_data <- Re(Signal_data)\n+ } else {\n+ # The signal is numeric Im(Signal_data) is zero anyway so let\'s avoid\n+ # using complex(real=...,imaginary=0) which would give a complex signal\n+ # in endTreatment()\n+ force.real <- FALSE\n+ }\n+ }\n+ \n+ \n+ # Return the formatted data and metadata entries --------------------\n+ \n+ return(list(start = proc.time(), vec = vec, force.real = force.real, \n+ Original_data = Original_data, Signal_data = Signal_data, Signal_info = Signal_info))\n+}\n+\n+# endTreatment \n+endTreatment <- function(name, begin_info, Signal_data) {\n+ \n+ # begin_info: object outputted from beginTreatment\n+ \n+ \n+ # Formatting the entries and printing process time -----------------------\n+ end_time <- proc.time() # record it as soon as possible\n+ start_time <- begin_info[["start"]]\n+ delta_time <- end_time - start_time\n+ delta <- delta_time[]\n+ cat("End", name, "\\n")\n+ cat("It lasted", round(delta["user.self"], 3), "s user time,", round(delta["sys.self"],3),\n+ "s system time and", round(delta["elapsed"], 3), "s elapsed time.\\n")\n+ \n+ \n+ if (begin_info[["force.real"]]) {\n+ # The imaginary part is left untouched\n+ i <- complex(real = 0, imaginary = 1)\n+ Signal_data <- Signal_data + i * Im(begin_info[["Original_data"]])\n+ }\n+ \n+ if (begin_info[["vec"]]) {\n+ Signal_data <- Signal_data[1, ]\n+ }\n+ \n+ # Return the formatted data and metadata entries --------------------\n+ return(Signal_data)\n+}\n+\n+# checkArg \n+checkArg <- function(arg, checks, can.be.null=FALSE) {\n+ check.list <- list(bool=c(is.logical, "a boolean"),\n+ int =c(function(x){x%%1==0}, "an integer"),\n+ num =c(is.numeric, "a numeric"),\n+ str =c(is.character, "a string"),\n+ pos =c(function(x){x>0}, "positive"),\n+ pos0=c(function(x){x>=0}, "positive or zero"),\n+ l1 =c(function(x){length(x)==1}, "of length 1")\n+ )\n+ if (is.null(arg)) {\n+ if (!can.be.null) {\n+ stop(deparse(substitute(arg)), " is null.")\n+ }\n+ } else {\n+ if (is.matrix(arg)) {\n+ stop(deparse(substitute(arg)), " is not scalar.")\n+ }\n+ for (c in checks) {\n+ if (!check.list[[c]][[1]](arg)) {\n+ stop(deparse(substitute(arg)), " is not ", check.list[[c]][[2]], ".")\n+ }\n+ }\n+ }\n+}\n+\n+# getArg \n+getArg <- function(arg, info, argname, can.be.absent=FALSE) {\n+ if (is.null(arg)) {\n+ start <- paste("impossible to get argument", argname, "it was not given directly and");\n+ if (!is.matrix(info)) {\n+ stop(paste(start, "the info matrix was not given"))\n+ }\n+ if (!(argname %in% colnames(info))) {\n+ if (can.be.absent) {\n+ return(NULL)\n+ } else {\n+ stop(paste(start, "is not in the info matrix"))\n+ }\n+ }\n+ if (nrow(info) < 1) {\n+ stop(paste(start, "the info matr'..b' returnBaseline = F) {\n+ \n+ # Data initialisation ----------------------------------------------\n+ begin_info <- beginTreatment("BaselineCorrection", Spectrum_data, force.real = T)\n+ Spectrum_data <- begin_info[["Signal_data"]]\n+ p <- p.bc\n+ lambda <- lambda.bc\n+ n <- dim(Spectrum_data)[1]\n+ \n+ # Data check\n+ checkArg(ptw.bc, c("bool"))\n+ checkArg(maxIter, c("int", "pos"))\n+ checkArg(lambda, c("num", "pos0"))\n+ checkArg(p.bc, c("num", "pos0"))\n+ checkArg(eps, c("num", "pos0"))\n+ checkArg(returnBaseline, c("bool"))\n+ \n+ \n+ \n+ # Baseline Correction implementation definition ----------------------\n+ \n+ # 2 Ways: either use the function asysm from the ptw package or by \n+ # built-in functions \n+ if (ptw.bc) {\n+ asysm <- ptw::asysm\n+ } else {\n+ difsmw <- function(y, lambda, w, d) {\n+ # Weighted smoothing with a finite difference penalty cf Eilers, 2003.\n+ # (A perfect smoother) \n+ # y: signal to be smoothed \n+ # lambda: smoothing parameter \n+ # w: weights (use0 zeros for missing values) \n+ # d: order of differences in penalty (generally 2)\n+ m <- length(y)\n+ W <- Matrix::Diagonal(x=w)\n+ E <- Matrix::Diagonal(m)\n+ D <- Matrix::diff(E, differences = d)\n+ C <- Matrix::chol(W + lambda * t(D) %*% D)\n+ x <- Matrix::solve(C, Matrix::solve(t(C), w * y))\n+ return(as.numeric(x))\n+ \n+ }\n+ asysm <- function(y, lambda, p, eps) {\n+ # Baseline estimation with asymmetric least squares\n+ # y: signal\n+ # lambda: smoothing parameter (generally 1e5 to 1e8)\n+ # p: asymmetry parameter (generally 0.001)\n+ # d: order of differences in penalty (generally 2)\n+ # eps: 1e-8 in ptw package\n+ m <- length(y)\n+ w <- rep(1, m)\n+ i <- 1\n+ repeat {\n+ z <- difsmw(y, lambda, w, d = 2)\n+ w0 <- w\n+ w <- p * (y > z + eps | y < 0) + (1 - p) * (y <= z + eps)\n+ if (sum(abs(w - w0)) == 0) {\n+ break\n+ }\n+ i <- i + 1\n+ if (i > maxIter) {\n+ warning("cannot find Baseline estimation in asysm")\n+ break\n+ }\n+ }\n+ return(z)\n+ }\n+ }\n+ \n+ # Baseline estimation ----------------------------------------------\n+ Baseline <- matrix(NA, nrow = nrow(Spectrum_data), ncol = ncol(Spectrum_data))\n+ \n+ for (k in 1:n) {\n+ Baseline[k, ] <- asysm(y = Spectrum_data[k, ], lambda = lambda, p = p, eps = eps)\n+ if (F & k == 1) {\n+ m <- ncol(Spectrum_data)\n+ graphics::plot(1:m, Spectrum_data[k, ], type = "l", col = "red")\n+ graphics::lines(1:m, Baseline[k, ], type = "l", col = "blue")\n+ graphics::lines(1:m, Spectrum_data[k, ] - Baseline[k, ], type = "l",\n+ col = "green")\n+ }\n+ \n+ Spectrum_data[k, ] <- Spectrum_data[k, ] - Baseline[k, ]\n+ }\n+ \n+ # Data finalisation ----------------------------------------------\n+ Spectrum_data <- endTreatment("BaselineCorrection", begin_info, Spectrum_data) # FIXME create removeImaginary filter ??\n+ \n+ if (returnBaseline) {\n+ return(list(Spectrum_data = Spectrum_data, Baseline = Baseline))\n+ } else {\n+ return(Spectrum_data)\n+ }\n+}\n+\n+\n+\n+## ====================================================\n+# NegativeValuesZeroing \n+## ====================================================\n+\n+NegativeValuesZeroing <- function(Spectrum_data) {\n+ # Data initialisation and checks ----------------------------------------------\n+ begin_info <- beginTreatment("NegativeValuesZeroing", Spectrum_data, force.real = T)\n+ Spectrum_data <- begin_info[["Signal_data"]]\n+ \n+ # NegativeValuesZeroing ----------------------------------------------\n+ Spectrum_data[Spectrum_data < 0] <- 0\n+ \n+ # Data finalisation ----------------------------------------------\n+ return(endTreatment("NegativeValuesZeroing", begin_info, Spectrum_data))\n+}\n+\n+\n+## ====================================================\n+# Warping \n+## ====================================================\n' |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/NmrPreprocessing_wrapper.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nmr_preprocessing/NmrPreprocessing_wrapper.R Fri May 05 07:13:02 2017 -0400 |
| [ |
| b'@@ -0,0 +1,320 @@\n+#!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file\n+\n+## 170116_NmrPreprocessing.R\n+## Manon Martin and Marie Tremblay-Franco\n+\n+##======================================================\n+##======================================================\n+# Preamble\n+##======================================================\n+##======================================================\n+\n+runExampleL <- FALSE\n+\n+##------------------------------\n+## Options\n+##------------------------------\n+strAsFacL <- options()$stringsAsFactors\n+options(stringsAsFactors = FALSE)\n+\n+##------------------------------\n+## Libraries laoding\n+##------------------------------\n+library(batch)\n+library(ptw)\n+library(Matrix)\n+library(ggplot2)\n+library(gridExtra)\n+library(reshape2)\n+\n+\n+# R script call\n+source_local <- function(fname)\n+{\n+\targv <- commandArgs(trailingOnly = FALSE)\n+\tbase_dir <- dirname(substring(argv[grep("--file=", argv)], 8))\n+\tsource(paste(base_dir, fname, sep="/"))\n+}\n+#Import the different functions\n+source_local("NmrPreprocessing_script.R")\n+source_local("DrawFunctions.R")\n+\n+##------------------------------\n+## Errors ?????????????????????\n+##------------------------------\n+\n+\n+##------------------------------\n+## Constants\n+##------------------------------\n+topEnvC <- environment()\n+flagC <- "\\n"\n+\n+\n+##------------------------------\n+## Script\n+##------------------------------\n+if(!runExampleL)\n+ argLs <- parseCommandArgs(evaluate=FALSE)\n+\n+# log file\n+print(argLs[["logOut"]])\n+\n+## Starting\n+cat("\\nStart of \'Preprocessing\' Galaxy module call: ", as.character(Sys.time()), sep = "")\n+\n+\n+##======================================================\n+##======================================================\n+## Parameters Loading\n+##======================================================\n+##======================================================\n+\n+# 1rst order phase correction ------------------------\n+ # Inputs\n+\t## Data matrix\n+Fid_data0 <- read.table(argLs[["dataMatrixFid"]],header=TRUE, check.names=FALSE, sep=\'\\t\')\n+Fid_data0 <- Fid_data0[,-1]\n+Fid_data0 <- as.matrix(Fid_data0)\n+\n+\t## Samplemetadata\n+samplemetadataFid <- read.table(argLs[["sampleMetadataFid"]],check.names=FALSE,header=TRUE,sep="\\t")\n+samplemetadataFid <- as.matrix(samplemetadataFid)\n+\n+\n+# water and solvent(s) correction ------------------------\n+ # Inputs\n+lambda <- argLs[["lambda"]]\n+ptwSS1 <- argLs[["ptwSS"]]\n+ptwSS <- FALSE\n+if (ptwSS1=="YES")\n+ \tptwSS <- TRUE\n+\n+\n+# apodization -----------------------------------------\n+ # Inputs\n+phase=0\n+rectRatio=1/2\n+gaussLB=1\n+expLB=1\n+apodization <- argLs[["apodizationMethod"]]\n+\n+if (apodization==\'exp\'){\n+ expLB <- argLs[["expLB"]]\n+ } else if (apodization==\'cos2\'){\n+ phase <- argLs[["phase"]]\n+ } else if (apodization==\'hanning\'){\n+ phase <- argLs[["phase"]]\n+ } else if (apodization==\'hamming\'){\n+ phase <- argLs[["phase"]]\n+ } else if (apodization==\'blockexp\'){\n+ rectRatio <- argLs[["rectRatio"]]\n+ expLB <- argLs[["expLB"]]\n+ } else if (apodization==\'blockcos2\'){\n+ rectRatio <- argLs[["rectRatio"]]\n+ } else if (apodization==\'gauss\'){\n+ rectRatio <- argLs[["rectRatio"]]\n+ gaussLB <- argLs[["gaussLB"]]\n+ }\t\t\n+\n+\n+# Fourier transform ----------------------------------\n+ # Inputs\n+FTGraph <- argLs[["FTGraph"]]\n+\n+# Internal referencering ----------------------------------\n+ # Inputs\n+shiftTreshold = 2 # c\n+ppm = TRUE\n+shiftReferencingRangeList = NULL # fromto.RC\n+pctNear0 = 0.02 # pc \n+rowindex_graph = NULL\n+ppm_ref = 0 # ppm.ref\n+\n+# \n+# shiftReferencing <- argLs[["shiftReferencing"]]\n+# print(shiftReferencing)\n+# \n+# if (shiftReferencing=="YES")\n+# {\n+# \n+\tshiftReferencingMethod <- argLs[["shiftReferencingMethod"]]\n+\n+\tif (shiftReferencingMethod == "thres")\t{\n+\t\tshiftTreshold <- argLs[["shiftTreshold"]]\n+\t}\n+\t\n+\tshiftReferencingRange <- argLs[["shiftReferencingRange"]]\n+\t\n+\tif (shiftReferencingRange == "near0"){\n+\t pctNear0 <- argLs[["pctNear0"]]\n+\t}\n+\t\n+\tif (shiftReferencingRa'..b'=============================================== \n+## Computation\n+##======================================================\n+##======================================================\n+\n+pdf(nomGraphe, onefile = TRUE, width = 13, height = 13)\n+\n+# FirstOrderPhaseCorrection ---------------------------------\n+Fid_data <- FirstOrderPhaseCorrection(Fid_data0, Fid_info = samplemetadataFid, group_delay = NULL)\n+\n+\n+# SolventSuppression ---------------------------------\n+Fid_data <- SolventSuppression(Fid_data, lambda.ss = lambda, ptw.ss = ptwSS, plotSolvent = F, returnSolvent = F)\n+\t\n+# Apodization ---------------------------------\t\n+Fid_data <- Apodization(Fid_data, Fid_info = samplemetadataFid, DT = NULL, \n+ type.apod = apodization, phase = phase, rectRatio = rectRatio, gaussLB = gaussLB, expLB = expLB, plotWindow = F, returnFactor = F)\n+\n+# FourierTransform ---------------------------------\n+Spectrum_data <- FourierTransform(Fid_data, Fid_info = samplemetadataFid, reverse.axis = TRUE)\n+\n+\n+if (FTGraph == "YES") {\n+ title = "Fourier transformed specta"\n+ DrawSignal(Spectrum_data, subtype = "stacked",\n+ ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, \n+ xlab = "Frequency", num.stacked = 4, \n+ main = title, createWindow=FALSE)\n+}\n+\n+# InternalReferencing ---------------------------------\n+# if (shiftReferencing=="YES") {\n+Spectrum_data <- InternalReferencing(Spectrum_data, samplemetadataFid, method = shiftReferencingMethod, range = shiftReferencingRange, \n+ ppm.ref = 0, shiftHandling = shiftHandling,ppm = TRUE,\n+\t\t\t\t\t\t\t\t\t\t\t c = shiftTreshold, fromto.RC = shiftReferencingRangeList, pc = pctNear0)\n+\n+# }\n+\n+# ZeroOrderPhaseCorrection ---------------------------------\n+Spectrum_data <- ZeroOrderPhaseCorrection(Spectrum_data, method = zeroOrderPhaseMethod, \n+ plot_rms = plot_rms, returnAngle = returnAngle, \n+ createWindow = createWindow,angle = angle, \n+ plot_spectra = plot_spectra, quant = quant, \n+ ppm = ppm, fromto.0OPC = fromto.0OPC)\n+\n+title = "Spectra after Zero Order Phase Correction"\n+DrawSignal(Spectrum_data, subtype = "stacked",\n+ ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, \n+ xlab = "Frequency", num.stacked = 4, \n+ main = title, createWindow=FALSE)\n+\n+# BaselineCorrection ---------------------------------\t\t\t\t\t\t\t\t\t \n+Spectrum_data <- BaselineCorrection(Spectrum_data, ptw.bc = TRUE, maxIter = maxIter, lambda.bc = lambdaBc, p.bc = pBc, eps = epsilon, returnBaseline = F) \n+\n+title = "Spectra after Baseline Correction"\n+DrawSignal(Spectrum_data, subtype = "stacked",\n+ ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, \n+ xlab = "Frequency", num.stacked = 4, \n+ main = title, createWindow=FALSE)\n+\n+# NegativeValuesZeroing ---------------------------------\n+if (NegativetoZero=="YES") {\n+\n+ Spectrum_data <- NegativeValuesZeroing(Spectrum_data)\n+\n+ title = "Spectra after Negative Values Zeroing"\n+ DrawSignal(Spectrum_data, subtype = "stacked",\n+ ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, \n+ xlab = "Frequency", num.stacked = 4, \n+ main = title, createWindow=FALSE)\n+ }\n+\n+invisible(dev.off())\n+\n+##======================================================\n+##======================================================\n+## Saving\n+##======================================================\n+##======================================================\n+\n+# Data Matrix\n+write.table(t(Re(Spectrum_data)),file=argLs$dataMatrixOut, quote=FALSE, row.names=TRUE, sep="\\t", col.names=TRUE)\n+\n+\n+## Ending\n+\n+cat("\\nEnd of \'Preprocessing\' Galaxy module call: ", as.character(Sys.time()), sep = "")\n+\n+sink()\n+\n+options(stringsAsFactors = strAsFacL)\n+\n+rm(list = ls())\n' |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/NmrPreprocessing_xml.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nmr_preprocessing/NmrPreprocessing_xml.xml Fri May 05 07:13:02 2017 -0400 |
| b |
| b'@@ -0,0 +1,366 @@\n+<tool id="NMR_Preprocessing" name="NMR_Preprocessing" version="3.0.1">\r\n+\t<description> Preprocessing of 1D NMR spectra </description>\r\n+\r\n+ <stdio>\r\n+ <exit_code range="1:" level="fatal" />\r\n+ </stdio>\r\n+\t\r\n+ \t<command>\r\n+ \t ## Wrapper\r\n+ Rscript $__tool_directory__/NmrPreprocessing_wrapper.R\r\n+\r\n+\t\t## First order phase correction\r\n+\t\t\t## Data matrix of FID spectra\r\n+\t\tdataMatrixFid $dataMatrixFid\r\n+\t\t\t## Sample metadata matrix\r\n+\t\tsampleMetadataFid $sampleMetadataFid\r\n+\r\n+\t\t\r\n+\t\t## Water and / or solvents suppression\r\n+\t\t ## Smoothing parameter\r\n+ lambda $lambda\r\n+\t\tptwSS $ptwSS\r\n+\t\t\r\n+\t\t\r\n+\t\t## Apodization\r\n+\t apodizationMethod $apodizationMethod.method\r\n+\t\t #if $apodizationMethod.method == "exp":\r\n+\t\t\t ## Line broadening for the exponential window\r\n+\t\t\t expLB $apodizationMethod.expLB \r\n+\t\t #end if\r\n+\t\t #if $apodizationMethod.method == "cos2":\r\n+\t\t\t ## Phase\r\n+\t\t\t phase $apodizationMethod.phase\r\n+\t\t #end if\r\n+\t\t #if $apodizationMethod.method == "hanning":\r\n+\t\t\t ## Phase\r\n+\t\t\t phase $apodizationMethod.phase\r\n+\t\t #end if\r\n+\t\t #if $apodizationMethod.method == "hamming":\r\n+\t\t\t ## Phase\r\n+\t\t\t phase $apodizationMethod.phase\r\n+\t\t #end if\r\n+\t\t #if $apodizationMethod.method == "blockexp":\r\n+\t\t\t ## Proportion of signal in the window\r\n+\t\t\t rectRatio $apodizationMethod.rectRatio\r\n+\t\t\t expLB $apodizationMethod.expLB \r\n+\t\t #end if\t\t\t\r\n+\t\t #if $apodizationMethod.method == "blockcos2":\r\n+\t\t\t ## Proportion of signal in the window\r\n+\t\t\t rectRatio $apodizationMethod.rectRatio\r\n+\t\t #end if\t\t\t\r\n+\t\t #if $apodizationMethod.method == "gauss":\r\n+\t\t\t ## Line broadening for the gaussian window\r\n+\t\t\t gaussLB $apodizationMethod.gaussLB \r\n+\t\t #end if\r\n+\t\r\n+\t\r\n+\t\t## Fourier transform\r\n+\t\t\t## Graphical display\r\n+\t\t\tFTGraph $FTGraph\r\n+\t\t\r\n+\t\t\r\n+\t\t## Shift referencing\r\n+\t\t\t## Method used to find the TMSP peaks in spectra\r\n+\t\t\tshiftReferencingMethod $shiftReferencingMethod.method\r\n+\t\t\t#if $shiftReferencingMethod.method == "thres":\r\n+\t\t\t\tshiftTreshold $shiftReferencingMethod.shiftTreshold\r\n+\t\t\t#end if\r\n+\t\t\t## Definition of the search zone\r\n+\t\t\tshiftReferencingRange $shiftReferencingRange.method\r\n+\t\t\t#if $shiftReferencingRange.method == "near0":\r\n+\t\t\t\tpctNear0 $shiftReferencingRange.pctNear0\r\n+\t\t\t#end if\r\n+\t\t\t#if $shiftReferencingRange.method == "window":\r\n+\t\t\t\t#for $i in $shiftReferencingRange.conditions:\r\n+ \t \tshiftReferencingRangeLeft ${i.shiftReferencingRangeLeft}\r\n+ \t \tshiftReferencingRangeRight ${i.shiftReferencingRangeRight}\r\n+ \t #end for\r\n+\t\t\t#end if\r\n+\t\t\tshiftHandling $shiftHandling\r\n+\t\t\r\n+\t\t\r\n+\t\t## Zero order phase correction\r\n+\t\tzeroOrderPhaseMethod $zeroOrderPhaseMethod.method\r\n+\t #if $zeroOrderPhaseMethod.method == "rms":\r\n+\t\t\t ## Line broadening for the exponential window\r\n+\t\t\t quant $zeroOrderPhaseMethod.quant \r\n+\t\t#end if\r\n+\t\tsearchZoneZeroPhase.choice ${searchZoneZeroPhase.choice}\r\n+ #if str($searchZoneZeroPhase.choice) == "YES":\r\n+ #for $i in $searchZoneZeroPhase.conditions:\r\n+ searchZoneZeroPhase_left ${i.searchZoneZeroPhase_left}\r\n+ searchZoneZeroPhase_right ${i.searchZoneZeroPhase_right}\r\n+ #end for\r\n+ #end if\t\t\r\n+\r\n+\r\n+\t\t## Baseline correction\r\n+\t\tptwBc $ptwBc\r\n+\t\tmaxIter $maxIter\r\n+\t\tlambdaBc $lambdaBc\r\n+\t\tpBc $pBc\r\n+\t\tepsilon $epsilon\r\n+\t\t\r\n+\t\t\r\n+\t\t## sets negative intensities to zero\r\n+\t\tNegativetoZero $NegativetoZero\r\n+\t\t\t\t\r\n+\t\t\r\n+ ## Outputs\r\n+ dataMatrixOut $dataMatrixOut\r\n+\t\tgraphOut $graphOut\r\n+\t\tlogOut $logOut\r\n+\t\t\r\n+\t</command>\r\n+\t\t\r\n+ \t<inputs>\r\n+\t\t<param name="dataMatrixFid" type="data" label="Data matrix of FIDs" help="" format="tabular" />\r\n+\t\t<param name="sampleMetadataFid" type="data" label="Sample metadata matrix" help="" format="tabular" />\r\n+\t\t\r\n+ \t<param name="lambda" label="Smoothing parameter for solvent suppression" type="float" value="1'..b'ption\r\n+-----------\r\n+\r\n+These steps correspond to the following steps in the PEPS-NMR R library (https://github.com/ManonMartin/PEPSNMR):\r\n+\r\n+* First order phase correction\r\n+* Removal of solvent residuals signal from the FID\r\n+* Apodization to increase the Signal-to-Noise ratio of the FID\r\n+* Fourier transformation\r\n+* Shift referencing to calibrate the spectra with internal compound referencing\r\n+* Zero order phase correction\r\n+* Baseline correction\r\n+* Setting of negatives values to 0\r\n+\r\n+----------------------\r\n+ Steps parameters \r\n+----------------------\r\n+\r\n+**Apodization**\r\n+----------------------\r\n+\r\n+The **types of apodization** are:\r\n+\r\n+\r\n+* exp: The signal is multiplied by a decreasing exponential exp(-t/LineBroadening).\r\n+\r\n+* cos2: The signal is multiplied by the value of a cosinus squared from 0 (where its value is 1) until pi/2 (where its value is 0).\r\n+\r\n+* blockexp: The first part of the signal (described by the proportion of signal in the window) is left unchanged and the second is multiplied by exp(-t/LineBroadening) starting at value 1.\r\n+\r\n+* blockcos2: The first part is left unchanged as with blockexp and the second part is multiplied by a cosinus squared where its value starts at 1 at the end of the block and ends at 0 at the end of the signal.\r\n+\r\n+* gauss: The signal is multiplied by a gaussian window centered at the beginning of the FID and with sigma=1/LineBroadening.\r\n+\r\n+* hanning: The signal is multiplied by a hanning window : 0.5 + 0.5 cos.\r\n+\r\n+* hamming: The signal is multiplied by a hamming window : 0.54 + 0.46 cos.\r\n+\r\n+\r\n+**Shift referencing**\r\n+----------------------\r\n+\r\n+Different **methods for shift referencing**:\r\n+\r\n+* max: the maximum intensity in the search zone is defined as the reference compound peak.\r\n+\r\n+* thres: the reference compound peak is the first peak in the spectrum higher than a predefined threshold (c) which is computed as: c*(cumulated_mean/cumulated_sd).\r\n+\r\n+\r\n+The **searching window** can be adapted:\r\n+\r\n+* near0: the search concentrates around the 0ppm.\r\n+\r\n+* all: search accross the whole ppm axis.\r\n+\r\n+* window: the search is operated in the windows defined by the search_zone bounds.\r\n+\r\n+\r\n+**shiftHandling**: spectra can be shifted differently, we can handle misalignment of the left and right of the spectrum by different ways:\r\n+\r\n+* zerofilling: The extremities at which a spectrum is not defined are replaced by 0. It makes sense since in practice the spectrum is close to zero at the extremities.\r\n+\r\n+* NAfilling: The extremities at which a spectrum is not defined are replaced by NA. \r\n+\r\n+* circular: The spectra are shifted circularly which means that the end of a spectrum is reproduced at the beginning. \r\n+\r\n+* cut: The ppm values for which some spectra are not defined are removed.\r\n+\r\n+\r\n+**ppm value of the reference compound**: By default, the ppm value of the reference compound is set to 0, but any arbitrary value in the ppm interval can be used instead.\r\n+\r\n+\r\n+**Zero Order Phase correction**\r\n+-----------------------------------\r\n+\r\n+**Zero Order Phase correction method**:\r\n+\r\n+* rms: A positiveness criterion is applied on the spectrum with a quantile probability parameter to trim the values.\r\n+\r\n+* manual: A vector of angles are specified wth the angle argument to manually rotate the spectra.\r\n+\r\n+* max: Optimization of the maximal spectral intensity.\r\n+\r\n+\r\n+**Search zone for the Zero order phase correction**: enables to optimize the criterion only on the selected spectral window(s).\r\n+\r\n+\r\n+**Baseline computation**\r\n+----------------------------\r\n+**Baseline computation method**: either use ptw (YES) or an in-house R routine (NO).\r\n+\r\n+**Smoothing parameter**: the larger it is, the smoother the estimated baseline will be.\r\n+\r\n+**Asymmetry parameter**: the smaller it is, the less the estimated baseline will try to follow peaks when it is under the spectrum and the more it will try to be under the spectrum.\r\n+\r\n+\r\n+\r\n+ </help>\r\n+ \r\n+</tool>\r\n' |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/ReadFids.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nmr_preprocessing/ReadFids.R Fri May 05 07:13:02 2017 -0400 |
| [ |
| b'@@ -0,0 +1,382 @@\n+################################################################################################\r\n+#\r\n+# Read FIDs in Bruker format\r\n+#\r\n+#\r\n+################################################################################################\r\n+\r\n+# ReadFid ==============================================================================\r\n+ReadFid <- function(path) {\r\n+ \r\n+ # Read 1D FID using Bruker XWinNMR and TopSpin format. It is inspired of the\r\n+ # matNMR matlab library which deals with 2D FID and also other formats\r\n+ # Read also the parameters in the acqus file\r\n+ \r\n+ paramFile <- file.path(path, "acqus")\r\n+ # BYTEORDA: 0 -> Little Endian 1 -> Big Endian\r\n+ params <- readParams(paramFile, c("TD", "BYTORDA", "DIGMOD", "DECIM", "DSPFVS", \r\n+ "SW_h", "SW", "O1"))\r\n+ \r\n+ if (params[["DSPFVS"]] >= 20) {\r\n+ # The group delay first order phase correction is given directly from version 20\r\n+ grpdly <- readParams(paramFile, c("GRPDLY"))\r\n+ params[["GRPDLY"]] <- grpdly[["GRPDLY"]]\r\n+ }\r\n+ TD <- params[["TD"]]\r\n+ endianness <- if (params$BYTORDA) \r\n+ "big" else "little"\r\n+ if (TD%%2 != 0) {\r\n+ stop(paste("Only even numbers are allowed for size in TD because it is complex \r\n+ data with the real and imaginary part for each element.", \r\n+ "The TD value is in the", paramFile, "file"))\r\n+ }\r\n+ \r\n+ # Interpret params Dwell Time, time between 2 data points in the FID\r\n+ params[["DT"]] <- 1/(2 * params[["SW_h"]])\r\n+ \r\n+ # Read fid\r\n+ fidFile <- file.path(path, "fid")\r\n+ fidOnDisk <- readBin(fidFile, what = "int", n = TD, size = 4L, endian = endianness)\r\n+ \r\n+ # Real size that is on disk (it should be equal to TD2, except for TopSpin/Bruker\r\n+ # (which is our case) according to matNMR as just discussed\r\n+ TDOnDisk <- length(fidOnDisk)\r\n+ if (TDOnDisk < TD) {\r\n+ warning("Size is smaller than expected, the rest is filled with zero so the size is the same for every fid")\r\n+ fidGoodSize <- sapply(vector("list", length = TD), function(x) 0)\r\n+ fidGoodSize[1:TDOnDisk] <- fidOnDisk\r\n+ \r\n+ } else if (TDOnDisk > TD) {\r\n+ warning("Size is bigger than expected, the rest ignored so the size is the same for every fid")\r\n+ fidGoodSize <- fidOnDisk(1:TD)\r\n+ \r\n+ } else {\r\n+ fidGoodSize <- fidOnDisk\r\n+ }\r\n+ \r\n+ fidRePart <- fidGoodSize[seq(from = 1, to = TD, by = 2)]\r\n+ fidImPart <- fidGoodSize[seq(from = 2, to = TD, by = 2)]\r\n+ fid <- complex(real = fidRePart, imaginary = fidImPart)\r\n+ \r\n+ return(list(fid = fid, params = params))\r\n+}\r\n+\r\n+\r\n+\r\n+\r\n+# getDirsContainingFid ==============================================================================\r\n+getDirsContainingFid <- function(path) {\r\n+ subdirs <- dir(path, full.names = TRUE)\r\n+ if (length(subdirs) > 0) {\r\n+ cond <- sapply(subdirs, function(x) {\r\n+ content <- dir(x)\r\n+ # subdirs must contain fid, acqu and acqus files\r\n+ return("fid" %in% content && "acqu" %in% content && "acqus" %in% content)\r\n+ })\r\n+ subdirs <- subdirs[cond]\r\n+ }\r\n+ return(subdirs)\r\n+}\r\n+\r\n+\r\n+\r\n+\r\n+\r\n+\r\n+# beginTreatment ==============================================================================\r\n+\r\n+beginTreatment <- function(name, Signal_data = NULL, Signal_info = NULL, \r\n+ force.real = FALSE) {\r\n+ \r\n+ cat("Begin", name, "\\n")\r\n+ \r\n+ \r\n+ # Formatting the Signal_data and Signal_info -----------------------\r\n+ \r\n+ vec <- is.vector(Signal_data)\r\n+ if (vec) {\r\n+ Signal_data <- vec2mat(Signal_data)\r\n+ }\r\n+ if (is.vector(Signal_info)) {\r\n+ Signal_info <- vec2mat(Signal_info)\r\n+ }\r\n+ if (!is.null(Signal_data)) {\r\n+ if (!is.matrix(Signal_data)) {\r\n+ stop("Signal_data is not a matrix.")\r\n+ }\r\n+ if (!is.complex(Signal_data) && !is.numeric(Signal_data)) {\r\n+ stop("Signal_data contains non-numerical values.")\r\n+ }\r\n+ }\r\n+ if (!is.null(Signal_info) && !is.matrix(Signal_info)) {\r\n+ stop("S'..b'ttern <- paste("\\\\$", paramName, "=", sep = "")\r\n+ occurences <- grep(pattern, lines)\r\n+ if (length(occurences) == 0L) {\r\n+ stop(paste(file, "has no field", pattern))\r\n+ }\r\n+ if (length(occurences) > 1L) {\r\n+ warning(paste(file, "has more that one field", pattern, " I take the first one"))\r\n+ }\r\n+ line <- lines[occurences[1]]\r\n+ \r\n+ # Cut beginning and end of the line \'##$TD= 65536\' -> \'65536\'\r\n+ igual = as.numeric(regexpr("=", line))\r\n+ \r\n+ first <- igual\r\n+ while (first <= nchar(line) & !isDigit(substr(line, first, first))) {\r\n+ first <- first + 1\r\n+ }\r\n+ last <- nchar(line)\r\n+ while (last > 0 & !isDigit(substr(line, last, last))) {\r\n+ last <- last - 1\r\n+ }\r\n+ params[paramName] <- as.numeric(substr(line, first, last))\r\n+ }\r\n+ return(params)\r\n+}\r\n+\r\n+\r\n+\r\n+# ReadFids ==============================================================================\r\n+ReadFids <- function(path, l = 1, subdirs = FALSE) {\r\n+ \r\n+ # Data initialisation and checks ----------------------------------------------\r\n+ begin_info <- beginTreatment("ReadFids")\r\n+ checkArg(path, c("str"))\r\n+ checkArg(l, c("pos"))\r\n+ if (file.exists(path) == FALSE) {\r\n+ stop(paste("Invalid path:", path))\r\n+ }\r\n+ \r\n+ \r\n+ # Extract the FIDs and their info ----------------------------------------------\r\n+ \r\n+ if (subdirs == FALSE) {\r\n+ fidDirs <- getDirsContainingFid(path)\r\n+ n <- length(fidDirs)\r\n+ if (n == 0L) {\r\n+ stop(paste("No valid fid in", path))\r\n+ }\r\n+ fidNames <- sapply(X = fidDirs, FUN = getTitle, l = l, subdirs = subdirs, USE.NAMES = F)\r\n+ for (i in 1:n) {\r\n+ fidList <- ReadFid(fidDirs[i])\r\n+ fid <- fidList[["fid"]]\r\n+ info <- fidList[["params"]]\r\n+ m <- length(fid)\r\n+ if (i == 1) {\r\n+ Fid_data <- matrix(nrow = n, ncol = m, dimnames = list(fidNames, \r\n+ info[["DT"]] * (0:(m - 1))))\r\n+ Fid_info <- matrix(nrow = n, ncol = length(info), dimnames = list(fidNames, \r\n+ names(info)))\r\n+ }\r\n+ Fid_data[i, ] <- fid\r\n+ Fid_info[i, ] <- unlist(info)\r\n+ }\r\n+ \r\n+ } else {\r\n+ maindirs <- dir(path, full.names = TRUE)\r\n+ Fid_data <- numeric()\r\n+ Fid_info <- numeric()\r\n+ \r\n+ fidDirs <- c()\r\n+ for (j in maindirs) {\r\n+ fd <- getDirsContainingFid(j)\r\n+ n <- length(fd)\r\n+ if (n > 0L) {\r\n+ fidDirs <- c(fidDirs, fd)\r\n+ } else {warning(paste("No valid fid in",j ))}\r\n+ }\r\n+ \r\n+ fidNames <- sapply(X = fidDirs, FUN = getTitle, l = l, subdirs = subdirs, USE.NAMES = F)\r\n+ for (i in 1:length(fidNames)) {\r\n+ fidList <- ReadFid(fidDirs[i])\r\n+ fid <- fidList[["fid"]]\r\n+ info <- fidList[["params"]]\r\n+ m <- length(fid)\r\n+ if (i == 1) {\r\n+ Fid_data <- matrix(nrow = length(fidNames), ncol = m, dimnames = list(fidNames, \r\n+ info[["DT"]] * (0:(m - 1))))\r\n+ Fid_info <- matrix(nrow = length(fidNames), ncol = length(info), dimnames = list(fidNames, \r\n+ names(info)))\r\n+ }\r\n+ Fid_data[i, ] <- fid\r\n+ Fid_info[i, ] <- unlist(info)\r\n+ }\r\n+ \r\n+ \r\n+ }\r\n+ \r\n+ # Check for non-unique IDs ----------------------------------------------\r\n+ NonnuniqueIds <- sum(duplicated(row.names(Fid_data)))\r\n+ cat("dim Fid_data: ", dim(Fid_data), "\\n")\r\n+ cat("IDs: ", rownames(Fid_data), "\\n")\r\n+ cat("non-unique IDs?", NonnuniqueIds, "\\n")\r\n+ if (NonnuniqueIds > 0) {\r\n+ warning("There are duplicated IDs: ", Fid_data[duplicated(Fid_data)])\r\n+ }\r\n+ \r\n+ \r\n+ # Return the results ----------------------------------------------\r\n+ return(list(Fid_data = endTreatment("ReadFids", begin_info, Fid_data), Fid_info = Fid_info))\r\n+ \r\n+}\r\n+\r\n+\r\n' |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/ReadFids_Manon.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nmr_preprocessing/ReadFids_Manon.R Fri May 05 07:13:02 2017 -0400 |
| [ |
| b'@@ -0,0 +1,386 @@\n+################################################################################################\r\n+#\r\n+# Read FIDs in Bruker format\r\n+#\r\n+#\r\n+################################################################################################\r\n+\r\n+# ReadFid ==============================================================================\r\n+ReadFid <- function(path) {\r\n+ \r\n+ # Read 1D FID using Bruker XWinNMR and TopSpin format. It is inspired of the\r\n+ # matNMR matlab library which deals with 2D FID and also other formats\r\n+ # Read also the parameters in the acqus file\r\n+ \r\n+ paramFile <- file.path(path, "acqus")\r\n+ # BYTEORDA: 0 -> Little Endian 1 -> Big Endian\r\n+ params <- readParams(paramFile, c("TD", "BYTORDA", "DIGMOD", "DECIM", "DSPFVS", \r\n+ "SW_h", "SW", "O1"))\r\n+ \r\n+ if (params[["DSPFVS"]] >= 20) {\r\n+ # The group delay first order phase correction is given directly from version 20\r\n+ grpdly <- readParams(paramFile, c("GRPDLY"))\r\n+ params[["GRPDLY"]] <- grpdly[["GRPDLY"]]\r\n+ }\r\n+ TD <- params[["TD"]]\r\n+ endianness <- if (params$BYTORDA) \r\n+ "big" else "little"\r\n+ if (TD%%2 != 0) {\r\n+ stop(paste("Only even numbers are allowed for size in TD because it is complex \r\n+ data with the real and imaginary part for each element.", \r\n+ "The TD value is in the", paramFile, "file"))\r\n+ }\r\n+ \r\n+ # Interpret params Dwell Time, time between 2 data points in the FID\r\n+ params[["DT"]] <- 1/(2 * params[["SW_h"]])\r\n+ \r\n+ # Read fid\r\n+ fidFile <- file.path(path, "fid")\r\n+ fidOnDisk <- readBin(fidFile, what = "int", n = TD, size = 4L, endian = endianness)\r\n+ \r\n+ # Real size that is on disk (it should be equal to TD2, except for TopSpin/Bruker\r\n+ # (which is our case) according to matNMR as just discussed\r\n+ TDOnDisk <- length(fidOnDisk)\r\n+ if (TDOnDisk < TD) {\r\n+ warning("Size is smaller than expected, the rest is filled with zero so the size is the same for every fid")\r\n+ fidGoodSize <- sapply(vector("list", length = TD), function(x) 0)\r\n+ fidGoodSize[1:TDOnDisk] <- fidOnDisk\r\n+ \r\n+ } else if (TDOnDisk > TD) {\r\n+ warning("Size is bigger than expected, the rest ignored so the size is the same for every fid")\r\n+ fidGoodSize <- fidOnDisk(1:TD)\r\n+ \r\n+ } else {\r\n+ fidGoodSize <- fidOnDisk\r\n+ }\r\n+ \r\n+ fidRePart <- fidGoodSize[seq(from = 1, to = TD, by = 2)]\r\n+ fidImPart <- fidGoodSize[seq(from = 2, to = TD, by = 2)]\r\n+ fid <- complex(real = fidRePart, imaginary = fidImPart)\r\n+ \r\n+ return(list(fid = fid, params = params))\r\n+}\r\n+\r\n+\r\n+\r\n+\r\n+# getDirsContainingFid ==============================================================================\r\n+getDirsContainingFid <- function(path) {\r\n+ subdirs <- dir(path, full.names = TRUE)\r\n+ if (length(subdirs) > 0) {\r\n+ cond <- sapply(subdirs, function(x) {\r\n+ content <- dir(x)\r\n+ # subdirs must contain fid, acqu and acqus files\r\n+ return("fid" %in% content && "acqu" %in% content && "acqus" %in% content)\r\n+ })\r\n+ subdirs <- subdirs[cond]\r\n+ }\r\n+ return(subdirs)\r\n+}\r\n+\r\n+\r\n+\r\n+\r\n+\r\n+\r\n+# beginTreatment ==============================================================================\r\n+\r\n+beginTreatment <- function(name, Signal_data = NULL, Signal_info = NULL, \r\n+ force.real = FALSE) {\r\n+ \r\n+ cat("Begin", name, "\\n")\r\n+ \r\n+ \r\n+ # Formatting the Signal_data and Signal_info -----------------------\r\n+ \r\n+ vec <- is.vector(Signal_data)\r\n+ if (vec) {\r\n+ Signal_data <- vec2mat(Signal_data)\r\n+ }\r\n+ if (is.vector(Signal_info)) {\r\n+ Signal_info <- vec2mat(Signal_info)\r\n+ }\r\n+ if (!is.null(Signal_data)) {\r\n+ if (!is.matrix(Signal_data)) {\r\n+ stop("Signal_data is not a matrix.")\r\n+ }\r\n+ if (!is.complex(Signal_data) && !is.numeric(Signal_data)) {\r\n+ stop("Signal_data contains non-numerical values.")\r\n+ }\r\n+ }\r\n+ if (!is.null(Signal_info) && !is.matrix(Signal_info)) {\r\n+ stop("S'..b'as no field", pattern))\r\n+ }\r\n+ if (length(occurences) > 1L) {\r\n+ warning(paste(file, "has more that one field", pattern, " I take the first one"))\r\n+ }\r\n+ line <- lines[occurences[1]]\r\n+ \r\n+ # Cut beginning and end of the line \'##$TD= 65536\' -> \'65536\'\r\n+ igual = as.numeric(regexpr("=", line))\r\n+ \r\n+ first <- igual\r\n+ while (first <= nchar(line) & !isDigit(substr(line, first, first))) {\r\n+ first <- first + 1\r\n+ }\r\n+ last <- nchar(line)\r\n+ while (last > 0 & !isDigit(substr(line, last, last))) {\r\n+ last <- last - 1\r\n+ }\r\n+ params[paramName] <- as.numeric(substr(line, first, last))\r\n+ }\r\n+ return(params)\r\n+}\r\n+\r\n+\r\n+\r\n+# ReadFids ==============================================================================\r\n+ReadFids <- function(path, l = 1, subdirs = FALSE, subdirs.names = FALSE) {\r\n+\r\n+ # Data initialisation and checks ----------------------------------------------\r\n+ begin_info <- beginTreatment("ReadFids")\r\n+ checkArg(path, c("str"))\r\n+ checkArg(l, c("pos"))\r\n+ if (file.exists(path) == FALSE) {\r\n+ stop(paste("Invalid path:", path))\r\n+ }\r\n+ \r\n+ \r\n+ # Extract the FIDs and their info ----------------------------------------------\r\n+ \r\n+ if (subdirs == FALSE) {\r\n+ fidDirs <- getDirsContainingFid(path)\r\n+ n <- length(fidDirs)\r\n+ if (n == 0L) {\r\n+ stop(paste("No valid fid in", path))\r\n+ }\r\n+ fidNames <- sapply(X = fidDirs, FUN = getTitle, l = l, subdirs = subdirs, USE.NAMES = F)\r\n+ for (i in 1:n) {\r\n+ fidList <- ReadFid(fidDirs[i])\r\n+ fid <- fidList[["fid"]]\r\n+ info <- fidList[["params"]]\r\n+ m <- length(fid)\r\n+ if (i == 1) {\r\n+ Fid_data <- matrix(nrow = n, ncol = m, dimnames = list(fidNames, \r\n+ info[["DT"]] * (0:(m - 1))))\r\n+ Fid_info <- matrix(nrow = n, ncol = length(info), dimnames = list(fidNames, \r\n+ names(info)))\r\n+ }\r\n+ Fid_data[i, ] <- fid\r\n+ Fid_info[i, ] <- unlist(info)\r\n+ }\r\n+ \r\n+ } else {\r\n+ maindirs <- dir(path, full.names = TRUE) # subdirectories\r\n+ Fid_data <- numeric()\r\n+ Fid_info <- numeric()\r\n+ \r\n+ fidDirs <- c()\r\n+ for (j in maindirs) {\r\n+ fd <- getDirsContainingFid(j) # recoved FIDs from subdirectories\r\n+ n <- length(fd)\r\n+ if (n > 0L) {\r\n+ fidDirs <- c(fidDirs, fd)\r\n+ } else {warning(paste("No valid fid in",j ))}\r\n+ }\r\n+ \r\n+ if (subdirs.names==TRUE) {\r\n+ fidNames <- dir(path)\r\n+ } else {fidNames <- sapply(X = fidDirs, FUN = getTitle, l = l, subdirs = subdirs, USE.NAMES = F)}\r\n+ \r\n+ for (i in 1:length(fidNames)) {\r\n+ fidList <- ReadFid(fidDirs[i])\r\n+ fid <- fidList[["fid"]]\r\n+ info <- fidList[["params"]]\r\n+ m <- length(fid)\r\n+ if (i == 1) {\r\n+ Fid_data <- matrix(nrow = length(fidNames), ncol = m, dimnames = list(fidNames, \r\n+ info[["DT"]] * (0:(m - 1))))\r\n+ Fid_info <- matrix(nrow = length(fidNames), ncol = length(info), dimnames = list(fidNames, \r\n+ names(info)))\r\n+ }\r\n+ Fid_data[i, ] <- fid\r\n+ Fid_info[i, ] <- unlist(info)\r\n+ }\r\n+ \r\n+ \r\n+ }\r\n+ \r\n+ # Check for non-unique IDs ----------------------------------------------\r\n+ NonnuniqueIds <- sum(duplicated(row.names(Fid_data)))\r\n+ cat("dim Fid_data: ", dim(Fid_data), "\\n")\r\n+ cat("IDs: ", rownames(Fid_data), "\\n")\r\n+ cat("non-unique IDs?", NonnuniqueIds, "\\n")\r\n+ if (NonnuniqueIds > 0) {\r\n+ warning("There are duplicated IDs: ", Fid_data[duplicated(Fid_data)])\r\n+ }\r\n+ \r\n+ \r\n+ # Return the results ----------------------------------------------\r\n+ return(list(Fid_data = endTreatment("ReadFids", begin_info, Fid_data), Fid_info = Fid_info))\r\n+ \r\n+}\r\n+\r\n+\r\n+\r\n' |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/ReadFids_wrapper.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nmr_preprocessing/ReadFids_wrapper.R Fri May 05 07:13:02 2017 -0400 |
| [ |
| @@ -0,0 +1,148 @@ +#!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file + +## 08122016_ReadFids_wrapper.R +## Manon Martin +## manon.martin@uclouvain.be + +##====================================================== +##====================================================== +# Preamble +##====================================================== +##====================================================== + +runExampleL <- FALSE + + +##------------------------------ +## Options +##------------------------------ +strAsFacL <- options()$stringsAsFactors +options(stringsAsFactors = FALSE) +options(warn=1) + +##------------------------------ +## Libraries laoding +##------------------------------ +library(batch) +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("ReadFids_Manon.R") +source_local("DrawFunctions.R") +##------------------------------ +## Errors ????????????????????? +##------------------------------ + + +##------------------------------ +## Constants +##------------------------------ +topEnvC <- environment() +flagC <- "\n" + + +##------------------------------ +## Script +##------------------------------ +if(!runExampleL) + argLs <- parseCommandArgs(evaluate=FALSE) + + +##====================================================== +##====================================================== +## Parameters Loading +##====================================================== +##====================================================== + + ## Inputs + # Path + ## Bruker FIDs +fileType="Bruker" +zipfile= argLs[["fidzipfile"]] +directory=unzip(zipfile, list=F) +path=paste(getwd(),strsplit(directory[1],"/")[[1]][2],sep="/") + + +# other inputs from ReadFids +l = argLs[["title_line"]] +subdirs <- argLs[["subdirectories"]] +subdirs.names <- argLs[["subdirs_names"]] + + +# Outputs +dataMatrixOut <- argLs[["dataMatrixOut"]] +sampleMetadataOut <- argLs[["sampleOut"]] +logOut <- argLs[["logOut"]] +nomGraphe <- argLs[["graphOut"]] + + + + +## Checking arguments +##------------------- +error.stock <- "\n" + +if(length(error.stock) > 1) + stop(error.stock) + + + +##====================================================== +##====================================================== +## Computation +##====================================================== +##====================================================== +sink(logOut,append=TRUE) + +if(length(warnings())>0){ # or !is.null(warnings()) + print("something happened") +} + +## Starting +cat("\nStart of 'ReadFids' Galaxy module call: ", as.character(Sys.time()), "\n\n", sep = "") + +outputs <- ReadFids(path = path, l=l, subdirs = subdirs, subdirs.names = subdirs.names) + +data_matrix <- outputs[["Fid_data"]] # Data matrix +data_sample <- outputs[["Fid_info"]] # Sample metadata + +pdf(nomGraphe, onefile = TRUE, width = 13, height = 13) +title = "Raw FID data" +DrawSignal(data_matrix, subtype = "stacked", + ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, + xlab = "Frequency", num.stacked = 4, + main = title, createWindow=FALSE) +invisible(dev.off()) + +##====================================================== +##====================================================== +## Saving +##====================================================== +##====================================================== + +# Data matrix +write.table(data_matrix,file=argLs$dataMatrixOut, quote=FALSE, row.names=TRUE, sep="\t", col.names=TRUE) + +# Sample metadata +write.table(data_sample,file=argLs$sampleOut, quote=FALSE, row.names=TRUE, sep="\t", col.names=TRUE) + +## Ending + +cat("\nEnd of 'ReadFids' Galaxy module call: ", as.character(Sys.time()), sep = "") + + +sink() + +options(stringsAsFactors = strAsFacL) + +rm(list = ls()) \ No newline at end of file |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/ReadFids_xml.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nmr_preprocessing/ReadFids_xml.xml Fri May 05 07:13:02 2017 -0400 |
| b |
| @@ -0,0 +1,174 @@ +<tool id="NMR_Read" name="NMR_Read" version="3.0.1"> + <description> Read NMR raw files</description> + + <stdio> + <exit_code range="1:" level="fatal" /> + </stdio> + + <command> + ## Wrapper + Rscript $__tool_directory__/ReadFids_wrapper.R + + ## File input + fidzipfile $fidzipfile + + ## Title line + title_line $title_line + + ## Subdirectories + subdirectories $subdirectories + + ## Use subdirectories names as FID names? + subdirs_names $subdirs_names + + ## Outputs + dataMatrixOut $dataMatrixOut + sampleOut $sampleOut + logOut $logOut + graphOut $graphOut + + </command> + + <inputs> + <param name="fidzipfile" type="data" format="no_unzip.zip" label="Bruker FID file" /> + + <param name="title_line" label="FID title" type="integer" value="1" size="100" help="Default value is line 1"/> + + <param name="subdirectories" label="Presence of subdirectories?" type="select" help="Select 'FALSE' when there is no subdirectories, 'TRUE' if there are subdirectories"> + <option value="FALSE"> FALSE </option> + <option value="TRUE"> TRUE </option> + </param> + + <param name="subdirs_names" label="Use subdirectories names as FID names?" type="select" help="Select 'TRUE' to use the subdirectories names as the FID names"> + <option value="FALSE"> FALSE </option> + <option value="TRUE"> TRUE </option> + </param> + </inputs> + + + <outputs> + <data format="tabular" name="dataMatrixOut" label="${tool.name}_dataMatrixOut" /> + <data format="tabular" name="sampleOut" label="${tool.name}_sampleMetadataOut" /> + <data format="txt" name="logOut" label="${tool.name}_log" /> + <data format="pdf" name="graphOut" label="${tool.name}_graph" /> + </outputs> + + <help> + +.. class:: infomark + +**Authors** Manon Martin (manon.martin@uclouvain.be) and Marie Tremblay-Franco (marie.tremblay-franco@inra.fr; Galaxy integration) + +.. class:: infomark + +**Please cite** +Rousseau, R. (2011). Statistical contribution to the analysis of metabonomics data in 1H NMR spectroscopy (Doctoral dissertation, PhD thesis. +Institut de statistique, biostatistique et sciences actuarielles, Université catholique de Louvain, Belgium). + + +============= +NMR Read +============= + +----------- +Description +----------- + +Nuclear Magnetic Resonance files reading (from the R 'SOAP' package) + +----------------- +Workflow position +----------------- + +**Upstream tools** + +========================= ================= ======= +Name output file format +========================= ================= ======= +NA NA NA +========================= ================= ======= + + +**Downstream tools** + ++---------------------------+----------------------+--------+ +| Name | Output file | Format | ++===========================+======================+========+ +|NMR_PreProcessing | dataMatrix.tsv | Tabular| ++---------------------------+----------------------+--------+ +| | sampleMetadata.tsv | Tabular| ++---------------------------+----------------------+--------+ +| | variableMetadata.tsv | Tabular| ++---------------------------+----------------------+--------+ +|NMR_Alignement | dataMatrix.tsv | Tabular| ++---------------------------+----------------------+--------+ +|NMR_Bucketing | dataMatrix.tsv | Tabular| ++---------------------------+----------------------+--------+ +|Normalization | dataMatrix.tsv | Tabular| ++---------------------------+----------------------+--------+ +|Univariate | variableMetadata.tsv | Tabular| ++---------------------------+----------------------+--------+ +|Multivariate | sampleMetadata.tsv | Tabular| ++---------------------------+----------------------+--------+ +| | variableMetadata.tsv | Tabular| ++---------------------------+----------------------+--------+ + + +----------- +Input files +----------- + ++---------------------------+-----------------+ +| Parameter : num + label | Format | ++===========================+=================+ +| 1 : Choose your inputs | FID | ++---------------------------+-----------------+ +| 1 : Choose your inputs | Jcamp-dx | ++---------------------------+-----------------+ + +**Choose your inputs** + + | Zip file (recommended) of FID Bruker files: you can put a zip file containing your FID Bruker files: myinputs.zip (containing all your conditions as sub-directories). + | Zip file (recommended) of Jcamp-dx files: you can put a zip file containing your Jcamp-dx files: myinputs.zip (containing all your conditions as sub-directories). + + +---------- +Parameters +---------- + +FID Title + | Line in the acqus file to find the FID title (name) + | + +subdirectories + | Organization of individual's files + | TRUE: will search inside subdirectories for FIDs and will merge them to have unique FID and info matrices. + | + +------------ +Output files +------------ + +dataMatrixOut.tsv + | tabular output + | Data matrix with n rows (samples) and p columns (time) containing the raw FIDs + | + +sampleOut.pdf + | tabular output + | Data matrix with n rows (samples) containing the acquisition parameters + +logOut.tsv + | Text output + | Contains warnings + | +--------------------------------------------------- + +--------------- +Example +--------------- + + + </help> + +</tool> |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/ptw.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nmr_preprocessing/ptw.R Fri May 05 07:13:02 2017 -0400 |
| [ |
| @@ -0,0 +1,98 @@ +ptw <- function (ref, samp, selected.traces, init.coef = c(0, 1, 0), + try = FALSE, warp.type = c("individual", "global"), optim.crit = c("WCC", + "RMS"), mode = c("forward", "backward"), smooth.param = ifelse(try, + 0, 1e+05), trwdth = 20, trwdth.res = trwdth, verbose = FALSE, + ...) +{ + optim.crit <- match.arg(optim.crit) + warp.type <- match.arg(warp.type) + mode <- match.arg(mode) + if (is.vector(ref)) + ref <- matrix(ref, nrow = 1) + if (is.vector(samp)) + samp <- matrix(samp, nrow = 1) + if (nrow(ref) > 1 && nrow(ref) != nrow(samp)) + stop("The number of references does not equal the number of samples") + if (length(dim(ref)) > 2) + stop("Reference cannot be an array") + if (length(dim(samp)) > 2) + stop("Sample cannot be an array") + if (nrow(samp) == 1) + warp.type <- "individual" + r <- nrow(samp) + if (!missing(selected.traces)) { + samp <- samp[selected.traces, , drop = FALSE] + if (nrow(ref) > 1) + ref <- ref[selected.traces, , drop = FALSE] + } + if (is.vector(init.coef)) + init.coef <- matrix(init.coef, nrow = 1) + if (warp.type == "global") { + if (nrow(init.coef) != 1) + stop("Only one warping function is allowed with global alignment.") + } + else { + if (nrow(init.coef) != nrow(samp)) + if (nrow(init.coef) == 1) { + init.coef <- matrix(init.coef, byrow = TRUE, + nrow = nrow(samp), ncol = length(init.coef)) + } + else { + stop("The number of warping functions does not match the number of samples") + } + } + if (warp.type == "individual") { + w <- matrix(0, nrow(samp), ncol(ref)) + a <- matrix(0, nrow(samp), ncol(init.coef)) + v <- rep(0, nrow(samp)) + warped.sample <- matrix(NA, nrow = nrow(samp), ncol = ncol(samp)) + for (i in 1:nrow(samp)) { + if (verbose & nrow(samp) > 1) + cat(ifelse(nrow(ref) == 1, paste("Warping sample", + i, "with the reference \n"), paste("Warping sample", + i, "with reference \n", i))) + if (nrow(ref) == 1) { + rfrnc <- ref + } + else { + rfrnc <- ref[i, , drop = FALSE] + } + quad.res <- pmwarp(rfrnc, samp[i, , drop = FALSE], + optim.crit, init.coef[i, ], try = try, mode = mode, + smooth.param = smooth.param, trwdth = trwdth, + trwdth.res = trwdth.res, ...) + w[i, ] <- quad.res$w + a[i, ] <- quad.res$a + v[i] <- quad.res$v + warped.sample[i, ] <- c(warp.sample(samp[i, , drop = FALSE], + w[i, ], mode = mode)) + } + } + else { + if (nrow(ref) == 1) + ref <- matrix(ref, nrow = nrow(samp), ncol = ncol(ref), + byrow = TRUE) + if (verbose) { + if (nrow(ref) == 1) { + cat("Simultaneous warping of samples with reference... \n") + } + else { + cat("Simultaneous warping of samples with references... \n") + } + } + quad.res <- pmwarp(ref, samp, optim.crit, c(init.coef), + try = try, mode = mode, smooth.param = smooth.param, + trwdth = trwdth, trwdth.res = trwdth.res, ...) + w <- t(as.matrix(quad.res$w)) + a <- t(as.matrix(quad.res$a)) + v <- quad.res$v + warped.sample <- t(warp.sample(samp, w, mode)) + } + if (verbose) + cat("\nFinished.\n") + result <- list(reference = ref, sample = samp, warped.sample = warped.sample, + warp.coef = a, warp.fun = w, crit.value = v, optim.crit = optim.crit, + mode = mode, warp.type = warp.type) + class(result) <- "ptw" + result +} \ No newline at end of file |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/static/images/NmrPreprocessing.png |
| b |
| Binary file nmr_preprocessing/static/images/NmrPreprocessing.png has changed |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/test-data/MTBLS1.zip |
| b |
| Binary file nmr_preprocessing/test-data/MTBLS1.zip has changed |
| b |
| diff -r 000000000000 -r 68e2d63bece0 nmr_preprocessing/test-data/MTBLS1_alignedSpectra.tabular --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nmr_preprocessing/test-data/MTBLS1_alignedSpectra.tabular Fri May 05 07:13:02 2017 -0400 |
| b |
| b'@@ -0,0 +1,29802 @@\n+ADG10003u_007\tADG10003u_009\tADG10003u_015\tADG10003u_017\tADG10003u_021\tADG10003u_023\tADG10003u_051\tADG10003u_053\tADG10003u_067\tADG10003u_071\tADG10003u_073\tADG10003u_087\tADG10003u_089\tADG10003u_097\r\n+9.2997\t-0.00172384578396177\t-0.000133260203713335\t0.00188546634303473\t0.00145955595536551\t0.00266060722603239\t-0.00456406825398832\t-0.00131195747929342\t0.000299817194153485\t0.000241533956195763\t-0.00349445210212112\t0.00449073847864704\t-0.00150200410788306\t-0.00479292634495511\t-0.00845701750085149\r\n+9.2994\t-0.002312062118707\t-0.000821498405692276\t0.00530663314770368\t0.0031937690024524\t-0.00295334219789086\t-0.00562040636826083\t0.000544476077311362\t0.000328645770514397\t0.00269356222036034\t-3.40676966964021e-05\t0.00203230500493516\t-0.000176706365633302\t-0.00197312886982787\t-0.00571063951344053\r\n+9.2991\t0.00293819295016693\t0.000917760076187904\t0.00843046810963003\t-0.000269595302370741\t0.00164391778956326\t-0.00116020450231351\t0.00180051068287967\t0.00190620953248653\t0.0037736553777559\t0.0093018774151955\t-0.000394479670034688\t-0.003086471186395\t0.00198994473348899\t-0.00300222608200099\r\n+9.2988\t-0.00452049127972714\t-0.00104938317666152\t0.00739844285463343\t-0.00723747846021565\t0.00239897800809775\t0.00247439274074376\t0.00110713566571126\t0.00128511387066643\t0.00345596457000934\t0.00696433078367433\t0.00224254345372155\t-0.00566245731651603\t0.000398110507158\t-0.00679345467506024\r\n+9.2985\t-0.00679644435394394\t0.00108834623376689\t0.00259733552669646\t-0.000234823010309861\t-0.00287387338820875\t0.00360433314397773\t0.00131195747929342\t0.000999710964693405\t0.000689024945235094\t0.00650078671387083\t0.0031892946618048\t-0.00307469076201945\t-0.00265731165999202\t-0.0097797246393341\r\n+9.2982\t-0.00277642443085087\t0.00566634173332424\t-0.00691481595689443\t-0.0028308607137917\t-0.00725637584172303\t-0.0016260295153459\t0.00171405323393041\t0.000538453742918813\t0.00428868693222166\t0.0079838484938265\t0.00181890167627456\t0.00292768088283286\t-0.00247963078733168\t-0.00681381248043622\r\n+9.2979\t0.0018982524688494\t0.00520336893713099\t-0.00477951684143076\t-0.00662236101521466\t-0.00714026189728136\t-0.00179792959811773\t0.00189382983412648\t-0.00159037646257698\t0.00846508420011296\t0.00553739693024306\t0.000491234569046061\t0.0065864843534735\t-0.000949387194169153\t-0.000196700416808458\r\n+9.2976\t0.00679969416794806\t0.00526852564607191\t0.000441679666812651\t-0.00120910742140804\t-0.00213479215285431\t-0.000118151042806551\t0.000152672876120703\t-0.00308433735310247\t0.00592967252184167\t-0.000913405212245339\t0.00251223643157335\t0.00205224809642454\t-0.00322256779992628\t0.00248695352160624\r\n+9.2973\t0.0142226304437968\t0.00623769076398707\t0.000127692305908127\t0.00213321409459558\t-0.00336112584864566\t0.000102171598492551\t-0.00158299472004702\t-0.00277298872840462\t0.00392207785486584\t0.00206109565013233\t0.00669824219090419\t-0.00475364666104364\t-0.000766236100799773\t0.00362478977883671\r\n+9.297\t0.0095699800612356\t0.0038963057105373\t0.00351092154143045\t-0.00115056672717896\t-0.00313273291572552\t0.00197273503440119\t0.00124299737120294\t-0.000204042257132233\t0.00372445916343294\t-0.00169081888112061\t0.00880017455307183\t-0.00277600791899761\t0.00325660472878254\t-0.00132985988091201\r\n+9.2967\t0.00329603358106473\t0.00164692585915988\t-0.00187004456695887\t0.00352938764417924\t-0.00634216495071527\t0.00176403380714863\t0.00254534846727977\t-0.000827380141558176\t0.00593828880796603\t0.000819579424622461\t0.00804286401174449\t-0.00315985341323439\t0.00436989594345592\t-0.00604406735284172\r\n+9.2964\t0.00459126500692792\t0.00211219059988875\t-0.0050096097404826\t0.00286563300585258\t-0.00620346737671246\t0.00252983657025749\t0.000727340443541328\t-0.00440212361031127\t0.00549691260262785\t0.00261260139583228\t0.00773247972122436\t0.0021570447882654\t0.000121155258666913\t-0.000861080146308356\r\n+9.2961\t-0.000111938037919595\t0.00239344493647291\t-0.00924072822465575\t0.00409058602870609\t-0.00303387627311829\t-0.00573516783197047\t0.00190686706849184\t-0.00350107044027521\t0.00336757815105623\t0.00231464867144653\t0.0069998100303461'..b'77632683724\t-0.00376535913578589\t-0.00189051397424559\t-0.00680797781709899\t0.00220797768670862\t0.000364413311229752\t0.00462651624799079\t0.00499289330198562\t-0.0113395726998991\r\n+0.2037\t-0.00173576176864353\t0.000437761406301544\t0.000580167215973881\t0.00201899371750953\t0.00560137929048828\t-0.0079577632683724\t-0.00376535913578589\t-0.00127838720284889\t-0.00680797781709899\t-0.00129066306664566\t-0.00171378245188136\t-0.000434157723340709\t0.00830663144706259\t-0.00234747504423297\r\n+0.2034\t-0.00247527499980266\t0.000416806484833108\t-0.00300200293092705\t0.00609923608971649\t0.00340480173600446\t-0.0079577632683724\t-0.00376535913578589\t0.00338543581731644\t-0.00680797781709899\t-7.53957221969555e-06\t-0.00361949226184836\t0.00424659756187919\t0.00552958273377266\t0.00552741926505671\r\n+0.2031\t0.00054344111957739\t-0.00068660109873922\t0.00441371231291134\t0.00608933258881308\t0.00585129058742422\t-0.0079577632683724\t-0.00376535913578589\t0.00275184777240662\t-0.00680797781709899\t0.00633435763820644\t-0.00334233921000737\t0.00861222649505304\t0.00292555507550204\t0.00768149515821787\r\n+0.2028\t8.88282494458724e-05\t-0.00183519273172786\t0.0110009697459545\t0.00410885248592794\t0.0110917579111525\t-0.0079577632683724\t-0.0011345824749015\t0.00165508060063147\t-0.00738054392729838\t0.00380245758946646\t-0.00781250596947208\t0.00850792065422783\t0.00171359728729891\t0.000210180585233094\r\n+0.2025\t0.00254316000344422\t-0.0030761169874368\t0.00813837967075322\t0.00410885248592794\t0.0103098530008755\t-0.0079577632683724\t0.00219162771384057\t0.00228354356529936\t-0.00133941557526741\t-0.00210186518880179\t-0.00742752093906001\t0.00598911866742981\t0.0043536878820958\t0.004245152631848\r\n+0.2022\t0.00234564353008287\t-0.0045118565286726\t-0.00333850608490233\t0.00410885248592794\t0.0111906145537597\t-0.0079577632683724\t0.00279820219694485\t0.00318267483024469\t-0.0019044771781972\t-0.00718297837804551\t-0.00844050856165338\t0.00561705359756858\t0.00229810050009832\t0.00193399151071813\r\n+0.2019\t-0.000822564033486575\t-0.0040515030976629\t-0.00796658108526813\t0.00410885248592794\t0.0111910406599779\t-0.0038357929737381\t0.00141728460956096\t1.28127006048498e-05\t0.000142307693408781\t-0.0043810499464742\t-0.00639260522019536\t0.00376237303494238\t0.00180395722938159\t-0.00401076276454617\r\n+0.2016\t-0.00322453767297415\t-0.00410487266327783\t-0.00458643620496098\t0.00410885248592794\t0.011761383832951\t-0.000322736352584289\t0.00112291758099565\t-0.000356513394329946\t0.00264936901088379\t-0.00804193112425971\t-0.00513162665317331\t0.00391821823241064\t-0.00198893172965399\t-0.00319700076046307\r\n+0.2013\t-0.00509715272023551\t0.00140921846875231\t0.00132164620970126\t-0.00839068612096886\t0.0126843299014305\t0.000420792033602021\t0.0003677872431492\t0.000191870191557626\t0.00253318812056178\t-0.00941497099849093\t-0.00805778415972288\t0.00426009596480952\t-0.00233902585503262\t0.000656676776114392\r\n+0.201\t0.000761178657853248\t0.00155000934736837\t0.00579704562691605\t-0.00447396155256717\t0.00786613383987363\t-0.00118828958625933\t0.00191750270705306\t-0.000602517245943062\t0.00493796578339922\t-0.00887212179867285\t-0.00843078785918261\t-0.0009743392660614\t-0.000338140680125548\t0.00136672401497122\r\n+0.2007\t0.00609592889083422\t0.00516571556261739\t0.00594756216141645\t-0.00749474940590659\t0.00803891991132721\t0.00248480359082712\t-0.0030404202880487\t-0.00130433292157371\t0.00596802889233075\t-0.00572141834108675\t-0.00758938193834072\t-0.00445790892411569\t0.00262611114187379\t0.00898632044062698\r\n+0.2004\t0.00381564273127859\t0.00426530878077054\t0.00435325895069397\t-0.00396602199512091\t0.00999751714298286\t0.00384741438778465\t-0.00201699739036767\t-0.00225087117875699\t0.00396905051147928\t-0.00392727939621475\t-0.00596308580869628\t-0.00598543728481245\t0.00439076382245708\t0.000253097039809485\r\n+0.2001\t0.00347693989396059\t0.00341008604834\t0.0123071941795799\t0.00155969135338893\t0.00947255428224105\t0.00357649017282454\t0.0018327606836782\t-0.0032220738846046\t0.00134108324354955\t0.00179022953705438\t-0.00319313773022348\t-0.00350001316707849\t0.0046579942341321\t-0.00175132147329041\r\n' |
