Mercurial > repos > melpetera > batchcorrection

comparison BC/batch_correction_all_loess_script.R @ 4:23314e1192d4 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author melpetera
date Thu, 14 Jan 2021 09:56:58 +0000
parents
children
comparison
equal deleted inserted replaced
3:73892ef177e3 4:23314e1192d4
1 loessF <- function(datVn, qcaVi, preVi, spnN) {
2
3 if(length(qcaVi) < 5) {
4
5 return(predict(lm(datVn[qcaVi] ~ qcaVi),
6 newdata = data.frame(qcaVi = preVi)))
7
8 } else {
9
10 return(predict(loess(datVn[qcaVi] ~ qcaVi,
11 control = loess.control(surface = "direct"),
12 span = spnN),
13 newdata = data.frame(qcaVi = preVi)))
14
15 }
16
17 ## Note:
18 ## the surface = 'direct' argument allows extrapolation
19
20 } ## loessF
21
22 plotBatchF <- function(datMN, samDF.arg, spnN.arg) {
23
24 maiC <- switch(gsub("MN", "", deparse(substitute(datMN))),
25 raw = "Raw",
26 nrm = "Normalized")
27
28 colVc <- c(samp = "green4",
29 biol = "green4",
30 pool = "red",
31 blan = "black",
32 other = "yellow")
33
34 par(font = 2, font.axis = 2, font.lab = 2, lwd = 2, pch = 18)
35
36 layout(matrix(c(1, 1, 2, 3), nrow = 2),
37 widths = c(0.7, 0.3))
38
39 obsNamVc <- rownames(datMN)
40
41 obsColVc <- sapply(substr(samDF.arg[, args$sample_type_col_name], 1, 4),
42 function(typC)
43 ifelse(typC %in% names(colVc), colVc[typC], colVc["other"]))
44
45 ## Graphic 1: Sum of intensities for each sample
46
47 par(mar = c(3.6, 3.6, 3.1, 0.6))
48
49 batTab <- table(samDF.arg[, args$batch_col_name])
50
51 sumVn <- rowSums(datMN, na.rm = TRUE)
52
53 plot(sumVn,
54 cex = 1.2,
55 col = obsColVc,
56 pch = 18,
57 xaxs = "i",
58 xlab = "",
59 ylab = "")
60
61 mtext("Injection order",
62 line = 2.2,
63 side = 1)
64 mtext("Sum of variable intensities",
65 line = 2.2,
66 side = 2)
67
68 mtext(maiC, cex = 1.2, line = 1.5, side = 3)
69
70 abline(v = cumsum(batTab) + 0.5,
71 col = "red")
72
73 mtext(names(batTab),
74 at = batTab / 2 + c(0, cumsum(batTab[-length(batTab)])))
75
76 obsColVuc <- obsColVc[sort(unique(names(obsColVc)))]
77
78 text(rep(batTab[1], times = length(obsColVuc)),
79 par("usr")[3] + (0.97 - length(obsColVuc) * 0.03 + 1:length(obsColVuc) * 0.03) * diff(par("usr")[3:4]),
80 col = obsColVuc,
81 font = 2,
82 labels = names(obsColVuc),
83 pos = 2)
84
85 for(batC in names(batTab)) {
86
87 batSeqVi <- which(samDF.arg[, args$batch_col_name] == batC)
88 batPooVi <- intersect(batSeqVi,
89 which(samDF.arg[, args$sample_type_col_name] == "pool"))
90 batSamVi <- intersect(batSeqVi,
91 which(samDF.arg[, args$sample_type_col_name] == "sample"))
92 if(length(batPooVi))
93 lines(batSeqVi,
94 loessF(sumVn, batPooVi, batSeqVi, spnN=spnN.arg),
95 col = colVc["pool"])
96 if (length(batSamVi))
97 lines(batSeqVi,
98 loessF(sumVn, batSamVi, batSeqVi, spnN=spnN.arg),
99 col = colVc["samp"])
100
101 }
102
103 ## Graphics 2 and 3 (right): PCA score plots of components 1-4
104
105 radVn <- seq(0, 2 * pi, length.out = 100)
106 epsN <- .Machine[["double.eps"]] ## [1] 2.22e-16
107
108 pcaMN <- datMN
109
110 if(any(is.na(pcaMN))) {
111 minN <- min(pcaMN, na.rm = TRUE)
112 pcaMN[is.na(pcaMN)] <- minN
113 }
114
115 pcaLs <- opls(pcaMN, predI = 4, algoC = "svd", printL = FALSE, plotL = FALSE)
116 tMN <- getScoreMN(pcaLs)
117 vRelVn <- pcaLs@modelDF[, "R2X"]
118
119 n <- nrow(tMN)
120 hotN <- 2 * (n - 1) * (n^2 - 1) / (n^2 * (n - 2))
121
122 hotFisN <- hotN * qf(0.95, 2, n - 2)
123
124 pcsLs <- list(c(1, 2), c(3, 4))
125
126 par(mar = c(3.6, 3.6, 0.6, 1.1))
127
128 for(pcsN in 1:length(pcsLs)) {
129
130 pcsVn <- pcsLs[[pcsN]]
131
132 tcsMN <- tMN[, pcsVn]
133
134 micMN <- solve(cov(tcsMN))
135
136 n <- nrow(tMN)
137 hotN <- 2 * (n - 1) * (n^2 - 1) / (n^2 * (n - 2))
138
139 hotFisN <- hotN * qf(0.95, 2, n - 2)
140
141 hotVn <- apply(tcsMN,
142 1,
143 function(x) 1 - pf(1 / hotN * t(as.matrix(x)) %*% micMN %*% as.matrix(x), 2, n - 2))
144
145 obsHotVi <- which(hotVn < 0.05)
146
147 xLabC <- paste("t",
148 pcsVn[1],
149 "(",
150 round(vRelVn[pcsVn[1]] * 100),
151 "%)",
152 sep = "")
153
154 yLabC <- paste("t",
155 pcsVn[2],
156 "(",
157 round(vRelVn[pcsVn[2]] * 100),
158 "%)",
159 sep = "")
160
161 xLimVn <- c(-1, 1) * max(sqrt(var(tcsMN[, 1]) * hotFisN), max(abs(tcsMN[, 1])))
162 yLimVn <- c(-1, 1) * max(sqrt(var(tcsMN[, 2]) * hotFisN), max(abs(tcsMN[, 2])))
163
164 plot(tcsMN,
165 main = "",
166 type = "n",
167 xlab = "",
168 ylab = "",
169 xlim = xLimVn,
170 ylim = yLimVn)
171
172 mtext(xLabC,
173 line = 2.2,
174 side = 1)
175 mtext(yLabC,
176 line = 2.2,
177 side = 2)
178
179 par(lwd = 1)
180
181 abline(v = axTicks(1),
182 col = "grey")
183
184 abline(h = axTicks(2),
185 col = "grey")
186
187 abline(v = 0)
188 abline(h = 0)
189
190 lines(sqrt(var(tcsMN[, 1]) * hotFisN) * cos(radVn),
191 sqrt(var(tcsMN[, 2]) * hotFisN) * sin(radVn))
192
193 points(tcsMN,
194 col = obsColVc,
195 pch = 18)
196
197 if(length(obsHotVi))
198 text(tcsMN[obsHotVi, 1],
199 tcsMN[obsHotVi, 2],
200 col = obsColVc[obsHotVi],
201 labels = obsNamVc[obsHotVi],
202 pos = 3)
203
204 } ## for(pcsN in 1:length(pcsLs)) {
205
206 return(invisible(list(sumVn = sumVn,
207 tcsMN = tcsMN)))
208
209 } ## plotBatchF
210
211 shiftBatchCorrectF <- function(rawMN.arg,
212 samDF.arg,
213 refC.arg,
214 spnN.arg) {
215
216 cat("\nReference observations are: ", refC.arg, "\n")
217
218 ## computing median off all pools (or samples) for each variable
219
220 refMeaVn <- apply(rawMN.arg[samDF.arg[, args$sample_type_col_name] == refC.arg, ],
221 2,
222 function(feaRefVn) mean(feaRefVn, na.rm = TRUE))
223
224 ## splitting data and sample metadata from each batch
225
226 batRawLs <- split(as.data.frame(rawMN.arg),
227 f = samDF.arg[, args$batch_col_name])
228 batRawLs <- lapply(batRawLs, function(inpDF) as.matrix(inpDF))
229
230 batSamLs <- split(as.data.frame(samDF.arg),
231 f = samDF.arg[, args$batch_col_name])
232
233 ## checking extrapolation: are there pools at the first and last observations of each batch
234
235 if(refC.arg == args$sample_type_tags$pool) {
236 pooExtML <- matrix(FALSE, nrow = 2, ncol = length(batRawLs),
237 dimnames = list(c("first", "last"), names(batRawLs)))
238 for(batC in names(batSamLs)) {
239 batSamTypVc <- batSamLs[[batC]][, args$sample_type_col_name]
240 pooExtML["first", batC] <- head(batSamTypVc, 1) == args$sample_type_tags$pool
241 pooExtML["last", batC] <- tail(batSamTypVc, 1) == args$sample_type_tags$pool
242 }
243 if(!all(c(pooExtML))) {
244 cat("\nWarning: Pools are missing at the first and/or last position of the following batches:\n")
245 pooExtBatVi <- which(!apply(pooExtML, 2, all))
246 for(i in 1:length(pooExtBatVi))
247 cat(names(pooExtBatVi)[i], ": ",
248 paste(rownames(pooExtML)[!pooExtML[, pooExtBatVi[i]]], collapse = ", "), "\n", sep = "")
249 cat("Extrapolating loess fits for these batches may result in inaccurate modeling!\n")
250 }
251 }
252
253 ## normalizing
254
255 nrmMN <- NULL ## normalized data matrix to be computed
256
257 cat("\nProcessing batch:")
258
259 for(batC in names(batRawLs)) { ## processing each batch individually
260
261 cat("\n", batC)
262
263 batRawMN <- batRawLs[[batC]]
264 batSamDF <- batSamLs[[batC]]
265
266 batAllVi <- 1:nrow(batRawMN)
267
268 batRefVi <- which(batSamDF[, args$sample_type_col_name] == refC.arg)
269
270 if(length(batRefVi) < 5)
271 cat("\nWarning: less than 5 '", refC.arg, "'; linear regression will be performed instead of loess regression for this batch\n", sep="")
272
273 ## prediction of the loess fit
274
275 batLoeMN <- apply(batRawMN,
276 2,
277 function(rawVn) loessF(rawVn, batRefVi, batAllVi, spnN=spnN.arg))
278
279 ## normalization
280
281 batLoeMN[batLoeMN <= 0] <- NA
282
283 batNrmMN <- batRawMN / batLoeMN
284
285 nrmMN <- rbind(nrmMN,
286 batNrmMN)
287
288 }
289
290 cat("\n")
291
292 nrmMN <- sweep(nrmMN, MARGIN = 2, STATS = refMeaVn, FUN = "*")
293
294 return(nrmMN)
295
296 } ## shiftBatchCorrectF