Raw File
---
title: "Explant scRNAseq datasets analysis"
subtitle: "Scanpy PAGA"
author: "Renaud Mevel"
output:
  html_document:
    self_contained: yes
    toc: true
    toc_float: true
    df_print: paged
    number_sections: false
editor_options: 
  chunk_output_type: console
---
```{r setup, echo=FALSE, message=FALSE, results='hide'}
library(knitr)
knitr::opts_chunk$set(cache=TRUE, error=FALSE, cache.lazy = TRUE)
```

## Objective

Import and plot data generated in scanpy, including ForceAtlas2/Neighborhood graph and PAGA.

## Prepare the environment

```{r , warning=FALSE, message=FALSE}
# Data wrangling
library(plyr)
library(dplyr)
library(tidyverse)
library(data.table)

# Plots
library(gridExtra)
library(ggpubr)
library(viridis)
# sc
library(Seurat)
library(sctransform)
library(MAST)
library(org.Mm.eg.db)
library(DoubletFinder)
library(slingshot)

# GO
library(gprofiler2)

# Palettes
pal.sa <- c("#DC2B19", "#F7A400", "#72ABF8", "#434CA7")
pal.sa.vi <- viridis(4,begin = 0.05, end = 0.95)

pal.cl <- c(
  "#F09900", #C0
  "#F1BA88", #C1
  "#D63A69", #C2
  "#9DDED6", #C3
  "#8DD593", #C4
  "#428654", #C5
  "#B2B8E0", #C6
  "#4A6FE3", #C7
  "#1037AA"  #C8 
)

# Directories
setwd(dir = "~/set-directory/")
pdf.dir <- "~/set-directory/"
fig.dir <- "~/set-directory/"

# Functions
source("Explant_functions.R")

# Seed
set.seed(1)
```

## Functions 1

```{r}
plotPAGAClustGraph <- function(embedding, edges, thresh = 0, colour = "in_silico_clusters") {

    is_discrete <- is.factor(embedding[[colour]])

    gg <- ggplot(embedding, aes(x = X, y = Y))

    if (is_discrete) {
        gg <- gg +
            geom_segment(data = filter(edges, Connectivity > thresh),
                         aes(x = FromX, y = FromY, xend = ToX, yend = ToY, #alpha = Connectivity
                             colour = Connectivity), size = (edges$Connectivity)*10 ) +
            scale_colour_viridis_c(direction = 1, begin = 0.2, end = 0.6, option="inferno")
    } else {
        gg <- gg +
            geom_segment(data = filter(edges, Connectivity > thresh),
                         aes(x = FromX, y = FromY, xend = ToX, yend = ToY,
                             alpha = Connectivity), colour = "grey30") +
            scale_alpha(limits = c(0, 1)) +
            scale_fill_viridis_c()
    }

    gg <- gg +
        geom_point(aes(fill = !!ensym(colour), size = Size), shape = 21) +
        scale_fill_manual(values = pal.cl) +
        geom_text(aes(label = in_silico_clusters)) +
        scale_size(range = c(5, 20)) +
        theme_void() +
        theme(legend.position = "none")

    return(gg)
}


###  -----------------------------------------

plotPAGACellGraph <- function(embedding, edges, thresh = 0, colour = "in_silico_clusters", label = FALSE) {

    is_discrete <- is.factor(embedding[[colour]])

    gg <- ggplot(embedding, aes(x = X, y = Y, colour = !!ensym(colour))) +
        geom_segment(data = filter(edges, Connectivity > thresh),
                     aes(x = FromX, y = FromY, xend = ToX, yend = ToY),
                     size = 0.1, colour = "grey50") +
        geom_point(size = 0.5, aes(colour = in_silico_clusters)) +
        scale_colour_manual(values = pal.cl) +
        theme_void() +
        theme(legend.position = "none")

    if (!is_discrete) {
        gg <- gg + scale_color_viridis_c()
    }

    if (label) {
        clust_data <- embedding %>%
            group_by(in_silico_clusters) %>%
            summarise(X = mean(X),
                      Y = mean(Y))

        gg <- gg +
          #geom_point(data = clust_data, aes(fill = in_silico_clusters),
          #             size = 10, shape = 21, colour = "white") +
          scale_fill_manual(values = pal.cl) +
          geom_text(data = clust_data, aes(label = in_silico_clusters),
                      colour = "black")
    }

    return(gg)
}

###  -----------------------------------------

plotPAGACompare <- function(clust_embedding, clust_edges, clust_thresh = 0,
                            cell_embedding, cell_edges, cell_thresh = 0,
                            colour = NA, label = FALSE) {

    clusts <- plotPAGAClustGraph(clust_embedding, clust_edges, clust_thresh,
                                 colour)

    cells <- plotPAGACellGraph(cell_embedding, cell_edges, cell_thresh,
                               colour, label)

    cowplot::plot_grid(clusts, cells, nrow = 1)
}
```


