--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dexseq-hts_1.0/src/run_DEXseq.R	Tue Oct 08 08:22:45 2013 -0400
@@ -0,0 +1,73 @@
+### load DEXSeq package
+suppressMessages(require("DEXSeq"))
+
+### get arguments 1: INFILE, 2: OUTFILE 3:SIZE
+args <- commandArgs()
+INFILE<-args[4]
+EXTRAPATH<-args[5]
+annodb<-args[6]
+OUTFILE<-args[7]
+
+INFILE_CONDS=c(paste(INFILE, "_CONDITIONS.tab", sep=""))
+
+### read count data from file
+condsTable <- read.delim( INFILE_CONDS, header=TRUE, stringsAsFactors=TRUE, row.names=1 )
+#head(condsTable)
+
+conditions<-factor( condsTable[ , 1] )
+#print(conditions)
+
+## unique condition to define the pair of tests 
+uniq_conds <- unique(conditions)
+#print(uniq_conds)
+
+## all possible pairs of conditions
+pw_tests <- list()
+for(i in 1:(length(uniq_conds)-1)) {
+    for(j in (i+1):length(uniq_conds)) {
+        pw_tests[[length(pw_tests)+1]] <- c(uniq_conds[i], uniq_conds[j])
+    }
+}
+#print(pw_tests)
+
+tab <- NULL
+## testing all possible pairs of conditions
+for(i in 1:length(pw_tests)) {
+    ## header name 
+    test_pair_name <- c(paste(pw_tests[[i]][1], "__vs__", pw_tests[[i]][2], sep=""))
+    print(test_pair_name)
+
+    sub.data <- subset(condsTable, (conditions %in% c(pw_tests[[i]][1],pw_tests[[i]][2])))
+    sub.data[[1]]<-as.factor(sub.data[[1]])
+
+    ecs = read.HTSeqCounts(countfiles=file.path(EXTRAPATH, row.names(sub.data)), design=sub.data, flattenedfile=annodb)
+    
+    ## Normalisation
+    ecs <- estimateSizeFactors(ecs)
+    print(sizeFactors(ecs))
+
+    suppressMessages(require("parallel"))
+    ## Dispersion estimation 
+    ecs <- estimateDispersions(ecs, nCores=2)
+    ecs <- fitDispersionFunction(ecs)
+
+    ## Testing for differential exon usage
+    ecs <- testForDEU(ecs, nCores=2)
+
+    ## estimate the fold change 
+    ecs <- estimatelog2FoldChanges(ecs, nCores=2)
+
+    ## Result table
+    resultTable <- DEUresultTable(ecs)
+    print(head(resultTable))
+
+    colnames(resultTable) <- paste(test_pair_name, colnames(resultTable), sep=":")
+    #print(colnames(resultTable))
+
+    if(is.null(tab)) {
+        tab<- resultTable
+    }
+    else tab<- cbind(tab, resultTable)
+}
+## printing the result to out file 
+write.table(tab, OUTFILE, quote=F, sep="\t", row.names=F)