#!/usr/bin/R
# ------------------------------------
# Load fusion calls for the 10x human:
# --follows P301S
# Updated: 01/11/23 
# ------------------------------------
library(cbrbase)
set_proj('DEVTRAJ', 'mosaicism')
options(width=170)


# Fusion call files:
# ------------------
fsdir = paste0(dbdir, 'fusions_AD430/')
prefix = 'AD430.aggregated.fusion_predictions'
predfile = paste0(fsdir, prefix, '.tsv')
bcfile = paste0(fsdir, prefix, '.byread.withBC.tsv')
procfile = paste0(fsdir, prefix, '.processed.tsv')
countfile = paste0(fsdir, prefix, '.processed_counts.tsv')

# Read files:
fsdf = read.delim(predfile, header=T)
bcdf = read.delim(bcfile, header=T)
merged.fsdf = read.delim(procfile, header=T)
ctdf = read.delim(countfile, header=T)

# Order variables:
ctdf$bd = factor(ctdf$bd, levels=c('CTRL','AD'))
ctdf$bd2 = factor(ctdf$bd2, levels=c('CTRL','AD'))
ctdf$nrad = factor(ctdf$nrad, levels=c('CTRL','AD'))
ctdf$cogdx.ad = factor(ctdf$cogdx.ad, levels=c('Low/Mild','AD'))

# Alternative regression vars:
ctdf$has_fusion = ctdf$count > 0
ctdf$capped_count = ctdf$count
ctdf$capped_count[ctdf$capped_count > 10] = 10
ctdf$rate = ctdf$count / ctdf$n_counts
ctdf$capped_rate = ctdf$capped_count / ctdf$n_counts
reg.vars = c('has_fusion', 'count', 'capped_count', 'rate', 'capped_rate')

ctdf$age_rescaled = ctdf$age_death / 10
ctdf$pmi_rescaled = ctdf$pmi / 10


# Subset data:
# ------------
# countcutoff = 25000
# countcutoff = 0
# sub.ctdf = ctdf[ctdf$n_counts >= countcutoff,]
sub.ctdf = ctdf


# Regression estimates for a number of selected covariates:
# ------------------------------------------------
advars = c('niareagansc','nrad','braaksc','bd', 'bd2','cogdx.ad')
covars = c('age_rescaled','msex','pmi_rescaled', 'Apoe_e4')
test.vars= c('cogdx.ad','bd','nrad', covars)
ext.covar = '+ pmi_rescaled + msex + age_rescaled'
# count.covar = '+ n_counts + n_genes + count_per_gene'
count.covar = '+ log(n_genes) + count_per_gene'
# count.covar = '+ batch * (n_counts + n_genes + count_per_gene)'


# Variable mappings:
varmap = c(
    'msex' = 'Sex (M)',
    'nradAD' = 'NIA-Reagan score (1-2)',
    'bd2AD' = 'Braak Stage (3-6)',
    'bdAD' = 'Braak Stage (5-6)',
    'Apoe_e4Yes' = 'Has APOE-e4',
    'cogdx.adAD' = 'AD Cog. Imp.',
    'pmi_rescaled' = 'PMI', 'age_rescaled' = 'Age of death')