## Functions 2 

Function to plot genes in PAGA and FA
```{r}

##### In PAGA space ---------------------------------------------------------
plotPAGAgene <- function(scdata, cl.embedding, cl.edges, sc.embedding, gene = "Runx1", legend=TRUE, thresh=0, show.clusters=FALSE, rm.cl = NA) {
  
  # remove cluster if asked
  if (!is.na(rm.cl)) { 
    cl.embedding <- dplyr::filter(cl.embedding, in_silico_clusters != rm.cl)
    cl.edges <- dplyr::filter(cl.edges, From != rm.cl, To != rm.cl)
    }
  
  sc.embedding[[gene]] <- logcounts(scdata)[gene, ]
  
  paga_gene <- sc.embedding %>%
    dplyr::select(-Cell, -X, -Y) %>%
    group_by(in_silico_clusters) %>%
    summarise_all(mean)
  
  paga_gene[, gene][paga_gene[, gene] == 0] <- NA 
  
  max.scale = max(paga_gene[[gene]])
  min.scale = min(paga_gene[[gene]], na.rm = TRUE)
  
  paga_embedding <- left_join(cl.embedding, paga_gene, by = "in_silico_clusters")
  
  cl.edges <- cl.edges %>% mutate(Connectivity = ifelse(Connectivity < thresh, NA, Connectivity))
  
  g <- ggplot(paga_embedding, aes(x = X, y = Y)) +
  geom_segment(data = cl.edges, aes(x = FromX, y = FromY, xend = ToX, yend = ToY), size = (cl.edges$Connectivity)*10, alpha=0.5) +
  geom_point(aes(fill = !!ensym(gene), size = Size), shape = 21, colour="grey30") +
  scale_fill_viridis(option = "plasma", na.value = "grey80", begin = 0, end = 0.9, limits = c(min.scale, max.scale)) +
  scale_size(range = c(6, 12), guide = FALSE) +
  theme_void()
  
  if (legend==FALSE) { g <- g + theme(legend.position = "none") }
  if (show.clusters==TRUE) { g <- g + geom_text(aes(label = in_silico_clusters)) }
  g
}

##### In ForceAtlas space ---------------------------------------------------------

plotFAgene <- function(scdata, cell.embedding, cell.edges, gene, legend=TRUE, pt.size=0.8, show.connections=TRUE, thresh=0) {
  
  cell.embedding[[gene]] <- logcounts(scdata)[gene, ]
  
  cell.embedding[, 5][cell.embedding[, 5] == 0] <- NA # replace 0 with NA
  
  max.scale = max(cell.embedding[[gene]])
  min.scale = min(cell.embedding[[gene]], na.rm = TRUE)
  
  cell.edges <- cell.edges %>% mutate(Connectivity = ifelse(Connectivity < thresh, NA, Connectivity))
  
  g <- ggplot(cell.embedding, aes(x = X, y = Y))
  
    if (show.connections==TRUE) { 
  g <- g + geom_segment(data = cell.edges, aes(x = FromX, y = FromY, xend = ToX, yend = ToY, size = Connectivity, alpha = Connectivity), colour="grey50")
  g <- g + scale_alpha_continuous(limits = c(0, 1), range = c(0.2, 1), guide = FALSE)
  g <- g + scale_size(range = c(0.1, 0.5), guide = FALSE)
  }
  
  g <- g + geom_point(aes(colour = !!ensym(gene)), size=pt.size)
  g <- g + scale_colour_viridis_c(option = "plasma", na.value = "grey80", end = 1, limits = c(min.scale, max.scale))
  g <- g + theme_void()

  if (legend==FALSE) { g <- g + theme(legend.position = "none") }
  
  g
}
```


