Mercurial > repos > marie-tremblay-metatoul > nmr_alignment

comparison nmr_alignement/NmrAlignment_wrapper.R @ 4:265ee538e120 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author marie-tremblay-metatoul
date Wed, 09 Aug 2017 06:28:10 -0400
parents f3ec6799c435
children
comparison
equal deleted inserted replaced
3:f3ec6799c435 4:265ee538e120
69 69
70 ## Spectral width 70 ## Spectral width
71 leftBorder <- argLs[["left_border"]] 71 leftBorder <- argLs[["left_border"]]
72 rightBorder <- argLs[["right_border"]] 72 rightBorder <- argLs[["right_border"]]
73 73
74 ##Exclusion zone(s) 74 ## Exclusion zone(s)
75 exclusionZones <- argLs[["zone_exclusion_choices.choice"]] 75 exclusionZones <- argLs[["zone_exclusion_choices.choice"]]
76 exclusionZonesBorders <- NULL 76 exclusionZonesBorders <- NULL
77 if (!is.null(argLs$zone_exclusion_left)) 77 if (!is.null(argLs$zone_exclusion_left))
78 { 78 {
79 for(i in which(names(argLs)=="zone_exclusion_left")) 79 for(i in which(names(argLs)=="zone_exclusion_left"))
89 nDivRange <- argLs[["nDivRange"]] 89 nDivRange <- argLs[["nDivRange"]]
90 90
91 ## Intensity threshold for peak removal 91 ## Intensity threshold for peak removal
92 baselineThresh <- argLs[["baselineThresh"]] 92 baselineThresh <- argLs[["baselineThresh"]]
93 93
94 ## Graphical zone(s)
95 graphicalZones <- argLs[["zone_graphical_choices.choice"]]
96 graphicalZonesBorders <- NULL
97 if (!is.null(argLs$zone_exclusion_left))
98 {
99 for(i in which(names(argLs)=="zone_graphical_left"))
100 {
101 graphicalZonesBorders <- c(graphicalZonesBorders,list(c(argLs[[i]],argLs[[i+1]])))
102 }
103 }
104
94 105
95 # Outputs 106 # Outputs
96 logOut <- argLs[["logOut"]] 107 logOut <- argLs[["logOut"]]
97 alignedSpectra <- argLs[["alignedSpectra"]] 108 alignedSpectra <- argLs[["alignedSpectra"]]
98 graphOut <- argLs[["graphOut"]] 109 graphOut <- argLs[["graphOut"]]
99 110
111
112 ## Checking R packages
113 ##--------------------
114 sink(logOut)
115 cat("\tPACKAGE INFO\n")
116 #pkgs=c("xcms","batch")
117 pkgs=c("batch","speaq")
118 for(pkg in pkgs) {
119 suppressPackageStartupMessages( stopifnot( library(pkg, quietly=TRUE, logical.return=TRUE, character.only=TRUE)))
120 cat(pkg,"\t",as.character(packageVersion(pkg)),"\n",sep="")
121 }
100 122
101 ## Checking arguments 123 ## Checking arguments
102 ##------------------- 124 ##-------------------
103 error.stock <- "\n" 125 error.stock <- "\n"
104 if(length(error.stock) > 1) 126 if(length(error.stock) > 1)
120 rownames(t.directory.aligned) <- colnames(directory.aligned) 142 rownames(t.directory.aligned) <- colnames(directory.aligned)
121 # colnames(t.directory.aligned) <- c("Bucket",colnames(t.directory.aligned)) 143 # colnames(t.directory.aligned) <- c("Bucket",colnames(t.directory.aligned))
122 write.table(t.directory.aligned,file=alignedSpectra,row.names=TRUE,quote=FALSE,sep="\t") 144 write.table(t.directory.aligned,file=alignedSpectra,row.names=TRUE,quote=FALSE,sep="\t")
123 145
124 146
125 excludedZone <- NULL
126 for (c in 1:length(exclusionZonesBorders))
127 {
128 excludedZone <- c(excludedZone,exclusionZonesBorders[[c]])
129 excludedZone <- sort(excludedZone)
130 }
131
132 ## Graphical output: overlay of raw and estimated spectra 147 ## Graphical output: overlay of raw and estimated spectra
133 pdf(graphOut,onefile=TRUE) 148 pdf(graphOut,onefile=TRUE)
149
150 graphical.zone.length <- length(graphicalZonesBorders)
134 par(mfrow=c(2,1)) 151 par(mfrow=c(2,1))
135 152
136 raw.spectra <- data.frame(directory.raw) 153 raw.spectra <- data.frame(directory.raw)
137 colnames(raw.spectra) <- substr(colnames(raw.spectra),2,7) 154 colnames(raw.spectra) <- substr(colnames(raw.spectra),2,7)
138 155
140 colnames(aligned.spectra) <- substr(colnames(aligned.spectra),2,7) 157 colnames(aligned.spectra) <- substr(colnames(aligned.spectra),2,7)
141 158
142 drawSpec(raw.spectra,xlab="", ylab="Raw spectra", main="") 159 drawSpec(raw.spectra,xlab="", ylab="Raw spectra", main="")
143 drawSpec(aligned.spectra,xlab="", ylab="Aligned spectra", main="") 160 drawSpec(aligned.spectra,xlab="", ylab="Aligned spectra", main="")
144 161
145 nbZones <- length(excludedZone)/2 162 if (graphical.zone.length != 0)
146 if (nbZones != 0)
147 { 163 {
148 n <- length(excludedZone) 164 # par(mfrow=c(2,1))
149 drawSpec(raw.spectra[,1:which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[n])[1]],xlab="", ylab="Raw spectra", main="") 165 for (g in 1:graphical.zone.length)
150 drawSpec(aligned.spectra[,1:which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[n])[1]],xlab="", ylab="Aligned spectra", main="")
151
152 n <- n - 1
153 while (n >= nbZones & nbZones > 1)
154 { 166 {
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="") 167 print(paste(g, graphicalZonesBorders[[g]][1], graphicalZonesBorders[[g]][2]))
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="") 168 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == graphicalZonesBorders[[g]][1])[1]):(which(round(as.numeric(colnames(raw.spectra)),2) == graphicalZonesBorders[[g]][2])[1])],xlab="", main="")
157 n <- n - 2 169 drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == graphicalZonesBorders[[g]][1])[1]):(which(round(as.numeric(colnames(aligned.spectra)),2) == graphicalZonesBorders[[g]][2])[1])],xlab="", main="")
158 } 170 }
159
160 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[1])[1]):ncol(raw.spectra)],xlab="", ylab="Raw 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="")
162 } 171 }
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="")
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="")
165
166 dev.off() 172 dev.off()
167 173
168 174
169 ## Ending 175 ## Ending
170 ##--------------------- 176 ##---------------------
171 cat("\nEnd of 'NMR alignment' Galaxy module call: ", as.character(Sys.time()), sep = "") 177 cat("\nEnd of 'NMR alignment' Galaxy module call: ", as.character(Sys.time()), sep = "")
178 sink()
172 options(stringsAsFactors = strAsFacL) 179 options(stringsAsFactors = strAsFacL)
173 rm(list = ls()) 180 rm(list = ls())