#!/usr/bin/R
# --------------------------------------------------
# Aggregate the gene fusions at the individual level
# for prioritizing 10x ONT
# Updated: 01/10/23
# --------------------------------------------------
library(cbrbase)
set_proj('DEVTRAJ')

# Libraries:
library(tidyr)
library(ggplot2)
library(ggpubr)
library(RColorBrewer)
options(width=170)

# Load in and process data (saves to matrices):
# project = 'ALZ_repl'
project = 'ALZ'
commandArgs <- function(){ c(project)}
source(paste0(bindir, 'mosaicism/load_metadata.R'))


# Make directories for analysis:
analysis = 'fusion'
imgdir = paste0(img, project, '/', analysis , '/')
imgpref = paste0(imgdir, analysis, "_mspred_", project, '_')
cmd = paste0('mkdir -p ', imgdir)
system(cmd)


# Load in the two cohorts of candidate gene fusion counts:
# --------------------------------------------------------
source(paste0(bindir, 'mosaicism/genefusions/load_fusion_data.R'))
source(paste0(bindir, 'mosaicism/genefusions/auxiliary_fusion_plotting.R'))


fulldf = c()
for (only.exc in c(TRUE, FALSE)){
    if (only.exc){
        df = sub.ctdf[sub.ctdf$relabel == 'Excitatory',]
    } else { df = sub.ctdf }
    aggdf = aggregate(cbind(count, n_counts) ~ projid + platform + cohort + bd + nrad, df, sum)
    aggdf$rate = with(aggdf, count / n_counts)
    aggdf = aggdf[order(aggdf$rate),]
    aggdf$only.exc = only.exc
    fulldf = rbind(fulldf, aggdf)
}


# Save this table
outpref = paste0(dbdir,'variants/ALZ_joint_fusion_datasets_')
filesuff = 'maxsensitivity_prediction_aggregated'
agg.tsv = paste0(outpref, filesuff, '.tsv')
write.table(fulldf, agg.tsv, quote=F, row.names=F, sep="\t")