## Load data

```{r}
seu <- readRDS(file.path("r_save/sce_final.diffmap.rds"))

# Drop C9 level
relevel <- c("C0", "C1", "C2", "C3", "C4", "C5", "C6", "C7", "C8")
Idents(seu) <- factor(Idents(seu), levels = relevel)
seu$in_silico_clusters <- Idents(seu)

sce <- as.SingleCellExperiment(seu, assay = "RNA")
```

## Load PAGA

### Load embeddings

```{r}
# cluster embeddings--------------------------------------------------------------------------
clust_embedding <- read_csv("scanpy_paga/cluster_embedding.csv",
                            col_types = cols(.default = col_double())) %>%
  mutate(Size = as.numeric(table(colData(sce)$in_silico_clusters))) %>%
  dplyr::rename(Cluster = in_silico_clusters) %>% 
  mutate(Cluster = factor(Cluster)) %>% 
  mutate(in_silico_clusters = factor(names(table(colData(sce)$in_silico_clusters)), levels = levels(colData(sce)$in_silico_clusters)))

# cluster edges--------------------------------------------------------------------------------
clust_edges <- read_csv("scanpy_paga/cluster_edges.csv",
                        col_types = cols(.default = col_double()))  %>%
  mutate(To = factor(To, levels = levels(clust_embedding$Cluster)),
         From = factor(From, levels = levels(clust_embedding$Cluster))) %>%
    left_join(clust_embedding, by = c("From" = "Cluster")) %>%
    dplyr::rename(FromX = X, 
                  FromY = Y,
                  From_in_silico_clusters = in_silico_clusters) %>%  
    dplyr::select(-Size) %>%
    left_join(clust_embedding, by = c("To" = "Cluster")) %>%
    dplyr::rename(ToX = X, 
                  ToY = Y, 
                  To_in_silico_clusters = in_silico_clusters) %>%
    dplyr::select(-Size)

# remove 'Cluster'----------------------------------------------------------------------------
clust_embedding <- dplyr::select(clust_embedding, -Cluster)
clust_edges <- clust_edges %>% 
  dplyr::mutate(From = From_in_silico_clusters,
                To = To_in_silico_clusters) %>% 
  dplyr::select(-From_in_silico_clusters, -To_in_silico_clusters)

# cell embeddings----------------------------------------------------------------------------
cell_embedding <- read_csv("scanpy_paga/cell_embedding.csv",
                           col_types = cols(.default = col_double(), Cell = col_character())
                           ) %>%
  mutate(in_silico_clusters = colData(sce)$in_silico_clusters)

# cell edges ---------------------------------------------------------------------------------
cell_edges <- read_csv("scanpy_paga/cell_edges.csv",
                        col_types = cols(.default = col_double(), From = col_character(), To = col_character())
                       )  %>%
  dplyr::left_join(cell_embedding, by = c("From" = "Cell")) %>%
  dplyr::rename(FromX = X, FromY = Y) %>%
  dplyr::select(-in_silico_clusters) %>%
  dplyr::left_join(cell_embedding, by = c("To" = "Cell")) %>%
  dplyr::rename(ToX = X, ToY = Y) %>%
  dplyr::select(-in_silico_clusters)

# pseudotime ----------------------------------------------------------------------------
dpt_pseudotime <- read_csv("scanpy_paga/dpt_pseudotime.csv", 
                           col_types = cols(.default = col_character(), dpt_pseudotime = col_double())) %>% 
  tibble::column_to_rownames(var = "X1")
```

