Mercurial > repos > marie-tremblay-metatoul > nmr_alignment

comparison nmr_alignement/NmrAlignment_wrapper.R @ 3:f3ec6799c435 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author marie-tremblay-metatoul
date Wed, 12 Apr 2017 04:32:23 -0400
parents 908e1345d7ca
children 265ee538e120
comparison
equal deleted inserted replaced
2:908e1345d7ca 3:f3ec6799c435
15 15
16 ##------------------------------ 16 ##------------------------------
17 ## Libraries loading 17 ## Libraries loading
18 ##------------------------------ 18 ##------------------------------
19 # ParseCommandArgs function 19 # ParseCommandArgs function
20 library(batch) 20 library(batch)
21 # Alignment 21 # Alignment
22 library(speaq) 22 library(speaq)
23 23
24 24
25 # R script call 25 # R script call
49 ## Script 49 ## Script
50 ##------------------------------ 50 ##------------------------------
51 if(!runExampleL) 51 if(!runExampleL)
52 argLs <- parseCommandArgs(evaluate=FALSE) 52 argLs <- parseCommandArgs(evaluate=FALSE)
53 53
54 54
55 ## Parameters Loading 55 ## Parameters Loading
56 ##------------------- 56 ##-------------------
57 # Inputs 57 # Inputs
58 ## Library of spectra to align 58 ## Library of spectra to align
59 if (!is.null(argLs[["zipfile"]])){ 59 if (!is.null(argLs[["zipfile"]])){
60 fileType="zip" 60 fileType="zip"
61 zipfile= argLs[["zipfile"]] 61 zipfile= argLs[["zipfile"]]
62 directory=unzip(zipfile, list=F) 62 directory=unzip(zipfile, list=F)
63 directory=paste(getwd(),strsplit(directory[1],"/")[[1]][2],sep="/") 63 directory=paste(getwd(),strsplit(directory[1],"/")[[1]][2],sep="/")
64 } else if (!is.null(argLs[["library"]])){
65 fileType="zip"
66 directory=argLs[["library"]]
67 if(!file.exists(directory)){
68 error_message=paste("Cannot access the directory :",directory,".Please verify if the directory exists or not.")
69 print(error_message)
70 stop(error_message)
71 }
72 } else if (!is.null(argLs[["tsvfile"]])){ 64 } else if (!is.null(argLs[["tsvfile"]])){
73 fileType="tsv" 65 fileType="tsv"
74 directory <- read.table(argLs[["tsvfile"]],check.names=FALSE,header=TRUE,sep="\t") 66 directory <- read.table(argLs[["tsvfile"]],check.names=FALSE,header=TRUE,sep="\t")
75 } 67 }
76 68
77 69
78 ## Spectral width 70 ## Spectral width
79 leftBorder <- argLs[["left_border"]] 71 leftBorder <- argLs[["left_border"]]
80 rightBorder <- argLs[["right_border"]] 72 rightBorder <- argLs[["right_border"]]
81 73
82 ##Exclusion zone(s) 74 ##Exclusion zone(s)
83 exclusionZones <- argLs[["zone_exclusion_choices.choice"]] 75 exclusionZones <- argLs[["zone_exclusion_choices.choice"]]
84 exclusionZonesBorders <- NULL 76 exclusionZonesBorders <- NULL
85 if (!is.null(argLs$zone_exclusion_left)) 77 if (!is.null(argLs$zone_exclusion_left))
86 { 78 {
87 for(i in which(names(argLs)=="zone_exclusion_left")) 79 for(i in which(names(argLs)=="zone_exclusion_left"))
88 { 80 {
89 exclusionZonesBorders <- c(exclusionZonesBorders,list(c(argLs[[i]],argLs[[i+1]]))) 81 exclusionZonesBorders <- c(exclusionZonesBorders,list(c(argLs[[i]],argLs[[i+1]])))
90 } 82 }
91 } 83 }
92 84
93 ## Reference spectrum 85 ## Reference spectrum
94 reference <- argLs[["reference"]] 86 reference <- argLs[["reference"]]
95 87
96 ## Size of a small nDivRange 88 ## Size of a small nDivRange
97 nDivRange <- argLs[["nDivRange"]] 89 nDivRange <- argLs[["nDivRange"]]
98 90
99 ## Intensity threshold for peak removal 91 ## Intensity threshold for peak removal
100 baselineThresh <- argLs[["baselineThresh"]] 92 baselineThresh <- argLs[["baselineThresh"]]
101 93
102 94
103 # Outputs 95 # Outputs
109 ## Checking arguments 101 ## Checking arguments
110 ##------------------- 102 ##-------------------
111 error.stock <- "\n" 103 error.stock <- "\n"
112 if(length(error.stock) > 1) 104 if(length(error.stock) > 1)
113 stop(error.stock) 105 stop(error.stock)
114 106
115 107
116 ## Computation 108 ## Computation
117 ##------------ 109 ##------------
118 directory.alignement <- nmr.alignment(fileType=fileType,directory=directory,leftBorder=leftBorder,rightBorder=rightBorder,exclusionZones=exclusionZones, 110 directory.alignement <- nmr.alignment(fileType=fileType,directory=directory,leftBorder=leftBorder,rightBorder=rightBorder,exclusionZones=exclusionZones,
119 exclusionZonesBorders=exclusionZonesBorders, reference=reference, nDivRange=nDivRange, 111 exclusionZonesBorders=exclusionZonesBorders, reference=reference, nDivRange=nDivRange,
120 baselineThresh=baselineThresh, maxshift=50, verbose=FALSE) 112 baselineThresh=baselineThresh, maxshift=50, verbose=FALSE)
121 directory.raw <- directory.alignement[[1]] 113 directory.raw <- directory.alignement[[1]]
122 directory.aligned <- directory.alignement[[2]] 114 directory.aligned <- directory.alignement[[2]]
123 115
124 ## Saving 116 ## Saving
125 ##------- 117 ##-------
126 # Aligned spectra 118 # Aligned spectra
127 t.directory.aligned <- t(directory.aligned) 119 t.directory.aligned <- t(directory.aligned)
128 rownames(t.directory.aligned) <- colnames(directory.aligned) 120 rownames(t.directory.aligned) <- colnames(directory.aligned)
129 # colnames(t.directory.aligned) <- c("Bucket",colnames(t.directory.aligned)) 121 # colnames(t.directory.aligned) <- c("Bucket",colnames(t.directory.aligned))
130 write.table(t.directory.aligned,file=alignedSpectra,row.names=TRUE,quote=FALSE,sep="\t") 122 write.table(t.directory.aligned,file=alignedSpectra,row.names=TRUE,quote=FALSE,sep="\t")
131 123
132 124
133 excludedZone <- NULL 125 excludedZone <- NULL
134 for (c in 1:length(exclusionZonesBorders)) 126 for (c in 1:length(exclusionZonesBorders))
135 { 127 {
162 { 154 {
163 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[n])[1]):(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[n-1])[1])],xlab="", ylab="Raw spectra", main="") 155 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[n])[1]):(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[n-1])[1])],xlab="", ylab="Raw spectra", main="")
164 drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[n])[1]):(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[n-1])[1])],xlab="", ylab="Aligned spectra", main="") 156 drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[n])[1]):(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[n-1])[1])],xlab="", ylab="Aligned spectra", main="")
165 n <- n - 2 157 n <- n - 2
166 } 158 }
167 159
168 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[1])[1]):ncol(raw.spectra)],xlab="", ylab="Raw spectra", main="") 160 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[1])[1]):ncol(raw.spectra)],xlab="", ylab="Raw spectra", main="")
169 drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[1])[1]):ncol(aligned.spectra)],xlab="", ylab="Aligned spectra", main="") 161 drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[1])[1]):ncol(aligned.spectra)],xlab="", ylab="Aligned spectra", main="")
170 } 162 }
171 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == 2.4)[1]):(which(round(as.numeric(colnames(raw.spectra)),2) == 2.8)[1])],xlab="", ylab="Raw spectra", main="") 163 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == 2.4)[1]):(which(round(as.numeric(colnames(raw.spectra)),2) == 2.8)[1])],xlab="", ylab="Raw spectra", main="")
172 drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == 2.4)[1]):(which(round(as.numeric(colnames(aligned.spectra)),2) == 2.8)[1])],xlab="", ylab="Aligned spectra", main="") 164 drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == 2.4)[1]):(which(round(as.numeric(colnames(aligned.spectra)),2) == 2.8)[1])],xlab="", ylab="Aligned spectra", main="")