ErinPoster.Rmd

---
title: "Heatmap: All stages modENCODE ChIP-seq overlapping promoters of cell sorted intestine RNA-seq"
author: "David C. King"
date: "03/27/2023"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache=TRUE)
library(tidyverse)
library(dplyr)
library(GenomicRanges)
library(ComplexHeatmap)
library(RMariaDB)
library(R.cache)
library(magrittr)
library(stringr)
library('randomForest')

sort_hc_by_svd = function(hc, mat, which.vec=1){
    sv = svd(t(mat))$v[,which.vec]
    dend = reorder(as.dendrogram(hc), wts = sv^2)
    sorted.hc = as.hclust(dend)
    sorted.hc$distance.method = hc$distance.method
    sorted.hc$method = hc$method
    sorted.hc
}

# for distances between rows/columns in a binary matrix
cosine.dist = function(df) { 
  M = as.matrix(df)
  sim = M / sqrt(rowSums(M*M))
  sim = sim %*% t(sim)
  as.dist(1-sim)
}
cosine.sim = function(df) { 
  M = as.matrix(df)
  sim = M / sqrt(rowSums(M*M))
  sim = sim %*% t(sim)
  sim
}
source('OnishRDB/R/DBPeakOverlaps.R', echo = T, chdir=T) # sources DB_connect.R
```

```{r intestine present genes}
intestinePresentDB = dbReadTableCached("intestine_present", "NishimuraLab")
cat(nrow(intestinePresentDB), "genes enriched or equal in any stage")
head(intestinePresentDB,10)
```

```{r }
#DBPeakOverlaps %<>% column_to_rownames("name")
cat(nrow(DBPeakOverlaps), "genes with ChIP-seq overlap information")
DBPeakOverlaps[1:10,1:10]
```

```{r merge overlaps with intestine-present}
DBPeakOverlapsIntestinePresent = intestinePresentDB %>% inner_join(DBPeakOverlaps, by = "name")
enrichedmentPattern = DBPeakOverlapsIntestinePresent %>% select(name, WBID, embryo, L1, L3)
DBPeakOverlapsIntestinePresent %<>% column_to_rownames("name") %>% select(-WBID, -embryo, -L1, -L3)


embryo.DBPeakOverlapsIntestinePresent = intestinePresentDB %>% 
  filter( is.na(intestinePresentDB$L1) & 
          is.na(intestinePresentDB$L3)) %>% 
  inner_join(DBPeakOverlaps, by = "name")

embryo.DBPeakOverlapsIntestinePresent %<>% column_to_rownames("name") %>% 
  select(-WBID, -embryo, -L1, -L3)

sum(!is.na(DBPeakOverlapsIntestinePresent))
naOrZero = function(x) { is.na(x) | (x == 0)}
mx = ifelse(naOrZero(DBPeakOverlapsIntestinePresent), 0, 1)
embryo.mx = ifelse(naOrZero(embryo.DBPeakOverlapsIntestinePresent), 0, 1)
```


```{r count zero fraction}

fractionZero = sum(mx == 0) / length(mx)
cat("fraction of zero cells:", fractionZero,"\n")

embryo.fractionZero = sum(embryo.mx == 0) / length(embryo.mx)
cat("fraction of zero embryo cells:", embryo.fractionZero,"\n")
```


```{r make-into-binary}


rs = rowSums(mx)
cs = colSums(mx)
cat("The following genes are bound in over 90% of the ChIP-seq experiments:")
rownames(mx)[rs > .9*ncol(mx)]

rs = rowSums(embryo.mx)
cs = colSums(embryo.mx)
cat("Embryo: The following genes are bound in over 90% of the ChIP-seq experiments:")
rownames(embryo.mx)[rs > .9*ncol(embryo.mx)]
```

```{r euclidean distance}
system.time({eu.mx.row = dist(mx)})
  #  user  system elapsed 
  # 16.157   0.057  16.228 
system.time({eu.mx.col = dist(t(mx))})
  #  user  system elapsed 
  # 2.311   0.010   2.322 
system.time({hc.mx.rows = hclust(eu.mx.row, method = "average")})
  #  user  system elapsed 
  # 0.536   0.037   0.577 
system.time({hc.mx.cols = hclust(eu.mx.col, method = "average")})
  #  user  system elapsed 
  # 0.005   0.001   0.006
```
```{r embryo euclidean distance}
system.time({embryo.eu.mx.row = dist(embryo.mx)})
  #  user  system elapsed 
  # 0.609   0.002   0.611
system.time({embryo.eu.mx.col = dist(t(embryo.mx))})
  #  user  system elapsed 
  # 0.230   0.002   0.231 
