--- a
+++ b/Fig5A_ligands_heatmap.R
@@ -0,0 +1,144 @@
+
+# Plot Figure 5A immunomodulatory molecule heatmap
+
+library(ComplexHeatmap)
+library(circlize)
+library(dendextend)
+library(readxl)
+library(RColorBrewer)
+library(matrixStats)
+
+# load Hemap gene expression data
+data = get(load("data9544_with_gene_symbols.RData"))
+
+# load Hemap annotations
+annot = get(load("Hemap_immunology_Annotations.Rdata"))
+
+annot = annot[!(annot$Category.specifying.lineage.tumor.origin=="MM"&annot$CELLS_SORTED==0),] # remove non-malignant myeloma samples
+
+data = data[annot[,1],]
+
+# modify annotations
+annot$Category.specifying.subtype.Heatmap <- annot$plotNormals
+
+annot$Category.specifying.subtype.Heatmap[annot$Sample.type%in%c("Cancer", "Prolif")] <- annot$Category.specifying.subtype[annot$Sample.type%in%c("Cancer", "Prolif")]
+
+annot$Category.specifying.subtype.Heatmap[annot$Category.specifying.lineage.tumor.origin=="MDS"] <- "MDS"
+annot$Category.specifying.subtype.Heatmap[annot$Category.specifying.lineage.tumor.origin=="AML"] <- "AML"
+annot$Category.specifying.subtype.Heatmap[annot$Category.specifying.subtype=="AILT"] <- "AITL"
+annot$Category.specifying.subtype.Heatmap[annot$Category.specifying.subtype=="LC"] <- "LCH"
+
+
+# -----------------------
+
+
+# list of immunomodulatory genes
+ligands = read_excel("ligands.xlsx")
+genelist = ligands$Gene
+
+# Select cancer samples and summarize by cancer type medians
+ligands_cancer <- data[annot$Sample.type%in%c("Cancer", "Prolif")&annot$Category.specifying.subtype.Heatmap!="na", colnames(data) %in% genelist]
+ligands_cancer_median <- aggregate(ligands_cancer, by = list(annot$Category.specifying.subtype.Heatmap[annot$Sample.type%in%c("Cancer", "Prolif")&annot$Category.specifying.subtype.Heatmap!="na"]), FUN = "median")
+rownames(ligands_cancer_median) <- ligands_cancer_median$Group.1
+ligands_cancer_median$Group.1 <- NULL
+
+# -----------------------
+
+# select genes from expression matrix
+mat_unscaled = t(ligands_cancer_median)
+mat_unscaled = mat_unscaled[genelist,]
+
+# scaale
+mat = t(apply(mat_unscaled, 1, scale))
+colnames(mat) <- colnames(mat_unscaled)
+
+# remove genes with max log2 expression <5 (not expressed)
+mat = mat[!rowMaxs(mat_unscaled)<5,]
+
+# make another data matrix for receptor names for later use
+mat_receptor <- mat
+
+# add more common names of ligands in parentheses
+ligands <- ligands[match(rownames(mat), ligands$Gene),]
+
+rownames(mat) <- paste(ligands$Gene, paste0("(", ligands$`Common name`, ")"))
+rownames(mat) <- gsub("\\ \\(NA\\)", "", rownames(mat))
+rownames(mat) <- gsub("PD-1H", "", rownames(mat))
+
+# make lineage annotation to heatmap
+lineage <- colnames(mat)
+lineage[grepl("AML|CML|MDS|JMML|Monocyte|Neutrophil|Myeloid|Dendritic|Macrophage|^LC|Langerhans|Erythroid", lineage)] <- "Myeloid"
+lineage[grepl("CLL|HCL|FL|MM|MALT|MCL|MZL|CHL|NLPHL|BL|B-ALL|DLBCL|B cell|Germinal|Plasma", lineage)] <- "B cell"
+lineage[grepl("AITL|ALCL|ATL|HSTCL|ENKTL|CTCL|PTCLNOS|T-ALL|T/NKcell", lineage)] <- "T/NK cell"
+
+lineage <- as.data.frame(lineage)
+
+mainclass <- colnames(mat)
+mainclass[grepl("AML|CML|JMML|T-ALL|CLL|pre-B-ALL|HCL", mainclass)] <- "Leukemia"
+mainclass[grepl("DLBCL|FL|MALT|MCL|MZL|CHL|NLPHL|BL|AI|ALCL|ATL|HSTCL|ENKTL|CTCL|PTCLNOS", mainclass)] <- "Lymphoma"
+mainclass[grepl("MM|LC|MDS", mainclass)] <- "Other"
+mainclass <- as.data.frame(mainclass)
+
+heatmap_annot <- cbind(lineage, mainclass)
+
+# heatmap annotations
+ha1 = HeatmapAnnotation(df = heatmap_annot,
+                        col = list(lineage = structure(brewer.pal(length(unique(heatmap_annot$lineage)), "Set2")[c(1,3,2)],
+                                                       names = as.character(unique(heatmap_annot$lineage))[c(2,1,3)]),
+                                   mainclass = structure(brewer.pal(10, "Set3")[c(10,7,9)],
+                                                         names = as.character(unique(heatmap_annot$mainclass))[c(2,1,3)])
+                        ),
+                        annotation_legend_param = list(lineage = list(title = "Lineage", title_gp = gpar(fontsize = 10, fontface = "plain"),
+                                                                      labels_gp = gpar(fontsize = 10)),
+                                                       mainclass = list(title = "Cancer type", title_gp = gpar(fontsize = 10, fontface = "plain"),
+                                                                        labels_gp = gpar(fontsize = 10))
+                        ))
+
+
+# ligand cell type information annotation
+ha2 = rowAnnotation(df = as.data.frame(ligands[5]),
+                    col = list(`Immune checkpoint function` = structure(brewer.pal(length(unique(ligands$`Immune checkpoint function`)), "Paired")[c(3,2,1,4)],
+                                                                        names = as.character(unique(ligands$`Immune checkpoint function`))[c(3,2,1,4)])),
+                    annotation_legend_param = list(`Immune checkpoint function` = list(title = "Immune checkpoint function", title_gp = gpar(fontsize = 10, fontface = "plain"),
+                                                                                       labels_gp = gpar(fontsize = 10)), width = unit(0.25, "cm")))
+
+
+# receptor names
+zero_col_mat = matrix(nrow = nrow(mat), ncol = 0)
+rownames(zero_col_mat) <- ligands$`Receptor common name`
+rownames(zero_col_mat) <- gsub("\\, TNFSF14", "", rownames(zero_col_mat))
+
+# dendrogram colors
+dend = hclust(as.dist(1-cor(mat, method="spearman")), method = "ward.D")
+dend = color_branches(dend, k = 3, col = c("#66C2A5", "#B3DE69", "#FC8D62"))
+
+
+# print heatmap
+pdf(file = "Figure5A_Hemap_ligands_heatmap.pdf", height = 8.5, width = 10)
+draw(Heatmap(mat,
+        split = ligands$Split,
+        gap = unit(2, "mm"),
+        name = "Expression (log2)\nZ-score",
+        col = colorRamp2(c(seq(-2, 2, length.out = 9), 4), rev(brewer.pal(9, "RdBu"))[c(1:9,9)]),
+        rect_gp = gpar(col = "white", lty = 1, lwd = 1),
+        row_names_side = "right",
+        row_names_gp = gpar(fontsize = 9),
+        column_names_gp = gpar(fontsize = 10),
+        cluster_rows = FALSE,
+        clustering_distance_rows = "spearman",
+        clustering_method_rows = "ward.D",
+        cluster_columns = dend,
+        clustering_distance_columns = "spearman",
+        clustering_method_columns = "ward.D",
+        top_annotation = ha1,
+        heatmap_legend_param = list(title_gp = gpar(fontsize = 10, fontface = "plain"))
+) +
+  Heatmap(zero_col_mat,
+          right_annotation = ha2,
+          row_names_gp = gpar(fontsize = 9)),
+auto_adjust = F
+)
+dev.off()
+
+
+