[Robast-commits] r511 - in pkg/RobLoxBioC: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Sep 12 13:50:05 CEST 2012
Author: stamats
Date: 2012-09-12 13:50:05 +0200 (Wed, 12 Sep 2012)
New Revision: 511
Added:
pkg/RobLoxBioC/R/00pre2160.R
Removed:
pkg/RobLoxBioC/R/sysdata.rda
Modified:
pkg/RobLoxBioC/DESCRIPTION
pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R
pkg/RobLoxBioC/R/robloxbiocMatrix.R
pkg/RobLoxBioC/inst/CITATION
pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd
pkg/RobLoxBioC/man/robloxbioc.Rd
Log:
several updates due to changes in R and package beadarray
Modified: pkg/RobLoxBioC/DESCRIPTION
===================================================================
--- pkg/RobLoxBioC/DESCRIPTION 2012-09-11 17:36:21 UTC (rev 510)
+++ pkg/RobLoxBioC/DESCRIPTION 2012-09-12 11:50:05 UTC (rev 511)
@@ -1,6 +1,6 @@
Package: RobLoxBioC
Version: 0.8.3
-Date: 2012-02-26
+Date: 2012-09-12
Title: Infinitesimally robust estimators for preprocessing omics data
Description: Functions for the determination of optimally robust influence curves and
estimators for preprocessing omics data, in particular gene expression data.
Added: pkg/RobLoxBioC/R/00pre2160.R
===================================================================
--- pkg/RobLoxBioC/R/00pre2160.R (rev 0)
+++ pkg/RobLoxBioC/R/00pre2160.R 2012-09-12 11:50:05 UTC (rev 511)
@@ -0,0 +1,43 @@
+## due to a change to .C and .Call in 2.16.0
+.getA1.locsc <- RobLox:::.getA1.locsc
+.getA2.locsc <- RobLox:::.getA2.locsc
+.geta.locsc <- RobLox:::.geta.locsc
+.getb.locsc <- RobLox:::.getb.locsc
+
+.getlsInterval <- function (r, rlo, rup){
+ if (r > 10) {
+ b <- 1.618128043
+ const <- 1.263094656
+ A2 <- b^2 * (1 + r^2)/(1 + const)
+ A1 <- const * A2
+ }
+ else {
+ A1 <- .getA1.locsc(r)
+ A2 <- .getA2.locsc(r)
+ b <- .getb.locsc(r)
+ }
+ if (rlo == 0) {
+ efflo <- (A1 + A2 - b^2 * r^2)/1.5
+ }
+ else {
+ A1lo <- .getA1.locsc(rlo)
+ A2lo <- .getA2.locsc(rlo)
+ efflo <- (A1 + A2 - b^2 * (r^2 - rlo^2))/(A1lo + A2lo)
+ }
+ if (rup > 10) {
+ bup <- 1.618128043
+ const.up <- 1.263094656
+ A2up <- bup^2 * (1 + rup^2)/(1 + const.up)
+ A1up <- const.up * A2up
+ effup <- (A1 + A2 - b^2 * (r^2 - rup^2))/(A1up + A2up)
+ }
+ else {
+ A1up <- .getA1.locsc(rup)
+ A2up <- .getA2.locsc(rup)
+ effup <- (A1 + A2 - b^2 * (r^2 - rup^2))/(A1up + A2up)
+ }
+ return(effup - efflo)
+}
+
+.onestep.locsc.matrix <- RobLox:::.onestep.locsc.matrix
+.kstep.locsc.matrix <- RobLox:::.kstep.locsc.matrix
\ No newline at end of file
Modified: pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R
===================================================================
--- pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R 2012-09-11 17:36:21 UTC (rev 510)
+++ pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R 2012-09-12 11:50:05 UTC (rev 511)
@@ -1,202 +1,252 @@
setMethod("robloxbioc", signature(x = "beadLevelData"),
- function(x, log = TRUE, imagesPerArray = 1, what = "G", probes = NULL, arrays = NULL,
- eps = NULL, eps.lower = 0, eps.upper = 0.05, steps = 3L,
- fsCor = TRUE, mad0 = 1e-4){
+ function(x, channelList = list(greenChannel), probeIDs = NULL, useSampleFac = FALSE,
+ sampleFac = NULL, weightNames = "wts", removeUnMappedProbes = TRUE,
+ eps = NULL, eps.lower = 0, eps.upper = 0.05, steps = 3L, fsCor = TRUE, mad0 = 1e-4){
BLData <- x
arraynms <- sectionNames(BLData)
- if(!is.null(arrays) && !is.character(arrays)) arraynms <- arraynms[arrays]
- if(is.character(arrays)) arraynms <- which(arraynms %in% arrays)
- len <- length(arraynms)
- what <- match.arg(what, c("G", "R", "RG", "M", "A", "beta"))
- whatelse <- ""
- if(what == "RG"){
- if(BLData at arrayInfo$channels == "two"){
- what <- "G"
- whatelse <- "R"
+ output <- vector("list", length(channelList))
+ if (useSampleFac) {
+ if (is.null(sampleFac)) {
+ if (!"SampleGroup" %in% names(BLData at sectionData)) {
+ cat("Could not determine sample factor from beadLevelData. Summarizing each section separately\n")
+ sList <- arraynms
+ sampleFac <- arraynms
+ newNames <- sList
+ }else{
+ sampleFac <- BLData at sectionData$SampleGroup[, 1]
+ sList <- unique(sampleFac)
+ dupList <- which(duplicated(sampleFac))
+ if (any(dupList)) {
+ newNames <- strtrim(arraynms[-dupList], 12)
+ }else{
+ newNames <- strtrim(arraynms, 12)
+ }
+ }
}else{
- stop("Need two-channel data to calculate summary R and G values")
+ if (length(sampleFac) != length(arraynms)) {
+ cat("Length of specified sample factor did not match number of sections\n")
+ cat("length of sample factor: ", length(sampleFac), sampleFac, "\n")
+ cat("number of sections: ", length(arraynms), arraynms, "\n")
+ stop("Aborting summarization\n")
+ }
+ else {
+ sList <- unique(sampleFac)
+ dupList <- which(duplicated(sampleFac))
+ if (any(dupList)) {
+ newNames <- strtrim(arraynms[-dupList], 12)
+ }else{
+ newNames <- strtrim(arraynms, 12)
+ }
+ }
}
}
- if(imagesPerArray == 1){
- pr <- getBeadData(BLData, what = "ProbeID", array = arraynms[1])
- sel <- pr != 0
- pr <- pr[sel]
- finten <- getBeadData(BLData, what = what, array = arraynms[1])[sel]
- if(log) finten <- log2(finten)
- nasinf <- !is.finite(finten) | is.na(finten)
- finten <- finten[!nasinf]
+ else {
+ cat("No sample factor specified. Summarizing each section separately\n")
+ sList <- arraynms
+ sampleFac <- arraynms
+ newNames <- arraynms
}
- else if(imagesPerArray == 2){
- if(length(arraynms)%%2 != 0)
- stop("Need an even number of arrays when 'imagesPerArray=2'")
- arrayord <- order(arraynms)
- arraynms <- arraynms[arrayord]
- tmp <- unlist(strsplit(arraynms, "_"))
- chipnums <- tmp[seq(1, length(tmp), by = 3)]
- pos <- tmp[seq(2, length(tmp), by = 3)]
- stripnum <- as.numeric(tmp[seq(3, length(tmp), by = 3)])
- check <- ((chipnums[seq(1, len, by = 2)] == chipnums[seq(2, len, by = 2)])
- & (pos[seq(1, len, by = 2)] == pos[seq(2, len, by = 2)])
- & (stripnum[seq(1, len, by = 2)] == stripnum[seq(2, len, by = 2)] - 1))
- if (sum(check) != length(check))
- stop("Missing arrays")
- sel1 <- getBeadData(BLData, what = "ProbeID", array = arraynms[1]) != 0
- sel2 <- getBeadData(BLData, what = "ProbeID", array = arraynms[2]) != 0
- pr <- append(getBeadData(BLData, what = "ProbeID", array = arraynms[1])[sel1],
- getBeadData(BLData, what = "ProbeID", array = arraynms[2])[sel2])
- finten <- append(getBeadData(BLData, what = what, array = arraynms[1])[sel1],
- getBeadData(BLData, what = what, array = arraynms[2])[sel2])
- if(log) finten <- log2(finten)
- nasinf <- !is.finite(finten) | is.na(finten)
- finten <- finten[!nasinf]
- }else{
- stop("You can only specify 1 or 2 images per array")
+ if (is.null(probeIDs)) {
+ cat("Finding list of unique probes in beadLevelData\n")
+ probeIDs <- beadarray:::uniqueProbeList(BLData)
+ cat(length(probeIDs), " unique probeIDs found\n")
}
- if(is.null(probes)) probes <- sort(unique(pr))
- probes <- probes[probes > 0 & !is.na(probes)]
- noprobes <- length(probes)
- pr <- pr[!nasinf]
- if (imagesPerArray == 1) {
- G <- GBeadStDev <- GNoBeads <- matrix(0, nrow = noprobes, ncol = len)
- colnames(G) <- colnames(GBeadStDev) <- colnames(GNoBeads) <- arraynms
- if (BLData at arrayInfo$channels == "two" && !is.null(BLData[[arraynms[1]]]$R) && whatelse == "R")
- R <- RBeadStDev <- RNoBeads <- G
- else R <- NULL
+ if (removeUnMappedProbes) {
+ annoName <- annotation(BLData)
+ if (!is.null(annoName)) {
+ annoLoaded <- require(paste("illumina", annoName, ".db", sep = ""), character.only = TRUE)
+ if (annoLoaded) {
+ mapEnv <- as.name(paste("illumina", annoName, "ARRAYADDRESS", sep = ""))
+ allMapped <- mappedkeys(revmap(eval(mapEnv)))
+ isMapped <- which(probeIDs %in% allMapped)
+ cat("Number of unmapped probes removed: ", length(probeIDs) - length(isMapped), "\n")
+ probeIDs <- probeIDs[isMapped]
+ }
+ }else{
+ cat("Could not determine annotation for this beadLevelData object.\n")
+ }
}
- else if (imagesPerArray == 2) {
- G <- GBeadStDev <- GNoBeads <- matrix(0, nrow = noprobes, ncol = (len/2))
- colnames(G) <- colnames(GBeadStDev) <- colnames(GNoBeads) <- arraynms[seq(1, len, by = 2)]
- if (BLData at arrayInfo$channels == "two" && !is.null(BLData[[arraynms[1]]]$R) && whatelse == "R")
- R <- RBeadStDev <- RNoBeads <- G
- else R <- NULL
+ cNames <- unlist(lapply(channelList, function(x) x at name))
+ if (any(duplicated(cNames))) {
+ uNames <- unique(cNames)
+ for (i in 1:length(uNames)) {
+ sPos <- grep(uNames[i], cNames)
+ if (length(sPos) > 1) {
+ for (j in 1:length(sPos)) {
+ cNames[sPos[j]] <- paste(cNames[sPos[j]], j, sep = ".")
+ warning("Duplicated channel names were found. Renaming...\n")
+ }
+ }
+ }
}
- i <- j <- 1
- while (j <= len) {
- probeIDs <- as.integer(pr)
- start <- 0
- blah <- rmxBeadSummary(x = finten, probeIDs = probeIDs, probes = probes,
- eps = eps, eps.lower = eps.lower, eps.upper = eps.upper,
- steps = steps, fsCor = fsCor, mad0 = mad0)
- G[, i] <- blah$fore
- GBeadStDev[, i] <- blah$sd
- GNoBeads[, i] <- blah$noBeads
- if (BLData at arrayInfo$channels == "two" && !is.null(BLData[[arraynms[i]]]$R) && whatelse == "R") {
- if (imagesPerArray == 1) {
- finten <- getBeadData(BLData, what = whatelse, array = arraynms[i])[sel]
- if(log) finten <- log2(finten)
- nasinf <- !is.finite(finten) | is.na(finten)
- finten <- finten[!nasinf]
+ for (cNum in 1:length(channelList)) {
+ template <- matrix(nrow = length(probeIDs), ncol = length(sList))
+ if (length(channelList) == 1) {
+ newCols <- newNames
+ }else{
+ newCols <- paste(cNames[cNum], newNames, sep = ":")
+ }
+ output[[cNum]][["eMat"]] <- template
+ colnames(output[[cNum]][["eMat"]]) <- newCols
+ rownames(output[[cNum]][["eMat"]]) <- probeIDs
+ output[[cNum]][["varMat"]] <- template
+ colnames(output[[cNum]][["varMat"]]) <- newCols
+ rownames(output[[cNum]][["varMat"]]) <- probeIDs
+ output[[cNum]][["nObs"]] <- template
+ colnames(output[[cNum]][["nObs"]]) <- newCols
+ rownames(output[[cNum]][["nObs"]]) <- probeIDs
+ }
+ for (s in 1:length(sList)) {
+ an <- which(sampleFac == sList[s])
+ pIDs <- wts <- NULL
+ values <- vector("list", length(channelList))
+ for (i in an) {
+ tmp <- BLData[[i]]
+ pidCol <- grep("ProbeID", colnames(tmp))
+ retainedBeads <- which(tmp[, pidCol] %in% probeIDs)
+ tmp <- tmp[retainedBeads, ]
+ wCol <- grep(weightNames, colnames(tmp))
+ pIDs <- c(pIDs, tmp[, pidCol])
+ if (length(wCol) == 0){
+ wts <- rep(1, length(pIDs))
+ }else{
+ wts <- c(wts, tmp[, wCol])
}
- else if (imagesPerArray == 2) {
- finten <- append(getBeadData(BLData, what = whatelse, array = arraynms[j])[sel1],
- getBeadData(BLData, what = whatelse, array = arraynms[j + 1])[sel2])
- if(log) finten <- log2(finten)
- nasinf <- !is.finite(finten) | is.na(finten)
- finten <- finten[!nasinf]
+ for (ch in 1:length(channelList)) {
+ chName <- channelList[[ch]]@name
+ transFun <- channelList[[ch]]@transFun[[1]]
+ cat("Summarizing ", chName, " channel\n")
+ cat("Processing Array", i, "\n")
+ newVals <- transFun(BLData, array = i)[retainedBeads]
+ if (length(newVals) != nrow(tmp))
+ stop("Transformation function did not return correct number of values")
+ values[[ch]] <- c(values[[ch]], newVals)
}
- blah <- rmxBeadSummary(x = finten, probeIDs = probeIDs, probes = probes,
- eps = eps, eps.lower = eps.lower, eps.upper = eps.upper,
- steps = steps, fsCor = fsCor, mad0 = mad0)
- R[, i] <- blah$fore
- RBeadStDev[, i] <- blah$sd
- RNoBeads[, i] <- blah$noBeads
}
- j <- j + imagesPerArray
- i <- i + 1
- rm(probeIDs, blah)
- gc()
- if ((imagesPerArray == 1) && (i <= len)) {
- sel <- getBeadData(BLData, what = "ProbeID", array = arraynms[i]) != 0
- pr <- getBeadData(BLData, what = "ProbeID", array = arraynms[i])[sel]
- finten <- getBeadData(BLData, what = what, array = arraynms[i])[sel]
- if(log) finten <- log2(finten)
- nasinf <- !is.finite(finten) | is.na(finten)
- pr <- pr[!nasinf]
- finten <- finten[!nasinf]
+ for (ch in 1:length(channelList)) {
+ exprFun <- channelList[[ch]]@exprFun[[1]] #!
+ varFun <- channelList[[ch]]@varFun[[1]] #!
+ values2 <- values[[ch]]
+ naVals <- which(is.na(values2) | is.infinite(values2))
+ pIDs2 <- pIDs
+ wts2 <- wts
+ if (length(naVals) > 0) {
+ values2 <- values2[-naVals]
+ pIDs2 <- pIDs[-naVals]
+ wts2 <- wts[-naVals]
+ }
+ pOrder <- order(pIDs2)
+ pIDs2 <- pIDs2[pOrder]
+ values2 <- values2[pOrder]
+ wts2 <- wts2[pOrder]
+ if (any(wts2 == 0)) {
+ values2 <- values2[-which(wts2 == 0)]
+ pIDs2 <- pIDs2[-which(wts2 == 0)]
+ wts2 <- wts2[-which(wts2 == 0)]
+ }
+ ## spezieller Teil für rmx
+ tmp <- split(wts2 * values2, pIDs2)
+ pMap <- match(names(tmp), probeIDs)
+ cat("Using rmx\n")
+ res.rmx <- rmxBeadSummary(x = tmp, eps = eps, eps.lower = eps.lower, eps.upper= eps.upper,
+ steps = steps, fsCor = fsCor, mad0 = mad0)
+ output[[ch]][["eMat"]][pMap, s] <- res.rmx$eMat
+ output[[ch]][["varMat"]][pMap, s] <- res.rmx$varMat
+ output[[ch]][["nObs"]][pMap, s] <- res.rmx$nObs
}
- else if ((imagesPerArray == 2) && (j < len)) {
- sel1 <- getBeadData(BLData, what = "ProbeID", array = arraynms[j]) != 0
- sel2 <- getBeadData(BLData, what = "ProbeID", array = arraynms[j + 1]) != 0
- pr <- append(getBeadData(BLData, what = "ProbeID", array = arraynms[j])[sel1],
- getBeadData(BLData, what = "ProbeID", array = arraynms[j + 1])[sel2])
- finten <- append(getBeadData(BLData, what = what, array = arraynms[j])[sel1],
- getBeadData(BLData, what = what, array = arraynms[j + 1])[sel2])
- if(log) finten <- log2(finten)
- nasinf <- !is.finite(finten) | is.na(finten)
- pr <- pr[!nasinf]
- finten <- finten[!nasinf]
+ }
+ cat("Making summary object\n")
+ eMat <- output[[1]][["eMat"]]
+ varMat <- output[[1]][["varMat"]]
+ nObs <- output[[1]][["nObs"]]
+ if (length(output) > 1) {
+ for (i in 2:length(output)) {
+ eMat <- cbind(eMat, output[[i]][["eMat"]])
+ varMat <- cbind(varMat, output[[i]][["varMat"]])
+ nObs <- cbind(nObs, output[[i]][["nObs"]])
}
}
- GBeadStDev <- GBeadStDev/sqrt(GNoBeads)
- if(!is.null(R)) RBeadStDev <- RBeadStDev/sqrt(RNoBeads)
- if (whatelse == "R") {
- rownames(G) <- rownames(R) <- rownames(GBeadStDev) <- rownames(RBeadStDev) <- rownames(GNoBeads) <- rownames(RNoBeads) <- probes
- BSData <- new("NChannelSet", R = R, G = G, GBeadStDev = GBeadStDev,
- RBeadStDev = RBeadStDev, GNoBeads = GNoBeads, RNoBeads = RNoBeads)
+ channelFac <- NULL
+ for (i in 1:length(channelList)) {
+ newfac <- cNames[i]
+ channelFac <- c(channelFac, rep(newfac, length(sList)))
}
+ BSData <- new("ExpressionSetIllumina")
+ annoName <- annotation(BLData)
+ if (!is.null(annoName)) {
+ annoLoaded <- require(paste("illumina", annoName, ".db", sep = ""), character.only = TRUE)
+ if (annoLoaded) {
+ mapEnv <- as.name(paste("illumina", annoName, "ARRAYADDRESS", sep = ""))
+ IlluminaIDs <- as.character(unlist(AnnotationDbi::mget(as.character(probeIDs), revmap(eval(mapEnv)), ifnotfound = NA)))
+ rownames(eMat) <- rownames(varMat) <- rownames(nObs) <- as.character(IlluminaIDs)
+ status <- rep("Unknown", length(probeIDs))
+ annoPkg <- paste("illumina", annoName, ".db", sep = "")
+ annoVers <- packageDescription(annoPkg, field = "Version")
+ message(paste("Annotating control probes using package ", annoPkg, " Version:", annoVers, "\n", sep = ""))
+ mapEnv <- as.name(paste("illumina", annoName, "REPORTERGROUPNAME", sep = ""))
+ t <- try(eval(mapEnv), silent = TRUE)
+ if (class(t) == "try-error") {
+ message(paste("Could not find a REPORTERGROUPNAME mapping in annotation package ", annoPkg,
+ ". Perhaps it needs updating?", sep = ""))
+ }
+ else {
+ status[which(!is.na(IlluminaIDs))] <- unlist(AnnotationDbi::mget(IlluminaIDs[which(!is.na(IlluminaIDs))],
+ eval(mapEnv), ifnotfound = NA))
+ status[which(is.na(status))] <- "regular"
+ }
+ featureData(BSData) <- new("AnnotatedDataFrame", data = data.frame(ArrayAddressID = probeIDs,
+ IlluminaID = IlluminaIDs, Status = status, row.names = IlluminaIDs))
+ BSData at annotation <- annoName
+ }
+ }
else {
- BSData <- new("ExpressionSetIllumina")
- assayData(BSData) <- assayDataNew(exprs = G, se.exprs = GBeadStDev,
- nObservations = GNoBeads, storage.mode = "list")
- rownames(exprs(BSData)) <- rownames(se.exprs(BSData)) <- rownames(nObservations(BSData)) <- probes
- featureData(BSData) <- new("AnnotatedDataFrame", data = data.frame(ProbeID = probes, row.names = probes))
+ cat("Could not map ArrayAddressIDs: No annotation specified\n")
+ featureData(BSData) <- new("AnnotatedDataFrame", data = data.frame(ProbeID = probeIDs,
+ row.names = probeIDs))
}
- if (nrow(pData(BLData)) == len) {
- if (imagesPerArray == 2)
- BSData at phenoData <- new("AnnotatedDataFrame", data = pData(BLData at phenoData)[arrayord, , drop = FALSE][seq(1, len, by = 2), , drop = FALSE])
- else BSData at phenoData <- BLData at phenoData
+ assayData(BSData) <- assayDataNew(exprs = eMat, se.exprs = varMat,
+ nObservations = nObs, storage.mode = "list")
+ sampInfo <- beadarray:::sampleSheet(BLData)
+ if (!is.null(sampInfo)) {
+ expIDs <- paste(sampInfo$Sentrix_ID, sampInfo$Sentrix_Position, sep = "_")
+ sampInfo <- sampInfo[sapply(newNames, function(x) grep(strtrim(x, 12), expIDs)), ]
+ p <- new("AnnotatedDataFrame", data.frame(sampInfo, row.names = newNames))
}
- else {
- BSData at phenoData <- new("AnnotatedDataFrame", data = data.frame(sampleName = colnames(G)))
+ else p <- new("AnnotatedDataFrame", data.frame(sampleID = newNames, SampleFac = unique(sampleFac),
+ row.names = newNames))
+ phenoData(BSData) <- p
+ qcNames <- names(BLData at sectionData)
+ qcNames <- setdiff(qcNames, "Targets")
+ if (length(qcNames) > 0) {
+ QC <- BLData at sectionData[[qcNames[[1]]]]
+ if (length(qcNames > 1)) {
+ for (i in 2:length(qcNames)) {
+ QC <- cbind(QC, BLData at sectionData[[qcNames[i]]])
+ }
+ }
+ QC <- new("AnnotatedDataFrame", data.frame(QC, row.names = arraynms))
+ BSData at QC <- QC
}
- if (!is.null(pData(BSData)$sampleName))
- sampleNames(BSData) <- pData(BSData)$sampleName
- else sampleNames(BSData) <- colnames(G)
- if (whatelse == "R") {
- varMetadata <- data.frame(labelDescription = colnames(BSData at phenoData@data), channel = "_ALL_")
- BSData at phenoData <- new("AnnotatedDataFrame", data = data.frame(BSData at phenoData@data), varMetadata = varMetadata)
- }
- BSData at annotation <- BLData at annotation
- if ("qcScores" %in% names(BLData at arrayInfo))
- t <- try(BSData at BeadLevelQC <- BLData at arrayInfo$qcScores, silent = TRUE)
+ BSData at channelData <- list(channelFac, channelList)
BSData
})
## computation of bead summaries via robloxbioc for "matrix"
-rmxBeadSummary <- function(x, probeIDs, probes, eps, eps.lower, eps.upper, steps, fsCor, mad0){
- comIDs <- intersect(probeIDs, probes)
- x <- x[probeIDs %in% comIDs]
- probeIDs <- probeIDs[probeIDs %in% comIDs]
- noBeads <- as.vector(table(probeIDs))
- noBeads.uni <- as.integer(names(table(noBeads)))
- probes1 <- comIDs
- len1 <- length(probes1)
- fore1 <- numeric(len1)
- SD1 <- numeric(len1)
- for(i in seq(along = noBeads.uni)){
- index <- noBeads == noBeads.uni[i]
- IDs <- probes1[index]
- if(noBeads.uni[i] == 1){
- fore1[index] <- x[probeIDs %in% IDs]
- SD1[index] <- mad0
+rmxBeadSummary <- function(x, eps, eps.lower, eps.upper, steps, fsCor, mad0){
+ nObs <- sapply(x, length)
+ noBeads <- as.integer(names(table(nObs)))
+ varMat <- eMat <- numeric(length(nObs))
+ for(i in seq(along = noBeads)){
+ index <- nObs == noBeads[i]
+ if(noBeads[i] == 1){
+ eMat[index] <- as.vector(unlist(x[index]))
+ varMat[index] <- mad0
}else{
- temp <- matrix(x[probeIDs %in% IDs], ncol = noBeads.uni[i], byrow = TRUE)
+ temp <- t(sapply(x[index], rbind))
rmx <- robloxbioc(temp, eps = eps, eps.lower = eps.lower, eps.upper = eps.upper,
steps = steps, fsCor = fsCor, mad0 = mad0)
- fore1[index] <- rmx[,"mean"]
- SD1[index] <- rmx[,"sd"]
+ eMat[index] <- rmx[,"mean"]
+ varMat[index] <- rmx[,"sd"]
}
- }
- len <- length(probes)
- fore <- numeric(len)
- SD <- numeric(len)
- noBeads1 <- numeric(len)
- nas <- !(probes %in% comIDs)
- fore[nas] <- NA
- fore[!nas] <- fore1
- SD[nas] <- NA
- SD[!nas] <- SD1
- noBeads1[nas] <- 0
- noBeads1[!nas] <- noBeads
-
- return(list(fore = fore, sd = SD, noBeads = noBeads1))
+ }
+ list(eMat = eMat, varMat = varMat, nObs = nObs)
}
Modified: pkg/RobLoxBioC/R/robloxbiocMatrix.R
===================================================================
--- pkg/RobLoxBioC/R/robloxbiocMatrix.R 2012-09-11 17:36:21 UTC (rev 510)
+++ pkg/RobLoxBioC/R/robloxbiocMatrix.R 2012-09-12 11:50:05 UTC (rev 511)
@@ -7,7 +7,6 @@
fsCor = TRUE, mad0 = 1e-4){
stopifnot(is.numeric(x))
if(ncol(x) <= 2){
- warning("Sample size <= 2! => Median and MAD are used for estimation.")
mean <- rowMedians(x, na.rm = TRUE)
sd <- rowMedians(abs(x-mean), na.rm = TRUE)/qnorm(0.75)
robEst <- cbind(mean, sd)
@@ -60,7 +59,7 @@
if(!is.null(eps)){
r <- sqrt(ncol(x))*eps
- if(fsCor) r <- .finiteSampleCorrection(r = r, n = ncol(x))
+ if(fsCor) r <- finiteSampleCorrection(r = r, n = ncol(x), model = "locsc")
if(r > 10){
b <- sd*1.618128043
const <- 1.263094656
@@ -76,7 +75,7 @@
mse <- A1 + A2
}
robEst <- .kstep.locsc.matrix(x = x, initial.est = cbind(mean, sd),
- A1 = A1, A2 = A2, a = a, b = b, k = steps)
+ A1 = A1, A2 = A2, a = a, b = b, k = steps)$est
colnames(robEst) <- c("mean", "sd")
}else{
sqrtn <- sqrt(ncol(x))
@@ -88,7 +87,7 @@
r <- uniroot(.getlsInterval, lower = rlo+1e-8, upper = rup,
tol = .Machine$double.eps^0.25, rlo = rlo, rup = rup)$root
}
- if(fsCor) r <- .finiteSampleCorrection(r = r, n = ncol(x))
+ if(fsCor) r <- finiteSampleCorrection(r = r, n = ncol(x), model = "locsc")
if(r > 10){
b <- sd*1.618128043
const <- 1.263094656
@@ -115,56 +114,9 @@
}
}
robEst <- .kstep.locsc.matrix(x = x, initial.est = cbind(mean, sd),
- A1 = A1, A2 = A2, a = a, b = b, k = steps)
+ A1 = A1, A2 = A2, a = a, b = b, k = steps)$est
colnames(robEst) <- c("mean", "sd")
}
return(robEst)
})
-###############################################################################
-## computation of k-step construction in case x is a matrix
-###############################################################################
-.onestep.locsc.matrix <- function(x, initial.est, A1, A2, a, b){
- mean <- initial.est[,1]
- sd <- initial.est[,2]
- u <- A1*(x-mean)/sd^2
- v <- A2*(((x-mean)/sd)^2-1)/sd - a
- ind <- b/sqrt(u^2 + v^2) <= 1
- IC1 <- rowMeans(u*(ind*b/sqrt(u^2 + v^2) + !ind), na.rm = TRUE)
- IC2 <- rowMeans(v*(ind*b/sqrt(u^2 + v^2) + !ind), na.rm = TRUE)
- IC <- cbind(IC1, IC2)
- return(initial.est + IC)
-}
-.kstep.locsc.matrix <- function(x, initial.est, A1, A2, a, b, mean, k){
- est <- .onestep.locsc.matrix(x = x, initial.est = initial.est, A1 = A1, A2 = A2, a = a, b = b)
- if(k > 1){
- for(i in 2:k){
- A1 <- est[,2]^2*A1/initial.est[,2]^2
- A2 <- est[,2]^2*A2/initial.est[,2]^2
- a <- est[,2]*a/initial.est[,2]
- b <- est[,2]*b/initial.est[,2]
- initial.est <- est
- est <- .onestep.locsc.matrix(x = x, initial.est = est, A1 = A1, A2 = A2, a = a, b = b)
- }
- }
- A1 <- est[,2]^2*A1/initial.est[,2]^2
- A2 <- est[,2]^2*A2/initial.est[,2]^2
- a <- est[,2]*a/initial.est[,2]
- b <- est[,2]*b/initial.est[,2]
-
- return(est)
-}
-.finiteSampleCorrection <- function(r, n){
- if(r >= 1.74) return(r)
-
- eps <- r/sqrt(n)
- ns <- c(3:50, seq(55, 100, by = 5), seq(110, 200, by = 10),
- seq(250, 500, by = 50))
- epss <- c(seq(0.001, 0.01, by = 0.001), seq(0.02, to = 0.5, by = 0.01))
- if(n %in% ns){
- ind <- ns == n
- }else{
- ind <- which.min(abs(ns-n))
- }
- return(max(r, approx(x = epss, y = .finiteSampleRadius.locsc[,ind], xout = eps, rule = 2)$y))
-}
Deleted: pkg/RobLoxBioC/R/sysdata.rda
===================================================================
(Binary files differ)
Modified: pkg/RobLoxBioC/inst/CITATION
===================================================================
--- pkg/RobLoxBioC/inst/CITATION 2012-09-11 17:36:21 UTC (rev 510)
+++ pkg/RobLoxBioC/inst/CITATION 2012-09-12 11:50:05 UTC (rev 511)
@@ -4,16 +4,16 @@
citHeader("To cite package RobLoxBioC in publications use:")
-citEntry(entry="Manual",
- title = "RobLoxBioC: Infinitesimally robust estimators for preprocessing omics data",
- author = personList(as.person("M. Kohl")),
+citEntry(entry="Article",
+ title = "Preprocessing of gene expression data by optimally robust estimators",
+ author = personList(as.person("M. Kohl"), as.person("H.P. Deigner")),
+ journal = "BMC Bioinformatics",
language = "English",
- year = year,
- note = note,
- type = "R package",
- url = "http://robast.r-forge.r-project.org/",
- textVersion = paste("Kohl, M.",
- sprintf("(%s).", year),
- "RobLoxBioC: Infinitesimally robust estimators for preprocessing omics data.",
- paste(note, ".", sep = ""),
- "URL http://robast.r-forge.r-project.org/"))
+ year = 2010,
+ Volume="11",
+ Pages="583",
+ url = "http://www.biomedcentral.com/1471-2105/11/583/abstract",
+ textVersion = paste("Kohl, M. and Deigner, H.P.",
+ "Preprocessing of gene expression data by optimally robust estimators.",
+ "BMC Bioinformatics 2010, 11:583",
+ "URL http://www.biomedcentral.com/1471-2105/11/583/abstract"))
Modified: pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd
===================================================================
--- pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd 2012-09-11 17:36:21 UTC (rev 510)
+++ pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd 2012-09-12 11:50:05 UTC (rev 511)
@@ -12,8 +12,8 @@
\details{
\tabular{ll}{
Package: \tab RobLoxBioC \cr
-Version: \tab 0.8.1 \cr
-Date: \tab 2012-01-16 \cr
+Version: \tab 0.8.3 \cr
+Date: \tab 2012-09-12 \cr
Depends: \tab R(>= 2.14.0), methods, Biobase, affy, beadarray, distr, RobLox,
lattice, RColorBrewer \cr
LazyLoad: \tab yes \cr
@@ -32,13 +32,18 @@
Kohl, M. (2005) \emph{Numerical Contributions to the Asymptotic Theory of Robustness}.
Bayreuth: Dissertation.
- Kohl M. and Deigner H.P. (2009). Using infinitesimally robust estimators for
- preprocessing gene expression data. In preparation.
+ Kohl M. and Deigner H.P. (2010). Preprocessing of gene expression data by optimally
+ robust estimators. \emph{BMC Bioinformatics}, 11:583.
+ M. Kohl, P. Ruckdeschel, and H. Rieder (2010). Infinitesimally Robust Estimation
+ in General Smoothly Parametrized Models. \emph{Statistical Methods and Application},
+ \bold{19}(3):333-354.
+
Rieder, H. (1994) \emph{Robust Asymptotic Statistics}. New York: Springer.
Rieder, H., Kohl, M. and Ruckdeschel, P. (2008) The Costs of not Knowing
- the Radius. Statistical Methods and Applications \emph{17}(1) 13-40.
+ the Radius. \emph{Statistical Methods and Applications} \bold{17}(1) 13-40.
+ Extended version: \url{http://www.stamats.de/RRlong.pdf}
}
\seealso{
\code{\link[RobLox]{roblox}}, \code{\link[RobLox]{rowRoblox}}
Modified: pkg/RobLoxBioC/man/robloxbioc.Rd
===================================================================
--- pkg/RobLoxBioC/man/robloxbioc.Rd 2012-09-11 17:36:21 UTC (rev 510)
+++ pkg/RobLoxBioC/man/robloxbioc.Rd 2012-09-12 11:50:05 UTC (rev 511)
@@ -22,9 +22,10 @@
mad0 = 1e-4, contrast.tau = 0.03, scale.tau = 10,
delta = 2^(-20), sc = 500)
-\S4method{robloxbioc}{beadLevelData}(x, log = TRUE, imagesPerArray = 1, what = "G", probes = NULL,
- arrays = NULL, eps = NULL, eps.lower = 0, eps.upper = 0.05,
- steps = 3L, fsCor = TRUE, mad0 = 1e-4)
+\S4method{robloxbioc}{beadLevelData}(x, channelList = list(greenChannel), probeIDs = NULL,
+ useSampleFac = FALSE, sampleFac = NULL, weightNames = "wts",
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 511
More information about the Robast-commits
mailing list