#!/usr/bin

library(DESeq)
library(hexbin)
library(latticeExtra)
library(gplots)
library(geneplotter)
library(Biobase)

##In a file called test_args.R
args <- commandArgs()


fileName <- args[4]
colNames <- as.integer(unlist(strsplit(args[5], ",")))
colCond1 <- as.integer(unlist(strsplit(args[6], ",")))
colCond2 <- as.integer(unlist(strsplit(args[7], ",")))
OUTPUTCSV <- args[8]
OUTPUTPNG <- args[9]

if(colNames[1]!=0){
	countsTable <- read.delim(fileName, row.names=1)
	conditions <- c((colNames[length(colNames)]+1):ncol(countsTable))
} else if(colNames[1]==0){
	countsTable <- read.delim(fileName)
	conditions <- c(1:ncol(countsTable))
	rownames(countsTable) <- paste( "Gene", 1:nrow(countsTable), sep="_" )}

for(i in colCond1){conditions[i] = "A"}
for(i in colCond2){conditions[i] = "B"}
conditions
#analysis with DESeq
cds <- newCountDataSet( countsTable, conditions )
cds <- estimateSizeFactors( cds )
cds <- estimateVarianceFunctions( cds )
result <- nbinomTest( cds, "A", "B" )
#stock the result dans un .tsv as output file
write.table(result, OUTPUTCSV, sep = "	", quote = FALSE, col.names = NA)

#figures for DE analysis
#pdf( OUTPUTPNG, width=4, height=4 )
png( filename=OUTPUTPNG, width=700, height=700 )
#png format is not as clear as pdf format!!!!!!!!!!!!!!!!!!!!!!!!!
print(xyplot(
	log2FoldChange ~ I(baseMean),
	result,
	pch=16, cex=.3,
	col=ifelse(result$padj < .1, "#FF000040","#00000040" ),
	panel = function( x, y, col, ...) {
		above <- (y > 5.8)
		below <- (y < -5.8)
		inside <- !( above | below )
		panel.xyplot( x=x[inside], y=y[inside], col=col[inside], ...)
		panel.arrows( x[above], 5.8, x[above], 5.95, col=col[above],length=".1", unit="native" )
		panel.arrows( x[below], -5.8, x[below], -5.95, col=col[below],length=".1", unit="native" ) },
	axis = function(side, ...) {
		if( side=="left") {
			panel.axis( side, outside=TRUE, at=seq(-14,14,by=1), labels=FALSE )
			panel.axis( side, outside=TRUE, at=seq(-10,10,by=5), labels=TRUE )
		}
		if( side=="bottom") {
			panel.axis( side, outside=TRUE, at=seq(-2,10,by=1), rot=0,
			labels = do.call( expression,
				lapply( seq(-2,10,by=1), function(a)
					substitute( 10^b, list(b=a) ) ) ) )
		} },
	xlab = "mean", ylab = "log2 fold change",
	scales = list(x = list( log=TRUE ),y = list( log=FALSE, limits=c( -6, 6 ) ) ) ))
dev.off()

#The volcano plot 
#pdf( "vulcano_fly.pdf", width=4, height=4 )
#print(xyplot( -log10( pval ) ~ log2FoldChange,
#	result,
#	pch=20, cex=.2,
#	col=ifelse( result$padj<.1, "#FF000050", "#00000050" ),
#	axis = function( side, ... ) {
#		if( side=="bottom") {
#			panel.axis( side, outside=TRUE, at=seq(-14,14,by=1), labels=FALSE )
#			panel.axis( side, outside=TRUE, at=seq(-10,10,by=5), labels=TRUE )
#			}
#		if( side=="left") {
#			panel.axis( side, outside=TRUE, at=seq(0,25,by=1), labels=FALSE )
#			panel.axis( side, outside=TRUE, at=seq(0,25,by=5),
#				labels = do.call( expression,
#					lapply( seq(0,25,by=5), function(a)
#						substitute( 10^-b, list(b=a) ) ) ) )
#		} },
#	xlab = "log2 fold change", ylab = "p value",
#	scales = list(
#		x = list( limits=c( -6, 6 ) ),
#		y = list( limits=c( 0, 25 ) ) ) ))
#dev.off()