library(monocle3)
library(Seurat)
library(SeuratData)
library(SeuratWrappers)
library(ggplot2)
library(dplyr)


wdir <- "Birocchi_SciTraMed2022/results"
macro_obj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
out_dir <- paste(wdir, "MacrophagesFiltered", "03-fastMNN_2-pseudotime", sep = "/")

# Load the data
expression_matrix <- macro_obj@assays$RNA@counts
cell_metadata <- macro_obj@meta.data
cell_metadata$orig.ident <- factor(cell_metadata$orig.ident, levels = unique(cell_metadata$orig.ident))
gene_annotation <- data.frame("gene_short_name" = rownames(macro_obj))
rownames(gene_annotation) <- gene_annotation$gene_short_name

# Make the CDS object
cds <- new_cell_data_set(expression_matrix,
                         cell_metadata = cell_metadata,
                         gene_metadata = gene_annotation)
cds <- preprocess_cds(cds, num_dim = 100)

cds <- align_cds(cds, num_dim = 50, alignment_group = "RNA_Group")
cds <- reduce_dimension(cds)

png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_Monocle_RNA_Group_umap_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds, color_cells_by = "RNA_Group", label_cell_groups = FALSE, group_label_size = 8, cell_size = 0.5)
dev.off()
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_Monocle_RNA_snn_res.0.6_umap_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds, color_cells_by = "RNA_snn_res.0.6", label_cell_groups = T, group_label_size = 8, cell_size = 0.5)
dev.off()

cds <- cluster_cells(cds = cds,resolution = 0.001)
plot_cells(cds, cell_size = 1, color_cells_by = "partition")
plot_cells(cds, genes=c("Pglyrp1", "Ace"))


marker_test_res <- top_markers(cds, group_cells_by="partition", reference_cells=1000, cores=8)
head(marker_test_res)
marker_test_res$cell_group
top_specific_markers <- marker_test_res %>%
    filter(fraction_expressing >= 0.10) %>%
    group_by(cell_group) %>% filter(cell_group == "6") %>%
    top_n(20, pseudo_R2)


marker_test_res %>%
    filter(fraction_expressing >= 0.10) %>%
    group_by(cell_group) %>% filter(cell_group == "3") %>%
    top_n(20, pseudo_R2) %>% select(gene_short_name) %>% list()


tt3 <- as.data.frame(top_specific_markers)
str(tt3)
tt3$cell_group
subset(tt3, tt3$cell_group %in% c("3", "4", "5", "6"))

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

cds <- learn_graph(cds)

png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_Monocle_RNA_snn_res.0.6_umap_trajectory_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds,
           color_cells_by = "RNA_snn_res.0.6",
           label_groups_by_cluster=FALSE,
           label_leaves=F,
           label_branch_points=F,
           graph_label_size=1.5,
           group_label_size = 6, cell_size = 0.5)
dev.off()

cds <- order_cells(cds)
png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_Monocle_pseudotime_umap_trajectory_plot3.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds,
           color_cells_by = "pseudotime",
           label_groups_by_cluster=FALSE,
           label_leaves=T,
           label_branch_points=T,
           graph_label_size=2.5,
           group_label_size = 8, cell_size = 0.5)
dev.off()


# a helper function to identify the root principal points:
get_earliest_principal_node <- function(cds, time_bin="7"){
    cell_ids <- which(colData(cds)[, "RNA_snn_res.0.6"] == time_bin)
    closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
    closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
    root_pr_nodes <- igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names(which.max(table(closest_vertex[cell_ids,]))))]
    root_pr_nodes
}
cds <- order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))

png(filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_Monocle_pseudotime_umap_trajectory_plot2.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds,
           color_cells_by = "pseudotime",
           label_groups_by_cluster=FALSE,
           label_leaves=T,
           label_branch_points=T,
           graph_label_size=2.5,
           group_label_size = 8, cell_size = 0.5)
dev.off()

cds_pr_test_res <- graph_test(cds, neighbor_graph = "principal_graph", cores = 4)
head(cds_pr_test_res, n = 10)
nrow(cds_pr_test_res)
pr_deg_ids <- subset(cds_pr_test_res, q_value < 0.05)
head(pr_deg_ids, n = 10)
nrow(pr_deg_ids)

