require(RUVSeq)
library(edgeR)
library(RColorBrewer)
library(gplots)

args <- commandArgs(trailingOnly = TRUE)
inputFile <- args[1]
outputFolder <- args[2]
minNumberOfReads <- 10

inputTable = read.table(
    inputFile, sep="\t", header=TRUE, stringsAsFactors=FALSE)
# Generate unique row name
row.names(inputTable) <- paste(inputTable$Locus_tag, inputTable$Gene_name,
                               inputTable$Start, inputTable$Feature,
                               inputTable$Attributes, sep="|")

countsTable = inputTable[,13:ncol(inputTable)-2]

colNames <- gsub("HeLa.S3.WT.", "", colnames(countsTable))
cond <- factor(c(gsub("...$", "", colNames)))
colnames(countsTable) <- colNames
colors <- brewer.pal(4, "Set1")
colorSelection = colors[c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4)]
set <- as.matrix(countsTable)
colnames(set) <- colNames
    
readFilterSelectionTable <- countsTable[,
        c("LB.00_h.R1", "LB.00_h.R2", "GFP_plus.24_h.R1", "GFP_plus.24_h.R1")]    
    
filter <- apply(readFilterSelectionTable, 1, function(x){
    length(x[x>minNumberOfReads])>=3})
filtered <- countsTable[filter,]
set <- as.matrix(filtered)
genes <- rownames(filtered)

pdf(file=sprintf("%s/edgeR-plots-all_comparisons.pdf", outputFolder))

plotRLE(set, ylim=c(-1,1), main="RLE, pre-norm", col=colorSelection)
plotPCA(set, main="PCA, pre-norm", col=colorSelection)
design <- model.matrix(~0+cond)

y <- DGEList(counts=set, group=cond)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)

plotRLE(cpm(y, normalized.lib.sizes=T), ylim=c(-1,1),
        main="RLE, post UQ norm", col=colorSelection)
plotPCA(cpm(y, normalized.lib.sizes=T),
        main="PCA, post UQ norm", col=colorSelection)

fit <- glmFit(y, design, prior.count=1)
    
contrasts <- makeContrasts(
    GFP_plus.04_h_vs_input.00_h=condGFP_plus.04_h-condinput.00_h,
    GFP_plus.24_h_vs_input.00_h=condGFP_plus.24_h-condinput.00_h,
    GFP_plus.04_h_vs_LB.00_h=condGFP_plus.04_h-condLB.00_h,
    GFP_plus.24_h_vs_LB.00_h=condGFP_plus.24_h-condLB.00_h,
    levels=design)

comparisons <- c(
    "GFP_plus.04_h_vs_input.00_h", "GFP_plus.24_h_vs_input.00_h",
    "GFP_plus.04_h_vs_LB.00_h", "GFP_plus.24_h_vs_LB.00_h")

for (i in seq(length(comparisons))){
    curLrt <- glmLRT(fit, contrast=contrasts[, comparisons[i]])
    top <- topTags(curLrt, n=nrow(set))$table
    write.table(top, file=sprintf("%s/edgeR-comparison-%s.csv",
                         outputFolder, gsub("\\.", "_", comparisons[i])),
                sep="\t", quote=FALSE, row.names=TRUE, col.names=NA)
}
