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