#!/usr/bin/R
# -----------------------------------------------
# Plot stats for fusions in the CKP25 mouse data:
# -----------------------------------------------
library(cbrbase)
set_proj('DEVTRAJ', 'mosaicism')

library(ggplot2)
library(ggpubr)
library(ggbeeswarm)
library(tidyr)
library(scales)
library(RColorBrewer)

source(paste0(sbindir, 'CKP25/load_metadata.R'))

# Directories:
imgdir = paste0(img, 'CKP25/')
imgdir2 = paste0(imgdir, 'fusions/')
imgpref = paste0(imgdir2, 'fusions_')
cmd = paste('mkdir -p',imgdir, imgdir2)
system(cmd)

# Load after defining imgpref
source(paste0(sbindir, 'CKP25/analysis/load_basic_fusions.R'))
source(paste0(sbindir, 'CKP25/analysis/load_umap_basis.R'))


# Plot fusion events on the UMAP:
# -------------------------------
# smeta$has_intra = smeta$cell %in% intradf$cell[intradf$close.fusion]
smeta$has_intra = smeta$cell %in% intradf$cell
smeta$has_inter = smeta$cell %in% interdf$cell
smeta$event_count = (smeta$has_intra + 2 * smeta$has_inter)
smeta = smeta[order(smeta$event_count, decreasing=F),]

smeta$event_type = 'None'
smeta$event_type[smeta$event_count == 1] = 'Intra'
smeta$event_type[smeta$event_count == 2] = 'Inter'
smeta$event_type[smeta$event_count == 3] = 'Both'

fusion.cols = c('grey80','indianred','slateblue2','purple')

png(paste0(imgpref, 'UMAP_intra_fusions.png'),res=300, width=3,height=3, units='in')
par(mar=c(0,0,0,0))
plot(smeta[[u1]], smeta[[u2]], pch=19, col=ifelse(smeta$has_intra, 'indianred', 'grey80'),
     axes=F, ylab='', xlab='', cex=.5)
dev.off()

png(paste0(imgpref, 'UMAP_inter_fusions.png'),res=300, width=3,height=3, units='in')
par(mar=c(0,0,0,0))
plot(smeta[[u1]], smeta[[u2]], pch=19, col=ifelse(smeta$has_inter, 'slateblue2', 'grey80'),
     axes=F, ylab='', xlab='', cex=.5)
dev.off()

png(paste0(imgpref, 'UMAP_event_count_fusions.png'),res=300, width=3,height=3, units='in')
par(mar=c(0,0,0,0))
plot(smeta[[u1]], smeta[[u2]], pch=19, col=fusion.cols[smeta$event_count + 1],
     axes=F, ylab='', xlab='', cex=.5)
dev.off()

png(paste0(imgpref, 'UMAP_event_has_fusion.png'),res=300, width=3,height=3, units='in')
par(mar=c(0,0,0,0))
plot(smeta[[u1]], smeta[[u2]], pch=19, col=ifelse(smeta$event_count > 0, 'slateblue3', 'grey80'),
     axes=F, ylab='', xlab='', cex=.5)
dev.off()

png(paste0(imgpref, 'UMAP_V1V2_event_has_fusion.png'),res=300, width=3,height=3, units='in')
par(mar=c(0,0,0,0))
plot(smeta$V1, smeta$V2, pch=19, col=ifelse(smeta$event_count > 0, 'slateblue3', 'grey80'),
     axes=F, ylab='', xlab='', cex=.5)
dev.off()



gplot = ggplot(smeta, aes(event_type, n_counts, fill=factor(event_count))) + 
    geom_violin() + 
    scale_fill_manual(values=fusion.cols, name='Sub-celltype:') + 
    scale_y_continuous(expand=c(0,0), labels=scales::comma) + 
    theme_pubr() + 
    labs(x='Fusion event type', y='Number of reads from the cell') + 
    theme(axis.text.x=element_text(angle=90,vjust=.5, hjust=1)) +
    theme(legend.position="none")
ggsave(paste0(imgpref, 'qcstats_cell_fusions.png'), gplot, width=4, height=5, units='in', dpi=450)
ggsave(paste0(imgpref, 'qcstats_cell_fusion.pdf'), gplot, width=4, height=5)


# Are the events co-localized to the same cells?
mat = table(smeta$has_inter, smeta$has_intra)
ft = fisher.test(mat)
sprintf('%0.2e', ft$p) # 6.72e-7 (highly colocalized)