Mercurial > repos > marie-tremblay-metatoul > nmr_preprocessing
changeset 1:cbea5e9fd0b4 draft
Uploaded
| author | marie-tremblay-metatoul |
|---|---|
| date | Tue, 23 May 2017 03:32:51 -0400 |
| parents | 68e2d63bece0 |
| children | 5e64657b4fe5 |
| files | nmr_preprocessing/NmrPreprocessing_script.R nmr_preprocessing/NmrPreprocessing_wrapper.R nmr_preprocessing/NmrPreprocessing_xml.xml |
| diffstat | 3 files changed, 255 insertions(+), 203 deletions(-) [+] |
line wrap: on
line diff
--- a/nmr_preprocessing/NmrPreprocessing_script.R Fri May 05 07:13:02 2017 -0400 +++ b/nmr_preprocessing/NmrPreprocessing_script.R Tue May 23 03:32:51 2017 -0400 @@ -550,7 +550,7 @@ # InternalReferencing ## ==================================================== -InternalReferencing <- function(RawSpect_data, RawSpect_info, method = c("max", "thres"), +InternalReferencing <- function(Spectrum_data, Fid_info, method = c("max", "thres"), range = c("near0", "all", "window"), ppm.ref = 0, shiftHandling = c("zerofilling", "cut", "NAfilling", "circular"), c = 2, pc = 0.02, fromto.RC = NULL, @@ -560,9 +560,9 @@ # Data initialisation and checks ---------------------------------------------- - begin_info <- beginTreatment("InternalReferencing", RawSpect_data, RawSpect_info) - RawSpect_data <- begin_info[["Signal_data"]] - RawSpect_info <- begin_info[["Signal_info"]] + begin_info <- beginTreatment("InternalReferencing", Spectrum_data, Fid_info) + Spectrum_data <- begin_info[["Signal_data"]] + Fid_info <- begin_info[["Signal_info"]] # Check input arguments @@ -582,9 +582,9 @@ if (!is.null(fromto.RC)) { diff <- diff(unlist(fromto.RC))[1:length(diff(unlist(fromto.RC)))%%2 !=0] for (i in 1:length(diff)) { - if (diff[i] <= 0) { - stop(paste("Invalid region removal because from > to")) - } + if (ppm == TRUE & diff[i] >= 0) { + stop(paste("Invalid region removal because from <= to in ppm")) + } else if (ppm == FALSE & diff[i] <= 0) {stop(paste("Invalid region removal because from >= to in column index"))} } } @@ -627,15 +627,15 @@ # Apply the method ('thres' or 'max') on spectra # ---------------------------------------------- - n <- nrow(RawSpect_data) - m <- ncol(RawSpect_data) + n <- nrow(Spectrum_data) + m <- ncol(Spectrum_data) # The Sweep Width has to be the same since the column names are the same - SW <- RawSpect_info[1, "SW"] # Sweep Width in ppm (semi frequency scale in ppm) + SW <- Fid_info[1, "SW"] # Sweep Width in ppm (semi frequency scale in ppm) ppmInterval <- SW/m # FIXME divide by two ?? if (range == "all") { - Data <- RawSpect_data + Data <- Spectrum_data } else { if (range == "near0") { fromto.RC <- list(c(-(SW * pc)/2, (SW * pc)/2)) # automatic fromto values in ppm @@ -644,7 +644,7 @@ # if ppm == TRUE, then fromto is in the colnames values, else, in the column # index if (ppm == TRUE) { - colindex <- as.numeric(colnames(RawSpect_data)) + colindex <- as.numeric(colnames(Spectrum_data)) } else { colindex <- 1:m } @@ -659,9 +659,9 @@ vector <- rep(0, m) vector[unlist(Int)] <- 1 if (n > 1) { - Data <- sweep(RawSpect_data, MARGIN = 2, FUN = "*", vector) # Cropped_Spectrum + Data <- sweep(Spectrum_data, MARGIN = 2, FUN = "*", vector) # Cropped_Spectrum } else { - Data <- RawSpect_data * vector + Data <- Spectrum_data * vector } # Cropped_Spectrum } @@ -701,17 +701,17 @@ ppmScale <- ppmScale + ppm.ref - Spectrum_data <- matrix(fill, nrow = n, ncol = -(end - start) + 1, - dimnames = list(rownames(RawSpect_data), ppmScale)) + Spectrum_data_calib <- matrix(fill, nrow = n, ncol = -(end - start) + 1, + dimnames = list(rownames(Spectrum_data), ppmScale)) for (i in 1:n) { shift <- (1 - TMSPpeaks[i]) + start - Spectrum_data[i, (1 + shift):(m + shift)] <- RawSpect_data[i, ] + Spectrum_data_calib[i, (1 + shift):(m + shift)] <- Spectrum_data[i, ] } if (shiftHandling == "cut") { - Spectrum_data = as.matrix(stats::na.omit(t(Spectrum_data))) - Spectrum_data = t(Spectrum_data) - base::attr(Spectrum_data, "na.action") <- NULL + Spectrum_data_calib = as.matrix(stats::na.omit(t(Spectrum_data_calib))) + Spectrum_data_calib = t(Spectrum_data_calib) + base::attr(Spectrum_data_calib, "na.action") <- NULL } @@ -730,13 +730,13 @@ } ppmScale <- ppmScale + ppm.ref - Spectrum_data <- matrix(nrow=n, ncol=end-start+1, - dimnames=list(rownames(RawSpect_data), ppmScale)) + Spectrum_data_calib <- matrix(nrow=n, ncol=end-start+1, + dimnames=list(rownames(Spectrum_data), ppmScale)) for (i in 1:n) { shift <- (maxpeak-TMSPpeaks[i]) - Spectrum_data[i,(1+shift):m] <- RawSpect_data[i,1:(m-shift)] + Spectrum_data_calib[i,(1+shift):m] <- Spectrum_data[i,1:(m-shift)] if (shift > 0) { - Spectrum_data[i,1:shift] <- RawSpect_data[i,(m-shift+1):m] + Spectrum_data_calib[i,1:shift] <- Spectrum_data[i,(m-shift+1):m] } } } @@ -757,9 +757,9 @@ fromto <- fromto.RC } else { fromto <- list() - idcol <- as.numeric(colnames(RawSpect_data)) + idcol <- as.numeric(colnames(Spectrum_data)) for (i in 1:length(fromto.RC)) { - fromto[[i]] <- as.numeric(colnames(RawSpect_data))[fromto.RC[[i]]] + fromto[[i]] <- as.numeric(colnames(Spectrum_data))[fromto.RC[[i]]] } } } else { @@ -767,7 +767,7 @@ } # TMSPloc in ppm - TMSPloc <- as.numeric(colnames(RawSpect_data))[TMSPpeaks[rowindex_graph]] + TMSPloc <- as.numeric(colnames(Spectrum_data))[TMSPpeaks[rowindex_graph]] # num plot per window num.stacked <- 6 @@ -778,7 +778,7 @@ Legend = "TMSP search zone and location") # vlines for TMSP peak - addlines <- data.frame(rowname = rownames(RawSpect_data)[rowindex_graph],TMSPloc) + addlines <- data.frame(rowname = rownames(Spectrum_data)[rowindex_graph],TMSPloc) nn <- length(rowindex_graph) i <- 1 @@ -789,7 +789,7 @@ last <- min(i + num.stacked - 1, nn) - melted <- reshape2::melt(Re(RawSpect_data[i:last, ]), + melted <- reshape2::melt(Re(Spectrum_data[i:last, ]), varnames = c("rowname", "ppm")) plots[[j]] <- ggplot2::ggplot() + ggplot2::theme_bw() + @@ -819,7 +819,7 @@ # Return the results ---------------------------------------------- - Spectrum_data <- endTreatment("InternalReferencing", begin_info, Spectrum_data) + Spectrum_data <- endTreatment("InternalReferencing", begin_info, Spectrum_data_calib) if (is.null(plots)) { return(Spectrum_data) @@ -837,8 +837,8 @@ ZeroOrderPhaseCorrection <- function(Spectrum_data, method = c("rms", "manual", "max"), plot_rms = NULL, returnAngle = FALSE, createWindow = TRUE, - angle = NULL, plot_spectra = FALSE, quant = 0.95, - ppm = TRUE, fromto.0OPC = NULL) { + angle = NULL, plot_spectra = FALSE, + ppm = TRUE, exclude = list(c(5.1,4.5))) { # Data initialisation and checks ---------------------------------------------- @@ -848,7 +848,7 @@ # optimal angles createWindow : for plot_rms plots angle : If angle is not NULL, # spectra are rotated according to the angle vector values # plot_spectra : if TRUE, plot rotated spectra - # quant: probability for sample quantile used to trim the spectral intensities + begin_info <- beginTreatment("ZeroOrderPhaseCorrection", Spectrum_data) @@ -861,7 +861,7 @@ # Check input arguments method <- match.arg(method) checkArg(ppm, c("bool")) - checkArg(unlist(fromto.0OPC), c("num"), can.be.null = TRUE) + checkArg(unlist(exclude), c("num"), can.be.null = TRUE) # method in c("max", "rms") ----------------------------------------- @@ -869,39 +869,25 @@ # angle is found by optimization # rms function to be optimised - rms <- function(ang, y, meth = c("max", "rms"), p = 0.95) { + rms <- function(ang, y, meth = c("max", "rms")) { # if (debug_plot) { graphics::abline(v=ang, col='gray60') } roty <- y * exp(complex(real = 0, imaginary = ang)) # spectrum rotation Rey <- Re(roty) if (meth == "rms") { - si <- sign(Rey) # sign of intensities - - Rey[abs(Rey) >= quantile(abs(Rey), p)] <- quantile(abs(Rey), p) # trim the values - Rey <- abs(Rey) * si # spectral trimmed values - - if ((sum(Rey==0)/length(Rey)) >= 0.99) { - stop("More than 99% of intensities are null, the rms criterion cannot work properly. \n - Either increase p or the window(s) in fromto.0OPC") - } ReyPos <- Rey[Rey >= 0] # select positive intensities - - # POSss = sum((ReyPos-mean(ReyPos))^2) # centred SS for positive intensities - POSss <- sum((ReyPos)^2) # SS for positive intensities - - # ss = sum((Rey - mean(Rey) )^2) # centred SS for all intensities - ss <- sum((Rey)^2) # SS for all intensities - + POSss <- sum((ReyPos)^2, na.rm = TRUE) # SS for positive intensities + ss <- sum((Rey)^2, na.rm = TRUE) # SS for all intensities return(POSss/ss) # criterion : SS for positive values / SS for all intensities - } else { - maxi <- max(Rey) - return(maxi) - } + } else { + maxi <- max(Rey, na.rm = TRUE) + return(maxi) + } } # Define the interval where to search for (by defining Data) - if (is.null(fromto.0OPC)) { + if (is.null(exclude)) { Data <- Spectrum_data } else { @@ -913,22 +899,23 @@ colindex <- 1:m } - # Second check for the argument fromto.0OPC - diff <- diff(unlist(fromto.0OPC))[1:length(diff(unlist(fromto.0OPC)))%%2 != 0] - for (i in 1:length(diff)) { - if (diff[i] < 0) { - stop(paste("Invalid region removal because from > to")) - } + # Second check for the argument exclude + diff <- diff(unlist(exclude))[1:length(diff(unlist(exclude)))%%2 !=0] + for (i in 1:length(diff)) { + if (ppm == TRUE & diff[i] >= 0) { + stop(paste("Invalid region removal because from <= to in ppm")) + } else if (ppm == FALSE & diff[i] <= 0) {stop(paste("Invalid region removal because from >= to in column index"))} } - Int <- vector("list", length(fromto.0OPC)) - for (i in 1:length(fromto.0OPC)) { - Int[[i]] <- indexInterval(colindex, from = fromto.0OPC[[i]][1], - to = fromto.0OPC[[i]][2], inclusive = TRUE) + + Int <- vector("list", length(exclude)) + for (i in 1:length(exclude)) { + Int[[i]] <- indexInterval(colindex, from = exclude[[i]][1], + to = exclude[[i]][2], inclusive = TRUE) } - vector <- rep(0, m) - vector[unlist(Int)] <- 1 + vector <- rep(1, m) + vector[unlist(Int)] <- 0 if (n > 1) { Data <- sweep(Spectrum_data, MARGIN = 2, FUN = "*", vector) # Cropped_Spectrum } else { @@ -953,8 +940,8 @@ # work in either [-pi;pi] or [0;2pi] (this is not easy to be convinced by that I # agree) and we can check which one it is simply by the following trick - f0 <- rms(0, Data[k, ], p = quant, meth = method) - fpi <- rms(pi, Data[k, ], p = quant, meth = method) + f0 <- rms(0, Data[k, ],meth = method) + fpi <- rms(pi, Data[k, ], meth = method) if (f0 < fpi) { interval <- c(-pi, pi) } else { @@ -967,7 +954,7 @@ x <- seq(min(interval), max(interval), length.out = 100) y <- rep(1, 100) for (K in (1:100)) { - y[K] <- rms(x[K], Data[k, ], p = quant, meth = method) + y[K] <- rms(x[K], Data[k, ], meth = method) } if (createWindow == TRUE) { grDevices::dev.new(noRStudioGD = FALSE) @@ -978,7 +965,7 @@ # Best angle best <- stats::optimize(rms, interval = interval, maximum = TRUE, - y = Data[k,], p = quant, meth = method) + y = Data[k,], meth = method) ang <- best[["maximum"]] @@ -1011,32 +998,12 @@ stop(paste("angle has length", length(angle), "and there are", n, "spectra to rotate.")) } for (k in 1:n) { - Spectrum_data[k, ] <- Spectrum_data[k, ] * exp(complex(real = 0, imaginary = ang)) + Spectrum_data[k, ] <- Spectrum_data[k, ] * exp(complex(real = 0, imaginary = - angle[k])) } } - # #================== Detect a 180° rotation due to the water signal MEAN_Q = c() - # for (i in 1:nrow(Spectrum_data)) { data = Re(Spectrum_data[i,]) data_p = - # data[data >= stats::quantile(data[data >=0 ], p.zo)] data_n = data[data <= - # stats::quantile(data[data <0 ], (1-p.zo))] mean_quant = (sum(data_p) + - # sum(data_n)) / (length(data_p) +length(data_n)) # mean(p.zo% higher pos and neg - # values) MEAN_Q = c(MEAN_Q, mean_quant) } vect = which(MEAN_Q < 0) if - # (length(vect)!=0) { warning('The mean of', p.zo,' positive and negative - # quantiles is negative for ', paste0(rownames(Spectrum_data)[vect],'; ')) - # if(rotation == TRUE) { warning(' An automatic 180 degree rotation is applied to - # these spectra') Angle[vect] = Angle[vect] + pi } } vect_risk = - # which(MEAN_Q<0.1*mean(MEAN_Q[MEAN_Q>0])) # is there any MEAN_Q with a very low - # value copared to mean of positive mean values? if (length(vect_risk)!=0) - # { warning('the rotation angle for spectra', - # paste0(rownames(Spectrum_data)[vect_risk],'; '), 'might not be optimal, you - # need to check visually for those spectra') } # result of automatic rotation for - # (k in vect_risk) { Spectrum_data[k,] <- Spectrum_data[k,] * exp(complex(real=0, - # imaginary=Angle[k])) } #================== - - - # Draw spectra if (plot_spectra == TRUE) { nn <- ceiling(n/4) @@ -1070,6 +1037,7 @@ } + ## ==================================================== # Baseline Correction ## ====================================================
--- a/nmr_preprocessing/NmrPreprocessing_wrapper.R Fri May 05 07:13:02 2017 -0400 +++ b/nmr_preprocessing/NmrPreprocessing_wrapper.R Tue May 23 03:32:51 2017 -0400 @@ -70,11 +70,22 @@ ##====================================================== ##====================================================== +# graphical inputs +FirstOPCGraph <- argLs[["FirstOPCGraph"]] +SSGraph <- argLs[["SSGraph"]] +ApodGraph <- argLs[["ApodGraph"]] +FTGraph <- argLs[["FTGraph"]] +SRGraph <- argLs[["SRGraph"]] +ZeroOPCGraph <- argLs[["ZeroOPCGraph"]] +BCGraph <- argLs[["BCGraph"]] +FinalGraph <- argLs[["FinalGraph"]] + + # 1rst order phase correction ------------------------ # Inputs ## Data matrix Fid_data0 <- read.table(argLs[["dataMatrixFid"]],header=TRUE, check.names=FALSE, sep='\t') -Fid_data0 <- Fid_data0[,-1] +# Fid_data0 <- Fid_data0[,-1] Fid_data0 <- as.matrix(Fid_data0) ## Samplemetadata @@ -120,7 +131,7 @@ # Fourier transform ---------------------------------- # Inputs -FTGraph <- argLs[["FTGraph"]] + # Internal referencering ---------------------------------- # Inputs @@ -138,11 +149,11 @@ # if (shiftReferencing=="YES") # { # - shiftReferencingMethod <- argLs[["shiftReferencingMethod"]] - - if (shiftReferencingMethod == "thres") { - shiftTreshold <- argLs[["shiftTreshold"]] - } + # shiftReferencingMethod <- argLs[["shiftReferencingMethod"]] + # + # if (shiftReferencingMethod == "thres") { + # shiftTreshold <- argLs[["shiftTreshold"]] + # } shiftReferencingRange <- argLs[["shiftReferencingRange"]] @@ -151,15 +162,18 @@ } if (shiftReferencingRange == "window"){ - shiftReferencingRangeList <- NULL - shiftReferencingRangeLeft <- argLs[["shiftReferencingRangeLeft"]] - shiftReferencingRangeRight <- argLs[["shiftReferencingRangeRight"]] - shiftReferencingRangeList <- list(shiftReferencingRangeList,c(shiftReferencingRangeLeft,shiftReferencingRangeRight)) + shiftReferencingRangeList <- list() + for(i in which(names(argLs)=="shiftReferencingRangeLeft")) + { + shiftReferencingRangeLeft <- argLs[[i]] + shiftReferencingRangeRight <- argLs[[i+1]] + shiftReferencingRangeList <- c(shiftReferencingRangeList,list(c(shiftReferencingRangeLeft,shiftReferencingRangeRight))) + } } shiftHandling <- argLs[["shiftHandling"]] - ppm_ref <- argLs[["ppm_ref"]] + # } @@ -171,27 +185,24 @@ createWindow = TRUE angle = NULL plot_spectra = FALSE -quant = 0.95 ppm = TRUE -fromto.0OPC = NULL +exclude = NULL zeroOrderPhaseMethod <- argLs[["zeroOrderPhaseMethod"]] -if (zeroOrderPhaseMethod=='rms') -{ - quant <- argLs[["quant"]] -} -if (zeroOrderPhaseMethod=='max') -{ + +if (zeroOrderPhaseMethod=='manual'){ angle <- argLs[["angle"]] } -searchZoneZeroPhase <- argLs[["searchZoneZeroPhase.choice"]] -if (searchZoneZeroPhase == "YES") { - searchZoneZeroPhaseList <- NULL - searchZoneZeroPhaseLeft <- argLs[["shiftReferencingRangeLeft"]] - searchZoneZeroPhaseRight <- argLs[["shiftReferencingRangeRight"]] - searchZoneZeroPhaseList <- list(searchZoneZeroPhaseList,c(searchZoneZeroPhaseLeft,searchZoneZeroPhaseRight)) +excludeZoneZeroPhase <- argLs[["excludeZoneZeroPhase.choice"]] +if (excludeZoneZeroPhase == 'YES') { + excludeZoneZeroPhaseList <- list() + for(i in which(names(argLs)=="excludeZoneZeroPhase_left")) { + excludeZoneZeroPhaseLeft <- argLs[[i]] + excludeZoneZeroPhaseRight <- argLs[[i+1]] + excludeZoneZeroPhaseList <- c(excludeZoneZeroPhaseList,list(c(excludeZoneZeroPhaseLeft,excludeZoneZeroPhaseRight))) } - + exclude <- excludeZoneZeroPhaseList +} # Baseline Correction ------------------------------- @@ -235,20 +246,45 @@ # FirstOrderPhaseCorrection --------------------------------- Fid_data <- FirstOrderPhaseCorrection(Fid_data0, Fid_info = samplemetadataFid, group_delay = NULL) +if (FirstOPCGraph == "YES") { + title = "FIDs after First Order Phase Correction" + DrawSignal(Fid_data, subtype = "stacked", + ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, + xlab = "Frequency", num.stacked = 4, + main = title, createWindow=FALSE) +} # SolventSuppression --------------------------------- Fid_data <- SolventSuppression(Fid_data, lambda.ss = lambda, ptw.ss = ptwSS, plotSolvent = F, returnSolvent = F) +if (SSGraph == "YES") { + title = "FIDs after Solvent Suppression " + DrawSignal(Fid_data, subtype = "stacked", + ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, + xlab = "Frequency", num.stacked = 4, + main = title, createWindow=FALSE) +} + + # Apodization --------------------------------- Fid_data <- Apodization(Fid_data, Fid_info = samplemetadataFid, DT = NULL, type.apod = apodization, phase = phase, rectRatio = rectRatio, gaussLB = gaussLB, expLB = expLB, plotWindow = F, returnFactor = F) +if (ApodGraph == "YES") { + title = "FIDs after Apodization" + DrawSignal(Fid_data, subtype = "stacked", + ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, + xlab = "Frequency", num.stacked = 4, + main = title, createWindow=FALSE) +} + + # FourierTransform --------------------------------- Spectrum_data <- FourierTransform(Fid_data, Fid_info = samplemetadataFid, reverse.axis = TRUE) if (FTGraph == "YES") { - title = "Fourier transformed specta" + title = "Fourier transformed spectra" DrawSignal(Spectrum_data, subtype = "stacked", ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, xlab = "Frequency", num.stacked = 4, @@ -257,45 +293,59 @@ # InternalReferencing --------------------------------- # if (shiftReferencing=="YES") { -Spectrum_data <- InternalReferencing(Spectrum_data, samplemetadataFid, method = shiftReferencingMethod, range = shiftReferencingRange, +Spectrum_data <- InternalReferencing(Spectrum_data, samplemetadataFid, method = "max", range = shiftReferencingRange, ppm.ref = 0, shiftHandling = shiftHandling,ppm = TRUE, - c = shiftTreshold, fromto.RC = shiftReferencingRangeList, pc = pctNear0) + c = shiftTreshold, fromto.RC = shiftReferencingRangeList, pc = pctNear0) + +if (SRGraph == "YES") { + title = "Spectra after Shift Referencing" + DrawSignal(Spectrum_data, subtype = "stacked", + ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, + xlab = "Frequency", num.stacked = 4, + main = title, createWindow=FALSE) +} # } # ZeroOrderPhaseCorrection --------------------------------- -Spectrum_data <- ZeroOrderPhaseCorrection(Spectrum_data, method = zeroOrderPhaseMethod, - plot_rms = plot_rms, returnAngle = returnAngle, - createWindow = createWindow,angle = angle, - plot_spectra = plot_spectra, quant = quant, - ppm = ppm, fromto.0OPC = fromto.0OPC) +Spectrum_data <- ZeroOrderPhaseCorrection(Spectrum_data, method = zeroOrderPhaseMethod, + plot_rms = plot_rms, returnAngle = returnAngle, + createWindow = createWindow,angle = angle, + plot_spectra = plot_spectra, + ppm = ppm, exclude = exclude) +if (ZeroOPCGraph == "YES") { title = "Spectra after Zero Order Phase Correction" DrawSignal(Spectrum_data, subtype = "stacked", ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, xlab = "Frequency", num.stacked = 4, main = title, createWindow=FALSE) +} # BaselineCorrection --------------------------------- -Spectrum_data <- BaselineCorrection(Spectrum_data, ptw.bc = TRUE, maxIter = maxIter, lambda.bc = lambdaBc, p.bc = pBc, eps = epsilon, returnBaseline = F) +Spectrum_data <- BaselineCorrection(Spectrum_data, ptw.bc = ptwBc, maxIter = maxIter, lambda.bc = lambdaBc, p.bc = pBc, eps = epsilon, returnBaseline = F) +if (BCGraph == "YES") { title = "Spectra after Baseline Correction" DrawSignal(Spectrum_data, subtype = "stacked", ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, xlab = "Frequency", num.stacked = 4, main = title, createWindow=FALSE) +} + # NegativeValuesZeroing --------------------------------- if (NegativetoZero=="YES") { + Spectrum_data <- NegativeValuesZeroing(Spectrum_data) +} - Spectrum_data <- NegativeValuesZeroing(Spectrum_data) - - title = "Spectra after Negative Values Zeroing" +if (FinalGraph == "YES") { + title = "Final preprocessed spectra" DrawSignal(Spectrum_data, subtype = "stacked", ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T, xlab = "Frequency", num.stacked = 4, main = title, createWindow=FALSE) - } +} invisible(dev.off())
--- a/nmr_preprocessing/NmrPreprocessing_xml.xml Fri May 05 07:13:02 2017 -0400 +++ b/nmr_preprocessing/NmrPreprocessing_xml.xml Tue May 23 03:32:51 2017 -0400 @@ -1,4 +1,4 @@ -<tool id="NMR_Preprocessing" name="NMR_Preprocessing" version="3.0.1"> +<tool id="NMR_Preprocessing" name="NMR_Preprocessing" version="3.1.0"> <description> Preprocessing of 1D NMR spectra </description> <stdio> @@ -10,6 +10,8 @@ Rscript $__tool_directory__/NmrPreprocessing_wrapper.R ## First order phase correction + ## Graphical display + FirstOPCGraph $FirstOPCGraph ## Data matrix of FID spectra dataMatrixFid $dataMatrixFid ## Sample metadata matrix @@ -20,9 +22,14 @@ ## Smoothing parameter lambda $lambda ptwSS $ptwSS + ## Graphical display + SSGraph $SSGraph ## Apodization + ## Graphical display + ApodGraph $ApodGraph + apodizationMethod $apodizationMethod.method #if $apodizationMethod.method == "exp": ## Line broadening for the exponential window @@ -54,18 +61,15 @@ gaussLB $apodizationMethod.gaussLB #end if - ## Fourier transform ## Graphical display FTGraph $FTGraph ## Shift referencing - ## Method used to find the TMSP peaks in spectra - shiftReferencingMethod $shiftReferencingMethod.method - #if $shiftReferencingMethod.method == "thres": - shiftTreshold $shiftReferencingMethod.shiftTreshold - #end if + ## Graphical display + SRGraph $SRGraph + ## Definition of the search zone shiftReferencingRange $shiftReferencingRange.method #if $shiftReferencingRange.method == "near0": @@ -78,24 +82,27 @@ #end for #end if shiftHandling $shiftHandling + ## Zero order phase correction - zeroOrderPhaseMethod $zeroOrderPhaseMethod.method - #if $zeroOrderPhaseMethod.method == "rms": - ## Line broadening for the exponential window - quant $zeroOrderPhaseMethod.quant - #end if - searchZoneZeroPhase.choice ${searchZoneZeroPhase.choice} - #if str($searchZoneZeroPhase.choice) == "YES": - #for $i in $searchZoneZeroPhase.conditions: - searchZoneZeroPhase_left ${i.searchZoneZeroPhase_left} - searchZoneZeroPhase_right ${i.searchZoneZeroPhase_right} + ## Graphical display + ZeroOPCGraph $ZeroOPCGraph + + zeroOrderPhaseMethod $zeroOrderPhaseMethod + + excludeZoneZeroPhase.choice ${excludeZoneZeroPhase.choice} + #if str($excludeZoneZeroPhase.choice) == "YES": + #for $i in $excludeZoneZeroPhase.conditions: + excludeZoneZeroPhase_left ${i.excludeZoneZeroPhase_left} + excludeZoneZeroPhase_right ${i.excludeZoneZeroPhase_right} #end for #end if ## Baseline correction + ## Graphical display + BCGraph $BCGraph ptwBc $ptwBc maxIter $maxIter lambdaBc $lambdaBc @@ -106,7 +113,11 @@ ## sets negative intensities to zero NegativetoZero $NegativetoZero - + ## final spectra + ## Graphical display + FinalGraph $FinalGraph + + ## Outputs dataMatrixOut $dataMatrixOut graphOut $graphOut @@ -118,14 +129,28 @@ <param name="dataMatrixFid" type="data" label="Data matrix of FIDs" help="" format="tabular" /> <param name="sampleMetadataFid" type="data" label="Sample metadata matrix" help="" format="tabular" /> - <param name="lambda" label="Smoothing parameter for solvent suppression" type="float" value="1000000" help="Default value is 1e6"/> + ## First order phase correction + <param name="FirstOPCGraph" label="Display the FIDs after 1st order phase correction?" type="select" help="Select 'YES' to display the spectra or 'NO' to not display them "> + <option value="NO"> NO </option> + <option value="YES"> YES </option> + </param> + + ## Water and / or solvents suppression + <param name="lambda" label="Solvent Suppression: Smoothing parameter" type="float" value="1000000" help="Default value is 1e6"/> <param name="ptwSS" label="Use of C library for computation of the solvent signal" type="select" help="Select 'Yes' to compute the solvent signal in C using the ptw package which is a lot faster"> <option value="NO"> NO </option> <option value="YES" selected="YES"> YES </option> </param> + + <param name="SSGraph" label="Display the FIDs after solvent suppression?" type="select" help="Select 'YES' to display the spectra or 'NO' to not display them "> + <option value="NO"> NO </option> + <option value="YES"> YES </option> + </param> + ## Apodization + <conditional name="apodizationMethod" > - <param name="method" label="Apodization method" type="select" help="Default method is Decreasing Exponential signal. See details below" > + <param name="method" label="Apodization: method" type="select" help="Default method is Decreasing Exponential signal. See details below" > <option value="exp" label="Decreasing Exponential window" /> <option value="cos2" label="Cosinus squared window"/> <option value="hanning" label="Hanning window"/> @@ -154,29 +179,25 @@ <param name="rectRatio" type="float" label="Proportion of signal in the window" value="0.5" help="Default value is 0.5" /> </when> <when value="gauss"> - <param name="gaussLB" type="float" label="Line broadening" value="1" help="Default value is 1" /> + <param name="gaussLB" type="float" label="Line broadening" value="0.3" help="Default value is 0.3" /> </when> - </conditional> + </conditional> - + <param name="ApodGraph" label="Display the FIDs after Apodization?" type="select" help="Select 'YES' to display the spectra or 'NO' to not display them "> + <option value="NO"> NO </option> + <option value="YES"> YES </option> + </param> + + ## Fourier transform <param name="FTGraph" label="Display the Fourier transformed spectra?" type="select" help="Select 'YES' to display the spectra or 'NO' to not display them "> <option value="NO"> NO </option> <option value="YES"> YES </option> </param> - - <conditional name="shiftReferencingMethod" > - <param name="method" label="Method of shift referencing used to find the reference compound peak" type="select" help="Method used to find the reference compound peak in spectra" > - <option value="max" label="Maximum real intensity location" /> - <option value="thres" label="Threshold"/> - </param> - <when value="max" /> - <when value="thres"> - <param name="shiftTreshold" type="float" label="Predefined threshold" value="2" help="Default treshold is 2" /> - </when> - </conditional> + ## shift Referencing + <conditional name="shiftReferencingRange" > - <param name="method" label="Definition of the search zone" type="select" help="Definition of the search zone" > + <param name="method" label="Shift Referencing: definition of the search zone" type="select" help="Definition of the search zone" > <option value="near0" label="Near the 0 ppm location" help="Near the 0 ppm location"/> <option value="all" label="Accross the whole ppm axis" help="Accross the whole ppm axis"/> <option value="window" label="Manually specified area of the ppm axis with the non-null parameter fromto" help="Manually specified area of the ppm axis with the non-null parameter fromto"/> @@ -187,56 +208,79 @@ <param name="pctNear0" type="float" label="percentage of the ppm axis around 0 ppm to look for the reference compound peak" value="0.02" help="Default treshold is 0.02" /> </when> <when value="window"> - <repeat name="conditions" title="Search_zone (in ppm)"> + <repeat name="conditions" title="Search_zone"> <param name="shiftReferencingRangeLeft" label="Search zone: left border" type="float" value="10.0" /> <param name="shiftReferencingRangeRight" label="Search zone: right border" type="float" value="10.0" /> </repeat> - </when> + </when> </conditional> - <param name="shiftHandling" type="select" label="shiftHandling" help="How to deal with shifts between spectra on their left and right sides" > + <param name="shiftHandling" type="select" label="Shift Referencing: shiftHandling" help="How to deal with shifts between spectra on their left and right sides" > <option value="zerofilling" selected="zerofilling"> zerofilling </option> <option value="cut" > cut </option> <option value="circular" > circular </option> </param> + + <param name="SRGraph" label="Display the spectra after Shift Referencing?" type="select" help="Select 'YES' to display the spectra or 'NO' to not display them "> + <option value="NO"> NO </option> + <option value="YES"> YES </option> + </param> - <conditional name="zeroOrderPhaseMethod" > - <param name="method" label="Zero Order Phase correction method" type="select" help="Zero Order Phase correction method" > - <option value="rms" selected="rms" label="A positiveness criterion is applied on the spectrum" /> - </param> - <when value="rms"> - <param name="quant" type="float" label="Quantile probability for the positiveness criterion" value="0.95" help="Default value is 0.95" /> - </when> - </conditional> + ## zero Order Phase Correction + + <param name="zeroOrderPhaseMethod" type="select" label="Zero Order Phase Correction: method" help="Method used to select the angles to rotate the spectra" > + <option value="rms" selected="rms"> rms </option> + <option value="max" > max </option> + </param> - <conditional name="searchZoneZeroPhase"> - <param name="choice" type="select" label="Search zone for Zero order phase correction" help="Choose if you want to exclude particular zone(s)" > - <option value="YES" > YES </option> - <option value="NO" selected="true"> NO </option> + <conditional name="excludeZoneZeroPhase"> + <param name="choice" type="select" label="Zero Order Phase Correction: exclusion area(s)" help="Choose if you want to exclude particular zone(s)" > + <option value="YES" selected="true" > YES </option> + <option value="NO" > NO </option> </param> <when value="YES"> - <repeat name="conditions" title="Search_zone"> - <param name="searchZoneZeroPhase_left" label="Search zone: left border" type="float" value="10.0" /> - <param name="searchZoneZeroPhase_right" label="Search zone: right border" type="float" value="10.0" /> + <repeat name="conditions" title="Exclusion_zone" min="1"> + <param name="excludeZoneZeroPhase_left" label="Excusion zone: left border" type="float" value="5.1" /> + <param name="excludeZoneZeroPhase_right" label="Excusion zone: right border" type="float" value="4.5" /> </repeat> </when> <when value="no"> </when> </conditional> - <param name="ptwBc" label="Use ptw for baseline computation?" type="select" help="If TRUE, calculates the baseline in C using the ptw library which is a lot faster. The R version is only kept because it is easier to understand than C and in case of problems with the installation of ptw.maxIter"> + <param name="ZeroOPCGraph" label="Display the spectra after the Zero Order Phase Correction?" type="select" help="Select 'YES' to display the spectra or 'NO' to not display them "> + <option value="NO"> NO </option> + <option value="YES"> YES </option> + </param> + + ## baseline correction + + <param name="ptwBc" label="Baseline Correction: use ptw for baseline computation?" type="select" help="If TRUE, calculates the baseline in C using the ptw library which is a lot faster. The R version is only kept because it is easier to understand than C and in case of problems with the installation of ptw.maxIter"> <option value="FALSE"> NO </option> <option value="TRUE" selected="true"> YES </option> </param> + <param name="maxIter" type="integer" label="Baseline Correction: maximum number of iterations if ptw.bc is set to FALSE" value="50" help="Maximum of iterations if ptw.bc is set to FALSE. Default value is 50" /> + <param name="lambdaBc" type="float" label="Baseline Correction: smoothing parameter" value="100000.0" help="Smoothing parameter, generally 1e5 – 1e8. Default value is 100000" /> + <param name="pBc" type="float" label="Baseline Correction: asymmetry parameter" value="0.05" help="Asymmetry parameter. Default value is 0.05" /> + <param name="epsilon" type="float" label="Baseline Correction: numerical precision for convergence when estimating the baseline" value="0.00000001" help="Numerical precision for convergence when estimating the baseline. Default value is 1e-8" /> - <param name="maxIter" type="integer" label="Maximum of iterations if ptw.bc is set to FALSE" value="50" help="Maximum of iterations if ptw.bc is set to FALSE. Default value is 50" /> - <param name="lambdaBc" type="float" label="Smoothing parameter" value="100000.0" help="Smoothing parameter, generally 1e5 – 1e8. Default value is 100000" /> - <param name="pBc" type="float" label="Asymmetry parameter" value="0.05" help="Asymmetry parameter. Default value is 0.05" /> - <param name="epsilon" type="float" label="Numerical precision for convergence when estimating the baseline" value="0.00000001" help="Numerical precision for convergence when estimating the baseline. Default value is 1e-8" /> + <param name="BCGraph" label="Display the spectra after Baseline Correction?" type="select" help="Select 'YES' to display the spectra or 'NO' to not display them "> + <option value="NO"> NO </option> + <option value="YES"> YES </option> + </param> + + ## NegativetoZero <param name="NegativetoZero" label="Set negative intensities to zero?" type="select" help="If YES, sets negative intensities to zero"> <option value="NO"> NO </option> <option value="YES"> YES </option> </param> + + ## final spectra + <param name="FinalGraph" label="Display the final spectra?" type="select" help="Select 'YES' to display the spectra or 'NO' to not display them "> + <option value="YES"> YES </option> + <option value="NO"> NO </option> + </param> + </inputs> <outputs> @@ -269,7 +313,7 @@ These steps correspond to the following steps in the PEPS-NMR R library (https://github.com/ManonMartin/PEPSNMR): -* First order phase correction +* Group Delay suppression (First order phase correction) * Removal of solvent residuals signal from the FID * Apodization to increase the Signal-to-Noise ratio of the FID * Fourier transformation @@ -287,7 +331,6 @@ The **types of apodization** are: - * exp: The signal is multiplied by a decreasing exponential exp(-t/LineBroadening). * 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). @@ -306,13 +349,6 @@ **Shift referencing** ---------------------- -Different **methods for shift referencing**: - -* max: the maximum intensity in the search zone is defined as the reference compound peak. - -* 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). - - The **searching window** can be adapted: * near0: the search concentrates around the 0ppm. @@ -333,8 +369,6 @@ * cut: The ppm values for which some spectra are not defined are removed. -**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. - **Zero Order Phase correction** ----------------------------------- @@ -348,7 +382,7 @@ * max: Optimization of the maximal spectral intensity. -**Search zone for the Zero order phase correction**: enables to optimize the criterion only on the selected spectral window(s). +** Exclusion area(s) for the Zero order phase correction**: enables to optimize the criterion with excluded spectral window(s), by default the water region is excluded. **Baseline computation**