plot_cells(cds, genes = rownames(head(pr_deg_ids[order(pr_deg_ids$q_value, decreasing = F),], n = 9)),
           show_trajectory_graph = T,
           label_cell_groups = FALSE,
           label_leaves = FALSE)

# plot_cells(cds, genes = rownames(pr_deg_ids)[1:6],
#            show_trajectory_graph = T,
#            label_cell_groups = FALSE,
#            label_leaves = FALSE)

gene_module_df <- find_gene_modules(cds[rownames(pr_deg_ids),], resolution=c(10^-10, 10^seq(-6,-1)))
cell_group_df <- tibble::tibble(cell=row.names(colData(cds)), 
                                cell_group=colData(cds)$RNA_snn_res.0.6)
agg_mat <- aggregate_gene_expression(cds, gene_module_df, cell_group_df)
row.names(agg_mat) <- stringr::str_c("Module ", row.names(agg_mat))
pheatmap::pheatmap(agg_mat,angle_col = 0, scale="column", clustering_method="ward.D2")

# mod <- 5 # Cluster 6
# mod <- 33 # Clsuter 10
# pr_deg_ids_mod <- pr_deg_ids[gene_module_df[gene_module_df$module == mod,][["id"]],]
# plot_cells(cds, genes = rownames(head(pr_deg_ids_mod[order(pr_deg_ids_mod$q_value, decreasing = F),], n = 9)),
#            show_trajectory_graph = T,
#            label_cell_groups = FALSE,
#            label_leaves = FALSE)

plot_cells(cds,
           genes = gene_module_df %>% filter(module %in% c(7, 25, 54, 14, 32, 12, 19, 38, 24, 58, 15, 5, 16, 43, 18, 29)),
           label_cell_groups = FALSE,
           show_trajectory_graph = FALSE)

# Analyzing branches in single-cell trajectories 
cds_subset <- choose_cells(cds)
subset_pr_test_res <- graph_test(cds_subset, neighbor_graph="principal_graph", cores=4)
pr_deg_ids <- row.names(subset(subset_pr_test_res, q_value < 0.05))
gene_module_df <- find_gene_modules(cds_subset[pr_deg_ids,], resolution=0.001)
agg_mat <- aggregate_gene_expression(cds_subset, gene_module_df)
module_dendro <- hclust(dist(agg_mat))
gene_module_df$module <- factor(gene_module_df$module, levels = row.names(agg_mat)[module_dendro$order])

plot_cells(cds_subset,
           genes=gene_module_df,
           label_cell_groups=FALSE,
           show_trajectory_graph=FALSE)

clu6_genes <- c("Mrc1", "Ccl8")
clu6_lineage_cds <- cds[rowData(cds)$gene_short_name %in% clu6_genes,
                        colData(cds)$RNA_snn_res.0.6 %in% c("6", "7", "10")]
plot_genes_in_pseudotime(clu6_lineage_cds, color_cells_by="RNA_snn_res.0.6", min_expr=0.5)

#####################################################
#####################################################
Idents(macro_obj) <- "RNA_Group"
# Ctrl
macro_obj_ctrl <- subset(macro_obj, idents = "AB_Ctrl")
expression_matrix_ctrl <- macro_obj_ctrl@assays$RNA@counts
cell_metadata_ctrl <- macro_obj_ctrl@meta.data
cell_metadata_ctrl$orig.ident <- factor(cell_metadata_ctrl$orig.ident, levels = unique(cell_metadata_ctrl$orig.ident))
gene_annotation_ctrl <- data.frame("gene_short_name" = rownames(macro_obj_ctrl))
rownames(gene_annotation_ctrl) <- gene_annotation_ctrl$gene_short_name
cds_ctrl <- new_cell_data_set(expression_data = expression_matrix_ctrl,
                              cell_metadata = cell_metadata_ctrl,
                              gene_metadata = gene_annotation_ctrl)
