[Robast-commits] r582 - in branches/robast-0.9/pkg/RobLoxBioC: . R inst inst/scripts man tests/Examples
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 1 05:12:24 CET 2013
Author: stamats
Date: 2013-02-01 05:12:19 +0100 (Fri, 01 Feb 2013)
New Revision: 582
Added:
branches/robast-0.9/pkg/RobLoxBioC/R/00pre2160.R
branches/robast-0.9/pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R
Removed:
branches/robast-0.9/pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R
branches/robast-0.9/pkg/RobLoxBioC/R/sysdata.rda
branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/AffySimStudy.R
Modified:
branches/robast-0.9/pkg/RobLoxBioC/DESCRIPTION
branches/robast-0.9/pkg/RobLoxBioC/R/0AllGeneric.R
branches/robast-0.9/pkg/RobLoxBioC/R/AffySimStudyFunction.R
branches/robast-0.9/pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R
branches/robast-0.9/pkg/RobLoxBioC/R/robloxbiocMatrix.R
branches/robast-0.9/pkg/RobLoxBioC/inst/CITATION
branches/robast-0.9/pkg/RobLoxBioC/inst/NEWS
branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R
branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/IlluminaExample.R
branches/robast-0.9/pkg/RobLoxBioC/inst/scripts/IlluminaSimStudy.R
branches/robast-0.9/pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd
branches/robast-0.9/pkg/RobLoxBioC/man/KolmogorovMinDist.Rd
branches/robast-0.9/pkg/RobLoxBioC/man/SimStudies.Rd
branches/robast-0.9/pkg/RobLoxBioC/man/robloxbioc.Rd
branches/robast-0.9/pkg/RobLoxBioC/tests/Examples/RobLoxBioC-Ex.Rout.save
Log:
merged trunk into branch, update of Rout.save, had to put some examples in \dontrun to reduce check time
Modified: branches/robast-0.9/pkg/RobLoxBioC/DESCRIPTION
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/DESCRIPTION 2013-01-31 19:09:59 UTC (rev 581)
+++ branches/robast-0.9/pkg/RobLoxBioC/DESCRIPTION 2013-02-01 04:12:19 UTC (rev 582)
@@ -1,10 +1,10 @@
Package: RobLoxBioC
Version: 0.9
-Date: 2010-12-03
-Title: Optimally robust estimation for omics data
+Date: 2013-01-31
+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.
-Depends: R(>= 2.8.1), methods, Biobase, affy, beadarray, distr, RobLox, lattice, RColorBrewer
+Depends: R(>= 2.14.0), methods, Biobase, affy, beadarray, distr, RobLox, lattice, RColorBrewer
Author: Matthias Kohl <Matthias.Kohl at stamats.de>
Maintainer: Matthias Kohl <Matthias.Kohl at stamats.de>
LazyLoad: yes
@@ -14,4 +14,4 @@
Encoding: latin1
LastChangedDate: {$LastChangedDate$}
LastChangedRevision: {$LastChangedRevision$}
-SVNRevision: 439
+SVNRevision: 511
Copied: branches/robast-0.9/pkg/RobLoxBioC/R/00pre2160.R (from rev 580, pkg/RobLoxBioC/R/00pre2160.R)
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/R/00pre2160.R (rev 0)
+++ branches/robast-0.9/pkg/RobLoxBioC/R/00pre2160.R 2013-02-01 04:12:19 UTC (rev 582)
@@ -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: branches/robast-0.9/pkg/RobLoxBioC/R/0AllGeneric.R
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/R/0AllGeneric.R 2013-01-31 19:09:59 UTC (rev 581)
+++ branches/robast-0.9/pkg/RobLoxBioC/R/0AllGeneric.R 2013-02-01 04:12:19 UTC (rev 582)
@@ -1,8 +1,4 @@
############# preparations ################
-.onLoad <- function(lib, pkg) {
- require("methods", character = TRUE, quietly = TRUE)
-}
-
if(!isGeneric("robloxbioc")){
setGeneric("robloxbioc",
function(x, ...) standardGeneric("robloxbioc"))
Modified: branches/robast-0.9/pkg/RobLoxBioC/R/AffySimStudyFunction.R
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/R/AffySimStudyFunction.R 2013-01-31 19:09:59 UTC (rev 581)
+++ branches/robast-0.9/pkg/RobLoxBioC/R/AffySimStudyFunction.R 2013-02-01 04:12:19 UTC (rev 582)
@@ -38,12 +38,13 @@
if(plot2){
- ind <- sample(1:M, min(M, 20))
+ ind <- if(M <= 20) 1:M else sample(1:M, 20)
if(plot1) dev.new()
+ M1 <- min(M, 20)
print(
- stripplot(rep(1:20, each = 20) ~ as.vector(Mre[ind,]),
+ stripplot(rep(1:M1, each = n) ~ as.vector(Mre[ind,]),
ylab = "samples", xlab = "x", pch = 20,
- main = "Randomly chosen samples")
+ main = ifelse(M <= 20, "Samples", "20 randomly chosen samples"))
)
}
@@ -72,11 +73,10 @@
abline(h = 0)
boxplot(Ergebnis2, col = myCol[c(1,2,4)], pch = 20, main = "Scale")
abline(h = 1)
- op <- par(mar = rep(2, 4), no.readonly = TRUE)
+ op <- par(mar = rep(2, 4))
plot(c(0,1), c(1, 0), type = "n", axes = FALSE)
legend("center", c("ML", "Med/MAD", "biweight", "rmx"),
fill = myCol, ncol = 4, cex = 1.5)
-# op$cin <- op$cra <- op$csi <- op$cxy <- op$din <- NULL
on.exit(par(op))
}
Modified: branches/robast-0.9/pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R 2013-01-31 19:09:59 UTC (rev 581)
+++ branches/robast-0.9/pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R 2013-02-01 04:12:19 UTC (rev 582)
@@ -39,12 +39,13 @@
if(plot2){
- ind <- sample(1:M, min(M, 20))
+ ind <- if(M <= 20) 1:M else sample(1:M, 20)
if(plot1) dev.new()
+ M1 <- min(M, 20)
print(
- stripplot(rep(1:20, each = 20) ~ as.vector(Mre[ind,]),
+ stripplot(rep(1:M1, each = n) ~ as.vector(Mre[ind,]),
ylab = "samples", xlab = "x", pch = 20,
- main = "Randomly chosen samples")
+ main = ifelse(M <= 20, "Samples", "20 randomly chosen samples"))
)
}
@@ -79,11 +80,11 @@
abline(h = 0)
boxplot(Ergebnis2, col = myCol, pch = 20, main = "Scale")
abline(h = 1)
- op <- par(mar = rep(2, 4), no.readonly = TRUE)
+ op <- par(mar = rep(2, 4))
plot(c(0,1), c(1, 0), type = "n", axes = FALSE)
legend("center", c("ML", "Med/MAD", "Illumina", "rmx"),
fill = myCol, ncol = 4, cex = 1.5)
-# op$cin <- op$cra <- op$csi <- op$cxy <- op$din <- NULL
+ op$cin <- op$cra <- op$csi <- op$cxy <- op$din <- NULL
on.exit(par(op))
}
Copied: branches/robast-0.9/pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R (from rev 580, pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R)
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R (rev 0)
+++ branches/robast-0.9/pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R 2013-02-01 04:12:19 UTC (rev 582)
@@ -0,0 +1,252 @@
+setMethod("robloxbioc", signature(x = "beadLevelData"),
+ 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)
+ 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{
+ 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)
+ }
+ }
+ }
+ }
+ else {
+ cat("No sample factor specified. Summarizing each section separately\n")
+ sList <- arraynms
+ sampleFac <- arraynms
+ newNames <- arraynms
+ }
+ 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 (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")
+ }
+ }
+ 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")
+ }
+ }
+ }
+ }
+ 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])
+ }
+ 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)
+ }
+ }
+ 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
+ }
+ }
+ 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"]])
+ }
+ }
+ 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 {
+ cat("Could not map ArrayAddressIDs: No annotation specified\n")
+ featureData(BSData) <- new("AnnotatedDataFrame", data = data.frame(ProbeID = probeIDs,
+ row.names = probeIDs))
+ }
+ 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 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
+ }
+ BSData at channelData <- list(channelFac, channelList)
+ BSData
+ })
+## computation of bead summaries via robloxbioc for "matrix"
+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 <- t(sapply(x[index], rbind))
+ rmx <- robloxbioc(temp, eps = eps, eps.lower = eps.lower, eps.upper = eps.upper,
+ steps = steps, fsCor = fsCor, mad0 = mad0)
+ eMat[index] <- rmx[,"mean"]
+ varMat[index] <- rmx[,"sd"]
+ }
+ }
+ list(eMat = eMat, varMat = varMat, nObs = nObs)
+}
Deleted: branches/robast-0.9/pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R 2013-01-31 19:09:59 UTC (rev 581)
+++ branches/robast-0.9/pkg/RobLoxBioC/R/robloxbiocBeadLevelList.R 2013-02-01 04:12:19 UTC (rev 582)
@@ -1,196 +0,0 @@
-setMethod("robloxbioc", signature(x = "BeadLevelList"),
- 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){
- BLData <- x
- arraynms <- arrayNames(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"
- }else{
- stop("Need two-channel data to calculate summary R and G values")
- }
- }
- if(imagesPerArray == 1){
- pr <- getArrayData(BLData, what = "ProbeID", array = arraynms[1])
- sel <- pr != 0
- pr <- pr[sel]
- finten <- getArrayData(BLData, what = what, log = log, array = arraynms[1])[sel]
- nasinf <- !is.finite(finten) | is.na(finten)
- finten <- finten[!nasinf]
- }
- 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 <- getArrayData(BLData, what = "ProbeID", array = arraynms[1]) != 0
- sel2 <- getArrayData(BLData, what = "ProbeID", array = arraynms[2]) != 0
- pr <- append(getArrayData(BLData, what = "ProbeID", array = arraynms[1])[sel1],
- getArrayData(BLData, what = "ProbeID", array = arraynms[2])[sel2])
- finten <- append(getArrayData(BLData, what = what, log = log, array = arraynms[1])[sel1],
- getArrayData(BLData, what = what, log = log, array = arraynms[2])[sel2])
- 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(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
- }
- 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
- }
- 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 <- getArrayData(BLData, what = whatelse, log = log, array = arraynms[i])[sel]
- nasinf <- !is.finite(finten) | is.na(finten)
- finten <- finten[!nasinf]
- }
- else if (imagesPerArray == 2) {
- finten <- append(getArrayData(BLData, what = whatelse, log = log, array = arraynms[j])[sel1],
- getArrayData(BLData, what = whatelse, log = log, array = arraynms[j + 1])[sel2])
- nasinf <- !is.finite(finten) | is.na(finten)
- finten <- finten[!nasinf]
- }
- 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 <- getArrayData(BLData, what = "ProbeID", array = arraynms[i]) != 0
- pr <- getArrayData(BLData, what = "ProbeID", array = arraynms[i])[sel]
- finten <- getArrayData(BLData, what = what, log = log, array = arraynms[i])[sel]
- nasinf <- !is.finite(finten) | is.na(finten)
- pr <- pr[!nasinf]
- finten <- finten[!nasinf]
- }
- else if ((imagesPerArray == 2) && (j < len)) {
- sel1 <- getArrayData(BLData, what = "ProbeID", array = arraynms[j]) != 0
- sel2 <- getArrayData(BLData, what = "ProbeID", array = arraynms[j + 1]) != 0
- pr <- append(getArrayData(BLData, what = "ProbeID", array = arraynms[j])[sel1],
- getArrayData(BLData, what = "ProbeID", array = arraynms[j + 1])[sel2])
- finten <- append(getArrayData(BLData, what = what, log = log, array = arraynms[j])[sel1],
- getArrayData(BLData, what = what, log = log, array = arraynms[j + 1])[sel2])
- nasinf <- !is.finite(finten) | is.na(finten)
- pr <- pr[!nasinf]
- finten <- finten[!nasinf]
- }
- }
- 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)
- }
- else {
- BSData <- new("ExpressionSetIllumina")
- assayData(BSData) <- assayDataNew(exprs = G, se.exprs = GBeadStDev,
- NoBeads = GNoBeads, storage.mode = "list")
- rownames(exprs(BSData)) <- rownames(se.exprs(BSData)) <- rownames(NoBeads(BSData)) <- probes
- featureData(BSData) <- new("AnnotatedDataFrame", data = data.frame(ProbeID = probes, row.names = probes))
- }
- 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
- }
- else {
- BSData at phenoData <- new("AnnotatedDataFrame", data = data.frame(sampleName = colnames(G)))
- }
- 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
- })
-## 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
- }else{
- temp <- matrix(x[probeIDs %in% IDs], ncol = noBeads.uni[i], byrow = TRUE)
- 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"]
- }
- }
- 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))
-}
Modified: branches/robast-0.9/pkg/RobLoxBioC/R/robloxbiocMatrix.R
===================================================================
--- branches/robast-0.9/pkg/RobLoxBioC/R/robloxbiocMatrix.R 2013-01-31 19:09:59 UTC (rev 581)
+++ branches/robast-0.9/pkg/RobLoxBioC/R/robloxbiocMatrix.R 2013-02-01 04:12:19 UTC (rev 582)
@@ -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
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 582
More information about the Robast-commits
mailing list