Mercurial > repos > melpetera > batchcorrection

comparison batch_correction_all_loess_script.R @ 1:f64656ae9ea4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/workflow4metabolomics/batchcorrection.git commit b1f8bd1260c1c4b73600fb3867ca3bc613f258a7
author melpetera
date Sun, 09 Oct 2016 09:54:00 -0400
parents 71d83d8920bf
children
comparison
equal deleted inserted replaced
0:71d83d8920bf 1:f64656ae9ea4
1 loessF <- function(datVn, qcaVi, preVi, spnN, vrbL=FALSE) { 1 loessF <- function(datVn, qcaVi, preVi, spnN) {
2 2
3 if(length(qcaVi) < 5) { 3 if(length(qcaVi) < 5) {
4
5 if(vrbL)
6 cat("\nWarning: less than 5 '", refC, "'; linear regression will be performed instead of loess regression for this batch\n", sep="")
7 4
8 return(predict(lm(datVn[qcaVi] ~ qcaVi), 5 return(predict(lm(datVn[qcaVi] ~ qcaVi),
9 newdata = data.frame(qcaVi = preVi))) 6 newdata = data.frame(qcaVi = preVi)))
10 7
11 } else { 8 } else {
106 radVn <- seq(0, 2 * pi, length.out = 100) 103 radVn <- seq(0, 2 * pi, length.out = 100)
107 epsN <- .Machine[["double.eps"]] ## [1] 2.22e-16 104 epsN <- .Machine[["double.eps"]] ## [1] 2.22e-16
108 105
109 pcaMN <- datMN 106 pcaMN <- datMN
110 107
111 pcaLs <- opls(pcaMN, predI = 4, printL = FALSE, plotL = FALSE) 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)
112 tMN <- getScoreMN(pcaLs) 114 tMN <- getScoreMN(pcaLs)
113 vRelVn <- pcaLs@modelDF[, "R2X"] 115 vRelVn <- pcaLs@modelDF[, "R2X"]
114 116
115 n <- nrow(tMN) 117 n <- nrow(tMN)
116 hotN <- 2 * (n - 1) * (n^2 - 1) / (n^2 * (n - 2)) 118 hotN <- 2 * (n - 1) * (n^2 - 1) / (n^2 * (n - 2))
261 263
262 batAllVi <- 1:nrow(batRawMN) 264 batAllVi <- 1:nrow(batRawMN)
263 265
264 batRefVi <- grep(refC.arg, batSamDF[, "sampleType"]) 266 batRefVi <- grep(refC.arg, batSamDF[, "sampleType"])
265 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
266 ## prediction of the loess fit 271 ## prediction of the loess fit
267 272
268 batLoeMN <- apply(batRawMN, 273 batLoeMN <- apply(batRawMN,
269 2, 274 2,
270 function(rawVn) loessF(rawVn, batRefVi, batAllVi, spnN=spnN.arg, vrbL=TRUE)) 275 function(rawVn) loessF(rawVn, batRefVi, batAllVi, spnN=spnN.arg))
271 276
272 ## normalization 277 ## normalization
273 278
274 batLoeMN[batLoeMN <= 0] <- NA 279 batLoeMN[batLoeMN <= 0] <- NA
275 280