library(data.table)
library(magrittr)
library(ggplot2)
library(Seurat)
library(RColorBrewer)
library(viridis)
library(compiler)
library(readr)
library(matrixStats)
library(Matrix)
library(kableExtra)
devtools::source_url("https://github.com/wbaopaul/scATAC-pro/blob/master/scripts/src/scATAC-pro_utils.R?raw=TRUE")

output_dir = 'pbmc10x_scATACpro_output/'
qc_dir = paste0(output_dir, 'summary/')
configure_user_file = 'pbmc10x_scATACpro_output/configure_user.txt'

## read configure_users.txt file
read_conf_file(configure_user_file)

## read mapping stats
mapping_qc_file = paste0(qc_dir, OUTPUT_PREFIX,  '.MappingStats')
mapping_qc = fread(mapping_qc_file, header = F)

## read fragments file
fragments_file = paste0(qc_dir, OUTPUT_PREFIX,  '.fragments.tsv.gz')


frags = fread(fragments_file)
names(frags) = c('chr', 'start', 'end', 'bc', 'ndup')

## read mapping stats for cell barcodes
cell_mapping_qc_file = paste0(qc_dir, 'cell_barcodes.MappingStats')
cell_mapping_qc = fread(cell_mapping_qc_file, header = F)


## read qc per barcode file
bc_stat_file = paste0(qc_dir, OUTPUT_PREFIX, '.', PEAK_CALLER, '.qc_per_barcode.txt')
bc_stat = fread(bc_stat_file)

## read cell barcodes
cell_barcodes = fread(paste0(output_dir, '/filtered_matrix/', PEAK_CALLER, '/', CELL_CALLER, '/barcodes_doubletsRemoved.txt'), header = F)$V1

lib_complx = nrow(frags)/sum(frags$ndup)
lib_complx = paste0(100 * round(lib_complx, 3), '%')

mapping_qc$frac = round(mapping_qc$V2/mapping_qc$V2[1], 3)
mapping_qc$frac = paste0(100*mapping_qc$frac, '%')

mapping_qc = rbind(mapping_qc, data.frame(V1 ='Library Complexity (nonredudant fraction)', V2 = '', frac = lib_complx))

kable(mapping_qc, col.names = NULL, format = 'html', caption = paste('Sample:', OUTPUT_PREFIX)) %>%
  kable_styling("striped", full_width = F, position = 'left', font_size = 15)

qc_sele = bc_stat[bc %in% cell_barcodes, ]
qc_nonsele = bc_stat[!bc %in% cell_barcodes, ]

peak_file = paste0(output_dir, 'peaks/', PEAK_CALLER, '/', OUTPUT_PREFIX, '_features_BlacklistRemoved.bed')
peaks = fread(peak_file)
ncells = length(cell_barcodes)
mapq30.frac.in.cell = paste0(round(sum(frags[bc %in% cell_barcodes]$ndup)/sum(frags$ndup), 3) * 100, '%')
#frac.in.cell = paste0(round(cell_mapping_qc$V2[2]/as.integer(as.character(mapping_qc$V2[2])), 3) * 100, '%')
med.frag.per.cell = round(median(qc_sele$total_frags))
frac.uniq = paste0(round(nrow(frags[bc %in% cell_barcodes])/sum(frags[bc %in% cell_barcodes]$ndup), 3) * 100, '%')
frac.in.peak = sum(bc_stat$total_frags * bc_stat$frac_peak)/sum(bc_stat$total_frags)
frac.in.peak = paste0(round(frac.in.peak, 3)*100, '%')
frac.in.cell = sum(qc_sele$total_frags)/sum(bc_stat$total_frags)
frac.in.cell = paste0(round(frac.in.cell, 3)*100, '%')
cell.table = data.frame(c(PEAK_CALLER, nrow(peaks),
                          frac.in.peak, 
                          CELL_CALLER,  paste0(ncells),
                          frac.in.cell,
                          paste0(med.frag.per.cell), 
                          frac.in.cell, mapq30.frac.in.cell)) 
rm(peaks)
rownames(cell.table) = c('Peak called by', 'Number of peaks', 
                         'Percentage of unique fragments in peaks',
                         'Cell called by', 'Estimated # of cells',
                         'Percentage of unique fragments in cells',
                         'Median # of unique fragments per cell', 
                         'Percentage of Mapped fragments in cells', 
                         'Percentage of MAPQ30 fragments in cells')  
kable(cell.table, row.names = T, col.names = NULL, format = 'html') %>%
  kable_styling("striped", full_width = F, position = 'left', font_size = 15)

if(CELL_MAP_QC){
    lib_complx = frac.uniq

    cell_mapping_qc$frac = round(cell_mapping_qc$V2/cell_mapping_qc$V2[1], 3)
    cell_mapping_qc$frac = paste0(100*cell_mapping_qc$frac, '%')

    cell_mapping_qc = rbind(cell_mapping_qc, data.frame(V1 ='Library Complexity (nonredudant fraction)', V2 = '', frac = lib_complx))

    kable(cell_mapping_qc, col.names = NULL, format = 'html') %>%
      kable_styling("striped", full_width = F, position = 'left', font_size = 15)
}