# Plot heatmap of cancer-germline antuigen expression and methylation in CCLE hematologic cell line data (Figure 6D)

library(CePa)

library(biomaRt)

library(ggplot2)

library(cowplot)

library(reshape2)

library(matrixStats)

library(data.table)

library(ggrepel)

library(dplyr)

library(gridExtra)

library(ggpubr)

library(ComplexHeatmap)

library(circlize)

library(RColorBrewer)

# load data

meth <- fread("CCLE_RRBS_TSS_1kb_20180614.txt", data.table = F)

meth_cpg <- fread("CCLE_RRBS_TSS_CpG_clusters_20180614.txt", data.table = F)

rpkm <- read.gct(file = "CCLE_DepMap_18q3_RNAseq_RPKM_20180718.gct") # read in CCLE RNA-seq rpkm table

annot <- fread(file = "DepMap-2018q3-celllines.csv", data.table = F)

annot2 <- fread(file = "CCLE_sample_info_file_2012-10-18_modified.csv", data.table = F, check.names = T) # read in CCLE annotation file (2012 publication) with lineage data added by Olli

# clean rpkm column names

colnames(rpkm) <- gsub("..ACH.*", "", colnames(rpkm))

# subset rpkm to methylation data

samples_rkpm_meth <- colnames(rpkm)[colnames(rpkm) %in% colnames(meth_cpg)]

# subset annotations to data

annot <- annot[annot$CCLE_Name%in%samples_rkpm_meth,]

# select hematopoietic cell lines from annot

annot_hem <- annot[grepl("HAEMATOPOIETIC", annot$CCLE_Name),]

colnames(meth) <- gsub("_name|TSS_|cluster_", "", colnames(meth_cpg))

colnames(meth_cpg) <- gsub("_name|TSS_|cluster_", "", colnames(meth_cpg))

# rpkm data frame with hematopoietic cell lines

rpkm_hem <- rpkm[,as.character(annot_hem$CCLE_Name)]

# combine TSS and CpG methylation data

colnames(meth[,c("id", "gene", as.character(annot_hem$CCLE_Name))])==colnames(meth_cpg[,c("id", "gene", as.character(annot_hem$CCLE_Name))]) # check that column names match before joining

meth_hem <- rbind(meth[,c("id", "gene", as.character(annot_hem$CCLE_Name))], meth_cpg[,c("id", "gene", as.character(annot_hem$CCLE_Name))])

meth_hem[,!colnames(meth_hem)%in%c("id", "gene")] <- sapply(meth_hem[,!colnames(meth_hem)%in%c("id", "gene")], as.numeric) # make numeric

# get gene symbols

ensembl_hs_mart <- useMart(biomart="ensembl", dataset="hsapiens_gene_ensembl")

ensembl_df <- getBM(attributes=c("ensembl_gene_id", "ensembl_gene_id_version",

                                 "hgnc_symbol", "description", "chromosome_name", "start_position"),

                    mart=ensembl_hs_mart)

gene_annot <- ensembl_df[match(gsub("\\..*", "", rownames(rpkm)), ensembl_df$ensembl_gene_id),] # match biomart data using Ensembl gene symbols

rownames(gene_annot) <- rownames(rpkm)

# get cga genes:

genelist <- as.character(read.table("cga.txt")[,1])

gene_annot_cga <- gene_annot[gene_annot$hgnc_symbol%in%gsub("C11orf85", "MAJIN", genelist),]

rpkm_hem_cga <- rpkm_hem

rownames(rpkm_hem_cga) <- gsub("\\..*", "", rownames(rpkm))

rpkm_hem_cga <- rpkm_hem_cga[gene_annot_cga$ensembl_gene_id,]

rownames(rpkm_hem_cga) <- gsub("MAJIN", "C11orf85", gene_annot_cga$hgnc_symbol)

rpkm_hem_cga_log2 <- log2(rpkm_hem_cga + 0.25)

meth_hem_cga <- as.data.frame(meth_hem[meth_hem$gene%in%gene_annot_cga$hgnc_symbol,])

rownames(meth_hem_cga) <- gsub("MAJIN", "C11orf85", meth_hem_cga$id)

meth_hem_cga$id <- NULL

meth_hem_cga$gene <- NULL

hem_cga <- t(rbind(rpkm_hem_cga_log2, meth_hem_cga))

hem_cga <- merge(hem_cga, annot_hem, by.x = "row.names", by.y = "CCLE_Name")

hem_cga <- merge(hem_cga, annot2, by.x = "Row.names", by.y = "CCLE.name")

hem_cga <- hem_cga[grepl("lymphoma|leukaemia|myeloma", hem_cga$Hist.Subtype1),]

colnames(hem_cga)[1] <- "CCLE_Name"

# add shortened annotations

