---

title: "my_monocle"

author: "TianyuChen"

date: "2022-10-08"

output: html_document

---

```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = TRUE)

```

```{r}

library(monocle3)

library(rhdf5)

library(viridis)

```

```{r}

f <- "C:/Users/10270/Desktop/硕一上/Wang/Trajectory/Dev/GitVersion/Monocle/Nature_Mouse_Merge/Batch_correct_fig/monocle_traj_batch_corrected.Rds"

m <- readRDS(f)

```

```{r}

umap <- m@reduce_dim_aux@listData[["UMAP"]]@listData[["model"]]@listData[["umap_model"]][["embedding"]]

```

```{r}

write.table(umap,"monocle3_umap.txt")

```

```{r}

meta <- colData(m)

write.table(meta,"monocle3_meta.txt")

```

```{r}

marker <- c("Hmga2","Sox2", "Pax6", "Hes5","Ube2c", "Id4","Olig1", "Olig2", "Pdgfra","Apoe",

            "Neurog2","Neurod1","Npy","Sst","Nxph1","Htr3a","Prox1","Cxcl14","Meis2","Etv1","Sp8",

            "Btg2","Neurog2","Hes6","Slc1a3","Dbi","Fabp7","Bcl11b")

m$Hmga2

```

```{r}

p1 <- plot_cells(m,

           color_cells_by = "Cluster2",

           label_cell_groups=FALSE,

           label_leaves=TRUE,

           label_branch_points=TRUE,

           graph_label_size=1.5)

p1

```

```{r}

png(file="traj_plot1.png",width=1000, height=700,res = 100)

p1

dev.off()

```

```{r}

ggsave(

  "traj_plot1.png",

  p1,

  width = 15,

  height = 10,

  dpi = 1200

)

```

```{r}

p2 <- plot_cells(m,

           color_cells_by = "Day",

           label_cell_groups=FALSE,

           label_leaves=TRUE,

           label_branch_points=TRUE,

           graph_label_size=1.5)

p2

```

```{r}

ggsave(

  "traj_plot2.png",

  p2,

  width = 15,

  height = 10,

  dpi = 1200

)

```

```{r}

m <- order_cells(m)

```

```{r}

p3 <- plot_cells(m,

           color_cells_by = "pseudotime",

           label_cell_groups=FALSE,

           label_leaves=FALSE,

           label_branch_points=FALSE,

           graph_label_size=1.5)

p3

```

```{r}

ggsave(

  "traj_plot3.png",

  p3,

  width = 15,

  height = 10,

  dpi = 1200

)

```

```{r}

ciliated_genes <- c(

                    "Hmga2")

p4 <- plot_cells(m,

           genes=ciliated_genes,

           label_cell_groups=FALSE,

           show_trajectory_graph=FALSE)

```

```{r}

ggsave(

  "NEC_marker.png",

  p4,

  width = 15,

  height = 10,

  dpi = 1200

)

```

```{r}

p4

```

```{r}

#saveRDS(cds, file = "monocle_traj_batch_corrected.Rds")

```

```{r}

ciliated_genes <- c("Sox2",

                    "Pax6","Hes5","Ube2c","Id4")

p5 <- plot_cells(m,

           genes=ciliated_genes,

           label_cell_groups=FALSE,

           show_trajectory_graph=FALSE)

p5

```

```{r}

ggsave(

  "RGC_marker.png",

  p5,

  width = 15,

  height = 10,

  dpi = 1200

)

```

```{r}

ciliated_genes <- c("Olig1",

                    "Olig2","Pdgfra","Apoe")

p6 <- plot_cells(m,

           genes=ciliated_genes,

           label_cell_groups=FALSE,

           show_trajectory_graph=FALSE)

p6

```

```{r}

ggsave(

  "OPC_marker.png",

  p6,

  width = 15,

  height = 10,

  dpi = 1200

)

```

```{r}

ciliated_genes <- c("Neurog2",

                    "Neurod1")

p7 <- plot_cells(m,

           genes=ciliated_genes,

           label_cell_groups=FALSE,

           show_trajectory_graph=FALSE)

p7

```

```{r}

ggsave(

  "inter_marker.png",

  p7,

  width = 15,

  height = 10,

  dpi = 1200

)

```

```{r}

ciliated_genes <- c("Bcl11b")

p8 <- plot_cells(m,

           genes=ciliated_genes,

           label_cell_groups=FALSE,

           show_trajectory_graph=FALSE)

p8

```

```{r}

ggsave(

  "layer6_marker.png",

  p8,

  width = 15,

  height = 10,

  dpi = 1200

)

```

```{r}

cd <- colData(m)

```

```{r}

ciliated_genes <- c("Npy","Sst","Nxph1","Htr3a","Prox1","Cxcl14","Meis2","Etv1","Sp8")

p9 <- plot_cells(m,

           genes=ciliated_genes,

           label_cell_groups=FALSE,

           show_trajectory_graph=FALSE)

```

```{r}

p9

```

```{r}

ggsave(

  "interneurons_marker.png",

  p9,

  width = 15,

  height = 10,

  dpi = 1200

)

```