### FA > Seurat

Add FA layout to the seurat object
```{r}
FA <- cell_embedding %>% 
  column_to_rownames("Cell") %>% 
  dplyr::select(-in_silico_clusters) %>% 
  dplyr::rename(FA1 = X, FA2 = Y) %>% 
  as.matrix()

seu[["FA"]] <- CreateDimReducObject(embeddings = FA, key = "FA_", assay = DefaultAssay(seu))
DimPlot(seu, reduction = "FA", cols = pal.cl)
```

Add dpt_pseudotime
```{r}
seu[["dpt_pseudotime"]] <- dpt_pseudotime
FeaturePlot(seu, reduction = "FA", features = "dpt_pseudotime") + scale_colour_viridis(option = "viridis")
```

Convert back to SingleCellExperiment
```{r}
sce <- as.SingleCellExperiment(seu, assay = "RNA")
```

## Save obkect

```{r}
saveRDS(seu, file =  "r_save/sce_final_paga_fa_dpt.rds")
```

## Visualise Scanpy output

### Cluster graph

```{r}
plotPAGAClustGraph(clust_embedding, clust_edges, thresh = 0)
```

### Edges threshold

Number of selected edges for different threshold connectivities.
```{r}
plot_data <- tibble(
    Threshold = seq(0, 1, 0.01)
) %>%
    mutate(Edges = map_int(Threshold, function(thresh) {
        sum(clust_edges$Connectivity > thresh)
    }))

con_thresh <- 0.7

ggplot(plot_data, aes(x = Threshold, y = Edges)) +
    geom_point() +
    geom_line() +
    geom_vline(xintercept = con_thresh, colour = "red") +
    xlab("Connectivity threshold") +
    ylab("Number of edges") +
    theme_minimal()
```

### Cell graph

```{r}
plotPAGACellGraph(cell_embedding, cell_edges, thresh = 0.0)
```

### Compare

```{r}
plotPAGACompare(clust_embedding, clust_edges, clust_thresh = 0,
                cell_embedding, cell_edges, cell_thresh = 0, colour = "in_silico_clusters")
```

### Gene expression

Bac a sable
In PAGA
```{r}
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Nkx3-1")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Runx1")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Tacstd2")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Mki67")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Psca")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Krt4")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Nupr1")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Meis2")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Shh")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Ar")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Foxa1")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Etv1")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Etv4")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Gata2")

plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Trp63")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Apoe")

(
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Krt5") /
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Krt14") /
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Krt7") /
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Krt19") 
) | (
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Krt8")/
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Krt18")/
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Krt17")/
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Krt4")
)

plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Gsto1")

plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Cd24a")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Itga6")

plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Dpp4")
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Krt20")

plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = "Sox2")
```

In ForceAtlas
```{r}
plotFAgene(sce, cell_embedding, cell_edges, gene = "Runx1", pt.size = 1, thresh=0, show.connections = TRUE, legend=FALSE)
plotFAgene(sce, cell_embedding, cell_edges, gene = "Nkx3-1", pt.size = 1, thresh=0, show.connections = TRUE, legend=FALSE)
plotFAgene(sce, cell_embedding, cell_edges, gene = "Shh", pt.size = 1, thresh=0, show.connections = TRUE)
```

#### Genes in PAGA