grep("lymphoma|leukaemia|myeloma", hem_cga$Hist.Subtype1, value=T)

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%"plasma_cell_myeloma"]="MM"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%"mantle_cell_lymphoma"]="MCL"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%"diffuse_large_B_cell_lymphoma"]="DLBCL"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%"chronic_lymphocytic_leukaemia-small_lymphocytic_lymphoma"]="CLL"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%"blast_phase_chronic_myeloid_leukaemia"]="CML"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%"anaplastic_large_cell_lymphoma"]="ALCL"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%"adult_T_cell_lymphoma-leukaemia"]="TCL other"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%"acute_myeloid_leukaemia"]="AML"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%"acute_lymphoblastic_T_cell_leukaemia"]="T-ALL"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%"acute_lymphoblastic_B_cell_leukaemia"]="pre-B-ALL"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%"Hodgkin_lymphoma"]="CHL"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%"Burkitt_lymphoma"]="BL"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%c("B_cell_lymphoma_unspecified")]="BCL other"

hem_cga$Hist.Subtype_hem[hem_cga$Hist.Subtype1%in%c("peripheral_T_cell_lymphoma_unspecified")]="TCL other"

# add lineage annotations

hem_cga$Lineage <- ifelse(grepl("ALCL|T-ALL|TCL", hem_cga$Hist.Subtype_hem), "T/NK cell", 

                          ifelse(grepl("AML|CML", hem_cga$Hist.Subtype_hem), "Myeloid", "B cell"))

# add methylation means to data frame for plotting

meth_mean <- function(gene){

  df <- data.frame(gene = rowMeans(data.frame(hem_cga[,colnames(hem_cga)[grepl(paste0(gene, "_"), colnames(hem_cga))]]), na.rm = T))

  return(df)

}

cga_meth_mean <- do.call(cbind, lapply(gsub("MAJIN", "C11orf85", genelist), meth_mean))

colnames(cga_meth_mean) <- paste0(gsub("MAJIN", "C11orf85", genelist), "_meth_mean")

hem_cga_meth_mean <- cbind(hem_cga, cga_meth_mean)

hem_cga_meth_mean <- hem_cga_meth_mean[order(hem_cga_meth_mean$Hist.Subtype_hem),]

# order cancers for plotting

hem_cga_meth_mean$Hist.Subtype_hem <- factor(hem_cga_meth_mean$Hist.Subtype_hem, levels = c("T-ALL", "pre-B-ALL", "AML", "CML", "CLL", "MM", "BL", "CHL", "DLBCL", "MCL", "BCL other", "ALCL", "TCL other"))

# data prepared

## --------------------------------------------------------------------------------------

## Test correlations for all methylation values averaged per gene

# function to correlate gexp to averaged same gene methylation features

corr_mean <- function(gene) {

  cor.result <- cor.test(hem_cga[,gene],

                         as.numeric( # numeric vector for cor.test

                           rowMeans( # average methylation values

                             data.frame( # data frame for rowMeans

                               hem_cga[,colnames(hem_cga)[grepl(paste0(gene, "_"), colnames(hem_cga))]]), na.rm = T)), method = "spearman") # grep columns with underscore to get methylation values only

  p <- cor.result$p.value

  r <- cor.result$estimate

  return(data.frame(gexp = gene, meth = gene, r = r, p = p))

}

# select genes that have methylation data

genes = genelist[genelist%in%gsub("_.*", "", colnames(hem_cga)[grepl("_", colnames(hem_cga))])]

# apply function on gene list

result_mean <- do.call(rbind, lapply(genes, corr_mean))

result_mean$q <- p.adjust(result_mean$p, method = "fdr")

# check results sorted by siginificance

result_mean %>%

  arrange(q)

# write sorted mean methylation correlations

result_mean_sorted <- result_mean %>%

  arrange(r)

# add significance codes to mean gexp to meth correlation for each gene

mean_corr <- result_mean %>%

  mutate(signifCode = ifelse(q<0.0001, "****",

                             ifelse(q<0.001, "***", 

                                    ifelse(q<0.01, "**",

                                           ifelse(q<0.05, "*", "")))))

## -------------------------------------------------------------------------------

# Heatmap 

mat <- data.matrix(rpkm_hem_cga_log2[,hem_cga_meth_mean$CCLE_Name])

mat <- t(apply(mat, 1, scale))

colnames(mat) <- gsub("_.*", "", colnames(rpkm_hem_cga_log2[,hem_cga_meth_mean$CCLE_Name]))

mat <- mat[gene_annot_cga$hgnc_symbol,]

mat 

# colors

cols <- read.table("colors_hemap_immunology.tsv", header = TRUE, sep = "\t", comment.char = " ")

cols <- cols %>% mutate(combinedcolor = ifelse(subtypecolor!="", as.character(subtypecolor), as.character(color)))

cols$sample <- gsub("^BCL", "BCL other", gsub("TCL", "TCL other", cols$sample))

# data frames for heatmap annotations

exprs_num <- colSums(rpkm_hem_cga[,hem_cga_meth_mean$CCLE_Name]>0.5)

meth_mat <- hem_cga_meth_mean[,grep("mean", colnames(hem_cga_meth_mean))]

colnames(meth_mat) <- gsub("_meth_mean", "", colnames(meth_mat))

meth_mat <- t(meth_mat)

colnames(meth_mat) <- gsub("_.*", "", hem_cga_meth_mean$CCLE_Name)

