Mercurial > repos > melpetera > batchcorrection
comparison BC/batch_correction_all_loess_script.R @ 4:23314e1192d4 draft default tip
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 |