Presto

Please refer https://github.com/immunogenomics/presto for details of presto #Run Preston Human Prefrontal Cortex

t.Presto.brain <- c()
for(i in 1:10){
  print(paste0("calculating time:",i))
  t1 <- Sys.time()
  wilcoxauc(snRNA.hpfc.sct, 'major_clust',logFC = 0.5)
  t2 <- Sys.time()
  t <- t2-t1
  t.Presto.brain <- append(t.Presto.brain,t)
}
t.test(t.Presto.brain,t.Star.brain)
mean(t.Presto.brain)
mean(t.Star.brain)
save(t.Presto.brain,file = "./out/rev_presto/t.Presto.brain.Rdata")

different of Presto and starTracer

time.1 <- as.numeric(t.Presto.brain)
time.2 <- as.numeric(t.Star.brain)
frame.time <- data.frame(Presto = time.1, starTracer = time.2)
frame.time$RunTime <- 1:10 %>% as.vector()
mat.use <- reshape2::melt(frame.time,id = "RunTime")
(p <- ggbarplot(mat.use, x = "variable", y = "value", add = c("mean_se"),color = "variable",fill = "white",palette = c("red","#123575"),stat = "identity",
          position = position_dodge(0.8)) + stat_compare_means(method = "t.test") + theme(axis.text.x = element_text(angle = 45,hjust = 1)))
ggsave(filename = "./out/rev_presto/time.brain.pdf",p,width = 5,height = 5)

Heart

t.Presto.heart <- c()
for(i in 1:10){
  print(paste0("calculating time:",i))
  t1 <- Sys.time()
  wilcoxauc(heart_big, 'cell_type__ontology_label',logFC = 0.5)
  t2 <- Sys.time()
  t <- t2-t1
  t.Presto.heart <- append(t.Presto.heart,t)
}

t.test(as.numeric(t.Presto.heart) * 60,as.numeric(t.Star.heart.new))$p.value
mean(t.Presto.heart)
mean(t.Star.heart.new)

time.1 <- as.numeric(t.Presto.heart) * 60
time.2 <- as.numeric(t.Star.heart.new)
frame.time <- data.frame(Presto = time.1, starTracer = time.2)
frame.time$RunTime <- 1:10 %>% as.vector()
mat.use <- reshape2::melt(frame.time,id = "RunTime")
(p <- ggbarplot(mat.use, x = "variable", y = "value", add = c("mean_se"),color = "variable",fill = "white",palette = c("red","#123575"),stat = "identity",
          position = position_dodge(0.8)) + stat_compare_means(method = "t.test") + theme(axis.text.x = element_text(angle = 45,hjust = 1)))
ggsave(filename = "./out/rev_presto/time.heart.pdf",p,width = 5,height = 5)

t.test(as.numeric(t.Presto.heart),as.numeric(t.Star.heart.new))$p.value

Kidney

t.Presto.kidney <- c()
for(i in 1:10){
  print(paste0("calculating time:",i))
  t1 <- Sys.time()
  wilcoxauc(sc_kidney, 'cell_type__ontology_label',logFC = 0.5)
  t2 <- Sys.time()
  t <- t2-t1
  t.Presto.kidney <- append(t.Presto.kidney,t)
}

t.test(t.Presto.kidney,t.Star.kidney.new)
mean(t.Presto.kidney)
mean(t.Star.kidney.new)

time.1 <- as.numeric(t.Presto.kidney)
time.2 <- as.numeric(t.Star.kidney.new)
frame.time <- data.frame(Presto = time.1, starTracer = time.2)
frame.time$RunTime <- 1:10 %>% as.vector()
mat.use <- reshape2::melt(frame.time,id = "RunTime")
(p <- ggbarplot(mat.use, x = "variable", y = "value", add = c("mean_se"),color = "variable",fill = "white",palette = c("red","#123575"),stat = "identity",
          position = position_dodge(0.8)) + stat_compare_means(method = "t.test") + theme(axis.text.x = element_text(angle = 45,hjust = 1)))
ggsave(filename = "./out/rev_presto/time.kidney.pdf",p,width = 5,height = 5)

t.test(as.numeric(t.Presto.kidney),as.numeric(t.Star.kidney.new))$p.value

PS A new method to get average expression mat

average.expression <- function(object, ident = "active.ident", assay = "default", slot = "counts", verbose = TRUE) {
  require(Matrix)

  ifelse(assay == "default",
    data <- GetAssayData(object, slot = slot),
    data <- GetAssayData(object, assay = assay, slot = slot))

  # convert data to a sparse matrix if it isn't already
  if (class(data)[1] != "dgCMatrix") data <- Matrix(data, sparse = TRUE)

  # get the idents class over which to average
  ifelse(ident == "active.ident", idents <- as.vector(Idents(object)), idents <- as.vector(object@meta.data[[ident]]))

  # loop through all idents, averaging them in data
  ident.names <- unique(idents)

  if (verbose > 0) pb <- txtProgressBar(char = "=", style = 3, max = length(ident.names), width = 50)
  m <- list()
  for (i in 1:length(ident.names)) {
    ifelse(object@meta.data[[ident]] %in% ident.names[i] %>% sum() > 1, 
           m[[i]] <- Matrix::rowMeans(data[, which(idents == ident.names[i])]),
           m[[i]] <- mean(data[, which(idents == ident.names[i])]))
    
    if (verbose > 0) setTxtProgressBar(pb = pb, value = i)
  }
  result <- do.call(cbind, m)
  colnames(result) <- ident.names
  return(result)
}

t1 <- Sys.time()
mat <- AverageExpression(snRNA.hpfc.sct)[[1]]
t2 <- Sys.time()
t2 - t1


t1 <- Sys.time()
mat2 <- average.expression(snRNA.hpfc.sct,assay = "RNA", ident = "major_clust")
starTracer::searchMarker(mat2,thresh.1 = 0.5,thresh.2 = 0.3,method = "del_MI",num = 2,gene.use = "HVG")
t2 <- Sys.time()
t2 - t1

#10s to 3.5 seconds

t1 <- Sys.time()
starTracer::searchMarker(heart_big,thresh.1 = 0.5,thresh.2 = 0.3,method = "del_MI",num = 3,gene.use = "HVG")
t2 <- Sys.time()
t2 - t1

t1 <- Sys.time()
mat2 <- average.expression(heart_big,assay = "RNA", ident = "cell_type__ontology_label")
starTracer::searchMarker(mat2,thresh.1 = 0.5,thresh.2 = 0.3,method = "del_MI",num = 3,gene.use = "HVG")
t2 <- Sys.time()
t2 - t1

now we are faster again!