system.time({embryo.hc.mx.rows = hclust(embryo.eu.mx.row, method = "average")})
  #  user  system elapsed 
  # 0.021   0.002   0.025 
system.time({embryo.hc.mx.cols = hclust(embryo.eu.mx.col, method = "average")})
  #  user  system elapsed 
  # 0.005   0.001   0.006
```
```{r embryo square heatmap}
embryo.eu.mx.col.mat = as.matrix(embryo.eu.mx.col)
hm=Heatmap(embryo.eu.mx.col.mat, 
           cluster_rows = embryo.hc.mx.cols, 
           cluster_columns = embryo.hc.mx.cols,
           heatmap_legend_param = list(title = "euclidean distance"))
```

```{r save embryo sq heatmap}
pdf(file=format(Sys.time(), "plots/embryo-square-heatmap-%Y-%m-%d-%H-%M.pdf"), 
    width=50, height=50)
hm
dev.off()
```

```{r cosine distance}
system.time({cos.mx.row = cosine.dist(mx)})
  #    user  system elapsed 
  # 7.337   0.214   7.590 
system.time({cos.mx.col = cosine.dist(t(mx))})
  #  user  system elapsed 
  # 0.845   0.010   0.860 

system.time({hc.mx.cosine.rows = hclust(cos.mx.row, method = "average")})
  # user  system elapsed 
  # 0.542   0.041   0.587 
system.time({hc.mx.cosine.cols = hclust(cos.mx.col, method = "average")})
  #  user  system elapsed 
  # 0.005   0.000   0.006 
```
```{r embryo cosine distance}
system.time({embryo.cos.mx.row = cosine.dist(embryo.mx)})
  #    user  system elapsed 
  # 1.099   0.016   1.119
system.time({embryo.cos.mx.col = cosine.dist(t(embryo.mx))})
  #  user  system elapsed 
  # 0.287   0.002   0.290 

system.time({embryo.hc.mx.cosine.rows = hclust(embryo.cos.mx.row, method = "average")})
  # user  system elapsed 
  # 0.019   0.000   0.020 
# system.time({embryo.hc.mx.cosine.cols = hclust(embryo.cos.mx.col, method = "average")})
# Error in hclust(embryo.cos.mx.col, method = "average") : 
#   NA/NaN/Inf in foreign function call (arg 10)
# Timing stopped at: 0.002 0 0.002
```

```{r scree}
plot(1:length(hc.mx.rows$height),rev(hc.mx.rows$height),  xlim=c(0,50), main="scree for euclidean rows")
plot(1:length(hc.mx.cols$height),rev(hc.mx.cols$height),  xlim=c(0,50), main="scree for euclidean columns")

plot(1:length(hc.mx.cosine.rows$height),rev(hc.mx.cosine.rows$height),  xlim=c(0,50), main="scree for cosine rows")
plot(1:length(hc.mx.cosine.cols$height),rev(hc.mx.cosine.cols$height),  xlim=c(0,50), main="scree for cosine columns")
```
```{r embryo scree}
plot(1:length(embryo.hc.mx.rows$height),rev(embryo.hc.mx.rows$height),  xlim=c(0,50), main="embryo scree for euclidean rows")
plot(1:length(embryo.hc.mx.cols$height),rev(embryo.hc.mx.cols$height),  xlim=c(0,50), main="embryo scree for euclidean columns")

plot(1:length(embryo.hc.mx.cosine.rows$height),rev(embryo.hc.mx.cosine.rows$height),  xlim=c(0,50), main="embryo scree for cosine rows")
#plot(1:length(hc.mx.cosine.cols$height),rev(hc.mx.cosine.cols$height),  xlim=c(0,50), main="scree for cosine columns")
```

```{r all cos and euclidean}

hmColors = c('#FFFFCC',  '#68001C')

hm.euclidean = Heatmap(mx, 
        cluster_rows = hc.mx.rows, 
        #cluster_columns = sort_hc_by_svd(hc.mx.cols, t(mx), 2), 
        cluster_columns = hc.mx.cols, 
        heatmap_legend_param = list(title = "euclidean\nTF bound/unbound"),
        show_heatmap_legend = FALSE,
        show_row_names = FALSE, 
        show_column_names = T, 
        show_row_dend = F,
        show_column_dend = T,
        col = hmColors,
        row_split = 4,
        column_split = 5,
        column_names_side = "top",
        use_raster = T)

