R/Prede.R

In libcell/deconvBench: A tool for approaches evaluation of cell fraction estimation

#' PREDE
#'
#' @param scdata single data with genes in rows and cells in columns.
#' @param bulkdata A matrix with genes in rows and samples in columns.
#' @param label The cell type label of single data.
#' @param celltypenum  Specify the range of underlying subpopulation number.
#' @param seed the random seed.
#' @importFrom dynwrap random_seed
#' @importFrom PREDE PREDE select_feature GetCelltypeNum
#' @return  A data frame of Mixed cellular fractions.
#' @export
#'
#' @examples
#'
#' Bulk <- Bulk_GSE60424
#' SC <- SC_GSE60424
#' Label <- Label_GSE60424$Label
#' res <- Prede(bulkdata = Bulk,
#'              scdata = SC,
#'              label = Label,
#'              celltypenum = 5,
#'              seed = 2022)
#'
#'

Prede <- function(scdata, bulkdata,label,celltypenum,seed=NULL) {

  gene <- intersect(rownames(scdata),rownames(bulkdata))

  scdata <- scdata[gene,]

  bulkdata <- bulkdata[gene,]

  factor_label <- label

  names(factor_label) <- colnames(scdata)

  all_cell_type <- unique(as.character(factor_label))

  pure_samples <- lapply(1:length(all_cell_type), function(i) {

    which(as.character(label) == all_cell_type[i])
  })


 names(pure_samples) = all_cell_type


  KT <- length(names(pure_samples))

  pure_expr <- sapply(1:KT, function(i){apply(scdata[,pure_samples[[i]]], 1, sum)})


  colnames(pure_expr) = all_cell_type


  bulk <- PREDE::generate_bulk(pure_expr,nSample =100,csd = 0.1)


  sigfeature <- PREDE::select_feature(mat = bulk$Y,
                                      method = "cv",
                                      nmarker = 1000,
                                      startn = 0)


  # Seed
  if(is.null(seed)){
    seed <- dynwrap::random_seed()
    cat(paste0("Seed is ", seed))
  }

  set.seed(seed)

  nn <- sample(c(1:ncol(pure_expr)),size=round(ncol(pure_expr)*3/4))

  ## determine the total number of cell types by AIC, by specifying only partial reference W1
  OptimalK <- PREDE::GetCelltypeNum(bulk$Y[sigfeature,],
                                    W1=pure_expr[sigfeature,nn],
                                    maxK = celltypenum)

  prede <- PREDE::PREDE(Y = as.matrix(bulkdata[sigfeature,]),  ###bulkdata
                       W1= pure_expr[sigfeature,nn],
                       type = "GE",
                       K = OptimalK$K,
                       iters = 100,
                       rssDiffStop = 1e-10)


  res_Prede <- prede$H

  return(res_Prede)

}