```{r}

ciliated_genes <- c("Htr3a","Prox1","Cxcl14")

plot_cells(m,

           genes=ciliated_genes,

           label_cell_groups=FALSE,

           show_trajectory_graph=FALSE)

```

```{r}

ciliated_genes <- c("Meis2","Etv1","Sp8")

plot_cells(m,

           genes=ciliated_genes,

           label_cell_groups=FALSE,

           show_trajectory_graph=FALSE)

```

```{r}

cl <- cd$Clusters

day <- cd$Day

source <- cd$Source

source <- as.character(source)

```

```{r}

cl <- as.character(cl)

cl[source == "2.0"] <- NA

```

```{r}

day <- as.character(day)

day[source == "2.0"] <- NA

```

```{r}

cd$mouse_day <- day

cd$mouse_clusters <- cl

```

```{r}

base_cl <- as.character(cd$Clusters)

base_cl[base_cl == "SCPN1"] <- "SCPN"

cd$base_cl <- base_cl

```

```{r}

colData(m) <- cd

```

```{r}

library(ggplot2)

```

```{r}

cols <- c("Endothelial Cell" = "#F8766D", "Immature Neuron" = "#E38900",

          "Interneurons" = "#C49A00", "IPC" = "#99A800",

          "Layer I" = "#53B400", "Layer II-IV" = "#00BC56",

          "Layer V-VI" = "#00C094", "Layer V-VI (Hippo)" = "#00BFC4",

          "Microglia" = "#00B6EB", "NEC" = "#06A4FF",

          "OPC" = "#A58AFF", "Pericyte" = "#DF70F8",

          "Pia" = "#FB61D7", "RGC" = "#FF66A8")

```

```{r}

p10 <- plot_cells(m,

           color_cells_by = "mouse_clusters",

           label_cell_groups=FALSE,cell_size = 0.7,alpha = 1,

           label_leaves=FALSE,

           label_branch_points=FALSE,

           graph_label_size=1.5)+scale_colour_manual(

             values = cols,

  aesthetics = "colour",

  breaks = waiver(),

  na.value = "#E5E8E8")

p10

```

```{r}

ggsave(

  "mouse_clusters.png",

  p10,

  width = 15,

  height = 10,

  dpi = 1200

)

```

```{r}

p11 <- plot_cells(m,

           color_cells_by = "mouse_day",

           label_cell_groups=FALSE,

           label_leaves=FALSE,

           label_branch_points=FALSE,

           graph_label_size=1.5)

p11

```

```{r}

ggsave(

  "mouse_days.png",

  p11,

  width = 15,

  height = 10,

  dpi = 1200

)

```

```{r}

p12 <- plot_cells(m,

           color_cells_by = "base_cl",

           label_cell_groups=FALSE,

           label_leaves=FALSE,

           label_branch_points=FALSE,

           graph_label_size=1.5)

p12

```

```{r}

ggsave(

  "traj_base_cl.png",

  p12,

  width = 15,

  height = 10,

  dpi = 1200

)

```

# Load model

```{r}

#cds <- readRDS("monocle_traj.Rds")

```

```{r}

plot_cells(cds,

           color_cells_by = "Clusters",

           label_cell_groups=FALSE,

           label_leaves=FALSE,

           label_branch_points=FALSE,

           graph_label_size=1.5)

```

```{r}

ciliated_genes <- c("Fezf2")

p9 <- plot_cells(cds,

           genes=ciliated_genes,

           label_cell_groups=FALSE,

           show_trajectory_graph=FALSE)

p9

```

```{r}

merge_18 <- rep(NA,dim(cds)[2])

day <- as.character(colData(cds)$Day)

merge_18[day=="E18"] <- 0

merge_18[day=="E18_S1"] <- 1

merge_18[day=="E18_S3"] <- 2

```

```{r}

colData(cds)$merge_18 <- merge_18

```

```{r}

p11 <- plot_cells(cds,

           color_cells_by = "merge_18",

           label_cell_groups=FALSE,

           label_leaves=FALSE,

           label_branch_points=FALSE,

           graph_label_size=1.5)

p11

```

```{r}

png(file="merge_18.png",width=1000, height=700)

p11

dev.off()

```

```{r}

merge_P1 <- rep(NA,dim(cds)[2])

day <- as.character(colData(cds)$Day)

merge_P1[day=="P1"] <- 0

merge_P1[day=="P1_S1"] <- 1

colData(cds)$merge_P1 <- merge_P1

```

```{r}

p13 <- plot_cells(cds,

           color_cells_by = "merge_P1",

           label_cell_groups=FALSE,

           label_leaves=FALSE,

           label_branch_points=FALSE,

           graph_label_size=1.5)

p13

```

```{r}

p14 <- plot_cells(cds,

           color_cells_by = "Source",

           label_cell_groups=FALSE,

           label_leaves=FALSE,

           label_branch_points=FALSE,

           graph_label_size=1.5)

p14

```

```{r}

png(file="Source.png",width=1000, height=700)

p14

dev.off()

```

```{r}

library(scales)

hex <- hue_pal()(14)

```