hm.cosine = Heatmap(mx, 
        cluster_rows = hc.mx.cosine.rows, 
        #cluster_columns = sort_hc_by_svd(hc.mx.cosine.cols,t(mx), 2), 
        cluster_columns = hc.mx.cosine.cols,
        heatmap_legend_param = list(title = "cosine\nTF bound/unbound"),
        show_heatmap_legend = FALSE,
        show_row_names = FALSE, 
        show_column_names = T, 
        show_row_dend = F,
        show_column_dend = T,
        col = hmColors,
        column_split = 5,
        column_names_side = "top",
        use_raster = T)


system.time(print(hm.euclidean))
 #   user  system elapsed 
 # 11.654   2.262  13.958 
system.time(print(hm.cosine)) 
 #   user  system elapsed 
 # 12.657   2.282  14.993 


pdf(file=format(Sys.time(), "plots/all-heatmap-euc-%Y-%m-%d-%H-%M.pdf"), 
    width=60, height=20)
hm.euclidean
dev.off()

pdf(file=format(Sys.time(), "plots/all-heatmap-cos-%Y-%m-%d-%H-%M.pdf"), 
    width=60, height=20)
hm.cosine
dev.off()
```

```{r embryo Heatmaps}
hmColors = c('#FFFFCC',  '#68001C')

embryo.hm.euclidean = Heatmap(embryo.mx, 
        cluster_rows = embryo.hc.mx.rows, 
        cluster_columns = embryo.hc.mx.cols, 
        heatmap_legend_param = list(title = "euclidean\nTF bound/unbound"),
        show_heatmap_legend = FALSE,
        show_row_names = FALSE, 
        show_column_names = T, 
        show_row_dend = F,
        show_column_dend = T,
        col = hmColors,
        row_split = 4,
        column_split = 5,
        column_names_side = "top",
        use_raster = T)
pdf(file=format(Sys.time(), "plots/embryo-heatmap-euc-%Y-%m-%d.pdf"), 
    width=60, height=20)
hm.euclidean
dev.off()

```

```{r embryo heatmap-ELT-2}
embryo.ELT2bound = embryo.mx[,'elt-2_LE_1']

embryo.mx.noELT2 = embryo.mx[,-which(colnames(embryo.mx) == 'elt-2_LE_1')]

dex = dist(t(embryo.mx.noELT2))
hc_cols = hclust(dex, method = "average")


num_binding = colSums(embryo.mx)[-which(colnames(embryo.mx) == 'elt-2_LE_1')]
q_num_binding = quantile(num_binding)
ecdf_num_binding = ecdf(num_binding)(num_binding)

#hc_cols_svd = sort_hc_by_svd(hc_cols, t(embryo.mx.noELT2), 2)
#min_binding_threshold = q_num_binding['50%']
min_binding_threshold = 0


system.time({
hm.embryo.elt2.binding = Heatmap(embryo.mx.noELT2[,num_binding>min_binding_threshold], 
        clustering_method_rows = "average",
        cluster_columns = sort_hc_by_svd(hclust(dist(t(embryo.mx.noELT2[,num_binding>min_binding_threshold])), 
                                                method = "average"), 
                                         t(embryo.mx.noELT2[,num_binding>min_binding_threshold]), 2),
        heatmap_legend_param = list(title = "euclidean\nTF bound/unbound"),
        show_heatmap_legend = FALSE,
        show_row_names = FALSE, 
        show_column_names = FALSE, 
        show_row_dend = F,
        show_column_dend = T,
        col = hmColors,
        row_split = ifelse(embryo.ELT2bound, "ELT-2 bound", "ELT-2 never bound"),
        column_split = 7,
        column_names_side = "top",
        use_raster = T,
        
        ### comment out this annotation to hide num_binding highlights (but restore show_column_names)
        # top_annotation = HeatmapAnnotation(
        #   labs = anno_text(colnames(embryo.mx.noELT2)[num_binding>min_binding_threshold],
        #                    gp = gpar(col=ifelse(num_binding>q_num_binding['50%'], 
        #                                         'black', 
        #                                         'grey'))),
        #   numGenesBound = anno_barplot(num_binding[num_binding>min_binding_threshold]))
        ###
        )
})

pdf(file=format(Sys.time(), "plots/heatmap-embryo-elt2-binding-%Y-%m-%d-%H-%M.pdf"), 
    width=60, height=20, family = "ArialMT")
hm.embryo.elt2.binding
dev.off()
```

```{r embryo heatmap subset}
embryo.ELT2bound = embryo.mx[,'elt-2_LE_1']
embryo.mx.noELT2 = embryo.mx[,-which(colnames(embryo.mx) == 'elt-2_LE_1')]

### deal with some long column names
cc = colnames(embryo.mx.noELT2)
cc[cc == 'W03F9.2_L4+yAd_1'] = 'W03F9.2_L4_1'
names_to_plot = sub('_[0-9]+$', '',cc) # remove last _[0-9]. 
                                       # The names are not unique
                                       # so must be plotted separately (can't be non-unique colnames)



