EBSeq/rsem-for-ebseq-find-DE

   1 #!/usr/bin/env Rscript
   2
   3 argv <- commandArgs(TRUE)
   4 if (length(argv) < 6) {
   5   cat("Usage: rsem-for-ebseq-find-DE path ngvector_file data_matrix_file output_file number_of_replicate_for_condition_1 number_of_replicate_for_condition_2 ...\n")
   6   q(status = 1)
   7 }
   8
   9 path <- argv[1]
  10 ngvector_file <- argv[2]
  11 data_matrix_file <- argv[3]
  12 output_file <- argv[4]
  13
  14 nc <- length(argv) - 4;
  15 num_reps <- as.numeric(argv[5:(5+nc-1)])
  16
  17 library(EBSeq, lib.loc = path)
  18
  19 DataMat <- data.matrix(read.table(data_matrix_file))
  20 n <- dim(DataMat)[2]
  21 if (sum(num_reps) != n) stop("Total number of replicates given does not match the number of columns from the data matrix!")
  22
  23 conditions <- as.factor(rep(paste("C", 1:nc, sep=""), times = num_reps))
  24 Sizes <- MedianNorm(DataMat)
  25 ngvector <- NULL
  26 if (ngvector_file != "#") {
  27   ngvector <- as.vector(data.matrix(read.table(ngvector_file)))
  28   stopifnot(!is.null(ngvector))
  29 }
  30
  31 if (nc == 2) {
  32   EBOut <- NULL
  33   EBOut <- EBTest(Data = DataMat, NgVector = ngvector, Conditions = conditions, sizeFactors = Sizes, maxround = 5)
  34   stopifnot(!is.null(EBOut))
  35
  36   PP <- as.data.frame(GetPPMat(EBOut))
  37   fc_res <- PostFC(EBOut)
  38
  39   results <- cbind(PP, fc_res$PostFC, fc_res$RealFC)
  40   colnames(results) <- c("PPEE", "PPDE", "PostFC", "RealFC")
  41   results <- results[order(results[,"PPDE"], decreasing = TRUE),]
  42   write.table(results, file = output_file, sep = "\t")
  43
  44 } else {
  45   patterns <- GetPatterns(conditions)
  46   eename <- rownames(patterns)[which(rowSums(patterns) == nc)]
  47   stopifnot(length(eename) == 1)
  48
  49   MultiOut <- NULL
  50   MultiOut <- EBMultiTest(Data = DataMat, NgVector = ngvector, Conditions = conditions, AllParti = patterns, sizeFactors = Sizes, maxround = 5)
  51   stopifnot(!is.null(MultiOut))
  52
  53   MultiPP <- GetMultiPP(MultiOut)
  54
  55   PP <- as.data.frame(MultiPP$PP)
  56   pos <- which(names(PP) == eename)
  57   probs <- rowSums(PP[,-pos])
  58
  59   results <- cbind(PP, MultiPP$MAP[rownames(PP)], probs)
  60   colnames(results) <- c(colnames(PP), "MAP", "PPDE")
  61   results <- results[order(results[,"PPDE"], decreasing = TRUE),]
  62   write.table(results, file = output_file, sep = "\t")
  63
  64   write.table(MultiPP$Patterns, file = paste(output_file, ".pattern", sep = ""), sep = "\t")
  65 }