Nothing
PoisMixClusWrapper <- function(y, gmin=1, gmax, conds, norm="TMM",
gmin.init.type = "small-em", init.runs=1, init.iter=10, split.init=TRUE,
alg.type="EM", cutoff = 10e-6, iter = 1000,
fixed.lambda=NA, equal.proportions = FALSE, verbose = FALSE,
interpretation = "sum", EM.verbose=FALSE, subset.index = NA)
{
if(is.matrix(y) == FALSE & is.data.frame(y) == FALSE)
stop(paste(sQuote("y"), "must be a matrix"))
if(min(y) < 0 | sum(as.numeric(round(y))) != sum(as.numeric(y)))
stop(paste(sQuote("y"), "must be a matrix made up of nonnegative counts"))
if(min(rowSums(y)) == 0)
stop(paste("at least one observation in", sQuote("y"),
"contains all 0's and must be removed from the data"))
if(is.vector(conds) == FALSE | length(conds) != ncol(y))
stop(paste(sQuote("conds"), "must be a vector the same length as the number of columns in", sQuote("y")))
if(gmin.init.type != "small-em" & gmin.init.type != "kmeans" & gmin.init.type != "split.small-em")
stop(paste(sQuote("gmin.init.type"), "must be one of", dQuote("small-em"), "or",
dQuote("kmeans"), "or", dQuote("split.small-em")))
if(alg.type != "EM" & alg.type != "CEM")
stop(paste(sQuote("alg.type"), "must be one of", dQuote("EM"), "or", dQuote("CEM")))
if(length(alg.type) > 1)
stop(paste(sQuote("alg.type"), "must be one of", dQuote("EM"), "or", dQuote("CEM")))
if(is.logical(verbose) == FALSE)
stop(paste(sQuote("verbose"), "must be", dQuote("TRUE"), "or", dQuote("FALSE")))
if(!inherits(fixed.lambda,"list") & is.na(fixed.lambda[1]) == FALSE)
stop(paste(sQuote("fixed.lambda"), "must be", dQuote("NA") , "or a list."))
## Grouping columns of y in order of condition (all replicates put together)
o.ycols <- order(conds)
y <- y[,o.ycols]
conds <- conds[o.ycols]
conds.names <- unique(conds)
d <- length(unique(conds))
r <- as.vector(table(conds))
if(length(rownames(y)) == 0) rn <- 1:nrow(y);
if(length(rownames(y)) > 0) rn <- rownames(y);
y <- as.matrix(y, nrow = nrow(y), ncol = ncol(y))
rownames(y) <- rn;
if(length(norm) != 1 & length(norm) != length(conds)) stop(paste(sQuote("norm"), "must be one of
the following: none, TC, UQ, Med, DESeq, TMM, or a vector of the same length as", sQuote("conds")))
## If estimated from data, all genes should be used
if(length(norm) == 1) {
if(norm == "none") s <- rep(1, cols);
if(norm == "TC") s <- colSums(y) / sum(as.numeric(y));
if(norm == "UQ") s <- apply(y, 2, quantile, 0.75) / sum(apply(y, 2, quantile, 0.75));
if(norm == "Med") s <- apply(y, 2, median) / sum(apply(y, 2, median));
if(norm == "DESeq") {
## Code from DESeq, v1.8.3
loggeomeans <- rowMeans(log(y))
s <- apply(y, 2, function(x) exp(median((log(x)-loggeomeans)[is.finite(loggeomeans)])))
s <- s / sum(s)
}
if(norm == "TMM") {
f <- calcNormFactors(as.matrix(y), method = "TMM")
s <- colSums(y)*f / sum(colSums(y)*f)
}
}
if(length(norm) == length(conds)) {
s <- norm / sum(norm)
}
## In case only a subset of data are to be used for analysis
if(is.na(subset.index)[1] == FALSE) {
y <- y[subset.index,]
n <- dim(y)[1];cols <- dim(y)[2]
w <- rowSums(y)
}
if(is.na(subset.index)[1] == TRUE) {
n <- dim(y)[1];cols <- dim(y)[2]
w <- rowSums(y)
}
all.results <- vector("list", length = gmax - gmin + 1)
names(all.results) <- paste("g=", seq(gmin,gmax, 1), sep = "")
## For gmin, run PoisMixClus with regular small-EM initialization
cat("Running g =", gmin, "...\n")
all.results[[1]] <- PoisMixClus(y = y, g = gmin, norm = s, conds = conds,
init.type = gmin.init.type, init.runs = init.runs, init.iter = init.iter,
alg.type = alg.type, cutoff = cutoff, iter = iter,
fixed.lambda = fixed.lambda, equal.proportions = equal.proportions,
prev.labels = NA, prev.probaPost = NA, verbose = verbose,
interpretation = interpretation,
EM.verbose = EM.verbose, wrapper=TRUE, subset.index=subset.index)
## For g > gmin, run PoisMixClus with Panos-like init using previous results
index <- 2
if(gmax > gmin) {
if(split.init == TRUE) {
for(K in seq((gmin+1),gmax,1)) {
cat("Running g =", K, "...\n")
prev.labels <- all.results[[index-1]]$labels
prev.probaPost <- all.results[[index-1]]$probaPost
all.results[[index]] <- PoisMixClus(y = y, g = K, norm=s,
conds = conds, init.type = "split.small-em",
alg.type = alg.type, cutoff = cutoff, iter = iter,
fixed.lambda = fixed.lambda,
equal.proportions = equal.proportions,
prev.labels = prev.labels, prev.probaPost = prev.probaPost,
init.runs = init.runs, init.iter = init.iter, verbose = verbose,
interpretation = interpretation, EM.verbose = EM.verbose,
wrapper=TRUE, subset.index=subset.index)
index <- index + 1
}
}
if(split.init == FALSE) {
for(K in seq((gmin+1),gmax,1)) {
cat("Running g =", K, "...\n")
all.results[[index]] <- PoisMixClus(y = y, g = K,
norm = s, conds = conds,
init.type = gmin.init.type,
alg.type = alg.type, cutoff = cutoff, iter = iter,
fixed.lambda = fixed.lambda,
equal.proportions = equal.proportions,
prev.labels = NA, prev.probaPost = NA, init.runs = init.runs,
init.iter = init.iter, verbose = verbose, interpretation = interpretation,
EM.verbose = EM.verbose, wrapper=TRUE, subset.index=subset.index)
index <- index + 1
}
}
}
logLike.all <- unlist(lapply(all.results, function(x) x$log.like))
ICL.all <- unlist(lapply(all.results, function(x) x$ICL))
ICL.choose <- which(ICL.all == max(ICL.all, na.rm = TRUE))
select.results <- all.results[[ICL.choose]]
select.results$model.selection <- "ICL"
BIC.all <- unlist(lapply(all.results, function(x) x$BIC))
BIC.choose <- which(BIC.all == max(BIC.all, na.rm = TRUE))
select.results2 <- all.results[[BIC.choose]]
select.results2$model.selection <- "BIC"
# Apply capushe: only if at least 10 models are considered
if(gmax-gmin+1 <= 10) {
message("Note: the slope heuristics approach for model selection may only be applied if more than 10 models are fit.")
DDSE.results <- NA
Djump.results <- NA
capushe <- NA
ResCapushe <- NA
}
if(c(gmax-gmin+1) > 10) {
message("Note: diagnostic plots for results corresponding to model selection via slope heuristics (Djump and DDSE) should be examined to ensure that sufficiently complex models have been considered.")
Kchoice <- gmin:gmax
np <- (Kchoice-1) + (length(unique(conds))-1)*(Kchoice)
mat <- cbind(Kchoice, np/n, np/n, -logLike.all/n)
ResCapushe <- capushe(mat, n)
DDSE <- ResCapushe@DDSE@model
Djump <- ResCapushe@Djump@model
DDSE.results <- all.results[[paste("g=", DDSE, sep="")]]
Djump.results <- all.results[[paste("g=", Djump, sep="")]]
DDSE.results$model.selection <- "DDSE"
Djump.results$model.selection <- "Djump"
}
RESULTS <- list(logLike.all = logLike.all, ICL.all = ICL.all,
capushe = ResCapushe,
all.results = all.results,
DDSE.results = DDSE.results,
Djump.results = Djump.results,
ICL.results = select.results,
BIC.results = select.results2)
class(RESULTS) <- "HTSClusterWrapper"
return(RESULTS)
}
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.