```{r}
gene <- "Runx1"
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Nkx3-1"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Tacstd2"
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Krt4"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Psca"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Psca"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Nupr1"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Nupr1"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Krt5"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Krt14"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Krt8"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Krt18"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Krt7"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Trp63"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Dcn"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Apoe"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Shh"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Notch1"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Meis2"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Mki67"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Top2a"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Hoxd13"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Hoxb13"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Sox9"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Ly6d"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Ly6a"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Ly6e"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Krt19"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Ar"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Foxa1"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Vcan"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Itga6"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

gene <- "Cd24a"
#plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene)
#ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, ".png")), device = "png", width = 5, height = 5, dpi = 300)
plotPAGAgene(sce, clust_embedding, clust_edges, cell_embedding, gene = gene, legend = FALSE)
ggsave(filename = file.path(fig.dir, paste0("6_paga/genes/PAGA_", gene, "_minimal",".png")), device = "png", width = 4, height = 4, dpi = 300)

```

#### Runx1/Nkx3-1

Plot cells which have co-expression of Runx1 and Nkx3-1

```{r}
known_genes <- c("Runx1", "Nkx3-1")

sc_embedding <- cell_embedding

for (gene in known_genes) {
    sc_embedding[[gene]] <- logcounts(sce)[gene, ]
}

sc_embedding <- sc_embedding %>% 
  mutate(cutoff = 0.1) %>% 
  mutate(coexpr = ifelse(Runx1 > cutoff & `Nkx3-1` > cutoff, "yes", "no"))

ggplot(sc_embedding, aes(x = X, y = Y)) +
  geom_point(aes(colour = coexpr), size=1, alpha=0.9) + 
  scale_colour_manual(values = c("grey80", "darkblue")) + 
  theme_void()


meta = seu@meta.data
known_genes <- c("Runx1", "Nkx3-1")
for (gene in known_genes) { meta[[gene]] <- logcounts(sce)[gene, ] }
meta <- meta %>% 
  mutate(cutoff = 0) %>% 
  mutate(coexpr = ifelse(Runx1 > cutoff & `Nkx3-1` > cutoff, "yes", "no"))

# Propotion by day
stat1 <- 
  meta %>% 
  dplyr::group_by(sample, coexpr) %>% 
  dplyr::summarise(count=n()) %>% 
  dplyr::mutate(perc=(count/sum(count)*100))

ggplot(stat1, aes(x = sample, y = perc, fill = coexpr, label = count)) +
  geom_bar(stat="identity", width = 0.5) +
  geom_text(size = 4, position = position_stack(vjust = 0.5)) +
  scale_fill_manual(values = c("grey80","#6D8FBA")) +
  labs(x="Days", y="Percentage", fill="Co-expression") +
   theme(
     panel.grid.minor = element_blank(),
     panel.grid.major = element_blank(),
     panel.background = element_blank(),
     panel.border = element_rect(colour = "grey30", fill=NA, size=0.5),
     strip.background = element_rect(fill="white", colour = "white")
  ) + rotate()

# Proportion by cluster
stat2 <- 
  meta %>% 
  dplyr::group_by(in_silico_clusters, coexpr) %>% 
  dplyr::summarise(count=n()) %>% 
  dplyr::mutate(perc=(count/sum(count)*100))

ggplot(stat2, aes(x = in_silico_clusters, y = perc, fill = coexpr, label = count)) +
  geom_bar(stat="identity", width = 0.5) +
  geom_text(size = 4, position = position_stack(vjust = 0.5)) +
  scale_fill_manual(values = c("grey80","#6D8FBA")) +
  labs(x="Clusters", y="Percentage", fill="Co-expression") +
   theme(
     panel.grid.minor = element_blank(),
     panel.grid.major = element_blank(),
     panel.background = element_blank(),
     panel.border = element_rect(colour = "grey30", fill=NA, size=0.5),
     strip.background = element_rect(fill="white", colour = "white")
  ) + rotate()
```

#### Genes in pseudotime