# get the spread of gene binding totals
num_binding = colSums(embryo.mx)[-which(colnames(embryo.mx) == 'elt-2_LE_1')]
q_num_binding = quantile(num_binding)
ecdf_num_binding = ecdf(num_binding)(num_binding)
# apply threshold
min_binding_threshold = q_num_binding['75%']
reduced.mx = embryo.mx.noELT2[,num_binding>min_binding_threshold]
# compute dist and clustering
dex = dist(t(reduced.mx))

hc_cols = sort_hc_by_svd(hclust(dex, method = "complete"), 
                         t(reduced.mx), 2)
```




```{r}

lvls=c("ELT-2 embryo bound", "ELT-2 not bound in embryo")
elt2_row_order = factor(ifelse(embryo.ELT2bound, lvls[1],lvls[2]),
                        levels=rev(lvls))

# plot primary reduced heatmap
hm.embryo.elt2.binding.subset = 
  Heatmap(name = "transcription factor binding profile", reduced.mx, 
        clustering_method_rows = "average",
        cluster_columns = hc_cols,
        heatmap_legend_param = list(title = "TF bound/unbound"),
        show_heatmap_legend = T,
        show_row_names = FALSE, 
        show_column_names = F, 
        show_row_dend = F,
        show_column_dend = T,
        col = hmColors,
        row_split = elt2_row_order,
        row_title="Intestine present genes",
        column_title="TFs binding >164 intestine present genes",
        column_split = 7,
        column_names_side = "top",
        height=unit(7, "inches"),
        top_annotation = HeatmapAnnotation(
          labs = anno_text(names_to_plot[num_binding>min_binding_threshold],
                           just=0, location=0, gp = gpar(fontfamily="Courier"))
          
          ),
        use_raster = T
        )
        # print(
        #   draw(hm.embryo.elt2.binding.subset, 
        #        column_title="TFs binding >164 intestine present genes",
        #        column_title_gp = gpar(fontsize = 18, fontfamily="ArialMT"))
        #   )
pdf(file=format(Sys.time(), "plots/subset-embryo-elt2-binding-%Y-%m-%d-%H-%M.pdf"), 
    width=20, height=20, family = "ArialMT")

hm.embryo.elt2.binding.subset
dev.off()
```


```{r}
# square
dex.asmx = as.matrix(dex)
hm_reduced_square_embryo = Heatmap(column_title = "TF similarity", 
                                   row_title = "TF similarity", 
                                   dex.asmx,
           cluster_rows = hc_cols, 
           cluster_columns = hc_cols,
           column_split = 7,
           row_split = 7,
           show_row_names = FALSE, 
           show_column_names = F,
           height=unit(10, "inches"),
            top_annotation = HeatmapAnnotation(
          labs = anno_text(names_to_plot[num_binding>min_binding_threshold],
                           just=0, location=0, gp = gpar(fontfamily="Courier"))
          
          ),
           heatmap_legend_param = list(title = "euclidean distance"))
pdf(file=format(Sys.time(), "plots/subset-embryo-square-%Y-%m-%d-%H-%M.pdf"), 
    width=20, height=25, family = "ArialMT")
draw(
hm.embryo.elt2.binding.subset %v% hm_reduced_square_embryo,
column_title="modENCODE TF binding of intestine present genes")
dev.off()
```

```{r}
euclidean.dex = dist(t(reduced.mx))
cos.dex = cosine.dist(t(reduced.mx))
hc_cosine = sort_hc_by_svd(hclust(cos.dex, method = "ward.D2"), 
                         t(reduced.mx), 2)
hc_euclidean = sort_hc_by_svd(hclust(euclidean.dex, method = "ward.D2"), 
                         t(reduced.mx), 2)
# square
cos.dmx = as.matrix(cos.dex)
scaled.euclidean.dex = ((euclidean.dex) - min(euclidean.dex))/diff(range(euclidean.dex))
#euclidean.dmx = as.matrix(scaled.euclidean.dex) # make euclidean.dmx from 0 to 1
euclidean.dmx = as.matrix(euclidean.dex)

# smashed: lower tri = euclidean, upper tri = cos
# smashed.dmx = cos.dmx
# smashed.dmx[lower.tri(smashed.dmx)] = euclidean.dmx[lower.tri(euclidean.dmx)]

hm_reduced_square_embryo = Heatmap(euclidean.dmx,
           cluster_rows = hc_euclidean,
           cluster_columns = hc_euclidean,
           heatmap_legend_param = list(title = "euclidean distance"),
           column_names_side = "top",
           row_names_side = "left",
           column_split = 7,
           row_split = 7)
