Monocle3

### Load data from Seurat format
p_load(monocle3, Seurat, tidyverse, patchwork)
sce <- readRDS("sce.rds")
data <- GetAssayData(sce, assay="RNA", layer = "counts")
cell_metadata <- sce@meta.data
gene_annotation <- data.frame(gene_short_name = rownames(data))
rownames(gene_annotation) <- rownames(data)
cds <- new_cell_data_set(data, cell_metadata = cell_metadata, 
  gene_metadata = gene_annotation)
rm(cell_metadata, data, gene_annotation)

### Normalization and PCA reduction
cds <- preprocess_cds(cds, num_dim = 30)
cds <- align_cds(cds, alignment_group = "orig.ident") ###Optional for merged samples
plot_pc_variance_explained(cds)

### UMAP reduction and clustering
cds <- reduce_dimension(cds, preprocess_method = "PCA", reduction_method = "UMAP")
cds <- cluster_cells(cds)
plot_cells(cds, reduction_method = "UMAP", show_trajectory_graph = F,
  graph_label_size = 5, group_label_size = 3)
# "colnames(colData(cds))" to show features for "color_cells_by ="

### Identification
# Use SingleR for identification
p_load(tidyverse, celldex, SingleR, reshape2)
load("D:/Saved_index_for_R/ref_Human_all.RData")
count_for_SingleR <- exprs(cds)
cds_idents <- SingleR(test = count_for_SingleR, ref = ref_Human_all, 
  labels = ref_Human_all$label.main)
table(cds_idents$labels, cds@clusters$UMAP$clusters)

# Use marker genes for identification
marker_test_res <- top_markers(cds, 
  group_cells_by = "cluster", reference_cells = 1000, cores = 8)
top_specific_markers <- marker_test_res %>%
    filter((fraction_expressing >= 0.1)& (specificity >= 0.3)) %>%
    group_by(cell_group) %>%
    top_n(3, pseudo_R2)
top_specific_marker_ids <- unique(top_specific_markers %>% pull(gene_id))
plot_genes_by_group(cds, top_specific_marker_ids,
  group_cells_by = "cluster", ordering_type = "maximal_on_diag",
  max.size = 3) + coord_flip()

#Use custom panel for identification
genes_to_check <- c("...")
plot_genes_by_group(cds, genes_to_check, 
  group_cells_by = "cluster", ordering_type = "maximal_on_diag", 
  max.size = 3) + coord_flip()

### Nomination
colData(cds)$assigned_cell_type = as.character(clusters(cds))
colData(cds)$assigned_cell_type = recode(colData(cds)$assigned_cell_type,
                                         "1"="...",
                                        )
plot_cells(cds, reduction_method = "UMAP", show_trajectory_graph = F, color_cells_by = "assigned_cell_type",
  graph_label_size = 5, group_label_size = 3, cell_size = 0.1) + facet_wrap(~orig.ident)

### Featureplot
plot_cells(cds, genes = c("..."), label_cell_groups = F, show_trajectory_graph = F)

### Trajactory identification
cds <- learn_graph(cds, use_partition = F)
plot_cells(cds, label_groups_by_cluster = F, label_leaves = F, label_branch_points = F)
cds <- order_cells(cds)
plot_cells(cds, color_cells_by = "pseudotime",
           label_cell_groups=FALSE, label_leaves=TRUE, label_branch_points=TRUE,
           graph_label_size=1.5, group_label_size=4,cell_size=1.5)