484 lines (400 with data), 7.8 kB

---
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)
```