```{r}
ExpUMAP_pseudo1 <- function(sceObje, gNAME){
  g <- gNAME
  #g.exp <- assay(sceObje, "logcounts")[g, ]
  g.exp <- sceObje[["RNA"]]@data[g,]
  df_plot <- data.frame(DPT = sceObje[["dpt_pseudotime"]]$dpt_pseudotime,
                        Gene = g.exp,
                        Population = sceObje$in_silico_clusters)
  #colp_fun <- df_col$Col[df_col$in_silico_clusters %in% sceObje$in_silico_clusters]
  pUMAPf <-
    ggplot(df_plot,
           aes(x=DPT, y=Gene, colour=Population)) +
    geom_point(size=1.5) +
    geom_smooth(colour="red", method="loess", se=TRUE, size=0.5) +
    scale_colour_manual(values=pal.cl) +
    theme_bw() +
    theme(legend.position = "none",
          panel.grid.major = element_line(linetype = "blank"), 
          panel.grid.minor = element_line(linetype = "blank")) +
    #geom_vline(xintercept = ps_ar_1, linetype = "dashed", size=0.5, colour="grey80") +
    labs(subtitle = gNAME) +
    ylab("Expression") +
    xlab("Pseudotime") #+ xlim(c(0.6,1))
  
  print(pUMAPf)
}

ExpUMAP_pseudo1(seu,"Runx1")
ExpUMAP_pseudo1(seu,"Nkx3-1")
ExpUMAP_pseudo1(seu,"Tacstd2")
```

#### Violin plots

```{r}
t = theme(axis.text.x = element_text(angle = 0))
VlnPlot(seu, features = c("Nkx3-1"), cols = pal.cl, pt.size = 0, group.by = "in_silico_clusters")# + NoLegend() + t
#ggsave(filename = file.path(fig.dir, "QC_before_genes.png"), device = "png", width = 4, height = 3, dpi = 300)
```

#### Heatmap

```{r}
geneList = c( 'Shh', "Foxa1", 'Notch1', 'Sox9', 'Hoxb13', 'Hoxd13', 'Ar', "Nkx3-1",  "Trp63", 'Wnt5a',
              'Krt5', 'Krt14', "Krt8", "Krt18", "Krt19", "Krt7",
              "Dcn", "Apoe", "Ly6a", "Ly6e",
              'Psca', 'Nupr1', 'Krt4', "Tacstd2", "Ly6d",
              'Mki67', 'Top2a')


DoHeatmap(object = seu, features = as.character(geneList), group.colors = pal.cl, draw.lines = TRUE, lines.width = 10,
          assay = "RNA", slot = "data", size = 4) + scale_fill_viridis(option = "viridis",na.value = "white")
```


## PAGA/FA Figures

### ggplot

Coloured connections
```{r, eval=FALSE}
ggplot(clust_embedding, aes(x = X, y = Y)) +
  geom_segment(data = clust_edges, aes(x = FromX, y = FromY, xend = ToX, yend = ToY, colour = Connectivity, alpha = Connectivity), 
               size = (clust_edges$Connectivity)*8) +
  scale_colour_viridis_c(direction = 1, begin = 0.2, end = 0.6, option="inferno") +
  scale_alpha_continuous(limits = c(0, 1), range = c(0.7, 0.9), guide = FALSE) +
  geom_point(aes(fill = in_silico_clusters, size = Size), shape = 21, colour="white") +
  scale_fill_manual(values = pal.cl) +
  scale_size(range = c(5, 10), guide = FALSE) +
  theme_void()  #+ scale_y_reverse() + coord_flip() 

ggsave(filename = file.path(fig.dir, "6_paga/PAGA_clusters_c.png"), device = "png", width = 5, height = 5, dpi = 300)


ggplot(clust_embedding, aes(x = X, y = Y)) +
  geom_segment(data = clust_edges, aes(x = FromX, y = FromY, xend = ToX, yend = ToY, colour = Connectivity, alpha = Connectivity), 
               size = (clust_edges$Connectivity)*8) +
  scale_colour_viridis_c(direction = 1, begin = 0.2, end = 0.6, option="inferno") +
  scale_alpha_continuous(limits = c(0, 1), range = c(0.7, 0.9), guide = FALSE) +
  geom_point(aes(fill = in_silico_clusters, size = Size), shape = 21, colour="white") +
  scale_fill_manual(values = pal.cl) +
  scale_size(range = c(5, 10), guide = FALSE) +
  theme_void() + theme(legend.position = "none")# + scale_y_reverse() + coord_flip() 

ggsave(filename = file.path(fig.dir, "6_paga/PAGA_clusters_c_nolegend.png"), device = "png", width = 5, height = 5, dpi = 300)

```

