3
|
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
|