roll back to original version plus fixed a small bug in calcCI.cpp

[rsem.git] / rsem-calculate-expression

diff --git a/rsem-calculate-expression b/rsem-calculate-expression
index bc976c32269a2938cb0607dabc36157e58509ffb..00f26cee381fdf03f6028748bda1caa1aa9190ae 100755 (executable)
--- a/rsem-calculate-expression
+++ b/rsem-calculate-expression
@@ -8,7 +8,7 @@ use strict;
 
 #const
 my $BURNIN = 200;
-my $CHAINLEN = 1000;
+my $NCV = 1000;
 my $SAMPLEGAP = 1;
 my $CONFIDENCE = 0.95;
 my $NSPC = 50;
@@ -191,16 +191,14 @@ if (!$is_sam && !$is_bam) {
     else { print "Oh, no!!!"; exit(2); }
     
     $command .= " -n $C -e $E -l $L";
-    
     if ($read_type == 2 || $read_type == 3) { $command .= " -I $minL -X $maxL"; }
+    if ($chunkMbs > 0) { $command .= " --chunkmbs $chunkMbs"; }
     
     if ($strand_specific || $probF == 1.0) { $command .= " --norc"; }
     elsif ($probF == 0.0) { $command .= " --nofw"; }
 
     $command .= " -p $nThreads -a -m $maxHits -S";
-    if ($quiet) { $command .= " --quiet"; }
-    
-    $command .= " --chunkmbs $chunkMbs" if $chunkMbs > 0;
+    if ($quiet) { $command .= " --quiet"; }    
 
     $command .= " $refName";
     if ($read_type == 0 || $read_type == 1) {
@@ -210,7 +208,7 @@ if (!$is_sam && !$is_bam) {
        $command .= " -1 $mate1_list -2 $mate2_list";
     }
 
-    $command .= " | gzip > $imdName.sam.gz";
+    $command .= " | gzip > $sampleName.sam.gz";
 
     if ($mTime) { $time_start = time(); }
 
@@ -218,7 +216,7 @@ if (!$is_sam && !$is_bam) {
 
     if ($mTime) { $time_end = time(); $time_alignment = $time_end - $time_start; }
 
-    $inpF = "$imdName.sam.gz";
+    $inpF = "$sampleName.sam.gz";
     $is_sam = 1; # output of bowtie is a sam file
 }
 
@@ -293,8 +291,8 @@ 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";
+    $command = $dir."rsem-run-gibbs $refName $sampleName $sampleToken $BURNIN $NCV $SAMPLEGAP";
+    $command .= " -p $nThreads";
     if ($quiet) { $command .= " -q"; }
     &runCommand($command);
 
@@ -303,7 +301,8 @@ if ($calcCI) {
     &collectResults("$imdName.iso_res", "$sampleName.isoforms.results"); # isoform level
     &collectResults("$imdName.gene_res", "$sampleName.genes.results"); # gene level
 
-    $command = $dir."rsem-calculate-credibility-intervals $refName $sampleName $sampleToken $CONFIDENCE $NSPC $NMB";
+    $command = $dir."rsem-calculate-credibility-intervals $refName $sampleName $sampleToken $CONFIDENCE $NCV $NSPC $NMB";
+    $command .= " -p $nThreads";
     if ($quiet) { $command .= " -q"; }
     &runCommand($command);
 
@@ -325,11 +324,11 @@ 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";
+    print OUTPUT "Aligning reads: $time_alignment s.\n";
+    print OUTPUT "Estimating expression levels: $time_rsem s.\n";
+    print OUTPUT "Calculating credibility intervals: $time_ci s.\n";
     my $time_del = $time_end - $time_start;
-    print OUTPUT "Delete: $time_del s.\n";
+#    print OUTPUT "Delete: $time_del s.\n";
     close(OUTPUT);
 }
 
@@ -471,7 +470,7 @@ Generate a BAM file, 'sample_name.genome.bam', with alignments mapped to genomic
 =item B<--sampling-for-bam>
 
 When RSEM generates a BAM file, instead of outputing all alignments a read has with their posterior probabilities, one alignment is sampled and outputed according to the posterior probabilities. If the sampling result is that the read comes from the "noise" transcript, nothing is outputed. (Default: off)
- 
+
 =item B<--calc-ci>
 
 Calculate 95% credibility intervals and posterior mean estimates.  (Default: off)
@@ -546,12 +545,16 @@ Number of bins in the RSPD. Only relevant when '--estimate-rspd' is specified.
 
 =item B<--ci-memory> <int>
 
-Amount of memory (in MB) RSEM is allowed to use for computing credibility intervals. (Default: 1024)
+Maximum size (in memory, MB) of the auxiliary buffer used for computing credibility intervals (CI). Set it larger for a faster CI calculation. However, leaving 2 GB memory free for other usage is recommended. (Default: 1024)
 
 =item B<--keep-intermediate-files>
 
 Keep temporary files generated by RSEM.  RSEM creates a temporary directory, 'sample_name.temp', into which it puts all intermediate output files. If this directory already exists, RSEM overwrites all files generated by previous RSEM runs inside of it. By default, after RSEM finishes, the temporary directory is deleted.  Set this option to prevent the deletion of this directory and the intermediate files inside of it. (Default: off)
 
+=item B<--time>
+
+Output time consumed by each step of RSEM to 'sample_name.time'. (Default: off)
+
 =item B<-q/--quiet>
 
 Suppress the output of logging information. (Default: off)
@@ -649,6 +652,18 @@ indicating the strand of the transcript it aligns to.
 'sample_name.genome.sorted.bam' and 'sample_name.genome.sorted.bam.bai' are the
 sorted BAM file and indices generated by samtools (included in RSEM package).
 
+=item B<sample_name.sam.gz>
+
+Only generated when the input files are raw reads instead of SAM/BAM format files
+
+It is the gzipped SAM output produced by bowtie aligner.
+
+=item B<sample_name.time>
+
+Only generated when --time is specified.
+
+It contains time (in seconds) consumed by aligning reads, estimating expression levels and calculating credibility intervals.
+
 =item B<sample_name.stat>
 
 This is a folder instead of a file. All model related statistics are stored in this folder. Use 'rsem-plot-model' can generate plots using this folder.