pdf(file=format(Sys.time(), "plots/subset-embryo-square-%Y-%m-%d-%H-%M.pdf"), 
    width=60, height=60, family = "ArialMT")
hm_reduced_square_embryo
dev.off()
```

```{r ELT-2 binding}
LE = mx[,'elt-2_LE_1']
L1 = mx[,'elt-2_L1_1']
L3 = mx[,'elt-2_L3_1']

elt2_anytime = LE | L1 | L3 
cat("ELT-2 binds", sum(elt2_anytime), "genes in any of the three stages.")
```


```{r heatmap-ELT-2}

hm.elt2.binding = Heatmap(mx, 
        clustering_method_rows = "average",
        clustering_method_columns = "average",
        # clustering_distance_columns =  function(u,v) {
        #   1 - ((u%*%v) / (norm(u,"2")*norm(v,"2")))
        # },
        # clustering_distance_rows =  function(u,v) {
        #   1 - ((u%*%v) / (norm(u,"2")*norm(v,"2")))
        # },
        heatmap_legend_param = list(title = "euclidean\nTF bound/unbound"),
        show_heatmap_legend = FALSE,
        show_row_names = FALSE, 
        show_column_names = T, 
        show_row_dend = F,
        show_column_dend = T,
        col = hmColors,
        row_split = ifelse(elt2_anytime, "ELT-2 bound", "ELT-2 never bound"),
        column_split = 5,
        column_names_side = "top",
        use_raster = T
)
# print(hm.elt2.binding)
pdf(file=format(Sys.time(), "plots/heatmap-elt2-binding-%Y-%m-%d-%H-%M.pdf"), 
    width=60, height=20)
hm.elt2.binding
dev.off()
```
```{r add evans chromatin-states}
source('OnishRDB/R/PromoterChromatinStatesOverlap.R', chdir=T)
```

```{r}
states.EE = PromoterChromatinStatesOverlap.df %>% 
  filter(stage == 'EE' & rank == 1) %>%
  select(geneName, stateDescriptionID) %>%
  rename(stateDescriptionID = "chromatinState.EE")

states.L3 = PromoterChromatinStatesOverlap.df %>% 
  filter(stage == 'L3' & rank == 1) %>%
  select(geneName, stateDescriptionID) %>%
  rename(stateDescriptionID = "chromatinState.L3")

states = full_join(states.EE, states.L3, by="geneName")

rownames(states) <- states$geneName

states_marked = states[ rownames(mx),]
states_marked$ELT2.bound = elt2_anytime
```

```{r EE.states-table}
EE.table = table(states_marked$chromatinState.EE, states_marked$ELT2.bound)
EE.table.df = as.data.frame.matrix(EE.table)
statedescripts = read.table("evans.states.list.txt", sep="\t", header=F,quote="")
rownames(EE.table.df) <- statedescripts$V2
knitr::kable(EE.table.df)
```
```{r L3.states-table}
L3.table = table(states_marked$chromatinState.L3, states_marked$ELT2.bound)
L3.table.df = as.data.frame.matrix(L3.table)
statedescripts = read.table("evans.states.list.txt", sep="\t", header=F,quote="")
rownames(L3.table.df) <- statedescripts$V2
colnames(L3.table.df) <- c("never","at least one stage")
knitr::kable(L3.table.df)
L3.table.mx = as.matrix(L3.table.df)
rownames(L3.table.mx) <- as.character(1:20)
crunk = hclust(dist(L3.table.mx))
plot(crunk, main="Reduce the number of states to show\nfor L3 chromatin states")
# 7 :Txn elongation V: introns/repeats
# 12,15,16: Repeats, intergenic, H3K9me2
# 17,20: Pc/H3K27me3: low/silent and pseudogenes
# 11, 6, 1, 3: border/low/promoter/Txn exon

mixed.table.mx = as.matrix(cbind(EE.table.df,L3.table.df))
rownames(mixed.table.mx) <- as.character(1:20)
crunk2 = hclust(dist(mixed.table.mx))
plot(crunk2)
```

```{r add states}
# 1 promoter
# 2 5-prime proximal gene body
# 3-7 Txn elongation
# 8-10 Enhancer
# 11 border
# 12,15 repeats
# 13 retro transposons, pseudogenes
# 14 mixed, tissue specific
# 17-20 silent/low expression
# 16 intergenic, silent genes, piRNA, repeats