cds_ctrl <- preprocess_cds(cds_ctrl, num_dim = 50)
cds_ctrl <- reduce_dimension(cds_ctrl)
### Construct and assign the made up partition
recreate.partition_ctrl <- c(rep(1, length(cds_ctrl@colData@rownames)))
names(recreate.partition_ctrl) <- cds_ctrl@colData@rownames
recreate.partition_ctrl <- as.factor(recreate.partition_ctrl)
cds_ctrl@clusters@listData[["UMAP"]][["partitions"]] <- recreate.partition_ctrl
### Assign the cluster info
list_cluster_ctrl <- macro_obj_ctrl@meta.data$RNA_snn_res.0.6
names(list_cluster_ctrl) <- macro_obj_ctrl@assays[["RNA"]]@data@Dimnames[[2]]
cds_ctrl@clusters@listData[["UMAP"]][["clusters"]] <- list_cluster_ctrl
### Could be a space-holder, but essentially fills out louvain parameters
cds_ctrl@clusters@listData[["UMAP"]][["louvain_res"]] <- "NA"
### Assign UMAP coordinate
reducedDims(cds_ctrl)[["UMAP"]] <- macro_obj_ctrl@reductions[["umap"]]@cell.embeddings
### Assign feature loading for downstream module analysis
##########cds_ctrl@preprocess_aux$gene_loadings <- macro_obj_ctrl@reductions[["mnn"]]@cell.embeddings
### Reset colnames and rownames (consequence of UMAP replacement)
rownames(cds_ctrl@principal_graph_aux[['UMAP']]$dp_mst) <- NULL
#########colnames(reducedDims(cds_ctrl)[["UMAP"]]) <- NULL
plot_cells(cds = cds_ctrl, color_cells_by = "RNA_snn_res.0.6", group_label_size = 6, cell_size = 1)
# Learn Graph
cds_ctrl <- learn_graph(cds_ctrl, use_partition = T)
png(filename = paste(out_dir, "MacrophagesFiltered_Ctrl_Monocle_RNA_snn_res.0.6_umap_trajectory_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_ctrl,
           color_cells_by = "RNA_snn_res.0.6",
           label_groups_by_cluster=FALSE,
           label_leaves=F,
           label_branch_points=F,
           graph_label_size=1.5,
           group_label_size = 6, cell_size = 1)
dev.off()
cds_ctrl <- order_cells(cds_ctrl)
png(filename = paste(out_dir, "MacrophagesFiltered_Ctrl_Monocle_pseudotime_umap_trajectory_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_ctrl,
           color_cells_by = "pseudotime",
           label_groups_by_cluster=FALSE,
           label_leaves=T,
           label_branch_points=T,
           graph_label_size=2.5,
           group_label_size = 8, cell_size = 1)
dev.off()
####
cds_ctrl_ps_res = graph_test(cds_ctrl, neighbor_graph = "principal_graph", cores = 4)
cds_ctrl_ps_res_sig = row.names(subset(cds_ctrl_ps_res, q_value == 0))
length(cds_ctrl_ps_res_sig)
png(filename = paste(out_dir, "MacrophagesFiltered_Ctrl_Monocle_pseudotime_umap_topGenes.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_ctrl, genes = cds_ctrl_ps_res_sig[0:25],
           show_trajectory_graph = T,
           label_cell_groups = FALSE,
           label_leaves = FALSE)
dev.off()
gene_module_df_ctrl <- find_gene_modules(cds_ctrl[cds_ctrl_ps_res_sig,], resolution=c(10^-100,10^seq(-6,-1)))
cell_group_df_ctrl <- tibble::tibble(cell=row.names(colData(cds_ctrl)), 
                                    cell_group=colData(cds_ctrl)$RNA_snn_res.0.6)
agg_mat_ctrl <- aggregate_gene_expression(cds_ctrl, gene_module_df_ctrl, cell_group_df_ctrl)
row.names(agg_mat_ctrl) <- stringr::str_c("Module ", row.names(agg_mat_ctrl))
pheatmap::pheatmap(agg_mat_ctrl, scale = "column", clustering_method = "ward.D2",
                   filename = paste(out_dir, "MacrophagesFiltered_Ctrl_Monocle_ModuleHM.png", sep = "/"))
png(filename = paste(out_dir, "MacrophagesFiltered_Ctrl_Monocle_pseudotime_umap_Modules.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds = cds_ctrl, cell_size = 1,
           genes = gene_module_df_ctrl,
           label_cell_groups = FALSE,
           show_trajectory_graph = FALSE)
dev.off()
######
head(gene_module_df_ctrl)
nrow(gene_module_df_ctrl)
for (mod in levels(gene_module_df_ctrl$module)) {
    tt <- as.data.frame(subset(gene_module_df_ctrl, module == mod))
    print(nrow(tt))
    lineage_cds_ctrl <- cds_ctrl[rowData(cds_ctrl)$gene_short_name %in% tt$id,]
    png(filename = paste(out_dir, paste0("MacrophagesFiltered_Ctrl_Monocle_pseudotime_umap_Module", mod, "_genes.png"), sep = "/"),
        width = 10, height = 12, units = "in", res = 300)
    print(plot_genes_in_pseudotime(cds_subset = lineage_cds_ctrl, ncol = 3, color_cells_by="RNA_snn_res.0.6", min_expr = 1))
    dev.off()
}

# High
macro_obj_high <- subset(macro_obj, idents = "CD_High")
expression_matrix_high <- macro_obj_high@assays$RNA@counts
cell_metadata_high <- macro_obj_high@meta.data
cell_metadata_high$orig.ident <- factor(cell_metadata_high$orig.ident, levels = unique(cell_metadata_high$orig.ident))
gene_annotation_high <- data.frame("gene_short_name" = rownames(macro_obj_high))
rownames(gene_annotation_high) <- gene_annotation_high$gene_short_name
cds_high <- new_cell_data_set(expression_data = expression_matrix_high,
                              cell_metadata = cell_metadata_high,
                              gene_metadata = gene_annotation_high)
cds_high <- preprocess_cds(cds_high, num_dim = 50)
cds_high <- reduce_dimension(cds_high)
### Construct and assign the made up partition
recreate.partition_high <- c(rep(1, length(cds_high@colData@rownames)))
names(recreate.partition_high) <- cds_high@colData@rownames
recreate.partition_high <- as.factor(recreate.partition_high)
cds_high@clusters@listData[["UMAP"]][["partitions"]] <- recreate.partition_high
### Assign the cluster info
list_cluster_high <- macro_obj_high@meta.data$RNA_snn_res.0.6
names(list_cluster_high) <- macro_obj_high@assays[["RNA"]]@data@Dimnames[[2]]
cds_high@clusters@listData[["UMAP"]][["clusters"]] <- list_cluster_high
### Could be a space-holder, but essentially fills out louvain parameters
cds_high@clusters@listData[["UMAP"]][["louvain_res"]] <- "NA"
### Assign UMAP coordinate
reducedDims(cds_high)[["UMAP"]] <- macro_obj_high@reductions[["umap"]]@cell.embeddings
### Assign feature loading for downstream module analysis
#########cds_high@preprocess_aux$gene_loadings <- macro_obj_high@reductions[["mnn"]]@cell.embeddings
### Reset colnames and rownames (consequence of UMAP replacement)
rownames(cds_high@principal_graph_aux[['UMAP']]$dp_mst) <- NULL
#########colnames(cds_high@reducedDims$UMAP) <- NULL
plot_cells(cds = cds_high, color_cells_by = "RNA_snn_res.0.6", group_label_size = 6, cell_size = 1)
# Learn Graph
cds_high <- learn_graph(cds_high, use_partition = T)
png(filename = paste(out_dir, "MacrophagesFiltered_High_Monocle_RNA_snn_res.0.6_umap_trajectory_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_high,
           color_cells_by = "RNA_snn_res.0.6",
           label_groups_by_cluster=FALSE,
           label_leaves=F,
           label_branch_points=F,
           graph_label_size=1.5,
           group_label_size = 6, cell_size = 1)
dev.off()
cds_high <- order_cells(cds_high)
png(filename = paste(out_dir, "MacrophagesFiltered_High_Monocle_pseudotime_umap_trajectory_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_high,
           color_cells_by = "pseudotime",
           label_groups_by_cluster=FALSE,
           label_leaves=T,
           label_branch_points=T,
           graph_label_size=2.5,
           group_label_size = 8, cell_size = 1)
dev.off()
####
cds_high_ps_res = graph_test(cds_high, neighbor_graph = "principal_graph", cores = 4)
cds_high_ps_res_sig = row.names(subset(cds_high_ps_res, q_value == 0))
length(cds_high_ps_res_sig)
png(filename = paste(out_dir, "MacrophagesFiltered_High_Monocle_pseudotime_umap_topGenes.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_high, genes = cds_high_ps_res_sig[0:25],
           show_trajectory_graph = T,
           label_cell_groups = FALSE,
           label_leaves = FALSE)
dev.off()
gene_module_df_high <- find_gene_modules(cds_high[cds_high_ps_res_sig,], resolution=c(10^-100,10^seq(-6,-1)))
cell_group_df_high <- tibble::tibble(cell=row.names(colData(cds_high)), 
                                    cell_group=colData(cds_high)$RNA_snn_res.0.6)
agg_mat_high <- aggregate_gene_expression(cds_high, gene_module_df_high, cell_group_df_high)
row.names(agg_mat_high) <- stringr::str_c("Module ", row.names(agg_mat_high))
pheatmap::pheatmap(agg_mat_high, scale = "column", clustering_method = "ward.D2",
                   filename = paste(out_dir, "MacrophagesFiltered_High_Monocle_ModuleHM.png", sep = "/"))
png(filename = paste(out_dir, "MacrophagesFiltered_High_Monocle_pseudotime_umap_Modules.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds = cds_high, cell_size = 1,
           genes = gene_module_df_high,
           label_cell_groups = FALSE,
           show_trajectory_graph = FALSE)
dev.off()
######
head(gene_module_df_high)
nrow(gene_module_df_high)
for (mod in levels(gene_module_df_high$module)) {
    tt <- as.data.frame(subset(gene_module_df_high, module == mod))
    print(nrow(tt))
    lineage_cds_high <- cds_high[rowData(cds_high)$gene_short_name %in% tt$id,]
    png(filename = paste(out_dir, paste0("MacrophagesFiltered_High_Monocle_pseudotime_umap_Module", mod, "_genes.png"), sep = "/"),
        width = 10, height = 12, units = "in", res = 300)
    print(plot_genes_in_pseudotime(cds_subset = lineage_cds_high, ncol = 3, color_cells_by="RNA_snn_res.0.6", min_expr = 1))
    dev.off()
}

# Low
macro_obj_low <- subset(macro_obj, idents = "EF_Low")
expression_matrix_low <- macro_obj_low@assays$RNA@counts
cell_metadata_low <- macro_obj_low@meta.data
cell_metadata_low$orig.ident <- factor(cell_metadata_low$orig.ident, levels = unique(cell_metadata_low$orig.ident))
gene_annotation_low <- data.frame("gene_short_name" = rownames(macro_obj_low))
rownames(gene_annotation_low) <- gene_annotation_low$gene_short_name
cds_low <- new_cell_data_set(expression_data = expression_matrix_low,
                             cell_metadata = cell_metadata_low,
                             gene_metadata = gene_annotation_low)
cds_low <- preprocess_cds(cds_low, num_dim = 50)
cds_low <- reduce_dimension(cds_low)
### Construct and assign the made up partition
recreate.partition_low <- c(rep(1, length(cds_low@colData@rownames)))
names(recreate.partition_low) <- cds_low@colData@rownames
recreate.partition_low <- as.factor(recreate.partition_low)
cds_low@clusters@listData[["UMAP"]][["partitions"]] <- recreate.partition_low
### Assign the cluster info
list_cluster_low <- macro_obj_low@meta.data$RNA_snn_res.0.6
names(list_cluster_low) <- macro_obj_low@assays[["RNA"]]@data@Dimnames[[2]]
cds_low@clusters@listData[["UMAP"]][["clusters"]] <- list_cluster_low
### Could be a space-holder, but essentially fills out louvain parameters
cds_low@clusters@listData[["UMAP"]][["louvain_res"]] <- "NA"
### Assign UMAP coordinate
reducedDims(cds_low)[["UMAP"]] <- macro_obj_low@reductions[["umap"]]@cell.embeddings
### Assign feature loading for downstream module analysis
############cds_low@preprocess_aux$gene_loadings <- macro_obj_low@reductions[["mnn"]]@cell.embeddings
### Reset colnames and rownames (consequence of UMAP replacement)
rownames(cds_low@principal_graph_aux[['UMAP']]$dp_mst) <- NULL
#colnames(cds_low@reducedDims$UMAP) <- NULL
#plot_cells(cds = cds_low, color_cells_by = "RNA_snn_res.0.6", group_label_size = 6, cell_size = 1)
# Learn Graph
cds_low <- learn_graph(cds_low, use_partition = T)
png(filename = paste(out_dir, "MacrophagesFiltered_Low_Monocle_RNA_snn_res.0.6_umap_trajectory_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_low,
           color_cells_by = "RNA_snn_res.0.6",
           label_groups_by_cluster=FALSE,
           label_leaves=F,
           label_branch_points=F,
           graph_label_size=1.5,
           group_label_size = 6, cell_size = 1)
dev.off()
cds_low <- order_cells(cds_low)
png(filename = paste(out_dir, "MacrophagesFiltered_Low_Monocle_pseudotime_umap_trajectory_plot.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_low,
           color_cells_by = "pseudotime",
           label_groups_by_cluster=FALSE,
           label_leaves=T,
           label_branch_points=T,
           graph_label_size=2.5,
           group_label_size = 8, cell_size = 1)
dev.off()
####
cds_low_ps_res = graph_test(cds_low, neighbor_graph = "principal_graph", cores = 4)
cds_low_ps_res_sig = row.names(subset(cds_low_ps_res, q_value == 0))
length(cds_low_ps_res_sig)
png(filename = paste(out_dir, "MacrophagesFiltered_Low_Monocle_pseudotime_umap_topGenes.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds_low, genes = cds_low_ps_res_sig[0:25],
           show_trajectory_graph = T,
           label_cell_groups = FALSE,
           label_leaves = FALSE)
dev.off()
gene_module_df_low <- find_gene_modules(cds_low[cds_low_ps_res_sig,], resolution=c(10^-100,10^seq(-6,-1)))
cell_group_df_low <- tibble::tibble(cell=row.names(colData(cds_low)), 
                                cell_group=colData(cds_low)$RNA_snn_res.0.6)
agg_mat_low <- aggregate_gene_expression(cds_low, gene_module_df_low, cell_group_df_low)
row.names(agg_mat_low) <- stringr::str_c("Module ", row.names(agg_mat_low))
pheatmap::pheatmap(agg_mat_low, scale = "column", clustering_method = "ward.D2",
                   filename = paste(out_dir, "MacrophagesFiltered_Low_Monocle_ModuleHM.png", sep = "/"))
png(filename = paste(out_dir, "MacrophagesFiltered_Low_Monocle_pseudotime_umap_Modules.png", sep = "/"),
    width = 10, height = 8, units = "in", res = 300)
plot_cells(cds = cds_low, cell_size = 1,
           genes = gene_module_df_low,
           label_cell_groups = FALSE,
           show_trajectory_graph = FALSE)
dev.off()
######
head(gene_module_df_low)
nrow(gene_module_df_low)
for (mod in levels(gene_module_df_low$module)) {
    tt <- as.data.frame(subset(gene_module_df_low, module == mod))
    print(nrow(tt))
    lineage_cds_low <- cds_low[rowData(cds_low)$gene_short_name %in% tt$id,]
    png(filename = paste(out_dir, paste0("MacrophagesFiltered_Low_Monocle_pseudotime_umap_Module", mod, "_genes.png"), sep = "/"),
        width = 10, height = 12, units = "in", res = 300)
    print(plot_genes_in_pseudotime(cds_subset = lineage_cds_low, ncol = 3, color_cells_by="RNA_snn_res.0.6", min_expr = 1))
    dev.off()
}

cds_obj <- list("Full" = cds, "AB_Ctrl" = cds_ctrl, "CD_High" = cds_high, "EF_Low" = cds_low)
saveRDS(object = cds_obj, file = paste(out_dir, "Monocle_cds_objects.rds", sep = "/"))