Switch to side-by-side view

--- a
+++ b/paper/Benchmarking with real and synthetic datasets/run_other_methods.R
@@ -0,0 +1,86 @@
+library(dyno)
+library(tidyverse)
+library(Matrix)
+library(hdf5r)
+library(Seurat)
+
+method_names <- c('paga','slingshot','monocle')
+container_names <- c('jaydu/ti_paga_tree:v2.0.0','dynverse/ti_slingshot:v2.0.1','jaydu/ti_monocle3:v3.0.0')
+
+path <- "../../data/*.h5"#"/Users/dujinhong/Downloads/new dataset/data_h5/*.h5"#
+save_path <- "/Users/dujinhong/Downloads/result/"#"result/other methods/"
+
+for(load_file in Sys.glob(path)){  
+    filename <- gsub(".h5", "", basename(load_file))
+    if(grepl('mouse_brain',filename)){
+        next
+    }
+
+    dir.create(file.path(save_path, filename))
+    
+    file.h5 <- H5File$new(load_file, mode="r")
+    counts <- t(as.matrix(file.h5[['count']][,]))
+    grouping <- data.frame(cell_id = paste0('c-',as.character(c(1:length(file.h5[['grouping']][])))), 
+                           group_id = as.character(as.vector(file.h5[['grouping']][])))    
+    milestone_net <- data.frame(file.h5[['milestone_network']][])
+    root_milestone_id <- file.h5[['root_milestone_id']][]    
+    if('w' %in% colnames(milestone_net)){
+        cell_ids <- which(milestone_net$from==root_milestone_id)
+        start_id <- paste0('c-',as.character(cell_ids[which.max(milestone_net[cell_ids,'w'])]))
+    }else{
+        start_id <- paste0('c-',sample(which(grouping$group_id==root_milestone_id), 1))
+    }    
+    file.h5$close_all()
+    rownames(counts) <- paste0('c-',as.character(c(1:dim(counts)[1])))
+    colnames(counts) <- paste0('g-',as.character(c(1:dim(counts)[2])))
+    
+    data <- CreateSeuratObject(counts = t(counts), project = "miller")
+    data <- NormalizeData(data, normalization.method = "LogNormalize", scale.factor = 10000)
+    data <- FindVariableFeatures(data, selection.method = "vst", nfeatures = 2000)
+    data <- ScaleData(data, features = rownames(data))
+  
+    counts <- Matrix(counts[,VariableFeatures(data)], sparse = TRUE)
+    expression <- Matrix(t(data[["RNA"]]@scale.data)[,VariableFeatures(data)], sparse = TRUE)
+    
+    
+    
+    dataset <- wrap_expression(
+      expression = expression,
+      counts = counts
+    )
+
+    for(i in c(1:3)){
+        method_name <- method_names[i]                
+
+        method <- create_ti_method_container(container_names[i])        
+        
+        if(i<3){
+            K <- length(unique(grouping$group_id))
+            if(i==1){
+                K <- K - 1
+            }
+            else
+                ndim <- 20
+            pca <- irlba::prcomp_irlba(expression, n = ndim)
+            rd <- pca$x
+            cluster_labels <- cluster::pam(rd, K)$clustering
+            grouping <- data.frame(cell_id = grouping$cell_id, 
+                                   group_id = cluster_labels)   
+        }
+         
+        if(i==2){
+            dataset <- dataset %>% add_prior_information(start_id = c(start_id), groups_id = grouping)
+        }else{
+            dataset <- dataset %>% add_prior_information(start_id = start_id, groups_id = grouping)
+        }
+        model <- infer_trajectory(dataset, method(), give_priors = c("start_id"), verbose = TRUE)
+        
+        
+        ix <- wrapr::match_order(model$progressions$cell_id, grouping$cell_id)
+        model$progressions <- model$progressions[ix,]
+        rownames(model$progressions) <- NULL
+        write.csv(x=model$progressions, file=file.path(save_path, filename, sprintf("%s_progressions.csv", method_name)))
+        
+        write.csv(x=model$milestone_network, file=file.path(save_path, filename, sprintf("%s_milestone_network.csv", method_name)))
+    }
+}
\ No newline at end of file