--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SMART/DiffExpAnal/testR.R	Tue Apr 30 14:33:21 2013 -0400
@@ -0,0 +1,93 @@
+#!/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()