diff deseq-hts_1.0/src/difftest_deseq.R @ 9:e27b4f7811c2 draft

Updated DESeq version 1.12

--- a/deseq-hts_1.0/src/difftest_deseq.R	Wed Jun 27 15:38:39 2012 -0400
+++ b/deseq-hts_1.0/src/difftest_deseq.R	Tue Oct 08 08:09:28 2013 -0400
@@ -1,4 +1,5 @@
-library( DESeq )
+### load DESeq package
+suppressMessages(require("DESeq"))
 
 ### get arguments 1: INFILE, 2: OUTFILE 3:SIZE
 args <- commandArgs()
@@ -31,26 +32,39 @@
 {
     cds <- estimateDispersions( cds )
 } else {
-    writeLines("\nYou did not enter any replicates! - The results may be less valuable without replicates!\n")
+    writeLines("\n***You did not enter any replicates! - The results may be less valuable without replicates!***\n")
     cds <- estimateDispersions( cds, method='blind', sharingMode='fit-only')
 }
 experiments <- levels(conds)
 
-res<-c()
+res_1<-c()
+res_2<-c()
+res_3<-c()
+res_4<-c()
+res_5<-c()
+res_6<-c()
+res_7<-c()
+res_8<-c()
 table_col_names<-c()
+
 for (i in 1:(length(experiments)-1))
 {
    for( j in (i+1):(length(experiments)))
    {
        print(c(i,j))
        tempres <- nbinomTest(cds,experiments[i],experiments[j])
-       res = cbind(res,tempres[,7])
-       #res = cbind(res,tempres[,8])
-       table_col_names = cbind(table_col_names,paste('cond_', experiments[i], '_vs._cond_', experiments[j], sep='')) 
+       res_1 = cbind(res_1,tempres[,1])
+       res_2 = cbind(res_2,tempres[,2])
+       res_3 = cbind(res_3,tempres[,3])
+       res_4 = cbind(res_4,tempres[,4])
+       res_5 = cbind(res_5,tempres[,5])
+       res_6 = cbind(res_6,tempres[,6])
+       res_7 = cbind(res_7,tempres[,7])
+       res_8 = cbind(res_8,tempres[,8])
+       table_col_names = cbind(table_col_names,paste('cond_', experiments[i], '_vs._cond_', experiments[j], sep='', 'test')) 
    }
 }
 
-DiffTable<-res
-rownames(DiffTable)<-rownames(countsTable)
-colnames(DiffTable)<-table_col_names
+DiffTable<-cbind(res_1,res_2,res_3,res_4,res_5,res_6,res_7,res_8)
+colnames(DiffTable)<-c('feature ID', 'base  mean', 'base mean A', 'base mean B', 'fold change', 'log2 fold change','p value', 'adjusted p value')
 write.table(DiffTable, file = OUTFILE, quote = FALSE, sep ="\t", eol ="\n", na = "1.000", dec = ".", row.names = TRUE,col.names =TRUE)