L3_states_labelled = case_match(
    states_marked$chromatinState.L3,
    3:7 ~ 'Txn Elongation',
    c(12,15) ~ 'repeats',
    1 ~ 'promoter',
    2 ~ '5-prime gene body',
    8:10 ~ 'Enhancer',
    16 ~ 'intergenic',
    11 ~ 'Border',
    13 ~ 'retro/pseudogenes',
    14 ~ 'mixed, tissue-specific',
    17:20 ~ 'silent/low expression'
)
# 7 :Txn elongation V: introns/repeats
# 12,15,16: Repeats, intergenic, H3K9me2
# 17,20: Pc/H3K27me3: low/silent and pseudogenes
# 11, 6, 1, 3: border/low/promoter/Txn exon
L3_states_grouped = case_match(
  states_marked$chromatinState.L3,
  7 ~ 'Txn elongation V: introns/repeats',
  c(12,15,16) ~ 'Repeats, intergenic, H3K9me2',
  c(17,20) ~ 'Pc/H3K27me3: low/silent and pseudogenes',
  c(11, 6, 1, 3,5) ~ 'border/low/promoter/Txn exon',
  .default = 'other'
)

grouped_names = c('Txn elongation V: introns/repeats',
                  'Repeats, intergenic, H3K9me2',
                  'Pc/H3K27me3: low/silent and pseudogenes',
                  'border/low/promoter/Txn exon',
                  'other')
L3_states_colors = RColorBrewer::brewer.pal(5,"Dark2")
names(L3_states_colors) <- grouped_names

hm.states = Heatmap(mx, 
        #cluster_rows = hc.mx.rows, 
        cluster_columns = hc.mx.cols, 
        heatmap_legend_param = list(title = "euclidean\nTF bound/unbound"),
        show_heatmap_legend = FALSE,
        show_row_names = FALSE, 
        right_annotation = rowAnnotation(
          statesL3 = L3_states_grouped,
          col = list(statesL3=L3_states_colors)
        ),
        show_column_names = T, 
        show_row_dend = F,
        show_column_dend = T,
        col = hmColors,
        row_split = ifelse(elt2_anytime, "ELT-2 bound", "ELT-2 never bound"),
        column_split = 5,
        column_names_side = "top",
        use_raster = T)

pdf(file=format(Sys.time(), "plots/heatmap-states-labeled-%Y-%m-%d-%H-%M.pdf"), 
    width=60, height=20)
hm.states
dev.off()
```

```{r states-grouped-numerically}

L3_states_renumbering = c(
  1.0,  # Promoter
  1.05, # 5' proximal, gene body
  1.1,  # Txn Elongation I: exon
  1.15, # Txn elongation II: exon/intron
  1.2,  # Txn elongation III: exon and gene end
  1.25, # Txn elongation IV: low expression/repeats
  1.3,  # Txn elongation V: introns/repeats
  2.0,  # Enhancer I: intronic
  2.05, # Enhancer II: intergenic
  2.1,  # Enhancer III: weak 
  3,    # border
  4,    # Repeats, low expr introns, AT rich
  5,    # Retrotransposons, pseudogenes, H3K9me3, H3K27me3
  6,    # mixed,tissue-specific
  7,    # Repeats, RNA pseudogenes, H3K9me2
  8,    # Intergenic, silent genes, piRNAs, repeats
  9,    # Pc/H3K27me3 I: low expr/silent and pseudogenes
  9.05, # Pc/H3K27me3 II:Low expr/silent
  9.1,  # Pc/H3K27me3 III:low expr/silent gene body
  9.15  # H3K9me3 & H3K27me3: silent genes and pseudogenes
)


L3_states_quantified = case_match(
    states_marked$chromatinState.L3,
  1 ~ 1.0,  # Promoter
  2 ~ 1.05, # 5' proximal, gene body
  3 ~ 1.1,  # Txn Elongation I: exon
  4 ~ 1.15, # Txn elongation II: exon/intron
  5 ~ 1.2,  # Txn elongation III: exon and gene end
  6 ~ 1.25, # Txn elongation IV: low expression/repeats
  7 ~ 1.3,  # Txn elongation V: introns/repeats
  8 ~ 2.0,  # Enhancer I: intronic
  9 ~ 2.05, # Enhancer II: intergenic
  10 ~ 2.1,  # Enhancer III: weak 
  11 ~ 3,    # border
  12 ~ 4,    # Repeats, low expr introns, AT rich
  14 ~ 5,    # mixed,tissue-specific
  15 ~ 6,    # Repeats, RNA pseudogenes, H3K9me2
  16 ~ 7,    # Intergenic, silent genes, piRNAs, repeats
  17 ~ 8,    # Pc/H3K27me3 I: low expr/silent and pseudogenes
  18 ~ 8.05, # Pc/H3K27me3 II:Low expr/silent
  19 ~ 8.1,  # Pc/H3K27me3 III:low expr/silent gene body
  20 ~ 8.15,  # H3K9me3 & H3K27me3: silent genes and pseudogenes
  13 ~ 8.2    # Retrotransposons, pseudogenes, H3K9me3, H3K27me3
)

