suppressPackageStartupMessages(library(Seurat))
suppressPackageStartupMessages(library(monocle3))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(dplyr))

# ED_Fig5B ----------------------------------------------------------------

## load scRNAseq data from public dataset GSE162950
load("Scarfo_HEC2023/scRNAseq/seurat_object.Rdata")
EHT_scorecard <- c("CXCR4","SOX17","GJA5","MECOM","HOXA9","SPN","CD44","ITGA2B","HLF","GFI1","MLLT3","KCNK17","MYB","STAT5A","SMIM24","RAB27B","SPINK2")
## Select HE cells
sample <- SetIdent(sample, value = sample@meta.data$seurat_clusters)
index1 = GetAssayData(sample, slot="counts")["RUNX1",]>0
index2 = GetAssayData(sample, slot="counts")["CDH5",]>0
index3 = GetAssayData(sample, slot="counts")["PTPRC",]==0
index4 = GetAssayData(sample, slot="counts")["FCGR2B",]==0
bars = colnames(sample)[Idents(sample) %in% c(0,3)]
bars1 = colnames(sample)[index1 & index2 & index3 & index4]
bars2 = intersect(bars, bars1)
## Establish identity for FCGR2B- HE cells
Idents(sample, cells=bars2) = "FCGR2B=0"
index5 = GetAssayData(sample, slot="counts")["FCGR2B",]>0
bars3 = colnames(sample)[index1 & index2 & index3 & index5]
bars4 = intersect(bars, bars3)
## Establish identity for FCGR2B+ HE cells
Idents(sample, cells=bars4) = "FCGR2B_exp"
sample$CellType <- Idents(sample)
sample_sub <- subset(sample, CellType=="FCGR2B_exp" | CellType=="FCGR2B=0")
DotPlot(sample_sub, 
             features=EHT_scorecard, 
             cols=c("grey90","red3"),
             dot.scale = 10) +
  coord_flip()

# ED_Fig5C ------------------------------------------------------------------

## load Seurat data
sample <- readRDS("Scarfo_HEC2023/scRNAseq/WNTd_HE.rds")
sample@meta.data[["Cell.Annotations"]] <- recode_factor(sample@meta.data[["RNA_snn_res.0.6"]],
                                                           "0" = "HECs", 
                                                           "1" = "HECs", 
                                                           "2" = "HECs", 
                                                           "3" = "Venous cells", 
                                                           "4" = "Venous cells", 
                                                           "5" = "Arterial ECs", 
                                                           "6" = "Arterial ECs", 
                                                           "7" = "Venous cells", 
                                                           "8" = "ECs (M phase)",
                                                           "9" = "Venous cells",
                                                           "10" = "EndoMT cells",
                                                           "11" = "HECs",
                                                           "12" = "Arterial ECs",
                                                           "13" = "Lymphatic ECs",
                                                           "14" = "Arterial ECs",
                                                           "15" = "Allantois/Placenta",
                                                           "16" = "HECs",
                                                           "17" = "HECs",
                                                           "18" = "EndoMT cells",
                                                           "19" = "SM cells",
                                                           "20" = "EndoMT cells",
                                                           "21" = "EndoMT cells")

DimPlot(object = sample,
        reduction = "umap",
        group.by = "Cell.Annotations",
        label = F,
        repel = T,
        pt.size = 0.8,
        label.size = 6,
        cols = c('HECs' = "#377EB8", 
                 'Lymphatic ECs' = "#4DAF4A",
                 'Arterial ECs' = "#E41A1C",
                 'Venous cells' = "#FF7F00",
                 'EndoMT cells' = "#E6AB02",
                 'Allantois/Placenta' = "#984EA3",
                 'SM cells' = "#A6761D",
                 "ECs (M phase)" = "#A50026"))


# ED_Fig5D - ED_Fig5E ------------------------------------------------------------------
FeaturePlot(object = sample,
            features = c("RUNX1", "FCGR2B"),
            cols = c("grey90","blue"),
            order = T,
            pt.size = 0.1)

# ED_Fig5F ------------------------------------------------------------------

runx1 <- readRDS("Scarfo_HEC2023/scRNAseq/RUNX1_clusters.rds")

## Monocle3
Idents(runx1) <- "RNA_snn_res.0.6"
start = WhichCells(runx1,idents = "0")
expression_matrix <- runx1@assays[["RNA"]]@counts
cell_metadata <- runx1@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(runx1))
rownames(gene_annotation) <- gene_annotation$gene_short_name
cds <- new_cell_data_set(expression_data = expression_matrix, 
                         cell_metadata = cell_metadata, 
                         gene_metadata = gene_annotation)
cds <- preprocess_cds(cds, num_dim = 50)
cds <- reduce_dimension(cds)
## Construct and assign the made up partition
recreate.partition <- c(rep(1, length(cds@colData@rownames)))
names(recreate.partition) <- cds@colData@rownames
recreate.partition <- as.factor(recreate.partition)
cds@clusters@listData[["UMAP"]][["partitions"]] <- recreate.partition
## Assign the cluster info
list_cluster <- runx1@meta.data[["RNA_snn_res.0.6"]]
names(list_cluster) <- runx1@assays[["RNA"]]@data@Dimnames[[2]]
cds@clusters@listData[["UMAP"]][["clusters"]] <- list_cluster
cds@clusters@listData[["UMAP"]][["louvain_res"]] <- "NA"
## Assign UMAP coordinate
reducedDims(cds)[["UMAP"]] <- runx1@reductions[["umap"]]@cell.embeddings
### Reset colnames and rownames (consequence of UMAP replacement)
rownames(cds@principal_graph_aux[['UMAP']]$dp_mst) <- NULL
## Learn Graph
cds <- learn_graph(cds = cds,use_partition = T,learn_graph_control=list(ncenter=220,minimal_branch_len=15),verbose = T)
cds <- order_cells(cds, root_cells = start)
plot_cells(cds,
           color_cells_by = "pseudotime",
           label_groups_by_cluster=FALSE,
           label_leaves=T,
           label_branch_points=T,
           group_label_size = 8,
           graph_label_size = 3,
           cell_size = 1,
           trajectory_graph_segment_size = 1)


# ED_Fig5G ------------------------------------------------------------------

cds <- estimate_size_factors(cds)
gi <- c("H19","KCNK17","RUNX1", "MYB", "SPN")
cds_gi <- cds[rowData(cds)$gene_short_name %in% gi,]
plot_genes_in_pseudotime(cds_subset = cds_gi, 
                         ncol = 3,
                         color_cells_by = "RNA_snn_res.0.6", 
                         min_expr = NULL,
                         cell_size = 1)