rsem-plot-model updated

diff --git a/PairedEndModel.h b/PairedEndModel.h
index 6e0f1a4fb385f7dbc53f904c85c95508f7f933a8..7c8737465c959da2e3089e39d0fcb091ab3e52ff 100644 (file)
--- a/PairedEndModel.h
+++ b/PairedEndModel.h
@@ -385,9 +385,8 @@ bool PairedEndModel::simulate(int rid, PairedEndRead& read, int& sid) {
                readseq2 = pro->simulate(sampler, mateL2, m2pos, m2dir, ref);
        }
 
-       std::ostringstream stdout;
-       stdout<<rid<<"_"<<dir<<"_"<<sid<<"_"<<pos<<"_"<<insertL;
-       name = stdout.str();
+       strout<<rid<<"_"<<dir<<"_"<<sid<<"_"<<pos<<"_"<<insertL;
+       name = strout.str();
 
        read = PairedEndRead(SingleRead(name + "/1", readseq1), SingleRead(name + "/2", readseq2));
 

diff --git a/PairedEndQModel.h b/PairedEndQModel.h
index 2b899c57fa6061e0974e56083da03f580766351d..dee535dda56ef89113b0ac8d06e470347f1e2249 100644 (file)
--- a/PairedEndQModel.h
+++ b/PairedEndQModel.h
@@ -199,8 +199,6 @@ public:
        }
 
        int getModelType() const { return model_type; }
-
-       bool simulate(int, PairedEndReadQ&, SingleReadQ&, int&);
  
 private:
        static const int model_type = 3;
@@ -405,9 +403,8 @@ bool PairedEndQModel::simulate(int rid, PairedEndReadQ& read, int& sid) {
                readseq2 = qpro->simulate(sampler, mateL2, m2pos, m2dir, qual2, ref);
        }
 
-       std::ostringstream stdout;
-       stdout<<rid<<"_"<<dir<<"_"<<sid<<"_"<<pos<<"_"<<insertL;
-       name = stdout.str();
+       strout<<rid<<"_"<<dir<<"_"<<sid<<"_"<<pos<<"_"<<insertL;
+       name = strout.str();
 
        read = PairedEndReadQ(SingleReadQ(name + "/1", readseq1, qual1), SingleReadQ(name + "/2", readseq2, qual2));
 

diff --git a/SingleModel.h b/SingleModel.h
index 49c103b956b9fa389edac2ad4a4cefcce000a2cc..756e3d5cdecf386339eb90c303e9aa6b908a26e5 100644 (file)
--- a/SingleModel.h
+++ b/SingleModel.h
@@ -423,9 +423,8 @@ bool SingleModel::simulate(int rid, SingleRead& read, int& sid) {
                }
        }
 
-       std::ostringstream stdout;
-       stdout<<rid<<"_"<<dir<<"_"<<sid<<"_"<<pos;
-       name = stdout.str();
+       strout<<rid<<"_"<<dir<<"_"<<sid<<"_"<<pos;
+       name = strout.str();
 
        read = SingleRead(name, readseq);
 

diff --git a/SingleQModel.h b/SingleQModel.h
index d2c8ca1fdb67f8075e699db44f097c5502c12419..f01920f2bfcf908e26aad049ca8bf8222733ec89 100644 (file)
--- a/SingleQModel.h
+++ b/SingleQModel.h
@@ -443,9 +443,8 @@ bool SingleQModel::simulate(int rid, SingleReadQ& read, int& sid) {
                }
        }
 
-       std::ostringstream stdout;
-       stdout<<rid<<"_"<<dir<<"_"<<sid<<"_"<<pos;
-       name = stdout.str();
+       strout<<rid<<"_"<<dir<<"_"<<sid<<"_"<<pos;
+       name = strout.str();
 
        read = SingleReadQ(name, readseq, qual);
 

diff --git a/rsem-calculate-expression b/rsem-calculate-expression
index 06604d96a8c415b4eaeaf9e74854d6bf772a8d17..93fabce4686c8d0b7e4fafeae279c6262e7d412e 100755 (executable)
--- a/rsem-calculate-expression
+++ b/rsem-calculate-expression
@@ -15,7 +15,7 @@ my $NSPC = 50;
 
 my $NMB = 1024; # default
 
-my $status;
+my $status = 0;
 
 my $read_type = 1; # default, single end with qual
 
@@ -55,6 +55,9 @@ my $keep_intermediate_files = 0;
 
 my $strand_specific = 0;
 
