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, row.names=10, stringsAsFactors=FALSE)
countsTable = inputTable[,11:ncol(inputTable)]
colNames <- gsub("HeLa.S3.", "", colnames(countsTable))
cond <- factor(c(gsub("...$", "", colNames)))
colnames(countsTable) <- colNames
colors <- brewer.pal(6, "Set1")
colorSelection = colors[c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6)]
set <- as.matrix(countsTable)
colnames(set) <- colNames

readFilterSelectionTable <- countsTable[,
   c("mock_GFP_minus.24_h.R1", "mock_GFP_minus.24_h.R2", "mock_GFP_minus.24_h.R3",
     "WT_GFP_plus.24_h.R1", "WT_GFP_plus.24_h.R1", "WT_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)

contrasts <- makeContrasts(
    WT_GFP_minus_04_h_vs_mock_GFP_minus_04_h=condWT_GFP_minus.04_h-condmock_GFP_minus.04_h,
    WT_GFP_minus_24_h_vs_mock_GFP_minus_24_h=condWT_GFP_minus.24_h-condmock_GFP_minus.24_h,
    WT_GFP_plus_04_h_vs_mock_GFP_minus_04_h=condWT_GFP_plus.04_h-condmock_GFP_minus.04_h,
    WT_GFP_plus_24_h_vs_mock_GFP_minus_24_h=condWT_GFP_plus.24_h-condmock_GFP_minus.24_h,        
    levels=design)
    
comparisons <- c(
    "WT_GFP_minus_04_h_vs_mock_GFP_minus_04_h",
    "WT_GFP_minus_24_h_vs_mock_GFP_minus_24_h",
    "WT_GFP_plus_04_h_vs_mock_GFP_minus_04_h",
    "WT_GFP_plus_24_h_vs_mock_GFP_minus_24_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)
}