PAGA - Black connections
```{r}
ggplot(clust_embedding, aes(x = X, y = Y)) +
  geom_segment(data = clust_edges, aes(x = FromX, y = FromY, xend = ToX, yend = ToY, alpha = Connectivity), 
               size = (clust_edges$Connectivity)*8, lineend="round", linejoin = "round") +
  scale_alpha_continuous(limits = c(0, 1), range = c(1, 1), guide = FALSE) +
  geom_point(aes(fill = in_silico_clusters, size = Size), shape = 21, colour="white") +
  scale_fill_manual(values = pal.cl) +
  scale_size(range = c(5, 10), guide = FALSE) +
  theme_void()  #+ scale_y_reverse() + coord_flip() 

ggsave(filename = file.path(fig.dir, "6_paga/PAGA_clusters.png"), device = "png", width = 5, height = 5, dpi = 300)



ggplot(clust_embedding, aes(x = X, y = Y)) +
  geom_segment(data = clust_edges, aes(x = FromX, y = FromY, xend = ToX, yend = ToY, alpha = 0.5), size = (clust_edges$Connectivity)*10) +
  geom_point(aes(fill = in_silico_clusters, size = Size), shape = 21, colour="white") +
  #geom_segment(data = clust_edges, aes(x = FromX, y = FromY, xend = ToX, yend = ToY, alpha = Connectivity),
  #             size = (clust_edges$Connectivity)*10) +
  #scale_alpha_continuous(limits = c(0, 1), range = c(0.7, 0.7), guide = FALSE) +
  #geom_point(aes(fill = in_silico_clusters, size = Size), shape = 21, colour="white") +
  scale_fill_manual(values = pal.cl) +
  scale_size(range = c(6, 12), guide = FALSE) +
  theme_void() + theme(legend.position = "none")  #+ scale_y_reverse() + coord_flip() 

ggsave(filename = file.path(fig.dir, "6_paga/PAGA_clusters_nolegend.png"), device = "png", width = 4, height = 4, dpi = 300)

```

FA
```{r}
ggplot(cell_embedding, aes(x = X, y = Y)) +
    geom_segment(data = cell_edges,
                 aes(x = FromX, y = FromY, xend = ToX, yend = ToY,
                     size = Connectivity, alpha = Connectivity), colour="grey50") +
    geom_point(aes(colour = in_silico_clusters), size = 0.8) +
    scale_colour_manual(values = pal.cl) +
    scale_alpha_continuous(limits = c(0, 1), range = c(0.05, 0.5), guide = FALSE) +
    scale_size(range = c(0.1, 0.5), guide = FALSE) +
    theme_void() # + scale_y_reverse() + coord_flip() 

ggsave(filename = file.path(fig.dir, "6_paga/PAGA_cells.png"), device = "png", width = 5, height = 5, dpi = 300)

ggplot(cell_embedding, aes(x = X, y = Y)) +
    geom_segment(data = cell_edges,
                 aes(x = FromX, y = FromY, xend = ToX, yend = ToY,
                     size = Connectivity, alpha = Connectivity), colour="grey50") +
    geom_point(aes(colour = in_silico_clusters), size = 0.8) +
    scale_colour_manual(values = pal.cl) +
    scale_alpha_continuous(limits = c(0, 1), range = c(0.05, 0.5), guide = FALSE) +
    scale_size(range = c(0.1, 0.5), guide = FALSE) +
    theme_void() + theme(legend.position = "none")  #+ scale_y_reverse() + coord_flip() 

ggsave(filename = file.path(fig.dir, "6_paga/PAGA_cells_nolegend.png"), device = "png", width = 5, height = 5, dpi = 300)

ggplot(cell_embedding, aes(x = X, y = Y)) +
    geom_point(aes(colour = in_silico_clusters), size = 0.8) +
    scale_colour_manual(values = pal.cl) +
    theme_void() + theme(legend.position = "none")   #+ scale_y_reverse() + coord_flip() 

ggsave(filename = file.path(fig.dir, "6_paga/PAGA_cells_noconnections.png"), device = "png", width = 5, height = 5, dpi = 300)
```