+my $mTime = 0;
+my ($time_start, $time_end, $time_alignment, $time_rsem, $time_ci) = (0, 0, 0, 0, 0);
+
 GetOptions("keep-intermediate-files" => \$keep_intermediate_files,
           "no-qualities" => \$no_qual,
           "paired-end" => \$paired_end,
@@ -82,6 +85,7 @@ GetOptions("keep-intermediate-files" => \$keep_intermediate_files,
           "out-bam" => \$genBamF,
           "calc-ci" => \$calcCI,
           "ci-memory=i" => \$NMB,
+          "time" => \$mTime,
           "q|quiet" => \$quiet,
           "h|help" => \$help) or pod2usage(-exitval => 2, -verbose => 2);
 
@@ -189,7 +193,13 @@ if (!$is_sam && !$is_bam) {
 
     $command .= " | gzip > $imdName.sam.gz";
     print "$command\n";
+
+    if ($mTime) { $time_start = time(); }
+
     $status = system($command);
+
+    if ($mTime) { $time_end = time(); $time_alignment = $time_end - $time_start; }
+
     if ($status != 0) {
        print "bowtie failed! Please check if you provide correct parameters/options for the pipeline!\n";
        exit(-1);
@@ -200,6 +210,8 @@ if (!$is_sam && !$is_bam) {
     $is_sam = 1; # output of bowtie is a sam file
 }
 
+if ($mTime) { $time_start = time(); }
+
 $command = $dir."rsem-parse-alignments $refName $sampleName $sampleToken";
 
 my $samInpType;
@@ -234,8 +246,8 @@ if ($status != 0) {
 }
 print "\n";
 
-$status = open(OUTPUT, ">$imdName.mparams");
-if ($status == 0) { print "Cannot generate $imdName.mparams!\n"; exit(-1); }
+my $doesOpen = open(OUTPUT, ">$imdName.mparams");
+if ($doesOpen == 0) { print "Cannot generate $imdName.mparams!\n"; exit(-1); }
 print OUTPUT "$minL $maxL\n";
 print OUTPUT "$probF\n";
 print OUTPUT "$estRSPD\n";
@@ -284,6 +296,10 @@ if ($genBamF) {
 &collectResults("$imdName.iso_res", "$sampleName.isoforms.results"); # isoform level
 &collectResults("$imdName.gene_res", "$sampleName.genes.results"); # gene level
 
+if ($mTime) { $time_end = time(); $time_rsem = $time_end - $time_start; }
+
+if ($mTime) { $time_start = time(); }
+
 if ($calcCI) {
     $command = $dir."rsem-run-gibbs $refName $sampleName $sampleToken $BURNIN $CHAINLEN $SAMPLEGAP";
 #    $command .= " -p $nThreads";
@@ -317,14 +333,30 @@ if ($calcCI) {
     &collectResults("$imdName.gene_res", "$sampleName.genes.results"); # gene level
 }
 
+if ($mTime) { $time_end = time(); $time_ci = $time_end - $time_start; }
+
+if ($mTime) { $time_start = time(); }
+
 if (!$keep_intermediate_files) {
-    $status = system ("rm -rf $temp_dir");
+    $status = system("rm -rf $temp_dir");
     if ($status != 0) {
        print "Fail to delete the temporary folder!\n";
        exit(-1);
     }
 }
 
+if ($mTime) { $time_end = time(); }
+
+if ($mTime) { 
+    open(OUTPUT, ">$sampleName.time");
+    print OUTPUT "Alignment: $time_alignment s.\n";
+    print OUTPUT "RSEM: $time_rsem s.\n";
+    print OUTPUT "CI: $time_ci s.\n";
+    my $time_del = $time_end - $time_start;
+    print OUTPUT "Delete: $time_del s.\n";
+    close(OUTPUT);
+}
+
 # inpF, outF
 sub collectResults {
     my $local_status;

diff --git a/rsem-plot-model b/rsem-plot-model
index c15a75f230c018da6ffc0cde6171f6d1f6fa0686..225d08f06bad85c9365fde9939632174d47122d5 100755 (executable)
--- a/rsem-plot-model
+++ b/rsem-plot-model
@@ -2,15 +2,21 @@
 
 argv <- commandArgs(TRUE)
 if (length(argv) != 2) {
-  cat("Usage: rsem-plot-model sample_name outF\n")
+  cat("Usage: rsem-plot-model sample_name output_plot_file\n")
   q(status = 1)
 }
 
 strvec <- strsplit(argv[1], split = "/")[[1]]
 token <- strvec[length(strvec)]
 
-modelF <- paste(argv[1], ".stat/", token, ".model", sep = "")
-cntF <- paste(argv[1], ".stat/", token, ".cnt", sep = "")
+stat.dir <- paste(argv[1], ".stat", sep = "")
+if (!file.exists(stat.dir)) {
+  cat("Error: directory does not exist: ", stat.dir, "\n", sep = "")
+  cat(strwrap("This version of rsem-plot-model only works with the output of RSEM versions >= 1.1.8"), sep="\n")
+  q(status = 1)
+}
+modelF <- paste(stat.dir, "/", token, ".model", sep = "")
+cntF <- paste(stat.dir, "/", token, ".cnt", sep = "")
 
 pdf(argv[2])
 
@@ -27,7 +33,11 @@ x <- (vec[1] + 1) : vec[2]
 y <- as.numeric(strsplit(strvec[2], split = " ")[[1]])
 mean <- weighted.mean(x, y)
 std <- sqrt(weighted.mean((x - mean)^2, y))
-plot(x, y, type = "h", main = "Fragment Length Distribution", sub = paste("Mean = ", mean, ", Std = ", std), xlab = "Fragment Length", ylab = "Probability")  
+plot(x, y, type = "h",
+     main = "Fragment Length Distribution",
+     sub = paste("Mean = ", round(mean, 1), ", Std = ", round(std, 1)),
+     xlab = "Fragment Length",
+     ylab = "Probability")
 
 # mate length distribution
 if (model_type == 0 || model_type == 1) bval <- as.numeric(readLines(con, n = 1)[1]) else bval <- 1
@@ -40,7 +50,11 @@ if (bval == 1) {
   y <- as.numeric(list[[2]])
   mean <- weighted.mean(x, y)
   std <- sqrt(weighted.mean((x - mean)^2, y))
-  plot(x, y, type = "h", main = "Mate Length Distribution", sub = paste("Mean = ", mean, ", Std = ", std), xlab = "Mate Length", ylab = "Probability")  
+  plot(x, y, type = "h",
+       main = "Read Length Distribution",
+       sub=paste("Mean = ", round(mean, 1), ", Std = ", round(std, 1)),
+       xlab = "Read Length",
+       ylab = "Probability")
 }
 strvec <- readLines(con, n = 1)
 
@@ -88,7 +102,10 @@ if (model_type == 1 || model_type == 3) {
     peT <- c(peT, ifelse(sum(vecT) < 1e-8, NA, -10 * log10(1.0 - vecT[4])))
   }
 
-  matplot(x, cbind(peA, peC, peG, peT), type = "b", lty = 1:4, pch = 0:3, col = 1:4, main = "Phred Quality Score vs. Observed Quality", xlab = "Quality Score", ylab = "Observed Quality")
+  matplot(x, cbind(peA, peC, peG, peT), type = "b", lty = 1:4, pch = 0:3, col = 1:4,
+          main = "Observed Quality vs. Phred Quality Score",
+          xlab = "Phred Quality Score",
+          ylab = "Observed Quality")
   legend("topleft", c("A", "C", "G", "T"), lty = 1:4, pch = 0:3, col = 1:4)
 } else {
   # Profile
@@ -124,6 +141,9 @@ if (model_type == 1 || model_type == 3) {
 close(con)
 
 pair <- read.table(file = cntF, skip = 3, sep = "\t")
-barplot(pair[,2], names.arg = pair[,1], xlab = "Number of Alignments", ylab = "Number of Reads", main = "Histogram of Reads with Different Number of  Alignments")
+barplot(pair[,2], names.arg = pair[,1],
+        xlab = "Alignments per fragment",
+        ylab = "Number of fragments",
+        main = "Histogram of Alignments per Fragment")
 
 dev.off()

diff --git a/simulation.cpp b/simulation.cpp
index 19f2b928b247f4557b24414f80d9eeda029743da..e516cc15c32848e143aaf4b3a56a1e738fbe3456 100644 (file)
--- a/simulation.cpp
+++ b/simulation.cpp
@@ -138,19 +138,19 @@ void simulate(char* modelF, char* resultsF) {
        fin.close();
        for (int i = 1; i <= M; i++) theta[i] = theta[i] / denom * (1.0 - theta[0]);
        
-       int tmp2140 = 0;
+       int resimulation_count = 0;
 
        //simulating...
        model.startSimulation(&sampler, theta);
        for (int i = 0; i < N; i++) {
-         while (!model.simulate(i, read, sid)) { ++tmp2140; }//;
+               while (!model.simulate(i, read, sid)) { ++resimulation_count; }
                read.write(n_os, os);
                ++counts[sid];
                if ((i + 1) % 1000000 == 0 && verbose) printf("GEN %d\n", i + 1);
        }
        model.finishSimulation();
 
-       printf("Total number of resimulation is %d\n", tmp2140);
+       printf("Total number of resimulation is %d\n", resimulation_count);
 }
 
 void writeResFiles(char* outFN) {