R/getiClusterBayes.R

In xlucpu/MOVICS: Multi-Omics integration and VIsualization in Cancer Subtyping

#' @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"))
}