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 <- colnames(countsTable)
    
cond <- factor(c(gsub("...$", "", colNames)))
colnames(countsTable) <- colNames
colors <- brewer.pal(12, "Set3")
colorSelection = colors[
    c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12,
      1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12,
      1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9,
      10, 10, 10, 11, 11, 11, 12, 12, 12, 1, 1, 1, 2, 2, 2,
      3, 3, 4, 4)]
set <- as.matrix(countsTable)
colnames(set) <- colNames
#readFilterSelectionTable <- countsTable
    readFilterSelectionTable <- countsTable[,
        c("AGS.WT.00_h.R1", "AGS.WT.00_h.R2",
          "AGS.WT.24_h.R1", "AGS.WT.24_h.R1",
          "CaCo.2.WT.00_h.R1", "CaCo.2.WT.00_h.R2",
          "CaCo.2.WT.24_h.R1", "CaCo.2.WT.24_h.R2",
          "HEK.WT.00_h.R1", "HEK.WT.00_h.R2",
          "HEK.WT.24_h.R1", "HEK.WT.24_h.R2",
          "HT29.WT.00_h.R1", "HT29.WT.00_h.R2",
          "HT29.WT.20_h.R1", "HT29.WT.20_h.R2",
          "LoVo.WT.00_h.R1", "LoVo.WT.00_h.R2",
          "LoVo.WT.20_h.R1", "LoVo.WT.20_h.R2",
          "THP.1_with_PMA.WT.00_h.R1", "THP.1_with_PMA.WT.00_h.R2",
          "THP.1_with_PMA.WT.24_h.R1", "THP.1_with_PMA.WT.24_h.R2",
          "THP.1_without_PMA.WT.00_h.R1", "THP.1_without_PMA.WT.00_h.R2",
          "THP.1_without_PMA.WT.24_h.R1", "THP.1_without_PMA.WT.24_h.R2",
          "BMDM.WT.00_h.R1", "BMDM.WT.00_h.R2",
          "BMDM.WT.20_h.R1", "BMDM.WT.20_h.R2",
          "L929.WT.00_h.R1", "L929.WT.00_h.R2",
          "L929.WT.06_h.R1", "L929.WT.06_h.R2",
          "MEF.WT.00_h.R1", "MEF.WT.00_h.R2",
          "MEF.WT.20_h.R1", "MEF.WT.20_h.R2",
          "RAW264.7.WT.00_h.R1", "RAW264.7.WT.00_h.R2",
          "RAW264.7.WT.24_h.R1", "RAW264.7.WT.24_h.R2",
           "X3D4_31.WT.00_h.R1",  "X3D4_31.WT.00_h.R2",  "X3D4_31.WT.00_h.R3",
          "X3D4_31.WT.16_h.R1",  "X3D4_31.WT.16_h.R2",  "X3D4_31.WT.16_h.R3",
          "IPEC.J2.WT.00_h.R1", "IPEC.J2.WT.00_h.R2",
          "IPEC.J2.WT.20_h.R1", "IPEC.J2.WT.20_h.R2")]
    
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_misc_cell_lines.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", p=0.85)
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(
    AGS_WT_04_h_vs_AGS_WT_00_h=condAGS.WT.04_h-condAGS.WT.00_h,
    AGS_WT_24_h_vs_AGS_WT_00_h=condAGS.WT.24_h-condAGS.WT.00_h,
    CaCo_2_WT_04_h_vs_CaCo_2_WT_00_h=condCaCo.2.WT.04_h-condCaCo.2.WT.00_h,
    CaCo_2_WT_24_h_vs_CaCo_2_WT_00_h=condCaCo.2.WT.24_h-condCaCo.2.WT.00_h,
    HEK_WT_04_h_vs_HEK_WT_00_h=condHEK.WT.04_h-condHEK.WT.00_h,
    HEK_WT_24_h_vs_HEK_WT_00_h=condHEK.WT.24_h-condHEK.WT.00_h,
    HT29_WT_02_h_vs_HT29_WT_00_h=condHT29.WT.02_h-condHT29.WT.00_h,
    HT29_WT_20_h_vs_HT29_WT_00_h=condHT29.WT.20_h-condHT29.WT.00_h,
    LoVo_WT_02_h_vs_LoVo_WT_00_h=condLoVo.WT.02_h-condLoVo.WT.00_h,
    LoVo_WT_20_h_vs_LoVo_WT_00_h=condLoVo.WT.20_h-condLoVo.WT.00_h,
    THP.1_with_PMA_WT_04_h_vs_THP.1_with_PMA_WT_00_h=condTHP.1_with_PMA.WT.04_h-condTHP.1_with_PMA.WT.00_h,
    THP.1_with_PMA_WT_24_h_vs_THP.1_with_PMA_WT_00_h=condTHP.1_with_PMA.WT.24_h-condTHP.1_with_PMA.WT.00_h,
    THP.1_without_PMA_WT_04_h_vs_THP.1_without_PMA_WT_00_h=condTHP.1_without_PMA.WT.04_h-condTHP.1_without_PMA.WT.00_h,
    THP.1_without_PMA_WT_24_h_vs_THP.1_without_PMA_WT_00_h=condTHP.1_without_PMA.WT.24_h-condTHP.1_without_PMA.WT.00_h,
    BMDM_WT_02_h_vs_BMDM_WT_00_h=condBMDM.WT.02_h-condBMDM.WT.00_h,
    BMDM_WT_06_h_vs_BMDM_WT_00_h=condBMDM.WT.06_h-condBMDM.WT.00_h,
    BMDM_WT_20_h_vs_BMDM_WT_00_h=condBMDM.WT.20_h-condBMDM.WT.00_h,    
    L929_WT_01_h_vs_L929_WT_00_h=condL929.WT.01_h-condL929.WT.00_h,
    L929_WT_06_h_vs_L929_WT_00_h=condL929.WT.06_h-condL929.WT.00_h,
    MEF_WT_04_h_vs_MEF_WT_00_h=condMEF.WT.04_h-condMEF.WT.00_h,
    MEF_WT_20_h_vs_MEF_WT_00_h=condMEF.WT.20_h-condMEF.WT.00_h,
    RAW264.7_WT_04_h_vs_RAW264.7_WT_00_h=condRAW264.7.WT.04_h-condRAW264.7.WT.00_h,
    RAW264.7_WT_24_h_vs_RAW264.7_WT_00_h=condRAW264.7.WT.24_h-condRAW264.7.WT.00_h,
    X3D4_31_WT_01_h_vs_X3D4_31_WT_00_h=condX3D4_31.WT.01_h-condX3D4_31.WT.00_h,
    X3D4_31_WT_02_h_vs_X3D4_31_WT_00_h=condX3D4_31.WT.02_h-condX3D4_31.WT.00_h,
    X3D4_31_WT_06_h_vs_X3D4_31_WT_00_h=condX3D4_31.WT.06_h-condX3D4_31.WT.00_h,
    X3D4_31_WT_16_h_vs_X3D4_31_WT_00_h=condX3D4_31.WT.16_h-condX3D4_31.WT.00_h,
    IPEC.J2_WT_02_h_vs_IPEC.J2_WT_00_h=condIPEC.J2.WT.02_h-condIPEC.J2.WT.00_h,
    IPEC.J2_WT_20_h_vs_IPEC.J2_WT_00_h=condIPEC.J2.WT.20_h-condIPEC.J2.WT.00_h,
    levels=design)