meth_mat <- meth_mat[,match(colnames(mat), colnames(meth_mat))]

meth_mat[meth_mat>0.5] = 1

meth_mat[meth_mat<0.5] = 0

meth_num <- colSums(meth_mat==0, na.rm = T)

disease_order <- data.frame(disease = unique(hem_cga_meth_mean$Hist.Subtype_hem), rank = c(13, 3, 11, 7, 8, 5, 4, 9, 10, 6, 2, 1, 12))

disease_order$disease <- as.character(disease_order$disease)

hem_cga_meth_mean <- merge(hem_cga_meth_mean, disease_order, by.x = "Hist.Subtype_hem", by.y = "disease")

sample_order <- order(hem_cga_meth_mean$rank, exprs_num)

exprs_num <- exprs_num[sample_order]

meth_num <- meth_num[sample_order]

annot_df <- data.frame(subtype = hem_cga_meth_mean$Hist.Subtype_hem[sample_order])

# heatmap annotations

ha1 = HeatmapAnnotation(exprs = anno_barplot(exprs_num, 

                                             bar_width = 0.75, 

                                             border = FALSE, 

                                             axis = TRUE,

                                             axis_param = list(gp = gpar(fontsize = 5, lwd = 0.5)),

                                             gp = gpar(col = NA, fill = "grey50")),

                        meth = anno_barplot(meth_num, 

                                            bar_width = 0.75, 

                                            border = FALSE, 

                                            axis = TRUE,

                                            axis_param = list(gp = gpar(fontsize = 5, lwd = 0.5)),

                                            gp = gpar(col = NA, fill = "#a52a2a")),

                        df = annot_df, 

                        col = list(subtype = structure(as.character(cols$combinedcolor[cols$sample %in% hem_cga_meth_mean$Hist.Subtype_hem]), 

                                                       names = as.character(cols$sample[cols$sample %in% hem_cga_meth_mean$Hist.Subtype_hem]))[as.character(disease_order$disease[order(disease_order$rank)])]

                        ),

                        annotation_legend_param = list(subtype = list(title = "Cancer type", title_gp = gpar(fontsize = 7, fontface = "bold"), 

                                                                      labels_gp = gpar(fontsize = 7), grid_height = unit(3, "mm"), grid_width = unit(3, "mm"))

                        ),

                        show_annotation_name = F,

                        height = unit(2, "cm"))

meth_corr_df <- data.frame(corr = mean_corr$r[match(rownames(mat), mean_corr$gexp)])

colnames(meth_corr_df) <- "Correlation to methylation"

mat <- mat[order(meth_corr_df$`Correlation to methylation`),sample_order]

meth_corr_df <- data.frame(corr = mean_corr$r[match(rownames(mat), mean_corr$gexp)])

colnames(meth_corr_df) <- "meth_corr"

ha2 = HeatmapAnnotation(df = meth_corr_df,

                        which = "row",

                        col = list(meth_corr = colorRamp2(c(-1, -0.5, -0.375, -0.25, -0.125, 0, 0.125, 0.25, 0.375, 0.5, 1), rev(brewer.pal(11, "PuOr")))),

                        annotation_legend_param = list(meth_corr = list(title = "Correlation to methylation", title_gp = gpar(fontsize = 7, fontface = "bold"), 

                                                                        labels_gp = gpar(fontsize = 7)), width = unit(0.1, "cm"), grid_height = unit(2, "mm"), grid_width = unit(2, "mm"), legend_direction = "horizontal", legend_width = unit(2, "cm"), title_position = "topcenter"))

# print heatmap

pdf(file = "Figure6D_CCLE_CGA_heatmap.pdf", height = 4, width = 5)

draw(

       ha2 +

       Heatmap(mat,

               name = "Expression Z-score", 

               col = colorRamp2(seq(-3, 3, length.out = 11), rev(brewer.pal(11, "RdBu"))),

               row_names_side = "right", 

               row_names_gp = gpar(fontsize = 6),

               show_column_names = FALSE,

               cluster_rows = FALSE,

               cluster_columns = FALSE,

               top_annotation = ha1, 

               heatmap_legend_param = list(title_gp = gpar(fontsize = 7, fontface = "bold"),

                                           labels_gp = gpar(fontsize = 7),

                                           legend_direction = "horizontal",

                                           legend_width = unit(2, "cm"), title_position = "topcenter",

                                           grid_height = unit(2, "mm"),

                                           grid_width = unit(2, "mm")

               )

       ),

     auto_adjust = F,

     padding = unit(c(2, 15, 2, 2), "mm"), heatmap_legend_side = "bottom")

decorate_annotation("exprs", {

  grid.text("# CGAs expressed", unit(0, "npc") + unit(1, "mm"), 0.5, 

            default.units = "npc", just = c("left", "bottom"), gp = gpar(fontsize = 7))

})

decorate_annotation("meth", {

  grid.text("# CGAs hypomethylated", unit(0, "npc") + unit(1, "mm"), 0.5, 

            default.units = "npc", just = c("left", "bottom"), gp = gpar(fontsize = 7))

})

dev.off()