### FA in Seurat

Clusters
```{r}
DimPlot(object = seu, reduction = "FA", group.by = "in_silico_clusters", label=TRUE, label.size = 5, cols = pal.cl, pt.size = 0.8)
ggsave(filename = file.path(fig.dir, "6_paga/FA_clusters.png"), device = "png", width = 6, height = 6, dpi = 300)
DimPlot(object = seu, reduction = "FA", group.by = "in_silico_clusters", cols = pal.cl, pt.size = 1) + NoLegend() + NoAxes()
ggsave(filename = file.path(fig.dir, "6_paga/FA_clusters_nolegend.png"), device = "png", width = 5, height = 5, dpi = 300)
```

Days
```{r}
DimPlot(object = seu, reduction = "FA", group.by = "sample", label=TRUE, label.size = 5, cols = pal.sa.vi, pt.size = 0.8)
ggsave(filename = file.path(fig.dir, "6_paga/FA_days.png"), device = "png", width = 6, height = 6, dpi = 300)
DimPlot(object = seu, reduction = "FA", group.by = "sample", cols = pal.sa.vi, pt.size = 0.8) + NoLegend() + NoAxes()
ggsave(filename = file.path(fig.dir, "6_paga/FA_days_nolegend.png"), device = "png", width = 5, height = 5, dpi = 300)
```

Pseudotime
```{r}
FeaturePlot(seu, reduction = "FA", features = "dpt_pseudotime", pt.size = 0.8) + scale_colour_viridis(option = "viridis")
ggsave(filename = file.path(fig.dir, "6_paga/FA_pseudotime.png"), device = "png", width = 6, height = 6, dpi = 300)

FeaturePlot(seu, reduction = "FA", features = "dpt_pseudotime", pt.size = 0.8) + scale_colour_viridis(option = "viridis") + NoLegend() + NoAxes() + ggtitle(label=NULL)
ggsave(filename = file.path(fig.dir, "6_paga/FA_pseudotime_nolegend.png"), device = "png", width = 5, height = 5, dpi = 300)
```

bac a sable
```{r}
FeaturePlot(seu, reduction = "FA", "Sox9", pt.size = 1) + scale_colour_viridis(option="plasma")
FeaturePlot(seu, reduction = "FA", "Wnt5a", pt.size = 1) + scale_colour_viridis(option="plasma")
FeaturePlot(seu, reduction = "FA", "Sfrp1", pt.size = 1) + scale_colour_viridis(option="plasma")
FeaturePlot(seu, reduction = "FA", "Fgfr2", pt.size = 1) + scale_colour_viridis(option="plasma")
FeaturePlot(seu, reduction = "FA", "Hoxb13", pt.size = 1) + scale_colour_viridis(option="plasma")
```

--------------------------------------------------------------------------------
## Session Information
```{r}
sessionInfo()
```
back to top