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")
10 ngvector_file <- argv[2]
11 data_matrix_file <- argv[3]
12 output_file <- argv[4]
14 nc <- length(argv) - 4;
15 num_reps <- as.numeric(argv[5:(5+nc-1)])
17 .libPaths(c(path, .libPaths()))
20 DataMat <- data.matrix(read.table(data_matrix_file))
22 if (sum(num_reps) != n) stop("Total number of replicates given does not match the number of columns from the data matrix!")
24 conditions <- as.factor(rep(paste("C", 1:nc, sep=""), times = num_reps))
25 Sizes <- MedianNorm(DataMat)
27 if (ngvector_file != "#") {
28 ngvector <- as.vector(data.matrix(read.table(ngvector_file)))
29 stopifnot(!is.null(ngvector))
34 EBOut <- EBTest(Data = DataMat, NgVector = ngvector, Conditions = conditions, sizeFactors = Sizes, maxround = 5)
35 stopifnot(!is.null(EBOut))
37 PP <- as.data.frame(GetPPMat(EBOut))
38 fc_res <- PostFC(EBOut)
40 results <- cbind(PP, fc_res$PostFC, fc_res$RealFC)
41 colnames(results) <- c("PPEE", "PPDE", "PostFC", "RealFC")
42 results <- results[order(results[,"PPDE"], decreasing = TRUE),]
43 write.table(results, file = output_file, sep = "\t")
46 patterns <- GetPatterns(conditions)
47 eename <- rownames(patterns)[which(rowSums(patterns) == nc)]
48 stopifnot(length(eename) == 1)
51 MultiOut <- EBMultiTest(Data = DataMat, NgVector = ngvector, Conditions = conditions, AllParti = patterns, sizeFactors = Sizes, maxround = 5)
52 stopifnot(!is.null(MultiOut))
54 MultiPP <- GetMultiPP(MultiOut)
56 PP <- as.data.frame(MultiPP$PP)
57 pos <- which(names(PP) == eename)
58 probs <- rowSums(PP[,-pos])
60 results <- cbind(PP, MultiPP$MAP[rownames(PP)], probs)
61 colnames(results) <- c(colnames(PP), "MAP", "PPDE")
62 ord <- order(results[,"PPDE"], decreasing = TRUE)
63 results <- results[ord,]
64 write.table(results, file = output_file, sep = "\t")
66 write.table(MultiPP$Patterns, file = paste(output_file, ".pattern", sep = ""), sep = "\t")
68 MultiFC <- GetMultiFC(MultiOut)
69 write.table(MultiFC$CondMeans[ord,], file = paste(output_file, ".condmeans", sep = ""), sep = "\t")