Heatmap displaying differential regions of combinatorial chromatin states in metastatic melanoma tumors

Christopher Terranova May 2019 based on code from Ayush Raman

Load necessary programs

library(circlize)
library(ComplexHeatmap)
library(pheatmap)
library(dplyr)
library(gplots)
library(ggplot2)
library(genefilter)
library(RColorBrewer)
library(reshape2)
library(vegan)

Identification of differential chromatin state regions

# set working directory to combined matrix state file from chromStatesMatrixforClustering_particularStates.pl output
  setwd("/TCGA_ChromHMM/tiled_segment/")

# set chromatin state for analysis -- can be changed for each chromatin state
  i = 6

# read in combined matrix state file from chromStatesMatrixforClustering_particularStates.pl output
  dat <- read.table(paste("CombinedMatrix-",i,"-10000bps.txt",sep =""),header = FALSE,sep = "\t",quote = "",
                    row.names = 1, na.strings = FALSE, stringsAsFactors = FALSE)
  colnames(dat) <- c("BRAF-M137","BRAF-M233","BRAF-M263","BRAF-M275",
                     "BRAF-M305","BRAF-M306","BRAF-M357","BRAF-M409",
                     "BRAF-M527","BRAF-M721","BRAF-M749","BRAF-M762",
                     "BRAF-M807","NRAS-M028","NRAS-M399","NRAS-M642",
                     "NRAS-M852","WT-M035","WT-M822","WT-M857")
  head(dat)

##                          BRAF-M137 BRAF-M233 BRAF-M263 BRAF-M275 BRAF-M305
## chr1:1000000-1009999             0         1         0         0         0
## chr1:100250000-100259999         0         0         0         0         0
## chr1:100260000-100269999         0         0         0         0         0
## chr1:100270000-100279999         0         0         0         0         0
## chr1:100990000-100999999         0         0         0         0         0
## chr1:1010000-1019999             0         0         0         0         0
##                          BRAF-M306 BRAF-M357 BRAF-M409 BRAF-M527 BRAF-M721
## chr1:1000000-1009999             0         0         0         3         0
## chr1:100250000-100259999         0         0         0         0         0
## chr1:100260000-100269999         0         0         0         0         0
## chr1:100270000-100279999         0         0         0         0         0
## chr1:100990000-100999999         0         0         0         0         0
## chr1:1010000-1019999             0         0         0         0         0
##                          BRAF-M749 BRAF-M762 BRAF-M807 NRAS-M028 NRAS-M399
## chr1:1000000-1009999             0         0         0         8         0
## chr1:100250000-100259999         0         0         0         0         0
## chr1:100260000-100269999         0         0         0         0         0
## chr1:100270000-100279999         0         0         0         0         0
## chr1:100990000-100999999         0         0         0         0         0
## chr1:1010000-1019999             0         0         0         2         0
##                          NRAS-M642 NRAS-M852 WT-M035 WT-M822 WT-M857
## chr1:1000000-1009999             0         0       4       0       0
## chr1:100250000-100259999         0         0       0       2       0
## chr1:100260000-100269999         0         0       0       2       0
## chr1:100270000-100279999         0         0       0       2       0
## chr1:100990000-100999999         0         0       0       2       0
## chr1:1010000-1019999             0         0       1       0       0

# sample type labels for plot
  mutation.type <- factor(c(rep("BRAF",13), rep("NRAS",4), rep("WT", 3)))
  
# filtering low variable regions -- Change be changed for other subtype comparisons 
  dat.state1 <- varFilter(as.matrix(dat), var.cutoff = 0.5)
  print(dim(dat.state1))

## [1] 3922   20

  ## Annotation Data frame
  df1 <- data.frame(Mutation = mutation.type)
  mut.gene.cols <- brewer.pal(7, "Set1")[1:3]
  mut.cols.assigned <- setNames(mut.gene.cols, unique(levels(df1$Mutation)))
  
# Feature Extraction by wilcox test
  idx1 <- which(mutation.type == "BRAF")
  idx2 <- which(mutation.type == "NRAS")
  idx3 <- apply(dat.state1, 1, function(r){pvalue = wilcox.test(r[idx1], r[idx2], exact = FALSE)$p.value}) < 0.05
  dat.state3 <- dat.state1[idx3,c(idx1,idx2)]
  dat.state3 <- dat.state3[!is.na(rownames(dat.state3)),]
  print(dim(dat.state3))

## [1] 703  17

  # Heatmap annotation  
  df1 <- data.frame(Mutation = df1[c(idx1,idx2),])
  annot2 <- HeatmapAnnotation(df = df1, col = list(Mutation = mut.cols.assigned))

Heatmap displaying differential regions based on bivalent polycomb chromatin state between NRAS and BRAF genotypes

print(Heatmap(log2(dat.state3+1), name = "log2(Frequency+1)", show_row_names = FALSE, show_column_names = TRUE,
                row_dend_reorder = TRUE, column_dend_reorder = TRUE, clustering_distance_rows = "pearson",
                clustering_distance_columns = "euclidean", clustering_method_rows = "complete",
                clustering_method_columns = "complete", top_annotation = annot2))

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.