4E. Construct trajectory
Related to Figure 2
- Subset data on glomerulus-associated and arterial cell populations
- Recluster and annotate cell types
- Construct single-cell trajectory and evaluate
celltype_names <- levels(obj_named)
trajectory_names <- grep("JG|tSMC|SMC_MC|SMC_RPC|RPC_MC|aLoH|MD|FB|glomEnC|EnC|arterialEnC|fEnC",
celltype_names, value = TRUE)
trajectory_obj <- subset(obj_named, idents = trajectory_names)
trajectory_obj[["RNA"]] <- split(trajectory_obj[["RNA"]], f = trajectory_obj$orig.ident)
trajectory_obj <- IntegrateLayers(object = trajectory_obj,
method = CCAIntegration, orig.reduction = "pca",
new.reduction = "integrated.trj_subset.cca", verbose = FALSE)
trajectory_obj <- IntegrateLayers(
object = trajectory_obj, method = RPCAIntegration,
orig.reduction = "pca", new.reduction = "integrated.trj_subset.rpca",
verbose = FALSE
)
trajectory_obj <- IntegrateLayers(
object = trajectory_obj, method = HarmonyIntegration,
orig.reduction = "pca", new.reduction = "trj_subset.harmony",
verbose = FALSE
)
trajectory_obj[["RNA"]] <- JoinLayers(trajectory_obj[["RNA"]])
DefaultAssay(trajectory_obj) <- "RNA"
trajectory_obj <- FindNeighbors(trajectory_obj, reduction = "integrated.trj_subset.cca", dims = 1:30)
trajectory_obj <- FindClusters(trajectory_obj, resolution = 1, cluster.name = "clusters.trj_subset.cca")
trajectory_obj <- RunUMAP(trajectory_obj, dims = 1:30, reduction = "integrated.trj_subset.cca",
reduction.name = "umap.trj_subset.cca",
min.dist = 0.01, seed.use = 99,
spread = 1)
DimPlot_scCustom(trajectory_obj, reduction = "umap.trj_subset.cca", label=T, pt.size=0.5,
group.by = c("cell_type", "orig.ident", "clusters.trj_subset.cca"))
Idents(trajectory_obj) <- trajectory_obj$clusters.trj_subset.cca
trajectory_obj <- FindSubCluster(trajectory_obj, cluster = "2", algorithm = 3,
graph.name = "RNA_snn", subcluster.name = "sub.clusters.trj_subset.cca",
resolution = 0.25)
Idents(trajectory_obj) <- trajectory_obj$sub.clusters.trj_subset.cca
new_cluster_ids <- c('tSMC_2', 'tSMC_2', 'SMC', 'JG', 'tSMC_1',
'SMC', 'SMC', 'SMC', 'SMC', 'fEnC', 'EnC', 'arterialEnC', 'aLoH', 'tSMC_1',
'MC', 'glomEnC', "MD", 'POD_EnC')
names(new_cluster_ids) <- c('0', '1', '2_0', '2_1', '2_2', '2_3', '2_4', '2_5',
'2_6', '3', '4', '5', '6', '7', '8', '9', '10', '11')
trajectory_obj <- RenameIdents(trajectory_obj, new_cluster_ids)
trajectory_obj[["cell_type_manual"]] <- Idents(trajectory_obj)
table(trajectory_obj$cell_type_manual)
trajectory_umap <- DimPlot_scCustom(trajectory_obj,
reduction = "umap.trj_subset.cca", label=T, pt.size=0.5,
group.by = c("clusters.trj_subset.cca", "cell_type_manual", "orig.ident"))
ggsave(paste0("figures/", sample_name,
"_UMAP_original-trajectory-pops-recluster_", Sys.Date(), ".png"),
plot = trajectory_umap, units = "mm", dpi = 600)
ggsave(paste0("figures/", sample_name,
"_UMAP_original-trajectory-pops-recluster_", Sys.Date(), ".svg"),
plot = trajectory_umap, units = "mm", dpi = 300)
saveRDS(trajectory_obj, paste0("trajectory_obj_", Sys.Date(), ".rds"))
Convert to cell_data_set object
DefaultAssay(trajectory_obj) <- "RNA"
tSMC_2_names <- rownames(trajectory_obj@meta.data[trajectory_obj$cell_type_manual == "tSMC_2",])
exclude_names <- names(trajectory_cds@clusters$UMAP$partitions[tSMC_2_names][trajectory_cds@clusters$UMAP$partitions[tSMC_2_names] == 2])
# Keep only the partition with the cells of interest - Partition 1
partition1_cellnames <- rownames(trajectory_obj@meta.data[trajectory_obj$cell_type_manual %in% c("SMC", "tSMC_1", "tSMC_2", "JG", "fEnC", "EnC", "arterialEnC", "glomEnC", "MC"),])
partition1_cellnames <- partition1_cellnames[partition1_cellnames %ni% exclude_names]
trajectory_cds <- as.cell_data_set(trajectory_obj[,partition1_cellnames],
default.reduction = "umap.trj_subset.cca",
group.by = "cell_type_manual")
# UMAP.TRJ_SUBSET.CCA == UMAP | UMAP == SEURAT_UMAP
reducedDimNames(trajectory_cds) <- c("PCA", "UMAP.UNINTEGRATED",
"INTEGRATED.CCA", "INTEGRATED.RPCA", "HARMONY", "INTEGRATED.MNN",
"SEURAT_UMAP", "UMAP.CCA", "UMAP.RPCA", "UMAP.HARMONY", "UMAP.MNN",
"UMAP", "INTEGRATED.TRJ_SUBSET.CCA",
"INTEGRATED.TRJ_SUBSET.RPCA", "TRJ_SUBSET.HARMONY",
"UMAP.TRJ_SUBSET.RPCA", "UMAP.TRJ_SUBSET.HARMONY")
trajectory_cds <- cluster_cells(trajectory_cds)
trajectory_cds <- preprocess_cds(trajectory_cds)
rowData(trajectory_cds)$gene_name <- rownames(trajectory_cds)
rowData(trajectory_cds)$gene_short_name <- rownames(rowData(trajectory_cds))
colData(trajectory_cds)[['Size_Factor']] <- NULL
trajectory_cds <- estimate_size_factors(trajectory_cds)
colData(trajectory_cds)$Size_Factor
trajectory_cds <- learn_graph(trajectory_cds,
learn_graph_control=list(ncenter=100, nn.k=NULL,
euclidean_distance_ratio = 0.5,
orthogonal_proj_tip = TRUE,
minimal_branch_len=4))
trajectory_plt <- plot_cells(trajectory_cds,
color_cells_by = "cell_type_manual",
label_groups_by_cluster=F,
label_leaves=T,
label_branch_points=T,
cell_size = 1,
cell_stroke = 0,
graph_label_size = 6, group_label_size = 4) +
scale_color_manual(values=pal_named)
saveRDS(trajectory_cds, paste0("trajectory-cds_", Sys.Date(), ".rds"))
ggsave(paste0("figures/monocle3-trajectory-plt_ReAnnotated_", Sys.Date(), ".png"),
plot = trajectory_plt + theme(aspect.ratio=1), width = 180, height = 180, units = "mm", dpi = 600)
ggsave(paste0("figures/monocle3-trajectory-plt_ReAnnotated_", Sys.Date(), ".svg"),
plot = trajectory_plt + theme(aspect.ratio=1), width = 180, height = 180, units = "mm", dpi = 300)
get_earliest_principal_node <- function(cds, starting_cell="SMC"){
cell_ids <- which(colData(cds)[, "cell_type_manual"] == starting_cell)
closest_vertex <-
cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
root_pr_nodes <-
igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names
(which.max(table(closest_vertex[cell_ids,]))))]
root_pr_nodes
}
get_earliest_principal_node(trajectory_cds, "SMC") # Y_5
trajectory_cds <- order_cells(trajectory_cds, root_pr_nodes=c("Y_5"))
pseudotime_plt <- plot_cells(trajectory_cds,
color_cells_by = "pseudotime",
label_groups_by_cluster = TRUE,
label_leaves=TRUE,
label_branch_points=TRUE,
cell_size = 1,
cell_stroke = 0,
graph_label_size=3,
group_label_size = 4)
ggsave(paste0("figures/", sample_name, "_SMC-root_Monocle3-Trajectory-Pseudotime_", Sys.Date(), ".png"),
plot = pseudotime_plt | trajectory_plt, width = 360, height = 180, units = "mm", dpi = 600)
ggsave(paste0("figures/", sample_name, "_SMC-root_Monocle3-Trajectory-Pseudotime_", Sys.Date(), ".svg"),
plot = pseudotime_plt | trajectory_plt, width = 360, height = 180, units = "mm", dpi = 300)
features <- rownames(trajectory_cds)
no_ribo <- features[!grepl("^Rp[sl][[:digit:]]|^Rplp[[:digit:]]|^Rpsa", features)]
no_ribo_mito <- sort(no_ribo[!grepl("^mt-(.*)", no_ribo)])
no_ribo_mito_hemo <- sort(no_ribo_mito[!grepl("^H2-(.*)", no_ribo_mito)])
filtered_features <- sort(no_ribo_mito_hemo[!grepl("^Hb[ab]-(.*)", no_ribo_mito_hemo)])
trajectory_cds_filtered <- trajectory_cds[rownames(trajectory_cds) %in% filtered_features, ]
saveRDS(trajectory_cds_filtered, paste0("trajectory-cds-filtered_", Sys.Date(), ".rds"))
trajectory_cds_results <- graph_test(trajectory_cds_filtered,
neighbor_graph="principal_graph", cores=1)
saveRDS(trajectory_cds_results,
paste0("trajectory-cds_results_", Sys.Date(), ".rds"))
cds_sig_genes <- row.names(subset(trajectory_cds_results, q_value < 0.05))
sig_dt <- as.data.table(subset(trajectory_cds_results,
morans_I > 0.01 & q_value < 0.05), keep.rownames=T)
sig_genes_I1Q01 <- row.names(subset(trajectory_cds_results,
morans_I > 0.1 & q_value < 0.01))
sig_genes_I15Q01 <- row.names(subset(trajectory_cds_results,
morans_I > 0.15 & q_value < 0.01))
sig_genes_I25Q01 <- row.names(subset(trajectory_cds_results,
morans_I > 0.25 & q_value < 0.01))
gene_modules_Q05 <- find_gene_modules(trajectory_cds_filtered[cds_sig_genes,],
resolution=c(10^seq(-6,-1)), k=20)
gene_modules_Q05_dt <- as.data.table(gene_modules_Q05)
fwrite(gene_modules_Q05_dt, paste0("figures/", sample_name, "_trajectory_GeneModules-Q05.csv"), sep=",")
gene_modules_I1Q01 <- find_gene_modules(trajectory_cds_filtered[sig_genes_I1Q01,],
resolution=c(10^seq(-6,-1)), k=20)
gene_modules_I1Q01_dt <- as.data.table(gene_modules_I1Q01)
fwrite(gene_modules_I1Q01_dt, paste0("figures/", sample_name, "_trajectory_GeneModules-I01-Q01.csv"), sep=",")
gene_modules_I15Q01 <- find_gene_modules(trajectory_cds_filtered[sig_genes_I15Q01,],
resolution=c(10^seq(-6,-1)), k=20)
gene_modules_I15Q01_dt <- as.data.table(gene_modules_I15Q01)
fwrite(gene_modules_I15Q01_dt, paste0("figures/", sample_name, "_trajectory_GeneModules-I015-Q01.csv"), sep=",")
gene_modules_I25Q01 <- find_gene_modules(trajectory_cds_filtered[sig_genes_I25Q01,],
resolution=c(10^seq(-6,-1)), k=20)
gene_modules_I25Q01_dt <- as.data.table(gene_modules_I25Q01)
fwrite(gene_modules_I25Q01_dt, paste0("figures/", sample_name, "_trajectory_GeneModules-I025-Q01.csv"), sep=",")
cellname_celltypes <- colData(trajectory_cds_filtered)$cell_type_manual
names(cellname_celltypes) <- row.names(colData(trajectory_cds_filtered))
trj_df <- tibble::tibble(cell=row.names(colData(trajectory_cds_filtered)),
cell_group=cellname_celltypes[colnames(trajectory_cds_filtered)])
agg_mattrj <- aggregate_gene_expression(trajectory_cds_filtered, gene_modules_Q05, trj_df)
row.names(agg_mattrj) <- stringr::str_c("Module ", row.names(agg_mattrj))
colnames(agg_mattrj) <- stringr::str_c("Partition ", colnames(agg_mattrj))
pheatmap_obj_ctrl <- pheatmap::pheatmap(agg_mattrj, cluster_rows=TRUE, cluster_cols=TRUE,
scale="column", clustering_method="ward.D2",
fontsize=6)[[4]]
heatmapplt <- grid.arrange(pheatmap_obj_ctrl, nrow=1, ncol=1)
ggsave(paste0("figures/", sample_name, "_trajectory_GeneModules-Q05_Heatmap_", Sys.Date(), ".png"),
plot = heatmapplt, width = 360, height = 180, units = "mm", dpi = 600)
ggsave(paste0("figures/", sample_name, "_trajectory_GeneModules-Q05_Heatmap_", Sys.Date(), ".svg"),
plot = heatmapplt, width = 360, height = 180, units = "mm", dpi = 300)
Run the hallmark pathways for samples split on condition and timepoint using SiPSiC
msigdb_path <- "/project/processed/msigdb/"
GO_hallmarks <- file.path(msigdb_path, "mh.all.v2024.1.Mm.symbols.gmt")
GO_cell_signature <- file.path(msigdb_path, "m8.all.v2024.1.Mm.symbols.gmt")
GO_regulatory_targets <- file.path(msigdb_path, "m3.all.v2024.1.Mm.symbols.gmt")
GO_canonical_pathways <- file.path(msigdb_path, "m2.cp.v2024.1.Mm.symbols.gmt")
GO_curated_gene_sets <- file.path(msigdb_path, "m2.all.v2024.1.Mm.symbols.gmt")
GO_biological_process <- file.path(msigdb_path, "m5.go.bp.v2024.1.Mm.symbols.gmt")
GO_cellular_component <- file.path(msigdb_path, "m5.go.cc.v2024.1.Mm.symbols.gmt")
GO_molecular_function <- file.path(msigdb_path, "m5.go.mf.v2024.1.Mm.symbols.gmt")
GO_list <- c(GO_hallmarks, GO_cell_signature, GO_regulatory_targets,
GO_canonical_pathways, GO_curated_gene_sets, GO_biological_process,
GO_cellular_component, GO_molecular_function)
names(GO_list) <- c("hallmarks", "cell-signature", "regulatory-targets",
"canonical-pathways", "curated-gene-sets", "BP", "CC", "MF")
Plot a heatmap of the significantly different enrichments of hallmark pathways
- Related to Figure 2
# Constants Definition
minimalClusterSize <- 10
minNumOfGenesExpressed <- 1000
cpm_obj <- NormalizeData(trj_object, method = "RC", scale.factor = 1e6)
data_mat <- as.matrix(SeuratObject::LayerData(cpm_obj, assay = "RNA", layer = "counts"))
data_mat_unfiltered <- as.matrix(SeuratObject::LayerData(cpm_obj, assay = "RNA", layer = "counts"))
# Filtering out cells which express less than the minimal number of genes
expressedGenesCounters <- apply(data_mat != 0, 2, sum)
cellsWithAThousandPlus <- expressedGenesCounters >= minNumOfGenesExpressed
data_mat <- data_mat[, cellsWithAThousandPlus]
expressedGenesCounters <- expressedGenesCounters[cellsWithAThousandPlus]
# Filtering out genes which are expressed by less than the minimal expected cluster size of cells
nonZeroCellCountsForGenes <- apply(data_mat != 0, 1, sum)
totalCellsCount <- ncol(data_mat)
minNumOfCellsInClust <- totalCellsCount * (minimalClusterSize / 100)
genesWithMinExpression <- (nonZeroCellCountsForGenes > minNumOfCellsInClust)
data_mat <- data_mat[genesWithMinExpression,]
# Check order of data_mat
data_colnames <- data.table(col_names=colnames(data_mat))
data_colnames[, c("condition", "timepoint", "barcode") := tstrsplit(col_names, "_", fixed=TRUE)]
data_colnames
table(data_colnames$condition, data_colnames$timepoint)
ctrl_1mo_names <- colnames(data_mat)[grep("ctrl_1mo", colnames(data_mat))]
ctrl_3mo_names <- colnames(data_mat)[grep("ctrl_3mo", colnames(data_mat))]
ctrl_6mo_names <- colnames(data_mat)[grep("ctrl_6mo", colnames(data_mat))]
trt_1mo_names <- colnames(data_mat)[grep("trt_1mo", colnames(data_mat))]
trt_3mo_names <- colnames(data_mat)[grep("trt_3mo", colnames(data_mat))]
trt_6mo_names <- colnames(data_mat)[grep("trt_6mo", colnames(data_mat))]
for (gene_set in "hallmarks") {
genesets <- getGmt(GO_list[[gene_set]])
allPathwayScores <- numeric()
all_effect_sizes_ctrl_1mo_vs_trt_1mo <- numeric()
all_effect_sizes_ctrl_1mo_vs_trt_3mo <- numeric()
all_effect_sizes_ctrl_1mo_vs_trt_6mo <- numeric()
all_effect_sizes_ctrl_1mo_vs_ctrl_3mo <- numeric()
all_effect_sizes_ctrl_1mo_vs_ctrl_6mo <- numeric()
all_effect_sizes_ctrl_3mo_vs_ctrl_6mo <- numeric()
all_effect_sizes_ctrl_3mo_vs_trt_1mo <- numeric()
all_effect_sizes_ctrl_3mo_vs_trt_3mo <- numeric()
all_effect_sizes_ctrl_3mo_vs_trt_6mo <- numeric()
all_effect_sizes_ctrl_6mo_vs_trt_1mo <- numeric()
all_effect_sizes_ctrl_6mo_vs_trt_3mo <- numeric()
all_effect_sizes_ctrl_6mo_vs_trt_6mo <- numeric()
all_effect_sizes_trt_1mo_vs_trt_3mo <- numeric()
all_effect_sizes_trt_1mo_vs_trt_6mo <- numeric()
all_effect_sizes_trt_3mo_vs_trt_6mo <- numeric()
pvals_ctrl_1mo_vs_ctrl_3mo <- numeric()
pvals_ctrl_1mo_vs_ctrl_6mo <- numeric()
pvals_ctrl_1mo_vs_trt_1mo <- numeric()
pvals_ctrl_1mo_vs_trt_3mo <- numeric()
pvals_ctrl_1mo_vs_trt_6mo <- numeric()
pvals_ctrl_3mo_vs_ctrl_6mo <- numeric()
pvals_ctrl_3mo_vs_trt_1mo <- numeric()
pvals_ctrl_3mo_vs_trt_3mo <- numeric()
pvals_ctrl_3mo_vs_trt_6mo <- numeric()
pvals_ctrl_6mo_vs_trt_1mo <- numeric()
pvals_ctrl_6mo_vs_trt_3mo <- numeric()
pvals_ctrl_6mo_vs_trt_6mo <- numeric()
pvals_trt_1mo_vs_trt_3mo <- numeric()
pvals_trt_1mo_vs_trt_6mo <- numeric()
pvals_trt_3mo_vs_trt_6mo <- numeric()
lapply(genesets@.Data, executePathwayCalculations, data_mat,
ctrl_1mo_names,
ctrl_3mo_names,
ctrl_6mo_names,
trt_1mo_names,
trt_3mo_names,
trt_6mo_names,
padj_method = "BH",
output_dir = paste0("figures/SiPSiC/", gene_set, "/"), debug=F)
saveRDS(allPathwayScores, file = paste0("figures/SiPSiC/", gene_set,
"/allPathwayScores_", gene_set,
".rds"))
}
scores_dt <- as.data.table(allPathwayScores, keep.rownames=T)
t_scores <- transpose(scores_dt, keep.names = "rn", make.names="rn")
t_scores$cell_name <- colnames(scores_dt)[-1]
t_scores[, c("condition", "timepoint", "barcode") := tstrsplit(cell_name, "_", fixed=TRUE)]
melted <- as.data.table(reshape2::melt(t_scores))
melted$condition_x_timepoint <- paste0(melted$condition, "_", melted$timepoint)
heat_mat <- as.data.frame.matrix(xtabs(value~variable + condition_x_timepoint, melted))
num_ctrl <- length(grep("ctrl", colnames(heat_mat)))
num_trt <- length(grep("trt", colnames(heat_mat)))
z_score <- t(scale(t(heat_mat)))
sig_pathways_1mo <- names(pvals_ctrl_1mo_vs_trt_1mo[pvals_ctrl_1mo_vs_trt_1mo < 0.01])
sig_pathways_3mo <- names(pvals_ctrl_3mo_vs_trt_3mo[pvals_ctrl_3mo_vs_trt_3mo < 0.01])
sig_pathways_6mo <- names(pvals_ctrl_6mo_vs_trt_6mo[pvals_ctrl_6mo_vs_trt_6mo < 0.01])
sig_pathways <- intersect(sig_pathways_1mo, intersect(sig_pathways_3mo, sig_pathways_6mo))
sig_pathways <- sort(c(sig_pathways, "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION"))
pheatmap_obj <- pheatmap::pheatmap(heat_mat[rownames(heat_mat) %in% sig_pathways, ],
cluster_rows=FALSE,
cluster_cols=FALSE,
scale="row",
clustering_method="ward.D2",
cellwidth = 5, cellheight = 5,
breaks = seq(-2,2, 0.04),
color = colorRampPalette(c("navy", "white", "red"))(100),
fontsize=6)[[4]]
heatmap_plt <- grid.arrange(pheatmap_obj, nrow=1, ncol=1)
ggsave(paste0("figures/SiPSiC/",
"captopril-multiome_significant-hallmark-pathways_heatmap_",
Sys.Date(), ".png"),
plot = heatmap_plt, width = 90, height = 90, units = "mm", dpi = 300)
ggsave(paste0("figures/SiPSiC/",
"captopril-multiome_significant-hallmark-pathways_heatmap_",
Sys.Date(), ".svg"),
plot = heatmap_plt, width = 90, height = 90, units = "mm", dpi = 300)
Heatmaps of treated vs untreated of each trajectory population's top differentially expressed genes
- Related to Figure 2.
trj_object$celltype_x_condition <- paste0(trj_object$cell_type_manual, "_", trj_object$condition)
Idents(trj_object) <- trj_object$celltype_x_condition
table(Idents(trj_object))
features <- rownames(trj_object)
no_ribo <- features[!grepl("^Rp[sl][[:digit:]]|^Rplp[[:digit:]]|^Rpsa", features)]
no_ribo_mito <- sort(no_ribo[!grepl("^mt-(.*)", no_ribo)])
no_ribo_mito_hemo <- sort(no_ribo_mito[!grepl("^H2-(.*)", no_ribo_mito)])
filtered_features <- sort(no_ribo_mito_hemo[!grepl("^Hb[ab]-(.*)", no_ribo_mito_hemo)])
for (pop in c("SMC", "tSMC_1", "tSMC_2", "MC", "JG", "EnC", "fEnC", "arterialEnC", "glomEnC")) {
diff_MAST <- as.data.table(FindMarkers(trj_object,
ident.1 = paste0(pop, "_trt"),
ident.2 = paste0(pop, "_ctrl"),
latent.vars = c("timepoint"),
features = filtered_features, test="MAST"), keep.rownames=T)
diff_MAST_sig <- diff_MAST[p_val_adj < 0.01,]
if (nrow(diff_MAST_sig) <= 20) {
diff_MAST_sig <- diff_MAST[p_val_adj < 0.10,]
}
trt_UP <- diff_MAST_sig[avg_log2FC > 0,][1:10,]
trt_UP <- trt_UP[complete.cases(trt_UP),]
ctrl_UP <- diff_MAST_sig[order(avg_log2FC),][avg_log2FC < 0,][1:10,]
ctrl_UP <- ctrl_UP[complete.cases(ctrl_UP),]
heat_mat <- as.matrix(xtabs(avg_log2FC~rn, rbind(trt_UP, ctrl_UP)))
break_max <- max(floor(heat_mat[,1]))
break_min <- min(floor(heat_mat[,1]))
break_step <- (break_max - break_min)/100
pheatmap_obj <- pheatmap::pheatmap(heat_mat,
cluster_rows=TRUE,
cluster_cols=FALSE,
scale="none",
clustering_method="ward.D2",
cellwidth = 5, cellheight = 5,
treeheight_row = 5,
breaks = seq(break_min, break_max, break_step),
color = colorRampPalette(c("navy", "white", "red"))(100),
fontsize=6)[[4]]
heatmap_plt <- grid.arrange(pheatmap_obj, nrow=1, ncol=1)
ggsave(paste0("figures/DEGs/trajectory-integrated_", pop, "-heatmap_Top20genes-clustered_", Sys.Date(), ".png"),
plot = heatmap_plt, units = "mm", dpi = 300, create.dir = TRUE)
ggsave(paste0("figures/DEGs/trajectory-integrated_", pop, "-heatmap_Top20genes-clustered_", Sys.Date(), ".svg"),
plot = heatmap_plt, units = "mm", dpi = 300, create.dir = TRUE)
}