source("/g/scb2/zaugg/daga/Universal_Scripts/Diffbind_visualization_functions.R")
setwd("/g/scb2/zaugg/daga/Autoimmune_all/Scripts")
library("DiffBind")


disease_meta = read.table("../Consensus_meta_table/Psoriasis_SLE_Graves_JIASF_ex_consensus_metadata_diffbind.txt",
                          header = T,sep = "\t")

metatable = disease_meta

####DBA object
data_dba =  dba(sampleSheet= as.data.frame(metatable))

############# Overlap Rate overall
#olap.rate <- dba.overlap(data_dba,mode=DBA_OLAP_RATE)
#plot(olap.rate,type='b',ylab='# peaks', xlab='Overlap at least this many peaksets')


output_file_name_path =  "../Consensus_peaks_raw_counts/Psoriais_SLE_Graves_JIASF_ex_ovelap3_with_nosummit.txt"

raw_reads_count = dba.count(DBA = data_dba ,minOverlap = 3,score = DBA_SCORE_READS)
raw_reads_df <- dba.peakset(raw_reads_count, bRetrieve=T, DataType=DBA_DATA_FRAME)
write.table(x = raw_reads_df,file = output_file_name_path,quote = F,sep = "\t",row.names = F,col.names = T)


to_filter_standard_chr_only = function (input_file,output_fn) {
  
  file_name = basename(input_file)
  
  df_input = read.table(file = input_file,header = T,sep = "\t")
  chr_name = paste("chr",c(seq(1:22),c("X","Y")),sep = "")
  
  select_rowids = which(df_input[,1]%in%chr_name == TRUE)
  df_mod = df_input[select_rowids,]
  
  write.table(x = df_mod ,file = output_fn,append = F,quote = F,sep = "\t",row.names = F,col.names = T)
  
}
to_filter_standard_chr_only(input_file = "../Consensus_peaks_raw_counts/Psoriais_SLE_Graves_JIASF_ex_ovelap3_with_nosummit.txt",
                            output_fn = "../Consensus_peaks_raw_counts/Psoriais_SLE_Graves_JIASF_ex_ovelap3_with_nosummit_with_nosummit_std_chr_only.txt")