#' @name getiClusterBayes
#' @title Get subtypes from iClusterBayes
#' @description This function wraps the iClusterBayes (Integrative clustering by Bayesian latent variable model) algorithm and provides standard output for `getMoHeatmap()` and `getConsensusMOIC()`.
#' @param data List of matrices with maximum of 6 subdatasets.
#' @param N.clust Number of clusters.
#' @param type Data type corresponding to the list of matrics, which can be gaussian, binomial or possion.
#' @param n.burnin An integer value to indicate the number of MCMC burnin.
#' @param n.draw An integer value to indicate the number of MCMC draw.
#' @param prior.gamma A numerical vector to indicate the prior probability for the indicator variable gamma of each subdataset.
#' @param sdev A numerical value to indicate the standard deviation of random walk proposal for the latent variable.
#' @param thin A numerical value to thin the MCMC chain in order to reduce autocorrelation.
#' @export
#' @return A list with the following components:
#'
#' \code{fit} an object returned by \link[iClusterPlus]{iClusterBayes}.
#'
#' \code{clust.res} a data.frame storing sample ID and corresponding clusters.
#'
#' \code{feat.res} the results of features selection process.
#'
#' \code{mo.method} a string value indicating the method used for multi-omics integrative clustering.
#' @import iClusterPlus
#' @importFrom dplyr %>%
#' @examples # There is no example and please refer to vignette.
#' @references Mo Q, Shen R, Guo C, Vannucci M, Chan KS, Hilsenbeck SG (2018). A fully Bayesian latent variable model for integrative clustering analysis of multi-type omics data. Biostatistics, 19(1):71-86.
getiClusterBayes <- function(data = NULL,
N.clust = NULL,
type = rep("gaussian", length(data)),
n.burnin = 18000,
n.draw = 12000,
prior.gamma = rep(0.5,length(data)),
sdev = 0.05,
thin = 3) {
# check data
n_dat <- length(data)
if(n_dat > 6){
stop('current verision of MOVICS can support up to 6 datasets.')
}
if(n_dat < 2){
stop('current verision of MOVICS needs at least 2 omics data.')
}
# remove features that made of categories not equal to 2
if(is.element("binomial",type)) {
bindex <- which(type == "binomial")
for (i in bindex) {
a <- which(rowSums(data[[i]]) == 0)
b <- which(rowSums(data[[i]]) == ncol(data[[i]]))
if(length(a) > 0) {
data[[i]] <- data[[i]][which(rowSums(data[[i]]) != 0),] # remove all zero
}
if(length(b) > 0) {
data[[i]] <- data[[i]][which(rowSums(data[[i]]) != ncol(data[[i]])),] # remove all one
}
if(length(a) + length(b) > 0) {
message(paste0("--", names(data)[i],": a total of ",length(a) + length(b), " features were removed due to the categories were not equal to 2!"))
}
}
}
data.backup <- data
data <- lapply(data, t)
if(n_dat == 1){
res <- iClusterBayes(dt1 = data[[1]],
type = type,
K = N.clust - 1,
n.burnin = n.burnin,
n.draw = n.draw,
prior.gamma = prior.gamma,
sdev = sdev,
thin = thin)
}
if(n_dat == 2){
res <- iClusterBayes(dt1 = data[[1]],
dt2 = data[[2]],
type = type,
K = N.clust - 1,
n.burnin = n.burnin,
n.draw = n.draw,
prior.gamma = prior.gamma,
sdev = sdev,
thin = thin)
}
if(n_dat == 3){
res <- iClusterBayes(dt1 = data[[1]],
dt2 = data[[2]],
dt3 = data[[3]],
type = type,
K = N.clust - 1,
n.burnin = n.burnin,
n.draw = n.draw,
prior.gamma = prior.gamma,
sdev = sdev,
thin = thin)
}
if(n_dat == 4){
res <- iClusterBayes(dt1 = data[[1]],
dt2 = data[[2]],
dt3 = data[[3]],
dt4 = data[[4]],
type = type,
K = N.clust - 1,
n.burnin = n.burnin,
n.draw = n.draw,
prior.gamma = prior.gamma,
sdev = sdev,
thin = thin)
}
if(n_dat == 5){
res <- iClusterBayes(dt1 = data[[1]],
dt2 = data[[2]],
dt3 = data[[3]],
dt4 = data[[4]],
dt5 = data[[5]],
type = type,
K = N.clust - 1,
n.burnin = n.burnin,
n.draw = n.draw,
prior.gamma = prior.gamma,
sdev = sdev,
thin = thin)
}
if(n_dat == 6){
res <- iClusterBayes(dt1 = data[[1]],
dt2 = data[[2]],
dt3 = data[[3]],
dt4 = data[[4]],
dt5 = data[[5]],
dt6 = data[[6]],
type = type,
K = N.clust - 1,
n.burnin = n.burnin,
n.draw = n.draw,
prior.gamma = prior.gamma,
sdev = sdev,
thin = thin)
}
message("clustering done...")
clustres <- data.frame(samID = rownames(data[[1]]),
clust = res$clusters,
row.names = rownames(data[[1]]),
stringsAsFactors = FALSE)
#clustres <- clustres[order(clustres$clust),]
featres <- data.frame(feature = as.character(unlist(lapply(data.backup, rownames))),
dataset = rep(names(data),as.numeric(sapply(data.backup, function(x) length(rownames(x))))),
post.prob = unlist(res$beta.pp),
stringsAsFactors = FALSE)
feat.res <- NULL
for (d in unique(featres$dataset)) {
tmp <- featres[which(featres$dataset == d),]
tmp <- tmp[order(tmp$post.prob, decreasing = TRUE),]
tmp$rank <- 1:nrow(tmp)
feat.res <- rbind.data.frame(feat.res,tmp)
}
message("feature selection done...")
return(list(fit = res, clust.res = clustres, feat.res = feat.res, mo.method = "iClusterBayes"))
}
Datasets
Models
