##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (segall, ratall = NULL, idcol = NULL, startcol = NULL,
endcol = NULL, medcol = NULL, madcol = NULL, errorcol = NULL,
chromcol = NULL, bpstartcol = NULL, bpendcol = NULL, annot = NULL,
annotstartcol = NULL, annotendcol = NULL, annotchromcol = NULL,
useend = F, blsize = NULL, minjoin = NULL, ntrial = 10, bestbic = -1e+07,
modelNames = "E", cweight = NULL, bstimes = NULL, chromrange = NULL,
myseed = NULL, distrib = c("vanilla", "Rparallel", "sge"),
njobs = 1, normalength = NULL, normalmedian = NULL, normalmad = NULL,
normalerror = NULL)
{
if (is.null(idcol)) {
cat("Found a single segmented profile with no ID", "")
if (!is.null(ratall)) {
if (sum(apply(ratall, 2, data.class) == "numeric") >
1)
stop("Ambiguity: more than 1 numeric column in raw data table
")
else {
idrat <- which(apply(ratall, 2, data.class) ==
"numeric")
segall <- data.frame(rep(as.character(idrat),
nrow(segall)), segall)
idcol <- "ID"
dimnames(segall)[[2]][1] <- idcol
}
}
}
if (is.null(ratall))
cat("No raw table, proceeding to comparison
")
else {
profnames <- unique(segall[, idcol])
if (!all(profnames %in% dimnames(ratall)[[2]]))
stop("Found unmatched segmented profile IDs
")
if (is.null(startcol) | is.null(endcol)) {
if (is.null(bpstartcol) | is.null(bpendcol) | is.null(chromcol))
stop("Unable to proceed: incomplete segment annotation
")
if (is.null(chromrange))
chromrange <- sort(unique(segall[, chromcol]))
if (is.null(annot))
stop("No annotation table; unable to determine boundary probes/bins
")
if (is.null(annotstartcol) | is.null(annotchromcol))
stop("No start and chrom column names provided for annotation table
")
if (useend & is.null(annotendcol))
stop("End column name required but not provided in annotation table
")
maxbpstart <- max(c(segall[, bpstartcol], annot[,
annotstartcol])) + 1
maxbpend <- ifelse(useend, max(c(segall[, bpendcol],
annot[, annotendcol])), max(c(segall[, bpendcol],
annot[, annotstartcol]))) + 1
startprobe <- match((segall[, chromcol] - 1) * maxbpstart +
segall[, bpstartcol], (annot[, annotchromcol] -
1) * maxbpstart + annot[, annotstartcol])
endprobe <- ifelse(useend, match((segall[, chromcol] -
1) * maxbpend + segall[, bpendcol], (annot[,
annotchromcol] - 1) * maxbpend + annot[, annotendcol]),
match((segall[, chromcol] - 1) * maxbpend + segall[,
bpendcol], (annot[, annotchromcol] - 1) * maxbpend +
annot[, annotstartcol]))
if (!all(!is.na(startprobe) & !is.na(endprobe)))
stop("Incomplete start and end annotation of segments
")
segall <- data.frame(segall, startprobe, endprobe)
startcol <- "StartProbe"
endcol <- "EndProbe"
}
profpack <- vector(mode = "list", length = length(profnames))
names(profpack) <- profnames
RNGkind("L'Ecuyer-CMRG")
if (!is.null(myseed))
set.seed(myseed)
s <- .Random.seed
distrib <- match.arg(distrib)
for (pn in profnames) {
profpack[[pn]] <- vector(mode = "list", length = 3)
names(profpack[[pn]]) <- c("seg", "rat", "rngseed")
profpack[[pn]]$seg <- segall[segall[, idcol] == pn,
c(startcol, endcol, chromcol), drop = F]
dimnames(profpack[[pn]]$seg)[[2]] <- c("StartProbe",
"EndProbe", "chrom")
profpack[[pn]]$rat <- ratall[, pn]
s <- nextRNGStream(s)
profpack[[pn]]$rngseed <- s
}
distrib <- match.arg(distrib)
if (distrib == "Rparallel") {
ncores <- min(njobs, length(profnames), detectCores())
cl <- parallel::makeCluster(getOption("cl.cores",
ncores))
parallel::clusterEvalQ(cl = cl, expr = library(parallel))
parallel::clusterEvalQ(cl = cl, expr = library(mclust))
}
if (distrib == "sge")
ncores <- min(njobs, length(profnames))
processed <- switch(distrib, vanilla = lapply(X = profpack,
FUN = CNclusterNcenter, blsize = blsize, minjoin = minjoin,
ntrial = ntrial, bestbic = bestbic, modelNames = modelNames,
cweight = cweight, bstimes = bstimes, chromrange = chromrange),
Rparallel = parLapply(cl, X = profpack, fun = CNclusterNcenter,
blsize = blsize, minjoin = minjoin, ntrial = ntrial,
bestbic = bestbic, modelNames = modelNames, cweight = cweight,
bstimes = bstimes, chromrange = chromrange),
sge = sge.parLapply(cl, X = profpack, FUN = CNclusterNcenter,
blsize = blsize, minjoin = minjoin, ntrial = ntrial,
bestbic = bestbic, modelNames = modelNames, cweight = cweight,
bstimes = bstimes, chromrange = chromrange, packages = c("parallel",
"mclust")))
segall <- cbind(segall, t(matrix(nrow = 6, data = unlist(lapply(processed,
t)))))
dimnames(segall)[[2]][(ncol(segall) - 5):ncol(segall)] <- c("segmedian",
"segmad", "mediandev", "segerr", "segz", "marginalprob")
medcol <- "mediandev"
madcol <- "segmad"
errorcol <- "segerr"
}
if (!(is.null(normalength) | is.null(normalmedian) | is.null(medcol))) {
if (is.null(bpstartcol) | is.null(bpendcol)) {
if (is.null(startcol) | is.null(endcol) | is.null(annot) |
is.null(annotstartcol) | is.null(annotendcol))
stop("Insufficient annotation for comparison")
tumorlength <- annot[segall[, endcol], annotendcol] -
annot[segall[, startcol], annotstartcol] + 1
}
else tumorlength <- segall[, bpendcol] - segall[, bpstartcol] +
1
tumormedian <- segall[, medcol]
if (!is.null(madcol))
tumormad <- segall[, madcol]
if (!is.null(errorcol))
tumorerror <- segall[, errorcol]
segall <- cbind(segall, normalComparison(normalmedian,
normalength, tumormedian, tumorlength, normalmad,
normalerror, tumormad, tumorerror))
}
return(segall)
}Run the code above in your browser using DataLab