hc_L3_states_quantified = hclust(dist(L3_states_quantified))

L3_states_colors = RColorBrewer::brewer.pal(8,"Dark2")
names(L3_states_colors) <- c('promoter/Txn elongation (7 states)',
                             'Enhancer (3 states)',
                             'border',
                             'Repeats, low expr introns, AT rich',
                             'mixed,tissue-specific',
                             'Repeats, RNA pseudogenes, H3K9me2',
                             'Intergenic, silent genes, piRNAs, repeats',
                             'Pc/H3K27me3, H3K9me3 (5 states)'
                  )

L3_states_legend = case_match(
    states_marked$chromatinState.L3,
  1:7 ~ 'promoter/Txn elongation (7 states)',
  8:10 ~ 'Enhancer (3 states)',
  3 ~  'border',
  4 ~ 'Repeats, low expr introns, AT rich',
  5 ~  'mixed,tissue-specific',
  6 ~ 'Repeats, RNA pseudogenes, H3K9me2',
  7 ~ 'Intergenic, silent genes, piRNAs, repeats',
  11 ~ 'border',
  12 ~ 'Repeats, low expr introns, AT rich',
  14 ~ 'mixed,tissue-specific',
  15 ~ 'Repeats, RNA pseudogenes, H3K9me2',
  16 ~ 'Intergenic, silent genes, piRNAs, repeats',
  c(17:20, 13) ~ 'Pc/H3K27me3, H3K9me3 (5 states)'
 
)

hm.states.clustered = Heatmap(mx, 
        row_order = hc_L3_states_quantified$order, 
        cluster_columns = hc.mx.cols, 
        heatmap_legend_param = list(title = "euclidean\nTF bound/unbound"),
        show_heatmap_legend = FALSE,
        show_row_names = FALSE, 
        right_annotation = rowAnnotation(
          statesL3 = L3_states_legend,
          col = list(statesL3=L3_states_colors)
        ),
        show_column_names = T, 
        show_row_dend = F,
        show_column_dend = T,
        col = hmColors,
        row_split = ifelse(elt2_anytime, "ELT-2 bound", "ELT-2 never bound"),
        column_split = 5,
        column_names_side = "top",
        use_raster = T)

pdf(file=format(Sys.time(), "plots/heatmap-states-clustered-%Y-%m-%d-%H-%M.pdf"), 
    width=50, height=20)
hm.states.clustered
dev.off()


```


```{r search-for-assoc}

# not significant
test.out = chisq.test(table(states_marked[,3:4]))
table(states_marked[,3:4])  %>% mosaicplot(shade = T, main = "No association between chromatin state in the\n promoters of intestine expressed genes and\n the status of ELT-2 binding", sub="p-value = .37") 


# visualize
states.ca = dudi.coa(t(as.data.frame.matrix(ex)), scannf=FALSE, nf=2)

	
get_ca_row(states.ca)

#fviz_ca_biplot(states.ca,repel=TRUE,col.col="brown", col.row="purple")
```

```{r ML}
dmx = data.frame(mx)
dmx$elt2.status = ifelse(elt2_anytime,1, 0)
dmx$elt.2_LE_1 = NULL
dmx$elt.2_L1_1 = NULL
dmx$elt.2_L3_1 = NULL
# Stratified sampling with the rsample package
set.seed(123)
split <- initial_split(dmx, prop = 0.7)
dmx_train  <- training(split)
dmx_test   <- testing(split)


set.seed(123)
cv_model_pcr <- train(
  elt2.status ~ ., 
  data = dmx_train, 
  method = "pcr",
  trControl = trainControl(method = "cv", number = 10),
  preProcess = c("center"),
  tuneLength = 100
  )

cv_model_pcr$bestTune
cv_model_pcr$results %>% dplyr::filter(ncomp == pull(cv_model_pcr$bestTune))
```

```{r diagnostics}
ggplot(cv_model_pcr)
vip(cv_model_pcr, num_features = 20, method = "model")

```
```{r logistic regression}
dmx$elt2.status = factor(ifelse(elt2_anytime,"bound", "never bound"))
split <- initial_split(dmx, prop = 0.7)
dmx_train  <- training(split)
dmx_test   <- testing(split)


set.seed(123)
cv_model1 <- train(
  elt2.status ~ fos.1_L2_1, 
  data = dmx_train, 
  method = "glm",
  family = "binomial",
  trControl = trainControl(method = "cv", number = 10)
)

