rm(list = ls())
options <- commandArgs(trailingOnly = TRUE)
print(sprintf("source code is %s", as.character(options[1])))
source(as.character(options[1]))
# source("/Volumes/yli/CRI-BIO-502-GMutulu-RNAseq/rna-seq/build/module-ve/DE-Analysis-AllComb-2.R")
countTabpath <- as.character(options[2])
metaTabpath <- as.character(options[3])
group1 <- as.character(options[4])
group2 <- as.character(options[5])
workDir <- as.character(options[6])
samp.cutoff <- as.numeric(options[7])
cpm.val <- as.numeric(options[8])
padjOn <- as.logical(options[9])

padj.val <- c(0.01, 0.05, 0.1)
fc.val <- c(0, 1.5, 2)
#######################################
###########################################
###DEseq2 analysis 
setwd(workDir)
print("============================")
print("Start DESeq2 DE analysis")
##This above path need to be kept the same in the source code, "/group/bioinformatics/Pipelines/Development/RNAseq_beta/RNAseq.v9.0.1/build/module-ve/DE-Analysis-AllComb.R"
print(sprintf("Raw count is from %s, low expression is removed at cpm = %.0f, group no = % .0f", 
              countTabpath,
              as.numeric(cpm.val),
              as.numeric(samp.cutoff)))
deseq2.res <- deseq2.de.analysis(countTabpath = as.character(countTabpath), 
                                 metaTabpath = as.character(metaTabpath), 
                                 cpmCutoff = as.numeric(cpm.val), sampCutoff=as.numeric(samp.cutoff), 
                                 group1=as.character(trim(group1)), group2=as.character(trim(group2)))
# save(deseq2.res, file = paste(deseq2.res$res.dir,"res.Rdata",sep = "/"))
GeneRegFilter="both"
print(sprintf("The filter is on %s", as.character(as.character(GeneRegFilter))))
print(sprintf("The comparision is between %s and %s", as.character(trim(group1)), as.character(trim(group2))))
print(sprintf("The results are saved in %s", as.character(deseq2.res$res.dir)))

##The FC is calculated on group2/group1
de.num.padj.fc.deseq2.both <- matrix(NA, nrow = 3, ncol = 3)
for (i in 1:3){
  for (j in 1:3){
    deseq2.filter.res <- get.filtered.deg(de.res = deseq2.res, 
                                         padjust = as.logical(padjOn),
                                         Pcutoff = as.numeric(padj.val[j]),
                                         FCcutoff = as.numeric(fc.val[i]),
                                         GeneRegFilter = GeneRegFilter,
                                         de.method = "DESeq2")
    
    de.num.padj.fc.deseq2.both[i,j] = length(deseq2.filter.res$degName)
  }
}
colnames(de.num.padj.fc.deseq2.both) <- paste("adj-p < ", as.numeric(padj.val), sep = "")
de.num.padj.fc.deseq2.both <- as.data.frame(de.num.padj.fc.deseq2.both)
de.num.padj.fc.deseq2.both$FC <- paste("|FC| > ", as.numeric(fc.val), sep = "")
de.num.padj.fc.deseq2.both$filter <- as.character(GeneRegFilter)

# fname <- paste(as.character(deseq2.res.path),"/DESeq2-summary-filter-", as.character(GeneRegFilter), ".txt", sep = "")
# write.table(de.num.padj.fc.deseq2, as.character(fname), quote = F, col.names = NA, row.names = T, sep = "\t")
print("**********")
GeneRegFilter="up"
print(sprintf("The filter is on %s", as.character(as.character(GeneRegFilter))))
de.num.padj.fc.deseq2.up <- matrix(NA, nrow = 3, ncol = 3)
for (i in 1:3){
  for (j in 1:3){
    deseq2.filter.res <- get.filtered.deg(de.res = deseq2.res, 
                                          padjust = as.logical(padjOn),
                                          Pcutoff = as.numeric(padj.val[j]),
                                          FCcutoff = as.numeric(fc.val[i]),
                                          GeneRegFilter = GeneRegFilter,
                                          de.method = "DESeq2")
    
    de.num.padj.fc.deseq2.up[i,j] = length(deseq2.filter.res$degName)
  }
}
colnames(de.num.padj.fc.deseq2.up) <- paste("adj-p < ", as.numeric(padj.val), sep = "")
de.num.padj.fc.deseq2.up <- as.data.frame(de.num.padj.fc.deseq2.up)
de.num.padj.fc.deseq2.up$FC <- paste("|FC| > ", as.numeric(fc.val), sep = "")
de.num.padj.fc.deseq2.up$filter <- as.character(GeneRegFilter)

# fname <- paste(as.character(deseq2.res.path),"/DESeq2-summary-filter-", as.character(GeneRegFilter), ".txt", sep = "")
# write.table(de.num.padj.fc.deseq2, as.character(fname), quote = F, col.names = NA, row.names = T, sep = "\t")
print("**********")
GeneRegFilter="down"
print(sprintf("The filter is on %s", as.character(as.character(GeneRegFilter))))
de.num.padj.fc.deseq2.down <- matrix(NA, nrow = 3, ncol = 3)
for (i in 1:3){
  for (j in 1:3){
    deseq2.filter.res <- get.filtered.deg(de.res = deseq2.res, 
                                          padjust = as.logical(padjOn),
                                          Pcutoff = as.numeric(padj.val[j]),
                                          FCcutoff = as.numeric(fc.val[i]),
                                          GeneRegFilter = GeneRegFilter,
                                          de.method = "DESeq2")
    
    de.num.padj.fc.deseq2.down[i,j] = length(deseq2.filter.res$degName)
  }
}
colnames(de.num.padj.fc.deseq2.down) <- paste("adj-p < ", as.numeric(padj.val), sep = "")
de.num.padj.fc.deseq2.down <- as.data.frame(de.num.padj.fc.deseq2.down)
de.num.padj.fc.deseq2.down$FC <- paste("|FC| > ", as.numeric(fc.val), sep = "")
de.num.padj.fc.deseq2.down$filter <- as.character(GeneRegFilter)

# fname <- paste(as.character(deseq2.res.path),"/DESeq2-summary-filter-", as.character(GeneRegFilter), ".txt", sep = "")
# write.table(de.num.padj.fc.deseq2, as.character(fname), quote = F, col.names = NA, row.names = T, sep = "\t")
print("**********")

de.num.padj.fc.deseq2 <- rbind(de.num.padj.fc.deseq2.both, de.num.padj.fc.deseq2.up, de.num.padj.fc.deseq2.down)
if ( grepl('batch.R$', as.character(options[1])) ) {
  fname <- paste(as.character(deseq2.res$res.dir), "/", "batch-DESeq2-summary-", as.character(trim(group1)), "_", as.character(trim(group2)), ".txt", sep = "")
  
} else {
  fname <- paste(as.character(deseq2.res$res.dir), "/", "DESeq2-summary-", as.character(trim(group1)), "_", as.character(trim(group2)), ".txt", sep = "")
  
}
print(sprintf("The NO. DEGs summary is saved in %s", fname))
write.table(de.num.padj.fc.deseq2, as.character(fname), quote = F, col.names = T, row.names = F, sep = "\t")

print("End DESeq2 DE analysis")
print("============================")
