Nothing
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@sectionData)) {
cat("Could not determine sample factor from beadLevelData. Summarizing each section separately\n")
sList <- arraynms
sampleFac <- arraynms
newNames <- sList
}else{
sampleFac <- BLData@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@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)]
}
## special part for 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, fields = "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 (inherits(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@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@sectionData)
qcNames <- setdiff(qcNames, "Targets")
if (length(qcNames) > 0) {
QC <- BLData@sectionData[[qcNames[[1]]]]
if (length(qcNames > 1)) {
for (i in 2:length(qcNames)) {
QC <- cbind(QC, BLData@sectionData[[qcNames[i]]])
}
}
QC <- new("AnnotatedDataFrame", data.frame(QC, row.names = arraynms))
BSData@QC <- QC
}
BSData@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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.