Fig5_plots.R

suppressPackageStartupMessages(library(Seurat))
suppressPackageStartupMessages(library(monocle3))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(reshape2))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(patchwork))
suppressPackageStartupMessages(library(clusterProfiler))

## load Seurat data
sample <- readRDS("Scarfo_HEC2023/scRNAseq/WNTd_HE.rds")

# Fig5A ------------------------------------------------------------------

DimPlot(object = sample,
        reduction = "umap",
        group.by = "RNA_snn_res.0.6",
        label = T,
        repel = T,
        pt.size = 0.8,
        label.size = 8,
        cols = c("0" = "#F8766D",
                 "1" = "#B79F00",
                 "2" = "#00BA38",
                 "3" = "#EC8239",
                 "4"= "#00B81B",
                 "5"= "#7C96FF",
                 "6"= "#00C085",
                 "7" = "#00C1A7",
                 "8" = "#00BFC4",
                 "9" = "#EF67EB",
                 "10" = "#00B2F3",
                 "11"= "#00BFC4",
                 "12"= "#00A6FF",
                 "13"= "#DB8E00",
                 "14"= "#FD61D3",
                 "15"= "#FF63B6",
                 "16"= "#619CFF",
                 "17"= "#F564E3",
                 "18"= "#00BADE",
                 "19"= "#B385FF",
                 "20"= "#D874FD",
                 "21"= "#FF6B94")) +
  ggtitle("") +
  xlab("UMAP1") +
  ylab("UMAP2") +
  theme(legend.text = element_text(size = 16*96/72),
        axis.text.x = element_text(color = "black", family = "Arial",size = 14*96/72),
        axis.text.y = element_text(color = "black", family = "Arial",size = 14*96/72),
        axis.title.x = element_text(color = "black", family = "Arial",size = 16*96/72),
        axis.title.y = element_text(color = "black", family = "Arial",size = 16*96/72),
        aspect.ratio = 1,
        panel.background = element_blank()) +
  guides(color = guide_legend(override.aes = list(size=4)))


## Seurat markers
markers <- FindAllMarkers(sample, 
                          only.pos = T, 
                          min.pct = 0.1, 
                          logfc.threshold = 0.25)
## define RUNX1 clusters
runx1.clusters <- subset(markers, markers$gene=="RUNX1")
runx1.clusters <- droplevels(runx1.clusters$cluster) 
## subset RUNX1 clusters into new object
runx1 <- subset(sample, idents = runx1.clusters)
DimPlot(object = runx1,
        reduction = "umap",
        group.by = "RNA_snn_res.0.6",
        label = T,
        repel = T,
        pt.size = 0.8,
        label.size = 8)


# Fig5B ----------------------------------------------------------------

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

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 = 1)


# Fig5C -------------------------------------------------------------------

## differential markers
degs_clu11_vs_clu0_1_2 <- FindMarkers(object = runx1, ident.1 = 11, ident.2 = c(0,1,2), only.pos = FALSE, min.pct = 0, test.use = "wilcox", logfc.threshold = 0, min.cells.group = 0)
degs_clu11_vs_clu16_17 <- FindMarkers(object = runx1, ident.1 = 11, ident.2 = c(16,17), only.pos = FALSE, min.pct = 0, test.use = "wilcox", logfc.threshold = 0, min.cells.group = 0)

## GSEA
adj.pval.thr <- 0.05
gmt.obj <- read.gmt("Scarfo_HEC2023/BulkRNAseq_1/c5.all.v7.2.symbols.gmt")
comp <- c("clu11_vs_clu0_1_2", "clu11_vs_clu16_17")
for (i in comp) {
  degs <- get(paste("degs",i,sep="_"))
  gene.res.sorted <- degs[order(degs$avg_logFC, decreasing = TRUE), ]
  # LogFC with Symbols
  gene.res.logfc.symbol <- as.vector(gene.res.sorted$avg_logFC)
  names(gene.res.logfc.symbol) <- row.names(gene.res.sorted)
  contr.gsea <- GSEA(gene.res.logfc.symbol, 
                     TERM2GENE = gmt.obj,
                     nPerm = 10000, 
                     pvalueCutoff = adj.pval.thr)
  contr.gsea.res <- as.data.frame(contr.gsea)
  assign(paste0("contr.gsea.res_",i), contr.gsea.res)
  if (i == "clu11_vs_clu0_1_2") {
    terms <- c("GO_CELL_CELL_JUNCTION_ORGANIZATION",
               "GO_ACTIN_CYTOSKELETON",
               "GO_EXTRACELLULAR_MATRIX_STRUCTURAL_CONSTITUENT",
               "GO_STRUCTURAL_CONSTITUENT_OF_RIBOSOME",
               "GO_RIBOSOME",
               "GO_TRANSLATIONAL_INITIATION")
  } else {
    terms <- c("GO_MITOTIC_NUCLEAR_DIVISION",
               "GO_CELL_DIVISION",
               "GO_MITOTIC_CELL_CYCLE",
               "GO_REGULATION_OF_NOTCH_SIGNALING_PATHWAY",
               "GO_CYTOSOLIC_RIBOSOME",
               "GO_ENDOTHELIAL_CELL_MIGRATION")
  }
  ti <- contr.gsea.res[c(contr.gsea.res$Description %in% terms), ]
  ti <- mutate(ti, color = case_when(ti$NES > 0 ~ "Upregulated",
                                     ti$NES < 0 ~ "Downregulated"))
  ti <- arrange(ti, NES)
  ti$Description <- factor(ti$Description, levels = ti$Description)
  p <- ggplot(data=ti, aes(x=NES,  y=Description, fill = color)) +
    geom_bar(stat="identity", width = 0.7) +
    scale_fill_manual(values = c("Upregulated" = "red","Downregulated" = "blue"),name = "") +
    ylab("Terms") +
    ggtitle(paste(i)) +
    theme_bw() +
    theme(
      axis.title.x = element_text(colour="black",family = "Arial", size=12*96/72),
      axis.text.x = element_text(colour="black",family = "Arial", size=10*96/72),
      axis.title.y = element_text(colour="black",family = "Arial", size=12*96/72),
      axis.text.y = element_text(colour="black",family = "Arial", size=8*96/72),
      aspect.ratio = 0.8)
  assign(paste0("p_",i), p)
}

