For genome BAM, modified MD tag accordingly

[rsem.git] / EBSeq / rsem-for-ebseq-find-DE

diff --git a/EBSeq/rsem-for-ebseq-find-DE b/EBSeq/rsem-for-ebseq-find-DE
index 07b989991cea46e63d064eacaa475c3731227398..f228bab752cc83b2b06bd12b9a71a7250c80b663 100755 (executable)
--- a/EBSeq/rsem-for-ebseq-find-DE
+++ b/EBSeq/rsem-for-ebseq-find-DE
@@ -14,7 +14,8 @@ output_file <- argv[4]
 nc <- length(argv) - 4;
 num_reps <- as.numeric(argv[5:(5+nc-1)])
 
-library(EBSeq, lib.loc = path)
+.libPaths(c(path, .libPaths()))
+library(EBSeq)
 
 DataMat <- data.matrix(read.table(data_matrix_file))
 n <- dim(DataMat)[2]
@@ -58,8 +59,12 @@ if (nc == 2) {
 
   results <- cbind(PP, MultiPP$MAP[rownames(PP)], probs)
   colnames(results) <- c(colnames(PP), "MAP", "PPDE")  
-  results <- results[order(results[,"PPDE"], decreasing = TRUE),]
+  ord <- order(results[,"PPDE"], decreasing = TRUE)
+  results <- results[ord,]
   write.table(results, file = output_file, sep = "\t")
 
   write.table(MultiPP$Patterns, file = paste(output_file, ".pattern", sep = ""), sep = "\t")
+
+  MultiFC <- GetMultiFC(MultiOut)
+  write.table(MultiFC$CondMeans[ord,], file = paste(output_file, ".condmeans", sep = ""), sep = "\t")
 }