Controls

We should add more controls such as: 1. Monocle 2. FindAllMarkers but with more methods for finding markers # Monocle3

calculation result

t1 <- Sys.time()
res.cds <- top_markers(cds, group_cells_by="major_clust",cores = 4)
t2 <- Sys.time()
t2 - t1 # around 40 seconds with 8 cores, 47.48599 with 4 cores

res.cds$cell_group <- factor(res.cds$cell_group,levels = levels(snRNA.hpfc.sct$major_clust))
res.cds <- arrange(res.cds,cell_group)
top_specific_markers <- res.cds %>%
                            filter(fraction_expressing >= 0.10) %>%
                            group_by(cell_group) %>%
                            top_n(3, wt = pseudo_R2)

top_specific_marker_ids <- unique(top_specific_markers %>% pull(gene_id))

p <- use_dotplot(snRNA.hpfc.sct,features = rev(top_specific_marker_ids), colors.gradient = c("#003ad0","#ffffff","#d30059"))
p

ggsave(filename = "./out/rev_controls/dotplot_monocle.pdf",p,width = 5,height = 7)

Calculation time

t.monocle <- c()
for (i in c(1,seq(2,16,2))){
  t1 <- Sys.time()
  top_markers(cds, group_cells_by="major_clust",cores = i)
  t2 <- Sys.time()
  t.1 <- t2 - t1
  
  t1 <- Sys.time()
  top_markers(cds, group_cells_by="major_clust",cores = i)
  t2 <- Sys.time()
  t.2 <- t2 - t1
  
  t1 <- Sys.time()
  top_markers(cds, group_cells_by="major_clust",cores = i)
  t2 <- Sys.time()
  t.3 <- t2 - t1
  
  t.i <- c(t.1,t.2,t.3)
  t.monocle <- append(t.monocle,t.i)
}
save(t.monocle,file = "./out/rev_controls/monocle_runtime.Rdata")

plot

t.monocle.frame <- matrix(data = t.monocle, ncol = 3, nrow = 9,byrow = T,
                    dimnames = list(c("core.1",paste0("core.",seq(2,16,2))),
                    c("run1","run2","run3"))) %>% as.data.frame()

t.monocle.frame <- t(t.monocle.frame) %>% reshape2::melt()
#save(t.monocle.frame,file = "./out/rev_controls/monocle_frame.Rdata")

t.monocle.frame <- Rmisc::summarySE(t.monocle.frame, measurevar = "value", groupvar = c("Var2"))

(p <- ggplot(data = t.monocle.frame) + 
  geom_errorbar(aes(x = Var2,y = value, ymin = value - se, ymax = value + se, color = Var2), width = 0.02) + 
  geom_line(aes(x = Var2,y = value, group = 1),color = "black",lwd = 0.5) +
  geom_point(aes(x = Var2,y = value, color = Var2)) + 
  scale_color_manual(values = viridis::magma(10)) + 
  theme_bw() + ylab("time") + ylim(0,150) + geom_hline(yintercept = 10.79,lty = 2))
ggsave("./out/rev_controls/line_time_monocle_brain.pdf",p,width = 5,height = 5)

#FindAllMarkers There are a lot of ways to perform FindAllMarkers, including: 1. test.use = “wilcox” 2. test.use = “bimod” 3. test.use = “roc” 4. test.use = “t” 5. test.use = “negbinom” 6. test.use = “poisson” 7. test.use = “LR” 8. test.use = “MAST” 9. test.use = “DESeq2” Here we use the same parameter as we did in “parameter_control_FindAllMarkers” FC will be set as 1.442839 pct.min will be set as 0.5

methods <- c("wilcox","bimod","roc","t","negbinom","poisson","LR","MAST","DESeq2")
t.methods <- c()
for(i in methods){
  t1 <- Sys.time()
  assign(paste0("res.seurat.method.",i),FindAllMarkers(snRNA.hpfc.sct,
                        logfc.threshold = mean.fc,
                        min.pct = 0.5,
                        only.pos = T,
                        test.use = i))
  t2 <- Sys.time()
  t.i1 <- t2 - t1
  
  t1 <- Sys.time()
  assign(paste0("res.seurat.method.",i),FindAllMarkers(snRNA.hpfc.sct,
                        logfc.threshold = mean.fc,
                        min.pct = 0.5,
                        only.pos = T,
                        test.use = i))
  t2 <- Sys.time()
  t.i2 <- t2 - t1
  
  t1 <- Sys.time()
  assign(paste0("res.seurat.method.",i),FindAllMarkers(snRNA.hpfc.sct,
                        logfc.threshold = mean.fc,
                        min.pct = 0.5,
                        only.pos = T,
                        test.use = i))
  t2 <- Sys.time()
  t.i3 <- t2 - t1
  
  t.i <- c(t.i1,t.i2,t.i3)
  t.methods <- append(t.methods,t.i)
}

done draw plot

t.methods.frame <- matrix(data = t.methods,nrow = 9,ncol = 3, byrow = T,
                          dimnames = list(methods,c("run1","run2","run3"))) %>% as.data.frame()
save(t.methods.frame,file = "./out/rev_controls/time.methods.frame.Rdata")

t.methods.frame <- t(t.methods.frame) %>% reshape2::melt()
t.methods.frame <- Rmisc::summarySE(t.methods.frame, measurevar = "value", groupvar = c("Var2"))

(p <- ggplot(data = t.methods.frame) + 
  geom_bar(aes(x = Var2,y = value, group = 1,fill = Var2),stat = "identity",color = "black",lwd = 0.5) + 
  scale_fill_manual(values = viridis::magma(10)) + 
  geom_errorbar(aes(x = Var2,y = value, ymin = value - se, ymax = value + se),color = "black", width = 0.02) + 
  geom_point(aes(x = Var2,y = value),color = "black") + 
  theme_bw() + ylab("time") + ylim(0,1100) + geom_hline(yintercept = 10.79,lty = 2))
ggsave(p,filename = "out/rev_controls/methods.bar.pdf",width = 5,height = 4)