#!/usr/bin/R
# ------------------------------------------------------
# Check if there is a subtype-specific depletion for HC:
# Updated: 01/13/23
# ------------------------------------------------------
library(cbrbase)
library(ggplot2)
library(ggpubr)
options(width=170)

set_proj('DEVTRAJ', 'mosaicism')
source(paste0(bindir, 'mosaicism/fusions_human10x/load_metadata.R'))
source(paste0(bindir, 'multiRegion/auxiliary_plotting_settings.R'))

# Directories:
plotdir = paste0(imgdir, 'metadata/')
imgpref = paste0(plotdir, 'human10x_')
cmd = paste('mkdir -p', plotdir)
system(cmd)


# Load the preliminary fusion file:
# ---------------------------------
# Subset cellmeta to HC data:
cellmeta$batch = sub("-.*", "", cellmeta$dataset)
hc.cellmeta = cellmeta[cellmeta$batch == 'HC',]
hc.cellmeta$CB = sub(".*:", "", hc.cellmeta$barcode)

ntotdf = agg.rename(CB ~ celltype + dataset, hc.cellmeta, length, 'total')
ctdf = agg.rename(CB ~ cthr + celltype + dataset, hc.cellmeta, length, 'count')
comb.df = unique(hc.cellmeta[,c('celltype','cthr')])
comb.df = merge(comb.df, data.frame(dataset=unique(hc.cellmeta$dataset)))
ctdf = merge(comb.df, ctdf, all.x=TRUE)
ctdf = merge(ctdf, ntotdf, all.x=TRUE)
ctdf$count[is.na(ctdf$count)] = 0
ctdf$total[is.na(ctdf$total)] = 0

ctdf = ctdf[ctdf$total > 25,]
ctdf$frac = ctdf$count / ctdf$total
ctdf = merge(ctdf, unique(fus.meta[,c('dataset', 'projid')]))
ctdf = merge(ctdf, unique(metadata[,c('projid','niareagansc', 'braaksc')]))
ctdf$nrad = ifelse(ctdf$niareagansc > 2, 'CTRL', 'AD')
ctdf$braaksc.ad = ifelse(ctdf$braaksc < 5, 'CTRL', 'AD')


gp = ggplot(ctdf, aes(cthr, frac, fill=nrad)) + 
    facet_wrap(~celltype, scales='free')+ 
    geom_boxplot() + 
    stat_compare_means() +
    theme_pubr() + coord_flip()


gp = ggplot(ctdf, aes(cthr, frac, fill=braaksc.ad)) + 
    facet_wrap(~celltype, scales='free')+ 
    geom_boxplot() + 
    stat_compare_means() +
    theme_pubr() + coord_flip()


# No change in HC proportions in this cohort