R/Cdseq_no.R

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

#' CDseq no reference
#'
#' @param bulkdata A matrix with genes in rows and samples in columns.
#' @param celltypenum  Specify the range of underlying subpopulation number, default is 5.
#' @importFrom CDSeq CDSeq
#' @return  A data frame of Mixed cellular fractions.
#'
#' @export
#'
#' @examples
#'
#' Bulk <- Bulk_GSE60424
#' res <- Cdseq_no(bulkdata = Bulk,
#'                 celltypenum = 5 )
#'
#'



Cdseq_no <- function(bulkdata,celltypenum=5) {

  cdseq.result <- CDSeq::CDSeq(bulk_data = bulkdata,
                               cell_type_number = celltypenum,
                               beta = 0.5,
                               alpha = 5,
                               mcmc_iterations = 200,
                               cpu_number = 1,
                               dilution_factor = 10)

  CDseq_compostion <- t(cdseq.result[["estProp"]])

  CDseq_sigmatrix <- cdseq.result[["estGEP"]]

  matrix <- as.data.frame(cdseq.result[["cellTypeAssignSplit"]])

  n =  dim( cdseq.result[["cellTypeAssignSplit"]])[3]

  type_matrix <- NULL

    for (c in 1:n ){

      cctype <- apply(matrix[,((c-1)*ncol(bulkdata)+1):(c*ncol(bulkdata))], 1, sum)

      type_matrix <- cbind(type_matrix,cctype)
    }



  marker_ann <- NULL
  for (i in 1:ncol(type_matrix)) {

    cellexp <-  type_matrix[,i]

    names(cellexp) <- rownames(type_matrix)

    cellexp <- sort(cellexp ,decreasing = T)[1:100]

    cellgene <- names(cellexp)

    marker_ann <- cbind(marker_ann ,cellgene)

  }

  ##########

  marker_ann <- as.data.frame(na.omit(marker_ann))

  anncell <- NULL

  probpath <- system.file("extdata", "PanglaoDB_markers_27_Mar_2020.tsv", package = "deconvBench")

  pl <- read.table(probpath,sep = "\t",header = T)

  Plhuman <- pl[pl$species!="Mm",]

  for (i in 1:ncol(marker_ann)) {

    marker <- marker_ann[,i]

    celltype <-NULL

    for (m in marker ) {

      ct <- Plhuman[which(Plhuman$official.gene.symbol==m),3]

      celltype <- c(celltype,ct)
    }

    uniquecell <- unique(celltype)

    times_cell <- tabulate(match(celltype, uniquecell))

    names(times_cell) <- uniquecell

    cell <- names(sort(times_cell,decreasing = T))[1]

    anncell <- c(anncell,cell)
  }

    colnames(CDseq_compostion) <- paste(anncell,1:n,sep = "_")

    res_Cdseq_no <- CDseq_compostion


    return(res_Cdseq_no)
}