set.seed(123)
cv_model2 <- train(
  elt2.status ~ fos.1_L2_1 + nhr.80_yAd_1, 
  data = dmx_train, 
  method = "glm",
  family = "binomial",
  trControl = trainControl(method = "cv", number = 10)
)

set.seed(123)
cv_model3 <- train(
  elt2.status ~ fos.1_L2_1 + nhr.80_yAd_1 + nhr.28_L4_1, 
  data = dmx_train, 
  method = "glm",
  family = "binomial",
  trControl = trainControl(method = "cv", number = 10)
)
```


```{r}
summary(
  resamples(
    list(
      model1 = cv_model1, 
      model2 = cv_model2, 
      model3 = cv_model3
    )
  )
)$statistics$Accuracy
pred_class2 <- predict(cv_model2, dmx_train)
pred_class3 <- predict(cv_model3, dmx_train)
confusionMatrix(
  data = relevel(pred_class2, ref = "bound"), 
  reference = relevel(dmx_train$elt2.status, ref = "bound")
)
```
```{r}
# Compute predicted probabilities
m1_prob <- predict(cv_model1, dmx_train, type = "prob")$bound
m2_prob <- predict(cv_model2, dmx_train, type = "prob")$bound

# Compute AUC metrics for cv_model1 and cv_model3
perf1 <- prediction(m1_prob, dmx_train$elt2.status) %>%
  performance(measure = "tpr", x.measure = "fpr")
perf2 <- prediction(m2_prob, dmx_train$elt2.status) %>%
  performance(measure = "tpr", x.measure = "fpr")

# Plot ROC curves for cv_model1 and cv_model2
plot(perf1, col = "black", lty = 2)
plot(perf2, add = TRUE, col = "blue")
legend(0.8, 0.2, legend = c("cv_model1", "cv_model2"),
       col = c("black", "blue"), lty = 2:1, cex = 0.6)
```

```{r random forests}

binaryClassifier = randomForest(y=dmx$elt2.status,x=dmx %>% select(-elt2.status), ntrees=2000, sampsize=c(60,60), importance = T)
binaryClassifier$confusion
as.data.frame(round(importance(binaryClassifier),2)) %>% arrange(-MeanDecreaseAccuracy) %>% head(20)
as.data.frame(round(importance(binaryClassifier),2)) %>% arrange(-MeanDecreaseAccuracy) %>% head(20) %>% rownames()

```

```{r heatmap-randomforest-importance}
dmxx = dmx %>% select(-elt2.status) %>% as.matrix()

# probability of each prediction
decision_p = apply(predict(binaryClassifier,type = "prob"), 1, max)

MeanDecreaseAccuracy = as.data.frame(importance(binaryClassifier))$MeanDecreaseAccuracy
dmxx.decr = dmxx[,order(-MeanDecreaseAccuracy)]
dmxx.decr[1:5,1:5]

MeanDecreaseAccuracy = MeanDecreaseAccuracy[order(-MeanDecreaseAccuracy)]

a = MeanDecreaseAccuracy > 14
b = MeanDecreaseAccuracy > 10
c = MeanDecreaseAccuracy > 5
labs = rep('[-2.3,5]',length(a))
labs[c] <- '(5,10]'
labs[b] <- '(10,14]'
labs[a] <- '(14,16]'

table(labs)
labs = factor(labs, levels = rev(c('[0,5]','(5,10]','(10,14]','(14,16]')))

# just put the values as a barplot on top of the graph
column_ha = HeatmapAnnotation(Importance=anno_barplot(MeanDecreaseAccuracy, baseline=-2.3))

row_ha = rowAnnotation(PredictionProbability=anno_barplot(decision_p))

hm.rf = Heatmap(dmxx.decr,
        cluster_rows = FALSE,
        row_order = order(decision_p,decreasing = T),
        column_title = "Random forest: which factors account for ELT-2 bound/not bound",
        cluster_columns = FALSE,
        top_annotation = column_ha,
        left_annotation = row_ha,
        heatmap_legend_param = list(title = "random forest classifier"),
        show_heatmap_legend = FALSE,
        show_row_names = FALSE, 
        show_column_names = T, 
        show_row_dend = F,
        show_column_dend = T,
        col = hmColors,
        row_split = ifelse(elt2_anytime, "ELT-2 bound", "ELT-2 never bound"),
        column_split = labs,
        column_names_side = "top",
        use_raster = T)

pdf(file=format(Sys.time(), "plots/%Y-%m-%d-%H-%M-heatmap-random-forests.pdf"), 
    width=50, height=20)
hm.rf
dev.off()
```