library(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 names
row.names(inputTable) <- paste(inputTable$Locus_tag, inputTable$Gene_name,
                               inputTable$Start, inputTable$Feature,
                               inputTable$Attributes, sep="|")
countsTable = inputTable[,31:ncol(inputTable)-2]
colNames <- gsub("HeLa.S3.", "", colnames(countsTable))
colnames(countsTable) <- colNames
cond <- factor(c(gsub("...$", "", colNames)))
colors <- brewer.pal(7, "Set1")
colorSelection = colors[c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5,
                          6, 6, 6, 7, 7, 7)]

readFilterSelectionTable <- countsTable[,
      c("WT.00_h.R1", "WT.00_h.R2", "WT.00_h.R3",
        "WT.24_h.R1", "WT.24_h.R2", "WT.24_h.R3")]

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_WT_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)

fit <- glmFit(y, design, prior.count=1.0)
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(
    WT.02_h_vs_WT.00_h=condWT.02_h-condWT.00_h,
    WT.04_h_vs_WT.00_h=condWT.04_h-condWT.00_h,
    WT.08_h_vs_WT.00_h=condWT.08_h-condWT.00_h,
    WT.16_h_vs_WT.00_h=condWT.16_h-condWT.00_h,
    WT.24_h_vs_WT.00_h=condWT.24_h-condWT.00_h,
    levels=design)

comparisons <- c(
    "WT.02_h_vs_WT.00_h",
    "WT.04_h_vs_WT.00_h",
    "WT.08_h_vs_WT.00_h",
    "WT.16_h_vs_WT.00_h",
    "WT.24_h_vs_WT.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)
    plot(top$logCPM, top$logFC, main=comparisons[i])
}

sessionInfo()