comparisons <- c(
    "AGS_WT_04_h_vs_AGS_WT_00_h",
    "AGS_WT_24_h_vs_AGS_WT_00_h",
    "CaCo_2_WT_04_h_vs_CaCo_2_WT_00_h",
    "CaCo_2_WT_24_h_vs_CaCo_2_WT_00_h",
    "HEK_WT_04_h_vs_HEK_WT_00_h",
    "HEK_WT_24_h_vs_HEK_WT_00_h",
    "HT29_WT_02_h_vs_HT29_WT_00_h",
    "HT29_WT_20_h_vs_HT29_WT_00_h",
    "LoVo_WT_02_h_vs_LoVo_WT_00_h",
    "LoVo_WT_20_h_vs_LoVo_WT_00_h",
    "THP.1_with_PMA_WT_04_h_vs_THP.1_with_PMA_WT_00_h",
    "THP.1_with_PMA_WT_24_h_vs_THP.1_with_PMA_WT_00_h",
    "THP.1_without_PMA_WT_04_h_vs_THP.1_without_PMA_WT_00_h",
    "THP.1_without_PMA_WT_24_h_vs_THP.1_without_PMA_WT_00_h",
    "BMDM_WT_02_h_vs_BMDM_WT_00_h",
    "BMDM_WT_06_h_vs_BMDM_WT_00_h",
    "BMDM_WT_20_h_vs_BMDM_WT_00_h",
    "L929_WT_01_h_vs_L929_WT_00_h",
    "L929_WT_06_h_vs_L929_WT_00_h",
    "MEF_WT_04_h_vs_MEF_WT_00_h",
    "MEF_WT_20_h_vs_MEF_WT_00_h",
    "RAW264.7_WT_04_h_vs_RAW264.7_WT_00_h",
    "RAW264.7_WT_24_h_vs_RAW264.7_WT_00_h",
    "X3D4_31_WT_01_h_vs_X3D4_31_WT_00_h",
    "X3D4_31_WT_02_h_vs_X3D4_31_WT_00_h",
    "X3D4_31_WT_06_h_vs_X3D4_31_WT_00_h",
    "X3D4_31_WT_16_h_vs_X3D4_31_WT_00_h",
    "IPEC.J2_WT_02_h_vs_IPEC.J2_WT_00_h",
    "IPEC.J2_WT_20_h_vs_IPEC.J2_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,
                             comparisons[i]), 
                sep="\t", quote=FALSE, row.names=TRUE, col.names=NA)
}
