Mercurial > repos > melpetera > batchcorrection

comparison batchcorrection/batch_correction_all_loess_script.R @ 2:016780b192a6 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author melpetera
date Thu, 02 Mar 2017 03:38:08 -0500
parents
children
comparison
equal deleted inserted replaced
1:f64656ae9ea4 2:016780b192a6
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