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

main <- function(){
    args <- commandArgs(trailingOnly = TRUE)
    inputFile <- args[1]
    outputFolder <- args[2]
    minNumberOfReads <- 50
    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]
    
    prefix <- "HeLa-S3-SPI-2_pinT_vector_05_min_ara_vs_HeLa-S3-SPI-2_empty_vector_05_min_ara"
    runEdgeR(countsTable[, 1:4], prefix, minNumberOfReads,
             outputFolder)

    prefix <- "in-vitro-SPI-1_pinT_vector_05_min_ara_vs_in-vitro-SPI-1_empty_vector_05_min_ara"
    runEdgeR(countsTable[, 5:8], prefix, minNumberOfReads,
             outputFolder)

    prefix <- "in-vitro-SPI-2_pinT_vector_05_min_ara_vs_in-vitro-SPI-2_empty_vector_05_min_ara"
    runEdgeR(countsTable[, 9:12], prefix, minNumberOfReads,
             outputFolder)
}

runEdgeR <- function(curCountsTable, prefix, minNumberOfReads,  outputFolder){
    print(colnames(curCountsTable))
    filter <- apply(curCountsTable, 1, function(x){
       length(x[x>minNumberOfReads])>=2})
    
    filtered <- curCountsTable[filter, ]
    set <- as.matrix(filtered)
    
    pdf(file=sprintf("%s/edgeR-plots-%s.pdf", outputFolder, prefix))    
    colors <- brewer.pal(6, "Set1")
    colorSelection = colors[c(1, 1, 2, 2)]
    plotRLE(set, ylim=c(-1,1), main="RLE, pre-norm", col=colorSelection)
    plotPCA(set, main="PCA, pre-norm", col=colorSelection)

    cond <- factor(c("empty", "empty", "pinT", "pinT"))

    # ~0 to remove the intercept
    design <- model.matrix(~0+cond)
    
    y <- DGEList(counts=set, group=cond)
    y <- calcNormFactors(y, method="upperquartile")
    y <- estimateGLMCommonDisp(y, design)
    y <- estimateGLMTagwiseDisp(y, design)
    fit <- glmFit(y, design, prior.count=1)
    
    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)
    
    contrasts <- makeContrasts(pinT.vs.empty=condpinT-condempty,
                               levels=design)

    lrt <- glmLRT(fit, contrast=contrasts[, "pinT.vs.empty"])
    top <- topTags(lrt, n=nrow(set))$table
    
    write.table(top,
                file=sprintf("%s/edgeR-comparison-%s.csv",
                             outputFolder, prefix),
                sep="\t", quote=FALSE, row.names=TRUE, col.names=NA)    
    dev.off()
}

main()
