## source("https://bioconductor.org/biocLite.R")
## biocLite("SummarizedExperiment")

file.prefix <- '~/work/TPON/fromMac/Nikolas_Final_Analysis_160922/'
se.file <- paste0(file.prefix,'output/data/summarisedExperiment.RData')
load(se.file)

library(nbHelpers)
library(SummarizedExperiment)
library(DESeq2)

## Get a hierarchy of the IgG treated cell types

colData(se)

## Wildtype IgG treated
se.b1 <- se[,se$batch == 'batch1' & se$treatment == 'IgG']
dim(se.b1)
colData(se.b1)

## get rlog transform of the data for the clustering
rld1 <- rlog(assay(se.b1))

## Fix the rld names
names1 <- paste0(se.b1$cell.population.summary,'-',se.b1$rna.prep.id)
names(names1) <- se.b1$rna.prep.id
colnames(rld1) <-names1[colnames(rld1)]

## Euclidian distance dendrogram on rlog
rld1.d.euclidian <- dist(t(rld1),method='euclidian')
hc1 <- hclust(rld1.d.euclidian)
png('hc.eucl.png')
plot(hc1,main='Euclidian Distance')
dev.off()

X11()

## Pearson correlation distance on rlog
rld1.d.corr.pearson <- as.dist(1 - cor(rld1))
hc2 <- hclust(rld1.d.corr.pearson)
png('hc.peason.png')
plot(hc2,main='Pearson Correlation Distance')
dev.off()

colData(se)

## Do a dendrogram with everything
rld1.all <- rlog(assay(se))

preserve.state()
##  "savepoint_2018-03-15_22:47:22_2734.RDataF"

library(nbHelpers)
load.image.fast("savepoint_2018-03-15_22:47:22_2734.RDataF")

###################################################################
## March 20th

## All Wildtype IgG treated
se.allb1 <- se[,se$batch == 'batch1' ]
dim(se.allb1)
colData(se.allb1)

## get rlog transform of the data for the clustering
rld1 <- rlog(assay(se.allb1))

preserve.state()
## "savepoint_2018-03-20_18:33:23_2609.RDataF"
load.image.fast("savepoint_2018-03-20_18:33:23_2609.RDataF")

colnames(rld1) <- colnames(assay(se.allb1))
## Fix the rld names
names1 <- paste0(se.allb1$cell.population.summary,'-',se.allb1$treatment,'-',se.allb1$rna.prep.id)
names(names1) <- se.allb1$rna.prep.id
colnames(rld1) <-names1[colnames(rld1)]
colnames(rld1)

## ## Euclidian distance dendrogram on rlog
## rld1.d.euclidian <- dist(t(rld1),method='euclidian')
## hc1 <- hclust(rld1.d.euclidian)
## png('hc.eucl.png')
## plot(hc1,main='Euclidian Distance')
## dev.off()

## Pearson correlation distance on rlog
rld1.d.corr.pearson <- as.dist(1 - cor(rld1))
hc2 <- hclust(rld1.d.corr.pearson, method='ward.D2')
hc2d <- as.dendrogram(hc2)

#install.packages('dendextend')
library(dendextend)

labs <- hc2d %>% get_leaves_attr('label')

pop <- strpart(labs, '-', 1)
trt <- strpart(labs, '-', 2)

library(RColorBrewer)

## Color labels by population
pop.pal <- rainbow(length(unique(pop)),1,0.8)
names(pop.pal) <- unique(pop)
pop.cols <- pop.pal[pop]

## 
trt.pal <- c('CD42b'='red','IgG'='blue')
trt.cols <- trt.pal[trt]
trt.cols

png('dend.png',width=1600,height=900)
par(mar=c(20,1,1,1),font=2)
hc2d %>% set('labels_col', pop.cols) %>% set('labels_cex',1.3) %>% set('leaves_pch',19) %>%  set('leaves_cex',2) %>%
    set('leaves_col',trt.cols) %>% plot(main='Batch 1 samples')
legend(x='topright',legend=names(trt.pal),fill=trt.pal)
legend(x='right', legend=names(pop.pal), fill=pop.pal)
dev.off()

new.labs <- lapply(lapply(strsplit(labs,'-'),'[',c(1,2)), function(x) {paste(x, collapse='-')} )

png('dend2.png',width=1600,height=900)
par(mar=c(20,4,4,4),font=2, cex.axis=2.5)
hc2d %>% set('labels_col', pop.cols) %>%
    set('labels_cex',1.3) %>%
    set('branches_lwd', 6) %>%
    set('leaves_pch',19) %>%
    set('leaves_cex',2) %>%
    set('labels',new.labs) %>%
    set('leaves_col',trt.cols) %>% plot(main='Batch 1 samples')
legend(x='topright',legend=names(trt.pal),fill=trt.pal,cex=1.3)
legend(x='right', legend=names(pop.pal), fill=pop.pal,cex=1.3)
dev.off()

pdf('dend2.pdf',width=16,height=9)
par(mar=c(15,4,4,4),font=2, cex.axis=2)
hc2d %>% set('labels_col', pop.cols) %>%
    set('labels_cex',1.3) %>%
    set('branches_lwd', 5) %>%
    set('leaves_pch',19) %>%
    set('leaves_cex',2) %>%
    set('labels',new.labs) %>%
    set('leaves_col',trt.cols) %>% plot(main='Batch 1 samples')
legend(x='right',legend=names(trt.pal),fill=trt.pal,cex=1)
legend(x='topright', legend=names(pop.pal), fill=pop.pal,cex=1)
dev.off()


pdf('dend2b.pdf',width=16,height=9)
par(mar=c(15,4,4,4),font=2, cex.axis=2)
hc2d %>% set('labels_col', pop.cols) %>%
    set('labels_cex',1.3) %>%
    set('branches_lwd', 4) %>%
    set('leaves_pch',19) %>%
    set('leaves_cex',2) %>%
    set('labels',new.labs) %>%
    set('leaves_col',trt.cols) %>% plot(main='Batch 1 samples')
legend(x='right',legend=names(trt.pal),fill=trt.pal,cex=1)
legend(x='topright', legend=names(pop.pal), fill=pop.pal,cex=1)
dev.off()

##

new.labs.2 <- lapply(lapply(strsplit(labs,'-'),'[',c(3)), function(x) {paste(x, collapse='-')} )

png('dend3.png',width=1600,height=900)
par(mar=c(10,1,1,1),font=2)
hc2d %>% set('labels_col', pop.cols) %>% set('labels_cex',1.3) %>% set('leaves_pch',19) %>%  set('leaves_cex',2) %>%
    set('labels',new.labs.2) %>%
    set('leaves_col',trt.cols) %>% plot(main='Batch 1 samples')
legend(x='topright',legend=names(trt.pal),fill=trt.pal)
legend(x='right', legend=names(pop.pal), fill=pop.pal)
dev.off()


pdf('dend3.pdf')
par(mar=c(10,1,1,1),font=2)
hc2d %>% set('labels_col', pop.cols) %>% set('labels_cex',1.3) %>% set('leaves_pch',19) %>%  set('leaves_cex',2) %>%
    set('labels',new.labs.2) %>%
    set('leaves_col',trt.cols) %>% plot(main='Batch 1 samples')
legend(x='topright',legend=names(trt.pal),fill=trt.pal)
legend(x='right', legend=names(pop.pal), fill=pop.pal)
dev.off()

setDisplay('10.0')

X11()

## png('hc.peason.png')
plot(hc2d,main='Pearson Correlation Distance')
## dev.off()