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