library(Seurat) library(openxlsx) library(readxl) library(ggplot2) library(dplyr) library(fgsea) library(pheatmap) library(scales) library(RColorBrewer) library(cetcolor) #username <- "/Users/Squadrito/" #username <- "/Users/bresesti.chiara/" #wdir <- paste0(username, "Dropbox (HSR Global)/90-935462466_scRNAseq_Bresesti/Analysis_MM/results2/CB1_CB3_CB4") #wdir_CB <- paste0(username, "Dropbox (HSR Global)/90-935462466_scRNAseq_Bresesti/CB_scripts/results2") #plot_dir <- paste0(username,"/Dropbox (HSR Global)/CancerGeneTherapy/Cancer Gene Therapy/MS/2024 Bresesti et al/Scripts/plots and tables used in figures") # fdir previously #dirCB <-paste0(username,"Dropbox (HSR Global)/90-935462466_scRNAseq_Bresesti/CB_scripts") #sig <- readRDS(paste0(wdir, "/miDB_sig5.MLS.rds")) #GO terms db wdir <- "/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_Bresesti/Analysis_MM/results2/CB1_CB3_CB4" wdir_CB <- "/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_Bresesti/Analysis_MM/CB_scripts/results2" plot_dir <- "/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_Bresesti/Analysis_MM/CB_scripts_final/plots and tables used in figures" dirCB <-"/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_Bresesti/Analysis_MM/CB_scripts" dir.create(plot_dir, showWarnings=F, recursive=T) # 16357 GO terms db. used for GSEA sig <- readRDS("/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_Bresesti/reference/miDB_sig5.MLS.rds") ################################################################################ # R script for single-cell RNA-seq data analysis and figure generation. # # This script performs cell annotation, visualization, differential gene expression # analysis, and gene set enrichment analysis (GSEA) on single-cell RNA-seq data. # # Figures generated/data produced: # - Figure 5A: UMAP visualization of cell clusters with final annotation # - Figure 5B: UMAP plot with overlayed density of mOrange+ cells # - Figure 5C: Dotplot of Slc7a11 gene expression across cell types and groups # - Figure 5D: Barplot of GSEA results for selected gene sets # - Figure 5E: Heatmap of cytokine signature GSEA results # - Supplementary Figure 5A: Dotplot of marker gene expression across cell types # - Supplementary Figure 5B: CSV tables for cell distribution per cluster and sample # - Supplementary Figure 5C: CSV tables for mOrange+ cell distribution per cluster and sample # # Input data files: # - RDS object: "CB1_CB3_CB4_final.rds" (Seurat object containing scRNA-seq data) # - RDS object: "miDB_sig5.MLS.rds" (GO terms database for GSEA) # ################################################################################ set.seed(42) obs1 <- readRDS(paste0(wdir, "/CB1_CB3_CB4_final.rds")) # Annotate each cell with a sample label based on its original identity (library) and hash ID obs1@meta.data$Sample <- "Undefined" obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB1"&obs1@meta.data$hash.ID=="Mouse1"] <- "Ctrl_01" obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB1"&obs1@meta.data$hash.ID=="Mouse2"] <- "Ctrl_02" obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB1"&obs1@meta.data$hash.ID=="Mouse3"] <- "Ctrl_03" obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB3"&obs1@meta.data$hash.ID=="Mouse1"] <- "miR_01" obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB3"&obs1@meta.data$hash.ID=="Mouse2"] <- "miR_02" obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB3"&obs1@meta.data$hash.ID=="Mouse3"] <- "miR_03" obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB4"&obs1@meta.data$hash.ID=="Mouse1"] <- "miR_04" obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB4"&obs1@meta.data$hash.ID=="Mouse2"] <- "miR_02" obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB4"&obs1@meta.data$hash.ID=="Mouse3"] <- "miR_05" obs1<-FindClusters(obs1, graph.name="RNA_snn_h.orig.ident", resolution=4) # Manually annotate the clusters obs1@meta.data$Annotation.final.CB<-"Undefined" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==0]<- "Monocytes" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==1]<- "TAMs_3" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==2]<- "TAMs_4" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==3]<- "LSECs" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==4]<- "TAMs_1" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==5]<- "Cancer cells" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==6]<- "cDC2_MoDCs" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==7]<- "TAMs_2" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==8]<- "T and NK cells" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==9]<- "Neutrophils" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==10]<- "pre DCs" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==11]<- "KCs" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==12]<- "TAMs_5" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==13]<- "B cells" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==14]<- "Hepatocytes" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==15]<- "Mig cDCs" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==16]<- "Cholangiocytes" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==17]<- "cDC1" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==18]<- "Macrophages" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==19]<- "pDCs" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.1.2==8]<-"KC-LSEC doublet" obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.4==51]<-"Basophils" # Visualize final annotation in Umap (Fig.5A) obs1 <- SetIdent(obs1,value="Annotation.final.CB") pdf(file = paste(plot_dir, "CB_Annotation.final.CB_Fig.5A.pdf", sep = "/"), width=9, height=7) DimPlot(obs1, reduction="umap.harmony.orig.ident", group.by="Annotation.final.CB", label=T, repel=T) dev.off() # Adding a simplified annotation obs1@meta.data$Shortlables<-obs1@meta.data$Annotation.final.CB DCs<-c("Mig cDCs","pDCs","pre DCs", "cDC1", "cDC2_MoDCs") obs1@meta.data$Shortlables[obs1@meta.data$Annotation.final.CB%in%DCs]<- "DCs" granu<-c("Neutrophils","Basophils") obs1@meta.data$Shortlables[obs1@meta.data$Annotation.final.CB%in%granu]<- "Granulocytes" TAMs<-c("TAMs_1", "TAMs_2", "TAMs_3", "TAMs_4", "TAMs_5") obs1@meta.data$Shortlables[obs1@meta.data$Annotation.final.CB%in%TAMs]<- "TAMs" # Dotplot with markers for population (Suppl.Fig.5A) obs1 <- SetIdent(obs1,value="Shortlables") obs1$Shortlables <- factor(x = obs1$Shortlables, levels = c("Cancer cells","Cholangiocytes", "Hepatocytes", "LSECs", "KC-LSEC doublet", "KCs", "Macrophages","TAMs", "Monocytes", "Granulocytes","DCs", "B cells", "T and NK cells")) gene <- c("Vil1", "Cdx2","Cldn4","Cldn7","Epcam", "Krt19", "Sox9", "Cftr", "Hnf1b", "Alb", "Ttr", "Hnf4a", "Cyp3a11", "Lyve1", "Stab2", "Fcgr2b", "Tek", "Vsig4", "Clec4f", "Timd4", "Cd5l", "Mertk", "Timd4", "Adgre1","Cd68","Marco", "Itgam", "Cd14","Arg1", "Ccl2", "Ccr2", "Itgam", "Cx3cr1", "Ly6g", "S100a8", "S100a9","Mcpt8", "Cpa3", "Cd209a","Itgax", "Flt3","H2-Aa","Xcr1", "Clec9a", "Cd19", "Cd22", "Cd79a", "Ms4a1", "Pax5", "Cd3e", "Cd4", "Cd8a", "Ncr1", "Gzmb") pdf(file = paste(plot_dir, "CB_DotPlot_Shortlables_Fig.S5A.pdf", sep="/"), width=12, height=6) DotPlot(obs1,features = unique(gene), cols = "RdBu") + theme(axis.text.x = element_text(angle=45, vjust=1, hjust=1)) dev.off() # Generate table with distribution clusters by Mouse # Table saved in excel and used to generate a barplot with GraphPad Prism (Suppl.Fig.5B) cells_labels <- table(obs1@meta.data$Annotation.final.CB, obs1@meta.data$Sample) cells_labels_p <- round(prop.table(table(obs1@meta.data$Annotation.final.CB, obs1@meta.data$Sample), 2) * 100, 1) write.csv(cells_labels, paste0(plot_dir, "/number_cells_sample_Annotation.final.CB_Fig.S5B.csv"), row.names=T) write.csv(cells_labels_p, paste0(plot_dir, "/number_cells_sample_Annotation.final.CB_percentage_Fig.S5B.csv"), row.names=T) #### Analysis of mOrange+ fraction (Fig.5B, Suppl.Fig.5C,D) #### # The metadata indicating mOrange positivity are: positive_cells (generated with the mapping), mOrange_label (generated with CellRanger), mOrange_union (union of the previous two labels) # Generate table with mOrange+ cells distribution by cluster # Table saved in excel and used to generate a barplot with GraphPad Prism (Suppl.Fig.5C) # Suppl.Fig.5D was subsequently generated dividing data from Suppl.Fig.5C and B obs2 <- subset(obs1, subset = mOrange_union == "mOrange+") cells_labels_orange <- table(obs2@meta.data$Annotation.final.CB, obs2@meta.data$Sample) cells_labels_orange_p <- round(prop.table(table(obs2@meta.data$Annotation.final.CB, obs2@meta.data$Sample), 2) * 100, 1) write.csv(cells_labels_orange, paste0(plot_dir, "/number_cells_sample_Annotation.final.CB_mOrange_fraction_Fig.S5C.csv"), row.names=T) write.csv(cells_labels_orange_p, paste0(plot_dir, "/number_cells_sample_Annotation.final.CB_percentage_mOrange_fraction_Fig.S5C.csv"), row.names=T) # DimPlot with overlayed mOrange density (Fig.5B) data <- Embeddings(object = obs2[["umap.harmony.orig.ident"]])[(colnames(obs2)), c(1,2)] data <- as.data.frame(data) obs1 <- SetIdent(obs1, value="cells") pdf(file = paste(plot_dir, "UMAP_mOrange_cells_Fig.5B.pdf", sep = "/"), width=9, height=7) DimPlot(obs1, reduction = "umap.harmony.orig.ident", pt.size = 0.1, cols = "azure3") + stat_density_2d(aes_string(x = "UMAPh_1", y = "UMAPh_2", fill = "after_stat(level)"), data=data, linewidth = 0.4, geom = "density_2d_filled", colour = "black", alpha = 0.2, n = 150, h = c(2,2))+ scale_fill_gradientn(colours = cet_pal(16, name = "fire")) dev.off() #### GSEA (Fig.5C-E) #### # Subset populations of interest (to plot) selected_pop <- c("DCs", "KCs", "TAMs", "Monocytes", "Cancer cells", "LSECs") obs3 <- subset(obs1, subset = Shortlables %in% selected_pop) # Dotplot of Slc7a11 expression by cluster and experimental group (Fig.5C, groups divided in Illustrator) obs3 <- SetIdent(obs3, value="Shortlables") obs3$Shortlables <- factor(x = obs3$Shortlables, levels = c("LSECs", "Cancer cells", "Monocytes", "TAMs", "KCs", "DCs")) pdf(file = paste(plot_dir, "CB_DotPlot_Shortlables_Slc7a11_Fig.5C.pdf", sep = "/"), width=9, height=5) DotPlot(obs3, features = "Slc7a11", dot.scale=20, cols = "RdBu", split.by = "PROP.Group") + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) + coord_flip() dev.off() # Run DEG between groups in each pop simpl_DEG <- createWorkbook() deg_results <- list() for(i in selected_pop){ obs3 <- SetIdent(obs3, value="Shortlables") # We want to calculate DGE between miR342 vs miRmut in each cell population in Shortlables markers <- FindMarkers(obs3, ident.1 = "miR342", group.by = "PROP.Group", subset.ident = i, min.pct=0.01, logfc.threshold=0, min.cells.feature=3, min.cells.group=3) # added min.pct=0.01, logfc.threshold=0, min.cells.feature=3, min.cells.group=3 markers$diff_pct <- abs(markers$pct.1-markers$pct.2) name <-paste("DEG_simpl", i, sep = "_") deg_results[[name]] <- markers assign(name, markers) addWorksheet(simpl_DEG, name) writeData(simpl_DEG, name, markers, rowNames = T) } saveWorkbook(simpl_DEG, file = paste0(plot_dir, "/DEG_simpl_Shortlables.xlsx"), overwrite = T) # Run GSEA for each population GSEA_simpl <- data.frame(matrix(nrow = 0, ncol = 9)) for(i in selected_pop){ df1 <- deg_results[[paste0("DEG_simpl_", i)]] df1 <- df1[order(df1$avg_log2FC),] ranks<-df1$avg_log2FC names(ranks)<-rownames(df1) test<-fgsea(pathways = sig, stats = ranks, minSize=5, maxSize=1000, nproc=1, nPermSimple=10000) test$cluster <- i colnames(GSEA_simpl) <- colnames(test) GSEA_simpl <- rbind(GSEA_simpl, test) } write.xlsx(GSEA_simpl, file = paste0(plot_dir, "/GSEA_simpl_Shortlables.xlsx")) # Plot barplot for selected gene sets (Fig.5D, names and colors edited on illustrator) gene_sets <- c("HALLMARK_TNFA_SIGNALING_VIA_NFKB", "HALLMARK_INFLAMMATORY_RESPONSE", "GOMF_CHEMOKINE_ACTIVITY", "GOBP_RESPONSE_TO_CYTOKINE", "GOBP_CYTOKINE_MEDIATED_SIGNALING_PATHWAY", "HALLMARK_INTERFERON_GAMMA_RESPONSE", "HALLMARK_INTERFERON_ALPHA_RESPONSE", "GOBP_RESPONSE_TO_TYPE_I_INTERFERON", "GOMF_ANTIOXIDANT_ACTIVITY", "GOBP_CELLULAR_OXIDANT_DETOXIFICATION", "GOMF_GLUTATHIONE_PEROXIDASE_ACTIVITY", "GOMF_GLUTATHIONE_TRANSFERASE_ACTIVITY", "GOBP_GLUTATHIONE_DERIVATIVE_METABOLIC_PROCESS", "GOBP_GLUTATHIONE_METABOLIC_PROCESS", "GOBP_ANTIGEN_PROCESSING_AND_PRESENTATION", "GOMF_ANTIGEN_BINDING", "GOMF_MHC_CLASS_II_PROTEIN_COMPLEX_BINDING", "GOMF_MHC_PROTEIN_COMPLEX_BINDING", "GOCC_MHC_CLASS_II_PROTEIN_COMPLEX", "GOCC_MHC_PROTEIN_COMPLEX", "GOBP_LEUKOCYTE_MIGRATION_INVOLVED_IN_INFLAMMATORY_RESPONSE", "GOBP_CYTOKINESIS", "GOBP_MICROTUBULE_CYTOSKELETON_ORGANIZATION", "GOBP_MICROTUBULE_BASED_PROCESS", "GOBP_POSITIVE_REGULATION_OF_CELL_CYCLE_PROCESS", "GOBP_CELL_DIVISION", "HALLMARK_G2M_CHECKPOINT", "GOBP_MITOTIC_CELL_CYCLE", "GOBP_CELL_CYCLE_PHASE_TRANSITION") plot <- filter(GSEA_simpl, pathway %in% gene_sets) plot$cluster <- factor(plot$cluster, levels = c("LSECs", "Cancer cells", "TAMs", "Monocytes","KCs","DCs")) # Remove the 'HALLMARK_' prefix and other prefixes from the 'pathway' column plot$pathway <- gsub(x=plot$pathway, pattern="HALLMARK_|GOBP_|GOMF_|GOCC_", replacement="") plot$pathway <- gsub(x=plot$pathway, pattern="_", replacement=" ") gene_sets <- gsub(x=gene_sets, pattern="HALLMARK_|GOBP_|GOMF_|GOCC_", replacement="") gene_sets <- gsub(x=gene_sets, pattern="_", replacement=" ") # Reorder the terms plot$pathway <- factor(plot$pathway, levels = gene_sets) write.xlsx(plot, file = paste0(plot_dir, "/GSEA_plotted.xlsx")) # Plot the GSEA results (Fig5D) p <- ggplot(data=plot, aes(x=NES, y=pathway, fill=log10(pval))) + geom_bar(stat="identity", color="black")+ scale_fill_gradientn( colors=c("red","white","blue"), values=rescale(c(-8,log10(0.05),0)), limits=c(-8,0), na.value = "red") + labs(y='Pathway', x='NES')+ theme(strip.text = element_text(face = "bold"), axis.text.x = element_text(colour = "black"), axis.text.y = element_text(colour = "black"))+ geom_vline(xintercept=0)+ facet_wrap(~cluster, ncol = length(unique(plot$cluster))) ggsave(filename = paste0(plot_dir, "/GSEA_Fig5D.pdf"), plot=p, width=12, height=7) # Create signature file with only signature cytokines for the more stringent GSEA in Fig.5E sig.cyto <- sig[15061:16357] # Define the populations we want to use for this GSEA selected_pop2 <- c("B cells","cDC1","cDC2_MoDCs","KCs","Macrophages","Mig cDCs","Monocytes","Neutrophils","pDCs","T and NK cells", "TAMs_1","TAMs_2","TAMs_3","TAMs_4","TAMs_5") # Run DEG between groups in each population simpl_DEG2 <- createWorkbook() deg_results <- list() for(i in selected_pop2){ obs1 <- SetIdent(obs1, value="Annotation.final.CB") # We want to calculate DGE between miR342 vs miRmut in each cell population in Annotation.final.CB markers <- FindMarkers(obs1, ident.1="miR342", group.by="PROP.Group", subset.ident=i, min.pct=0.01, logfc.threshold=0, min.cells.feature=3, min.cells.group=3) # added min.pct=0.01, logfc.threshold=0, min.cells.feature=3, min.cells.group=3 markers$diff_pct <- abs(markers$pct.1-markers$pct.2) name <-paste("DEG_simpl", i, sep = "_") deg_results[[name]] <- markers assign(name, markers) addWorksheet(simpl_DEG2, name) writeData(simpl_DEG2, name, markers, rowNames = T) } saveWorkbook(simpl_DEG2, file = paste0(plot_dir, "/DEG_Annotation.final.CB.xlsx"), overwrite = T) # Repeat GSEA for cytokines signatures only GSEA_cyto <- data.frame(matrix(nrow=0, ncol=9)) for(i in selected_pop2){ df1 <- deg_results[[paste0("DEG_simpl_", i)]] df1<-df1[order(df1$avg_log2FC),] ranks<-df1$avg_log2FC names(ranks)<-rownames(df1) test<-fgsea(pathways = sig.cyto, stats = ranks, minSize=2, maxSize=1000, nproc=1, nPermSimple=30000) test$cluster <- i colnames(GSEA_cyto) <- colnames(test) GSEA_cyto <- rbind(GSEA_cyto,test) } write.xlsx(GSEA_cyto, file = paste0(plot_dir, "/GSEA_cyto.xlsx")) # Subset the GSEA_cyto to match the cytokine signature calculated in the populations with the same populations that we have in our scRNAseq. # This is necessary because each cytokine was tested and the effect in each population was saved as a different signature. df <- GSEA_cyto df_Bcells <- df[grepl("B_cell", df$pathway) & df$cluster=="B cells",] df_cDC1 <- df[grepl("cDC1", df$pathway) & df$cluster=="cDC1",] df_cDC2 <- df[grepl("cDC2", df$pathway) & df$cluster=="cDC2_MoDCs",] df_macrophage <- df[grepl("Macrophage", df$pathway) & grepl("TAMs|KCs|Macrophages", df$cluster),] df_MigDC <- df[grepl("MigDC", df$pathway) & df$cluster=="Mig cDCs",] df_Monocyte <- df[grepl("Monocyte", df$pathway) & df$cluster=="Monocytes",] df_Neutrophil <- df[grepl("Neutrophil", df$pathway) & df$cluster=="Neutrophils",] df_pDC <- df[grepl("pDC", df$pathway) & df$cluster=="pDCs",] df_Tcells <- df[grepl("T_cell_CD8", df$pathway) & df$cluster=="T and NK cells",] # Combine again all the generated subsets df <- rbind(df_Bcells, df_cDC1, df_cDC2, df_macrophage, df_MigDC, df_Monocyte, df_Neutrophil, df_pDC, df_Tcells) # Remove the redundant data frames rm(df_Bcells, df_cDC1, df_cDC2, df_macrophage, df_MigDC, df_Monocyte, df_Neutrophil, df_pDC, df_Tcells) # Generate a column indicating the cytokine analysed. df$cytokine <- gsub(".*_","",df$pathway) df$padj[is.na(df$padj)] <- 1 # Convert the NA into padj=1 # Assing marks for the significance: *=padj<0.05; **=padj<0.01; ***=padj<0.001. df$label <- "" df <- df %>% mutate(label = cut( padj, breaks = c(0, 0.001, 0.01, 0.05, 1), label = c("***", "**", "*", ""))) paletteLength <- 100 myColor <- colorRampPalette(rev(brewer.pal(5, "RdBu")))(paletteLength) # Define breaks to set the zero to the center df$NES[is.na(df$NES)] <- 0 # Convert the NA into NAS=0 MIN <- min(df$NES) MAX <- max(df$NES) myBreaks <- c(seq(MIN, 0, length.out=ceiling(paletteLength/2) + 1), seq(MAX/paletteLength, MAX, length.out=floor(paletteLength/2))) # Plot heatmap cytokine signatures (Fig.5E, removed dendrograms in illustrator) mat.df <- with(df, tapply(NES, list(cluster, cytokine), FUN = mean)) mat.label <- with(df, tapply(label, list(cluster, cytokine), FUN = paste)) pheatmap(mat.df, scale = "none", display_numbers = mat.label, clustering_method = "ward.D2", cutree_rows = 2, cutree_cols = 3, col=myColor, breaks = myBreaks, na_col = "grey", filename = paste0(plot_dir, "/Heatmap_cytokine_sig_Fig.5E_new.pdf"), cellwidth=10, cellheight=10) # width=12, height=7