# Fig5D -------------------------------------------------------------------

runx1@meta.data[["Phase"]] <- recode_factor(runx1@meta.data[["Phase"]],
                                            "G1" = "G1", 
                                            "S" = "S", 
                                            "G2M" = "G2/M")
## grouping clusters 0,1,2 and cluster 16,17
runx1@meta.data[["RNA_snn_res.0.6"]] <- recode_factor(runx1@meta.data[["RNA_snn_res.0.6"]],
                                                          "0" = "0,1,2", 
                                                          "1" = "0,1,2", 
                                                          "2" = "0,1,2",
                                                          "11" = "11",
                                                          "16" = "16,17",
                                                          "17" = "16,17")
cc_umap <- DimPlot(object = runx1,
                   reduction = "umap",
                   group.by = "Phase",
                   split.by = "RNA_snn_res.0.6",
                   label = F,
                   repel = T,
                   pt.size = 0.8,
                   label.size = 8,
                   cols = c("G1" = "#E6AB02", 
                            "S" = "#7570B3", 
                            "G2/M" = "#66A61E")) +
  ggtitle("") +
  xlab("UMAP1") +
  ylab("UMAP2") +
  theme(legend.text = element_text(size = 12*96/72),
        axis.text.x = element_text(color = "black", family = "Arial",size = 10*96/72),
        axis.text.y = element_text(color = "black", family = "Arial",size = 10*96/72),
        axis.title.x = element_text(color = "black", family = "Arial",size = 12*96/72),
        axis.title.y = element_text(color = "black", family = "Arial",size = 12*96/72),
        aspect.ratio = 1,
        plot.margin = margin(0, 0, 0, 0, "pt"),
        panel.background = element_blank()) +
  guides(color = guide_legend(override.aes = list(size=4)))

## average cell counts per cluster
ggp_obj_clu <- melt(table(runx1$Phase, runx1$RNA_snn_res.0.6))
colnames(ggp_obj_clu) <- c("Phase","Cluster", "Count")
cc.clusters <- c("0,1,2","11","16,17")
ggp_obj_clu <- ggp_obj_clu[c(ggp_obj_clu$Cluster %in% cc.clusters), ]
for (i in cc.clusters) {
  totalB <- subset(ggp_obj_clu, ggp_obj_clu$Cluster==i)
  totalB <- totalB %>% 
    mutate(percentage = Count/sum(Count))
  totalB$percentage <- totalB$percentage*100
  totalB$fraction = totalB$Count / sum(totalB$Count)
  # Compute the cumulative percentages
  totalB$ymax <- cumsum(totalB$fraction)
  totalB$ymin <- c(0, head(totalB$ymax, n=-1))
  totalB$labelPosition <- (totalB$ymax + totalB$ymin) / 2
  # Compute a good label
  totalB$label <- paste0(totalB$Phase, ":","\n",round(totalB$percentage, 2), "%")
  assign(paste0("totalB_clu_", i), totalB)
  p <- ggplot(totalB, aes(ymax=ymax, ymin=ymin, xmax=4, xmin=3, fill=Phase)) +
    geom_rect() +
    geom_text(x=5, aes(y=labelPosition, label=label), size=5) +
    #scale_fill_brewer(palette = "Set1") +
    scale_fill_manual(values = c("G1" = "#E6AB02", 
                                 "S" = "#7570B3", 
                                 "G2/M" = "#66A61E")) +
    coord_polar(theta="y") +
    xlim(c(2, 5)) +
    theme_void() +
    theme(plot.margin = margin(0, 0, 0, 0, "pt")) +
    theme(legend.position = "none")
  assign(paste0("p_clu_",i), p)
}
p_final <- cc_umap / (`p_clu_0,1,2` | p_clu_11 | `p_clu_16,17`)

# Fig5E -------------------------------------------------------------------

cds <- order_cells(cds, root_cells = start)
cds <- estimate_size_factors(cds)
gi <- c("FCGR2B","HES1","HEY1", "HEY2")
cds_gi <- cds[rowData(cds)$gene_short_name %in% gi,]
plot_genes_in_pseudotime(cds_subset = cds_gi, 
                         ncol = 2,
                         color_cells_by = "RNA_snn_res.0.6", 
                         min_expr = NULL,
                         cell_size = 1)